Thalamic deep brain stimulation paradigm to reduce consciousness: Cortico-striatal dynamics implicated in mechanisms of consciousness

Anesthetic manipulations provide much-needed causal evidence for neural correlates of consciousness, but non-specific drug effects complicate their interpretation. Evidence suggests that thalamic deep brain stimulation (DBS) can either increase or decrease consciousness, depending on the stimulation target and parameters. The putative role of the central lateral thalamus (CL) in consciousness makes it an ideal DBS target to manipulate circuit-level mechanisms in cortico-striato-thalamic (CST) systems, thereby influencing consciousness and related processes. We used multi-microelectrode DBS targeted to CL in macaques while recording from frontal, parietal, and striatal regions. DBS induced episodes of abnormally long, vacant staring with low-frequency oscillations here termed vacant, perturbed consciousness (VPC). DBS modulated VPC likelihood in a frequency-specific manner. VPC events corresponded to decreases in measures of neural complexity (entropy) and integration (Φ*), proposed indices of consciousness, and substantial changes to communication in CST circuits. During VPC, power spectral density and coherence at low frequencies increased across CST circuits, especially in thalamo-parietal and cortico-striatal pathways. Decreased consciousness and neural integration corresponded to shifts in cortico-striatal network configurations that dissociated parietal and subcortical structures. Overall, the features of VPC and implicated networks were similar to those of absence epilepsy. As this same multi-microelectrode DBS method–but at different stimulation frequencies–can also increase consciousness in anesthetized macaques, it can be used to flexibly address questions of consciousness with limited confounds, as well as inform clinical investigations of other consciousness disorders.


Introduction
Deep brain stimulation (DBS) is a powerful tool to treat severe conditions impervious to other interventions, yet the mechanisms remain unclear [1]. DBS may lead to neural excitation/inhibition, regularized activity or altered neural communication. Further, effectiveness varies depending on stimulation location, parameters, and electrode type [2].
Intralaminar thalamic DBS is an effective tool to increase consciousness in minimally conscious state patients [3] and anesthetized primates [4,5]. Of the intralaminar nuclei, the central lateral nucleus (CL) seems the most promising target, with anatomical connectivity well-suited to modulate consciousness. In addition to thalamocortical interactions, theories of consciousness implicate interactions within and/or between frontal and parietal cortex [6][7][8][9]. CL projects to both areas, providing input to superficial layers, and sharing reciprocal projections with deep layers [10,11], which have been shown to be critical for consciousness [4,12,13]. CL receives input from the reticular activating system and projects directly to the basal ganglia [14,15], serving as both a key input and output of the structure [16]. While the role of the basal ganglia in consciousness is debated [17][18][19], the striatum contributes to integrated information [20], a measure of neural complexity associated with consciousness [6], and contributes strongly to decoding neural differences between conscious states [20]. Further, the basal ganglia are linked to altered consciousness with hallucinogens [21] and are suppressed during general anesthesia [22] and absence epilepsy [23,24].
Evidence suggests CL manipulations can bidirectionally influence consciousness. CL lesions are linked to consciousness disorders, like spatial neglect and coma [25]; and thalamic DBS in cats has been reported to induce sleep, produce a syndrome of vacant staring with behavioral arrest, or absence seizures, depending on the DBS method (e.g., stimulation frequency, current amplitude and electrode placement) [26]. In rodents, optogenetic stimulation of CL at 40Hz or 100Hz can drive arousal and activity in cortico-striatal-thalamic (CST) networks during sleep, while 10Hz leads to behavioral arrest [27]. Similarly, using a multi-microelectrode stimulation method (rather than applying current via 1-2 large electrode contacts in conventional DBS) in macaques (Fig 1), we have shown that DBS specific to CL, more than neighboring thalamic regions, awakens animals from anesthesia when applied at 50Hz, rather than 10Hz or 200Hz [4]. These frequency-specific effects may relate to natural rhythms of CL. In cat [28,29] and monkey [4], a subset of CL neurons exhibit high firing rate (around 50Hz) during wakefulness, but lower firing during sleep and anesthesia. Thus, frequency-specific effects with respect to consciousness may relate to stimulation parameters mimicking these firing patterns [4]. We hypothesize that stimulation frequencies mimicking CL activity patterns during sleep/anesthesia, or inducing abnormal CL activity, will disrupt consciousness in awake, healthy primates.
To test this hypothesis, we used a CL-tailored method of DBS in resting-state (wake) and rewarded (fixation task) conditions while simultaneously recording from frontal cortex, parietal cortex, and striatum of macaques; we also recorded from CL in a DBS-free condition. We applied DBS at different frequencies (10,50,200Hz) across a series of stimulation experiments, in addition to experimental series without DBS, to test effect specificity (Fig 2). We predicted that thalamic DBS would cause perturbations in consciousness in a frequency-dependent manner: Stimulations at 50 Hz, mimicking wake-state CL activity patterns, would have little effect on waking behavior; but 10/200Hz DBS (atypical frequencies for wakefulness) would perturb consciousness. Moreover, reduced consciousness would relate to aberrant communication along CST pathways, thought to play a role in consciousness disorders [18,30].

Ethics statement
All procedures conformed to the National Institutes of Health Guide for the Care and Use of Laboratory Animals and were approved by the University of Wisconsin-Madison Institutional Animal Care and Use Committee (G005777). All surgical and invasive procedures were performed using sterile techniques. Isoflurane was used for general anesthesia during head implant surgeries and buprenorphine was used for post-operative analgesia. Ketamine and dexmedetomidine were used for anesthesia during minor procedures and meloxicam for analgesia. Animal health and welfare was monitored multiple times daily by experimenters, veterinary and/or animal care staff at the Wisconsin National Primate Research Center.

Model
We acquired data from two male monkeys (Macaca mulatta, 4.3-5.5 years old, 7.63-10.30 kg body weight) housed at the Wisconsin National Primate Research Center.

Surgery and electrode placement
We performed a stereotaxic head implant and craniotomy surgery using aseptic techniques under isoflurane general anesthesia (1-2%) according to the procedures described previously [4,20]. We placed 2.5mm craniotomies over the frontal eye fields (FEF), lateral intraparietal area (LIP), CL, and caudate nucleus (CN) using stereotaxic measurements based on high-quality structural MRI acquired in advance of the surgery and comparisons to a stereotaxic atlas [31]. We also acquired MRI subsequent to surgery with electrodes in situ to verify electrode placements. We used a GE MR750 3T scanner to acquire 3D T1-weighted structural images from the anesthetized monkeys, with an inversion-recovery prepared gradient echo sequence (FOV = 128 mm 2 ; matrix = 256 x 256; no. of slices = 166; 0.5 mm isotropic; TR = 9.68 ms; TE = 4.192 ms; flip angle = 12˚; inversion time (TI) = 450 ms). We averaged 6-10 T1-weighted images for the pre-surgery high-quality structural image of each monkey, and we averaged 2 T1-weighted images to visualize electrodes in situ. Further imaging details have been described previously [4,20]. Electrode placements were further verified online using functional criteria. Local electrical stimulation at low current in FEF produced eye movements [32], many neurons in the LIP region of interest (ROI) exhibited the classic response characteristic of peri-saccadic activity, and the CL ROI had a subset of neurons with high firing rates around 40-50Hz as previously reported in CL [28]. Finally, we performed post-mortem histology in one monkey, and visualized Hematoxylin and Eosin-stained sections (8 μm) under a microscope, to verify electrode tracks in ROIs.
We separated cortical ROIs into superficial and deep layers using inverse current source density (iCSD) analysis of auditory-evoked ERPs [33]. Both FEF [34][35][36] and LIP [37,38] have been shown to respond to auditory stimuli. The demarcation between superficial and deep layers was the bottom of the earliest sink after auditory stimulation (approximate boundary between layers 4 and 5). We then assigned electrode contacts to iCSD-defined superficial and deep layers, and confirmed layer assignments using reconstructions of recording sites along each electrode track (computed from electrode depth measurements and MRI of electrodes in situ) and spiking activity (demarcating gray from white matter). We used auditory-evoked ERPs to allow comparison of wake data with anesthesia data. CSDs were computed using the CSDplotter toolbox (https://github.com/espenhgn/CSDplotter) for MATLAB (dt = 1 ms, cortical conductivity value = 0.4 S/m, diameter = 0.5 mm). Tones were played in a sequence (200 ms duration; 800 ± 100 ms jitter between tones) comprised of 80% standard tones (0.9 kHz) and 20% deviant/oddball (1 kHz). Further details have been previously described [4].

Electrophysiological recordings and deep brain stimulation
We recorded broadband neural activity (filtered 0.1-7,500Hz, amplified and sampled at 40kHz) using a preamplifier with a high input impedance headstage and OmniPlex data acquisition system controlled by PlexControl software. Probes were MRI-compatible linear microelectrode arrays (MicroProbes) with 16 (FEF and CN) or 24 (LIP and CL) contacts; all were platinum/iridium electrodes 12.5μm in diameter with 200μm spacing between contacts and impedance between 0.8 and 1.0 MO. We recorded EEG with titanium skull screws located above dorsal fronto-parietal cortex.
We performed DBS using electrode arrays (24-contacts, 200μm spacing) previously used as recording electrodes, which have reduced impedance. Stimulations (400μs biphasic pulses, negative phase first) occurred simultaneously across the deepest 16 electrode contacts for 60s at 10, 50 or 200Hz with 200μA of current. These parameters had been previously titrated and successfully used to influence consciousness [4,20]. We targeted electrodes to maximize the number of contacts across the dorsal-to-ventral extent of CL (Fig 1), which we have previously shown to optimize stimulation effects [4].

Experimental paradigm
There were 40 total 2-4 hour recording days (18 for monkey R; 22 for monkey W). Monkeys sat upright in a primate chair in a quiet dark room. The animal's head was immobilized using a head post and/or four rods that slid into fabricated hollow slots in the acrylic implant. Each day, we inserted electrodes into ROIs (CL, CN, FEF, and LIP) through chronically implanted guide tubes [4].
We exposed animals to a paradigm that included a series of recording days in which animals received thalamic DBS (Stim) and separate series of recording days in which animals did not receive DBS (No-Stim). The first two days after switching (as well as weekends, when recordings did not take place) between these series types were considered transitional, and all other days in the same series were considered established, as this allowed sufficient data for analysis. Fig 2 depicts the schematized paradigm based on Monkey W. The paradigm was similar for monkey R, but started in the No-Stim condition following previously described anesthetized recordings and DBS experiments [4]. During a non-stimulation series, recordings included 4 brain areas (CL, CN, FEF, LIP). During a stimulation series, we applied DBS through the CL probe, precluding acquisition of CL recordings. On each recording day, monkeys alternated between task-on (rewarded fixation task, "FX") and task-off (resting wake, "RW") blocks. For a stimulation series, each pair of task on/off blocks included thalamic DBS pseudorandomly assigned to one of four conditions (10/50/200 Hz stimulation frequency and No-Stim control). After completing a run of all conditions (8 blocks, i.e., task on/off at 10Hz, 50Hz, 200Hz and No-Stim control), we started another run. Each block of about seven minutes duration included three stimulation replications: typically, stimulations occurred 1-2, 3-4 and 5-6 minutes after the start of the block, with one minute allowed between for recovery. To characterize the timing of observed effects during the 60s stimulations, we considered specific time windows relative to the start of DBS: onset (within 2s of the start of DBS), early (2-10s since start), mid (10-40s since start), and late (40-60s since start). Onset thus described the region of stimulation very close to initiation. In its entirety, this paradigm allowed us to investigate effects at the level of individual stimulations, blocks, recording days, and across the entire history of Stim and No-Stim series including transitions between them.

Behavioral monitoring and video processing
During recordings, experimenters closely observed animal behavior. We monitored eye movements using an eye-tracker and videos focused on the animal's face recorded evidence of facial (mouth and nose) movements and eye openness. Offline, videos were processed using MATLAB and converted into continuous luminance traces based on the average luminance in small windows around the animal's eyes, mouth, and nose. Eye traces revealed periods of open (low luminance) and closed (high luminance) eyes. Derivatives of luminance around the mouth and nose revealed periods when the mouth and nose were moving, or when movement was reduced.

Detecting changes in consciousness
DBS often triggered long-duration periods of behavioral arrest where animals stared vacantly into space and seemed to be unresponsive to their environment. We refer to these as vacant, perturbed consciousness (VPC) events. For systematic detection of VPC, we used custom code identifying periods of time when eye movement was decreased (based on the derivative of the polar coordinates of the eye tracker data), eyes remained open and not partially lidded (based on video eye traces), and facial movement was reduced and/or small and repetitive (based on video mouth and nose traces). Onset times were defined based on when the changes in behavior occurred (stillness and little/no eye movement) and offset was defined by when normal behavior resumed (eye/face/limb movement resumed). Due to technical issues, eye tracker data were not recorded on a subset of VPC events (n = 4), and video data were not recorded on a separate subset (n = 3). Detection was then based on the remaining signals, as well as behavioral notes. Due to the high degree of correlation between these metrics, this did not substantially affect detection. We defined pre-and post-event conditions as the entire time series up to 30s before or after event onset or offset, respectively. We identified possible events as those lasting for longer than 4.2s -the 95% confidence interval of fixation durations in the awake data set, and if the location of the animal's gaze did not match the location of the central fixation target in the fixation task (within 3˚of visual angle).
Many states of reduced consciousness including NREM sleep [39,40], general anesthesia [41,42], absence epilepsy [43][44][45][46], temporal lobe epilepsy [30,[47][48][49], and coma [50,51] coincide with low-frequency activity in the EEG. Similarly, online, we observed VPC events to coincide with low-frequency EEG activity. To distinguish VPC events from natural long-duration stares (the top 5 th percentile), which also showed reduced movement of the face and eyes, we searched offline for evidence of low-frequency activity. Events were assigned as VPC if they involved significant increases at any frequencies 1-9Hz in the EEG (in addition to the behavioral stillness and little/no eye movement) relative to the pre-VPC condition. This frequency range was used to best cover the full range of low-frequency activity observed across different states of reduced consciousness without biasing results towards the assumptions of any one clinical/atypical condition. Automatically detected events that lacked evidence of low-frequency activity were considered to be natural stares and used as the primary control in this study (see below for additional controls), as they shared similarities in behavior to VPC.

Controls
To identify neural correlates of VPC, we compared neural activity during VPC to a variety of controls accounting for activity differences which could be attributed to behavior, consciousness, or the typical dynamics of wakefulness and attention. To account for time-specific differences, we compared VPC to data recorded up to 30 seconds prior to (Pre) or after (Post) event onset/offset. To account for differences related to consciousness but not behavior per se, we compared VPC to resting wake data, when the animals were conscious but not performing a specific task. Resting data derived from epochs when eye position was stable (for at least 1s), and eyes were open (other face or body movements may have occurred during these windows). The resting state controls were also case-matched (data taken from same recording day with same task, DBS type, and electrode positioning) to account for the experimental paradigm (Stim/No-Stim and stimulation frequency) and electrode positions.
VPC was associated with staring. To account for possible goal-directed fixations, VPC was compared to data from a fixation task, where animals received juice rewards the longer they fixated on a faint, central target. Facial movements were common as animals licked juice rewards. Again, controls were case-matched and we only analyzed data during which the eyes were open and the animal was fixating. The fixation task also allowed comparisons between neural mechanisms present during resting wakefulness and cognitive engagement, which both reflect consciousness. Finally, to account for differences in consciousness while controlling for behavior, we compared VPC to long-duration stares without atypical, low-frequency EEG activity (stare, STR). While these control stares could not be case-matched, as they were rare and occurred randomly in the data set, they could be temporally aligned, similarly to VPC, with a defined onset, offset, pre and post period. This allowed us to compare the temporal dynamics of VPC directly to these specific controls.

Local field potential (LFP) and spike preprocessing
For VPC and stare events, we extracted data up to 30 seconds before behavioral onset (Pre), and up to 30 seconds after offset (Post) when possible. For resting wake controls and the fixation task, we first determined stable eye epochs, to match eye behavior between different conscious states (wakefulness, VPC, general anesthesia). These epochs started 200ms after a saccade and ended 200ms before the next saccade. For all events and controls, we divided epochs of data into non-overlapping 1s windows for analysis, which can be considered similar to trials. For LFPs, we bandpass filtered data from 1-250Hz (Butterworth, order 6, zero-phase filter) and downsampled to 1000Hz sampling frequency. We linearly detrended LFPs, then removed powerline noise (significant sine waves at 60Hz) using the function rmlinesc from the Chronux data analysis toolbox for MATLAB (http://chonux.org/). When individual electrode contacts had signal amplitude greater than 5 standard deviations from the mean of other contacts on the same probe, they were excluded from the analysis. For power spectral density and coherence analyses, we calculated bipolar derivations of the LFPs (the difference between adjacent electrode contacts) to curtail possible effects of a common reference and volume conduction [52][53][54]. Prior to entropy (H) and F � analyses, we performed ICA (Independent Component Analysis) to minimize the linear dependencies among LFP contacts in each area, then we normalized each LFP to its mean across all epochs for that recording day. We then binarized each LFP with respect to its median amplitude over the 1s epoch, to remove potential biases related to amplitude differences across channels or conditions.
For spiking activity, we bandpass filtered data from 250-5,000 Hz (Butterworth, order 4, zero-phase filter) and sorted spikes using Plexon Offline Sorter software. Spikes were initially detected with thresholds greater than 3 standard deviations from the mean. We used principal components analysis to extract spike shape features and the t-distribution expectation maximization estimation algorithm to determine clusters of spikes with similar features.
During DBS, stimulations induced a brief artifact caused by the applied current. To remove this artifact, we excised a 1ms window around the artifact, linearly interpolated across this window, and removed any significant sine waves at the stimulation frequency using the Chronux function rmlinesc.

Spike rate
We calculated the average spike count across the same 1s windows as the LFP analyses. For epochs that included thalamic DBS, artifact removal decreased the viable window as a function of the stimulation frequency. Thus, the window size was adjusted accordingly and accounted for in the spike rate calculation. For each neuron, spike rate was normalized to data collected during the pre-event condition, to allow comparison with all possible controls. We only analyzed neurons if they were recorded during a VPC or control stare event (in addition to other controls).

Power and coherence spectrograms
For VPC and control stare events, we calculated power spectrograms using 1s sliding windows with 0.1s stride/overlap across the entire time course of the event aligned to the onset of the behavior. Spectrograms were calculated for every bipolar-derived LFP using multi-taper methods (5 Slepian taper functions, time bandwidth product of 3) with the Chronux function mtspecgramc and then averaged across all viable electrode contacts.
To measure the temporal relationship between LFPs, we used coherence within and between CL, CN, FEF and LIP. Coherence is given by: where S(f) is the estimated power spectral density with subscripts 1 and 2 referring to the simultaneously recorded LFPs at two different sites. We applied multi-taper methods (5 Slepian taper functions, time bandwidth product of 3) to yield coherence spectrograms, using the Chronux function cohgramc. For VPC and control stare events, we calculated coherence in 1s sliding windows, stepped 0.1s, for all bipolar-derived LFP pairs corresponding to anatomically plausible pathways. We then averaged coherence values for all LFP pairs representing each pathway.

F � Computation
Integrated Information (F) is a proposed index of consciousness measuring both the complexity and integrated structure of a given system [6,55]. While measuring F directly is intractable in larger systems [55], surrogate measures such as F � , based on mutual information, have been shown to adequately approximate integrated information and detect changes in consciousness [20,[55][56][57][58].
We calculated F � [20,[55][56][57] as a practicable estimate of integrated information (F). To compute F � based on LFPs, we formed the state of a subsystem (t) at time t (1 ms time bins considering 1 kHz sampling frequency) as: where its elements are the bipolar-derived LFP signals, ranging from 1 to 23 signals (resulting from 24 electrode contacts) in each brain area: CN, CL, FEF (F) and LIP (L). The component of X(t) for each area is N ch × T-sized, where N ch specifies the number of bipolar-derived channels for each of four areas and T = 1000 (1s windows, sampled at 1kHz). We next calculated the uncertainty about the states assuming a multivariate Gaussian distribution of the states across time as: where S(X) (t removed to be concise) is the covariance matrix of X estimated over a 1s data window, μ(X) is the mean of the state vector X over the window, |S(X)| is the determinant of the covariance matrix that can be considered a measure of uncertainty about the state X at any time point within the window, and N is the total number of channels across areas. We used the Henze-Zirkler Multivariate Normality test [59] to verify that our states had multivariate Gaussian distributions.
The entropy for the states X(t), given its probability density function (pdf) as p(X), is defined as: The entropy is maximized if (X) is multivariate Gaussian, and can be calculated in closed form as: which can be described as the uncertainty about state X(t) at time t. The mutual information, i.e., the reduction of uncertainty about state (t) at time t, given its past at time t − τ can be calculated as: Where H(X t |X t−τ ) is the conditional entropy of state (t) given its past (t − τ), and can be derived in canonical form with a Gaussian distribution assumption of states as: where S(X t |X t−τ ) is the covariance matrix of the conditional distribution p(X t |X t−τ ) (conditional covariance) that can be expressed analytically as: where S(X t |X t−τ ) is the cross-covariance matrix. Mutual information, I, can be regarded as a measure of the information the current state has about its past, and it is used to calculate F � , a measure of integrated information. F � of the subsystem (t) is information that cannot be partitioned into independent parts of X and can be defined as: Where I � (disconnected I) is called mismatched information. We calculated I � for every bipartition of system X. System parts of interest were defined CN, F s , F m , F d , L s , L m , L d ; where s, m and d subscripts correspond to superficial, middle and deep layers respectively. Within each cortical layer (F s , F d , etc.), there were multiple bipolar-derived channels. We did not partition these individual layers, and so each layer was effectively a subsystem part consisting of multiple parallel channels.
The partition P that minimizes the normalized F � is the minimum information partition (MIP), as defined in [60]: Here m is the number of partitions and M k is the k th part of subsystem X. N P counterbalances inevitable asymmetries introduced by computing F � across a number of partitions of unequal sizes. The MIP is the weakest link between the parts of X, where its partition into subsystems results in minimal information loss. We calculated the covariance and cross-covariance matrices for each 1s window using the shrinkage method for a more stable result and averaged them across all windows for each recording day, to calculate the MIP and its corresponding F � , using gradient ascent and L-BFGS optimization method. We incorporated 2 NVIDIA GTX 1080ti GPUs to accelerate the process of searching for the MIP to calculate F � . We calculated I and F � with timelag τ = 15 ms, as this lag maximized F � in our previous study of this same CST system [20].
F � was computed using all available epochs for each condition (pre, post, VPC, control stare, resting wake, and fixation). Because F � is a relative metric that cannot be easily compared between systems with different composition (i.e., different brain areas), we excluded recordings where a given brain region was absent (for example, when no electrode contacts were identified in the deep cortical layers). To account for small differences in F � related to electrode placement, F � was centered relative to the values obtained during recording tasks in the paired fixation experiments. We also computed F � across time for VPC and control stare conditions, aligned to the onset of the event. In this case, F � was normalized (Z score) across all available time and then averaged for VPC and control stare events. For both centered and Z-score values, higher values of F � indicate increased consciousness, while lower values indicate decreased consciousness.
There are other available approximations of integrated information, such as F G . To ensure that our estimates of F � were reliable estimates of integration, we also calculated our effects with F G [61,62]. As we have demonstrated in another paper [20], the effects computed with F � and F G were qualitatively similar. Further, as an additional control, we computed F � after extrapolating the finite sequences into an infinite time series [62][63][64] to ensure that the finite approximations were not contributing to inaccurate estimates of integration. The effects were qualitatively and statistically similar. This validates that our F � calculations based on 1s epochs were robust.

Machine learning
We computed the power spectrogram of bipolar derived LFPs in 1s epochs (0.5s sliding window) from 11 seconds prior to VPC or control stare onset up to 12 seconds after offset (5 Slepian tapers). Data were only further analyzed between -10.5 to 4.5 seconds (relative to event onset), as this represented the minimum consistent aligned segments across all events. Power values were averaged for each frequency band (delta 1-4 Hz; theta, 4-8 Hz, alpha, 8-15 Hz; beta, 15-30Hz; gamma 30-90Hz) and all contacts in each area of interest (CN, FEFs, FEFd, LIPs, LIPd) which resulted in 25 input features for each 1s time bin.
We trained a Bayesian classifier [65] for each 1s time bin to classify the data into two states: VPC and Control Stare. For the Bayesian classifier, the class-conditional densities were modeled as gaussian distributions with equal covariances. We used an uninformative prior for classes (Normal-Jeffreys prior [66]) to avoid any bias toward any of the classes. Given X as input feature vector at time t, the posterior probability that the data belongs to class c can be derived as follows: Where θ(X t ) will be calculated at each time bin t as: Where μ c (t) is the expected value of class c at time t; and ∑(t) is the common covariance matrix for all classes, which is calculated as the unweighted mean of two covariance matrices for two classes and subsequently regularized using shrinkage where the regularization parameter was derived during training and cross-validation. All the parameters in equations above were calculated using the training set at each fold of cross-validation. We used 10-fold cross validation with 60% of the data as training set and 40% as test set. To test feature importance, we computed Mean Decrease in Accuracy (MDA). One at a time for each feature, we randomly permuted the sample labels (1000 times) and compared average accuracy to that of the unpermuted model. Permutation effectively erases the relationship between the given feature, condition, and other features, and thus the drop in accuracy can be interpreted as the relative importance of the permuted feature [67].

VPC probability and causal power
To compare the relative effects of different stimulation frequencies, we computed VPC and control stare probability by recording block for each condition. To compare DBS effects across longer time scales (different experimental days) we computed VPC and control stare probability in sliding windows (step size = 1 block) consisting of 32 recording blocks (the typical number of blocks recorded in two days). Stimulation predominance (the number of blocks that included thalamic DBS over the 32 recording blocks) was also calculated in each window. Both metrics were normalized to the maximum for correlation analysis.
We computed causal power (CP) [68], for probabilities that represented a theoretical causal structure between a dependent outcome (VPC event) and an independent experimental variable (stimulation frequency). CP is a differential probability ratio indicating the degree to which a proposed cause predicts an outcome more so than alternative causes, and is calculated based on the conditional probability matrix of event (e) co-occurring with/out cause (c), as seen in Eq 15. Positive/negative CP indicates that the probabilistic relationship is causal, and the proposed cause is more/less likely to lead to an event than alternatives.

Statistical analyses
We performed statistical analyses using correlations and linear mixed effect models in R, the PyMC3 [69] package in Python, and MATLAB. We controlled all p-values (p) using the Holm-Bonferroni method, and report them as p c when applicable. Models are shown for all interaction tests. VPC probabilities (with respect to stimulation onset and frequency), causal power, and MIP association probabilities were compared to their generated null-distribution (permutation test, see Fig 3 for an example). We generated distributions (Fig 3A) to represent a given population in our study containing N elements with different member components (for example, N A , N B , N C and N D could represent the total number of experiments with each stimulation condition (10/50/200Hz and NS control), and n observations (n A , n B , n C and n D could represent the total number of VPC events observed with each condition). From that distribution, we drew n simulated random samples (Fig 3B), without replacement if the population was finite and unrepeatable, or with replacement if otherwise, from the population, where n = the total number of observations. The process was repeated Nboot times, yielding a distribution of plausible observations under the null hypothesis of the given metric for each member category. Bootstrapping was always performed with replacement. Actual observations were then tested against these distributions to yield Z statistics and p-values (Fig 3C), which were corrected for multiple comparisons. Comparisons of probability ( Fig 3D) and sample proportion (Fig 3E), for example, can then be represented across member conditions relative to the null hypothesis, with error bars derived from the standard deviation of null distributions. Note, with this approach, increasing Nboot can increase the goodness of fit between the discrete distributions ( Fig 3C, bar histograms) and the fitted distribution ( Fig 3C, colored curves). However, the shape of the distribution, and thus the estimated error depends on the dimensions of the source population (N) and the sample rate (n). Because those values match what was done in the study, the resulting error bars reflect a reasonable estimate of the true, expected variance under the null hypothesis and cannot be inflated by Nboot.
To compare MIPs across states, we used a kernel density approach to estimate the probability mass functions for each MIP distribution for each state. We then applied a two-sample Kolmogorov-Smirnov test using the Pymc3 package in Python to determine whether these probability mass functions were different. We also compared the probability of different brain areas associating on the same side of the MIP for each state. These probabilities were then compared against the null distribution with 10,000 bootstraps.
We compared EEG power, LFP power, and LFP coherence time-frequency spectrograms between VPC and control stare conditions using the permutest function in MATLAB [70], based on the Maris-Oostenveld method [71]. Tsum statistics, electrode counts, p-values, and bootstrap parameters are reported for significant clusters. We also performed t-tests with multiple comparisons correction at each frequency and time-point, yielding similar results.
To further compare EEG power spectral density between pre, during and post-events, for VPC and control stare conditions, we used paired t-tests at each frequency in MATLAB; p-values were corrected across frequencies. Similarly, we compared F � values across time, with paired t-tests between the bin with maximum F � , and all other time-points; p-values were corrected across time. We evaluated machine learning accuracy using a one-sample t-test, comparing results to chance level performance (50%) at each sample across time; p-values were corrected across time.
We compared other effects using linear models (LM in R), reporting t statistics and corrected p-values for all relevant main and interaction effects, or correlations (cor in R), reporting correlation coefficients and p-values. For behavioral changes, we compared data across condition (During versus Pre/Post) for VPC events. For F � controls, we compared data between task condition (resting wakefulness versus fixation task) for all non-stimulation PLOS COMPUTATIONAL BIOLOGY samples. Statistical analysis was then performed separately for each stimulation condition (10/ 50/200Hz stim-on compared to stim-off). We normalized spike rates to the pre-event condition, then tested with linear models to produce an intercept and p-value. P-values were corrected for all brain regions in the same condition.
For H and F � , we centered values relative to those obtained during fixation controls to account for day-to-day variations. Condition (C, During vs Pre/Post) and event (E, VPC vs control stare) were coded as centered, dichotomous variables. We included animal (A) as a covariate, though similar results were obtained without its inclusion.
For long-duration stimulation effects, we tested the relationship between local event probability (LEP) and local stimulation predominance (LSP), first using the cor function in R, followed by linear modeling to produce best-fit lines. Both variables were normalized to their maximum, thus ranging between 0 and 1. Comparisons between VPC and control stare conditions used the following interaction model, where event (E) was a centered, dichotomous variable (VPC versus control stare). Animal (A) was included as a covariate, though results were similar either way.
This effect could be simplified by comparing the event probability according to stimulation series type (SS, Stim versus No-Stim). Event (E) was a centered, dichotomous variable (VPC versus control stare); and animal (A) was included as a covariate, though results were similar either way.
To account for differences seen in the transition between stimulation and non-stimulation series, relative to the established phase, we updated the model to include a 3-way interaction between SS, E, and establishment (EST, Transition versus Established), as a centered, dichotomous variable. Simple effects were consistent between the two models.
We compared MDA results along the dimensions of frequency and brain area using a dummy coding approach. Patterns were compared between condition (C, during versus pre) with the following models. P-values were corrected across all frequency and area comparisons.

Results
Here we first report behavioral and neural measures of perturbed consciousness resulting from thalamic DBS. We call this state vacant perturbed consciousness (VPC). Second, we present evidence that the probability of VPC is modulated by the short-and long-term experimental history of DBS in a frequency-dependent manner (VPC probability higher for 200 and 10Hz cf. 50Hz DBS). Finally, we compare the neural signature of VPC with that in anesthetized states, as well as contrast this signature with that in typical wakeful states, to converge on neural correlates of consciousness.

Thalamic stimulation during wakefulness reduces behavioral and EEG signatures of consciousness
We simultaneously stimulated across the dorsal-ventral extent of CL in awake monkeys via 16 contacts of a linear electrode array (Fig 1). On a subset of recording blocks (17.75% of all nonfixation task blocks, including those with and without stimulation; stimulation frequency-and history-dependence detailed below), animals went from normal behaviors, including scanning the room and making movements with their face and limbs, to suddenly holding their gaze and staring vacantly (Fig 4A). These events occurred without any overt sensory cause, such as a noise that might draw the animal's attention, or visual cue. Animals became still except for occasional small eye flutters and lip movements commonly seen in certain kinds of epilepsies [43][44][45]72], and would hold their gaze for an abnormally long period of time (> 4.2s, M = 17.64s, SD = 5.96s, 99.93 percentile of all stable fixations which lasted at least 1s). Electrophysiological examination of these events also revealed transient increases in lower frequency (1-9Hz) activity, especially in the EEG. These slower rhythms are commonly associated with reduced consciousness [30,50,73], and below we show corroborating evidence of this.

VPC corresponds to selective reductions in neural complexity and integration
Measures of complexity are promising neural indicators of consciousness for both animal models and human subjects [20,57,74,75]. We computed entropy (H), a measure of the complexity of neural interactions and information-carrying capacity, and F � , which additionally measures integration of network parts, i.e., how partitioning the network influences the generation of information. To assess the relationship between VPC and consciousness, we compared H and F � before, during, and after all VPC events and control stares, as well as during resting wake epochs, the fixation task, and two types of general anesthesia (isoflurane and propofol). We have previously shown that both F � and H tend to be significantly reduced during anesthesia or sleep compared to wakefulness [20]. Relative to the levels recorded during the fixation task, both H ( Fig 5A) and F � (Fig 5B) were decreased during VPC, but not control stares, compared to pre and post conditions. The interaction of condition (During vs Pre/ Post) and event (VPC vs control stare) was significant for both H (F(295) = 14.578, n = 300, p = 1.64x10 -4 , ANCOVA) and F � (F(295) = 34.261, n = 300, p = 1.28x10 -8 , ANCOVA). Differences in F � cannot be explained by general effects of thalamic DBS. When excluding VPC events, there were no significant differences in F � between stim-on (S+) or stim-off (S-)  (Fig 5C and 5D) at any stimulation frequency (t 50Hz (62) = .29, n = 64, p = .773; t 10Hz (88) = .65, n = 90, p = .517; t 200Hz (114) = .732, n = 116, p = .466, paired t-tests). F � also proved statistically similar between resting wake and on-task fixation conditions (Fig 5C and  5D, t(118) = 0.925, n = 120, p = .357, paired t-test), meaning decreases in F � were selective to VPC. These results support the conclusion that decreases in F � during VPC events indicate an event-specific reduction in information capacity, complexity and integration associated with consciousness.

Changes in consciousness correspond to changes in cortico-striatal associations
Calculating F � requires computing the minimum information partition (MIP), the cut in the system resulting in the minimum amount of information loss. This indicates system parts which are most integrated. We compared the distribution of MIPs for the largest consistent subsystem (containing CN, layers of frontal and/or parietal cortex) across different states of consciousness (isoflurane/propofol anesthesia, resting wake, fixation, pre/post/during VPC and control stares) using kernel density estimates and two-sample Kolmogorov-Smirnov tests (Fig 5E, 5F and 5G). The results revealed that VPC events exhibit a substantially different MIP distribution relative to pre (KS (28)  To further interrogate this effect, we used the MIP distributions to determine the probability of different areas associating on the same side of the MIP as a measure of connection strength (Fig 6A and 6B). While cortical connections remained relatively consistent across conditions, CN commonly associated with parietal regions in pre/post-VPC, STR, resting wakefulness and fixation task (Fig 6A, red dashed border), but associated with frontal regions in VPC, isoflurane and propofol (Fig 6B, blue dashed border). We have previously shown a similar switch in configuration to reflect changes in consciousness [20]. These configuration differences were substantial relative to those expected by chance. Ld-CN connection differences were especially strong and significant (|Z| � 2.15, p c � .0312, Nboot = 10000, permutation test) compared to a null distribution (Fig 6C). To further link this effect to consciousness, we simulated a null distribution representing consciousness by randomly selecting combinations of MIPs found in the wake or fixation task conditions. For VPC, isoflurane and propofol, of each brain area associating with any other on the same side of the MIP for states assumed to be (A) conscious (pre-VPC, post-VPC, resting wake (RW), fixation task (FX) and control stares) or (B) less conscious (VPC, propofol (Prop), or isoflurane (Iso)). C, D, Results (Z statistics ± SD of the null distribution) of the permutation test comparing cortico-striatal MIP associations under the null hypotheses that (C) MIP types occur randomly, or (D) MIP types reflect the same patterns as conscious states (defined by the resting wake and fixation task samples). Both approaches separate more conscious from less conscious states. E, Normalized F � (±SE) for VPC and control stares (STR) calculated in 1s sliding windows (.1s steps) and aligned to event onset (black vertical line). Horizontal pink line shows regions of significance for the pairwise t-test across time comparing VPC F � at each sample to the maximum. No significant differences were found in the control stare condition. F, Cortico-striatal significance (Z) based on the sliding window approach in E for VPC and stare conditions. Computations used the same approach as in C, but now for each sample in the sliding window. the associations of CN with cortical areas Ld, Ls, Fd, and Fs were all substantially different (|Z| � 6.62, p � 3.62x10 -11 , Nboot = 10000, permutation test) relative to the expected distribution under consciousness (Fig 6D). These results imply that network-level activity during VPC events shifts into a configuration associated with a reduction in consciousness.
We computed F � in 1s sliding windows (.1s steps) leading up to, and beyond, the behavioral onset of VPC and control stare events (Fig 6E and 6F). We analyzed data 10.5 seconds prior to event onset, and up to 4.5 seconds after, as these limits represented the maximum consistent duration across all events. On average, F � began to decrease about 3.0s before VPC onset (Fig 6E), when the MIP configuration switched from what is commonly seen in conscious states to the MIP most favored under anesthesia (Fig 6F). This change in MIP configuration continued until 2.9s after VPC onset when it switched back to the wake-favoring MIP. Following the MIP switch, F � was significantly lower than the maximum value between -1.6 and 3.9s relative to event onset (|t(82)| � 3.72, n = 84, p c � 0.049, paired t-tests) for VPC, but was stable during control stares and never significantly lower than the maximum (Fig 6E). Similarly, the MIP configuration during control stares remained in the wake-favoring configuration throughout and never changed relative to event onset ( Fig 6F). These results imply that the decreases in consciousness associated with VPC predate the onset of abnormal behavior by 1-3s.

VPC, not control stares, are modulated by thalamic DBS in a frequencydependent manner across acute and long-term time scales
We argue below that thalamic DBS (specifically at 10 and 200Hz) induces atypical network dynamics biasing towards VPC. That is, not all VPC events occurred as the direct result of an individual DBS stimulation, but rather, VPC in general is modulated by the entire experimental history of DBS in a frequency-dependent manner. We performed 338 experiment blocks in the resting wake condition. VPC events were significantly more common during recording blocks that included thalamic DBS (23.13%, n = 37) than those without (12.92%, n = 23, Z = -2.45, p = .014). Each block involved periods of no stimulation with three replicated stimulation epochs (Fig 2). VPC was significantly more likely to occur within the bounds of the stimulations (64.86%, n = 24) than outside (35.14%, n = 13) compared to a null distribution representing random occurrence (Z = 2.36, p c = .036, n = 37, permutation test). Moreover, VPC events were overrepresented near the onset (Fig 7A) of a stimulation epoch (within 2s of stimulation onset, 25%, n = 6; Z = -2.58, p = .031, Nboot = 100,000, permutation test). These results imply that thalamic DBS serves as a trigger for VPC under certain conditions.
In line with our hypotheses, VPC probability was higher with 10Hz (23.9%) and 200Hz (30.0%) stimulations, than it was for 50Hz (14.8%) and no-stimulation controls (12.9%) ( Fig  7B). Permutation tests revealed that 200Hz stimulations gave rise to significantly increased probability of VPC (Z = 2.71, p c = .025, Nboot = 100,000, permutation test) compared to the average occurrence (17.8%), while experiments without stimulation (NS) had a significantly decreased probability (Z = -2.44, p c = .044, Nboot = 100,000, permutation test). We found similar results with causal power (CP, Fig 7D-7F), a probability ratio indicating the strength of causal relationship between two variables relative to alternatives (see METHODS). Positive or negative CP indicates that a given cause (stimulation condition) increases or decreases effect likelihood (VPC occurrence) more so than other known plausible causes [68]. Across the experimental history for both monkeys, CP was typically positive for 200Hz and 10Hz stimulations, and negative for no stimulation and 50Hz stimulation (Fig 7C and 7D). Simulated nulldistribution tests indicated causal power for VPC was significantly positive for 200 Hz stimulations (ΔCP = .175, Z = 2.66, p = .031, Nboot = 100,000, permutation test), and significantly negative for no stimulation (ΔCP = -.133, Z = -2.58, p = .031, Nboot = 100,000, permutation test) relative to expectations under a null distribution. These results imply that abnormally high frequency activity applied to CL promotes VPC with our DBS method.
The occurrence of VPC also appeared to depend on the stimulation order. In our study, we pseudorandomly assigned the order of stimulation experiments between no-stimulation control, 10Hz, 50Hz and 200Hz conditions (Fig 2) within each experimental run. Thus, for established days in a stimulation series, we could compare the probability of VPC occurring during blocks with a given stimulation frequency relative to that of the preceding stimulation ( Fig 7F). Interestingly, 50Hz stimulations, while the least likely to produce VPC, were more likely to do so if preceded by a 200Hz (CP = .091) rather than 10Hz stimulation (CP = -.147). 10Hz stimulations more commonly exhibited VPC after 200Hz (CP = .063) relative to 50Hz (CP = -.273) stimulations; and 200Hz stimulations more commonly exhibited VPC after 10Hz (CP = .515) relative to 50Hz (CP = .038) stimulations. Although none of these effects reached significance with permutation tests (due to the reduced sample sizes when considering all possible combinations of stimulation conditions), the synergistic effect of 10Hz on 200Hz was trending (Z = 1.63, p = .102, Nboot = 100,000, permutation test). Though the limited power of these analyses does not warrant strong claims, it is still interesting to note that abnormal stimulation frequencies relative to CL (10Hz and 200Hz) both tend to exacerbate VPC and have positive CP. In comparison, 50 Hz stimulations, which we have previously shown to increase consciousness in macaques [4], tend to have negative CP and decrease VPC occurrence.
Irrespective of stimulation frequency, VPC probability increased with the proportion of DBS blocks over consecutive days (Fig 7G). When computing the probability of VPC in sliding windows across 32 blocks (approximating 4 experimental runs and 2 recording days), the local probability of VPC was positively correlated (Fig 7G) with the predominance of stimulations within the same window (r = .5105, t(309) = 10.40, p = 6.54x10 -16 ). Control stares were not significantly correlated with stimulations (r = -.0466, t(309) = -.82, p = .4141). An interaction model verified that the effect of stimulation predominance on event probability was substantially different between VPC and control stares (F(613) = 68.043, n = 618, p = 9.79x10 -16 , ANCOVA). Our experimental paradigm consisted of shifting between series of recording days that either included stimulation or did not, and each series typically lasted for at least one week, with the first two days of a series being considered transitional (Fig 2). VPC on no-stimulation experimental runs was less common (38.33%, n = 23) than those with stimulation (61.67%, n = 37; Z = -2.45, p = .014, Nboot = 100,000, permutation test). Further, local event probability was significantly higher for experimental runs that included thalamic DBS (Fig 7H, "All" bars) compared to those that did not (t(613) = 6.639, p = 6.96x10 -11 , n = 618, t-test). This effect was substantially larger than the same effect computed for stare controls, yielding a significant interaction (F(618) = 18.787, n = 618, p = 1.71x10 -5 , ANCOVA). Interestingly, local VPC probability was higher in the transition period for no-stimulation series, i.e., when the previous few recording days had included thalamic DBS, and decreased markedly after the series was established (Fig 7H). In contrast, local VPC probability was reduced when transitioning from a no-stimulation into a stimulation series, and increased after the series was established. A 3-way ANCOVA verified these effects were significant and specific to VPC events rather than stare controls; the three-way interaction of experimental series type (Stim, No-Stim), establishment (Transition, Established), and event (VPC vs Control Stare) was significant (F(609) = 92.294, n = 618, p = 2.00x10 -16 , ANCOVA). These results imply that thalamic DBS influences the prevalence of VPC over longer periods of time. Thus, while DBS can directly initiate VPC events, it can also produce longer-duration vulnerabilities that gradually dissipate after stimulation experiments end.
We also included a fixation task in this study as a control (331 experimental blocks recorded), where DBS did not have the same effect. VPC only occurred on 2.72% of fixation task blocks compared to 17.75% of resting wake blocks; this difference was significant (Z = 6.389, p = 1.67x10 -10 , Nboot = 100,000, permutation test).

VPC events correspond to reduced neural firing
We compared the firing rate of the same neurons, normalized to the pre-condition, across VPC, control stare, post, resting wake and fixation task epochs (although this reduced the relatively small sample size further, it better enabled comparisons with all controls). Firing rates were very similar to the pre-condition for all conditions except VPC (Fig 8). While all areas had lower firing on average during VPC, this difference was only significant for Fs (t (22)

Reduced consciousness in VPC is associated with increases in lowfrequency power spectral density
To further evaluate neural signatures of decreased consciousness, we compared LFP spectrograms in different brain regions between all available VPC and control stare conditions relative to event onset. In general, power tended to be lower 3-5s prior to VPC onset (Fig 9A-9F). This was strongest in CL (Fig 9F) which had a significant cluster from θ-γ H frequencies from approximately -5 to -1.8s relative to event onset (Tsum = -650.76, p = 9.99x10 -5 , Nboot = 10,000). Power tended to increase in most areas at low frequencies just prior to event onset, and remain elevated at lower frequencies, especially in Ls ( Fig 9B) and Ld (Fig 9D), as the event continued. Power in CL tended to be much higher -1 to 4s around onset, especially at θ and α frequencies. Fig 9F shows the significant cluster (Tsum = 282.74, p = 0.0013, Nboot = 10,000). See Table 1 for full statistical results.
To further investigate the changes in power spectral density influencing VPC, we trained a Bayesian classifier to differentiate between VPC events and control stares, using the LFP power in different frequency bands (δ, θ, α, β, and γ L ) from CN, Ld, Ls, Fd, and Fs calculated in 1s sliding windows (0.5s stride) relative to event onset. CL was excluded because it was only available on a subset of samples, hindering convergence of the parameter estimators. Decoding was significant (Fig 10A) as early as 4s prior to event onset and remained significant thereafter (|t -4 to 4.5 (9)|�4.66, n = 10, p c �.0272). This window was very similar to the time course of changes in F � (Fig 6E and 6F). We used mean decrease in accuracy (MDA) analysis to estimate which of the LFP features (5 areas x 5 frequency bands) contributed most to decoding accuracy across time (Fig 10B). The MDA is computed by comparing accuracy when the machine learning model includes all possible features to an impoverished model where the information content of a feature has been negated by permuting the sample orders. The greater the MDA, the greater the contribution of that feature. Prior to event onset (-2.5 to -0.5s, Fig 10A, blue horizontal bar), exclusion of CN β , Ld θ,α,β , or Ls θ,α resulted in the largest decreases in accuracy. Overall, features from Ld (Fig 10C, |t(1998  Between .5 to 4.5s following VPC onset (Fig 10A, thick red horizontal bar), feature importance shifted to lower frequencies. Feature importance was significantly higher for δ, θ, and α frequencies, and lower for β and γ frequencies after the shift (|t(1998)| � 127.3, n = 2000, p c � 1.0x10 -15 ). The importance of CN significantly increased, and importance decreased for all other areas (|t(1998)| � 78.75, n = 2000, p c � 1.0x10 -15 ). Overall, features from CN and Ld ( Fig   Fig 9. VPC is associated with increases in low-frequency power relative to control stares. A-F Spectrograms of power differences (VPC-control stare) for each frequency band across time aligned to event onset (white vertical line). Results include all recorded VPC and control stares for (A) superficial FEF, (B) superficial LIP, (C) deep FEF, (D) deep LIP, (E) CN, and (F) only the subset of events without thalamic DBS for CL. Significant clusters are outlined in black with numerical labels referenced in Table 1 for statistical details.
These results indicate that VPC is associated with a substantial increase in low frequency activity between 1-15 Hz. As this activity is most substantial and informative when found in subcortical and parietal structures, these regions may play the largest role in inducing and promoting VPC.

Reduced consciousness in VPC is associated with changes in CST functional connectivity
To evaluate functional connectivity signatures associated with reduced consciousness, we compared LFP-LFP coherence between the different brain regions representing anatomicallymotivated pathways (Fig 11A) relative to event onset between VPC and STR conditions. Changes in coherence were substantial and dynamic across time, leading to a large number of significant clusters (Table 2) in time-frequency coherograms. Thalamo-cortical coherence ( Fig  11B and 11C) increased transiently at δ about 3.5s before onset, followed by θ and α about 2.5s before onset (CL-Fs cluster 3, CL-Fd clusters 4 and 5, CL-Ls cluster 2, and, CL-Ld cluster 4, Tsum � 43.23, p � .0084, Nboot = 10,000). This θ and α coherence only lasted about 1.5s in thalamo-frontal pathways (Fig 11B), but continued in thalamo-parietal pathways (Fig 11C) at least 1.5s after event onset, matching the time course of abnormal θ dynamics exhibited by CL.
To summarize the complex temporal dynamics across all relevant pathways, we compared the percentage of significant time-frequency tiles (hereon referred to as pixels) within different temporal windows, indicating increases or decreases in high (β,γ L ,γ H ) or low (δ,θ,α) frequency coherence (Fig 12). This can be considered the predominance of different coherence patterns. Window P1 (Fig 12A, -5 to -2.6s) represents the period of time prior to event onset where decoding accuracy (based on power spectral density) is low and differences in consciousness (F � ) are insignificant. Significant pixels are rare in this period, indicating similarities in coherence between VPC and control events. Most significant pixels that are present reflect decreases in low-frequency coherence between CL and other regions. Window P2 (Fig 12B, -2.5 to -0.6s) represents the time period prior to event onset where decoding accuracy increases, and decreases in F � become significant. This corresponds to an increase in the number of significant pixels indicating increases in low-frequency coherence. These changes span most pathways, but especially those involving CN and deep layers of cortical areas, where the effects are more broadband. At the same time, pixels representing the relationship between CL and other areas reflect high-frequency decreases in coherence, whereas cortico-cortical pathways have increased coherence at higher frequencies predominating.
Window EO (Fig 12C, -.5 to .5s) represents the second of data around event onset, where decoding has reached maximum accuracy and changes in F � hit the minimum F � value. In this window, the majority of pixels representing the relationship between cortex and CN reflect broadband coherence increases. At the same time, there are broadband decreases in coherence predominating for pathways connecting CL to other areas. Window A1 (Fig 12D, .6 to 2.5s) represents the first two seconds after event onset, where decoding remains high and decreases in F � remain significant. In this window, there is increased low-frequency coherence within areas (Fig 12D, third row, along diagonal), but not between most areas. Finally, in window A2 (Fig 12E, 2.5 to 4.5s) when behavioral signs of VPC are still present, significant decoding persists and F � is gradually returning towards the baseline level, there is predominantly increased low-frequency coherence in parietal cortex, in addition to CL.

Discussion
We used thalamic DBS to induce acute and spontaneous perturbations in consciousness (VPC events). These events include long-duration, vacant stares coinciding with lower-frequency EEG activity and decreases in neural measures of complexity and integration indicative of consciousness. VPC was best predicted by power changes in parietal cortex and striatum relative to frontal cortex. Just before VPC events, consciousness, as indexed by F � , decreased as lowerfrequency coherence increased between cortical and subcortical regions (Fig 13A). There was a burst of broadband coherence between the cortex and striatum around VPC onset (Fig 13B), which was followed by dominant lower-frequency coherence within each cortical and subcortical area during VPC (Fig 13C), indicative of increased isolation of individual areas. VPC also involved reduced firing rate in superficial FEF, which could result from reduced input, reflected in the decreased broadband coherence of the feedforward pathway during VPC.
We also saw reduced firing rates in LIP deep layers during VPC reminiscent of the reduced excitability of deep-layer cortical neurons previously shown under anesthesia [4,13], which for pairs of LFPs within/between brain area(s). Color scales with the strength of the resulting t-stat. Significant clusters are outlined in black with numerical labels referenced in Table 2   could result from perturbed intracolumnar communication. Overall, these patterns reflect disruption of typical CST network dynamics, resulting in reduced behavior and consciousness.

The nature of VPC
The VPC state that we describe in this manuscript proved inconsistent with many naturalistic states observed in macaques. Despite EEG signatures of reduced consciousness and reductions in integration and complexity, VPC was not associated with drowsiness or general decreases in arousal. Animals remained upright with their eyes wide open. VPC was also behaviorally inconsistent with goal-directed attention, despite stillness and decreased eye movement. During vacant stares the eye position would typically freeze and then drift at a very slow rate, unlike either visual search or the goal-oriented fixations observed during task performance. During rest, it is possible to engage in internally-directed activity such as mind-wandering, typically associated with reduced attention to the external environment [78] and shifts in functional network connectivity towards default mode dynamics [79]. Reduced external attention may explain the vacant staring observed during VPC. However, internally-directed states are not expected to involve substantial reductions in the brain's overall capacity for carrying information nor the complexity of information flow [80], but rather a different allocation of neural resources [78]. We found that VPC involved significant reductions in both H and F � , reflecting reduced complexity of neural dynamics and information-carrying capacity. Such reductions in H and F � were not observed when comparing on-task behavior (fixation task) relative to resting wakefulness, despite the conditions being distinguished by different levels of attention and the sensory experience of reward. But reductions in H and F � were observed in the same recorded network and animals during sleep or general anesthesia [20]. Thus, even if VPC does not involve observable decreases in arousal, it does appear to involve neural processes more typical of a loss of consciousness.
If VPC does not resemble typical wakeful states, perhaps it better resembles abnormal conditions induced by chemical intervention or disease. For example, VPC could be similar to dissociative states induced by psychedelics, which do often alter activity in cortico-striatothalamic networks [21,81] in some ways similar to VPC. However, behavioral changes during psychedelic use [82] are typically different to those observed during VPC. Further, prevailing evidence and theories about the effect of psychedelics suggest that they will cause increases in measures of complexity, such as entropy [83][84][85]. Rather, an intriguing possibility is that VPC is related to epileptic disorders, such as temporal lobe or absence epilepsy [30,43,73]. Generally speaking, these conditions are associated with a loss of consciousness despite sustained behavioral arousal with reduced but repetitive motor movements, i.e., automatisms. For absence  epilepsy, seizures coincide with reduced neural complexity [86], just like VPC. Further, absence epilepsy, which in humans is associated with the sudden onset of vacant staring during an episode, is a particularly close match to the behavioral signature of VPC [43][44][45]. Moreover, just as VPC appeared to be initiated and kindled by thalamic DBS, similar to mechanisms known to influence seizure rates in epilepsies [87], so too has thalamic DBS triggered absence epilepsy in a number of past studies [26,88,89].

Shared mechanisms connect VPC to absence epilepsy
Intralaminar thalamic DBS was originally reported to induce decreases in consciousness such as sleep or behavioral arrest in cats [88]. VPC matches the description of the arrest reaction as described by Hunter and Jasper (1949), which included sudden decreases in behavior with reduced responsiveness to sensory stimuli. The arrest reaction sometimes progressed into a state resembling an absence seizure [88], complete with 3 Hz spike-and-wave oscillations common to absence epilepsy [43][44][45]. Subsequent thalamic DBS studies have been able to replicate these effects in both cats [26] and monkeys [89]. Similar in effect, our study used a unique thalamic DBS protocol to produce VPC, which shares many behavioral correlates with absence epilepsy, including the sudden onset of vacant staring and the small mouth and eye movements consistent with automatisms [43][44][45]. Further, our results reveal considerable overlap between the neural mechanisms of VPC and absence. Absence epilepsy has long been associated with altered thalamocortical activity [90], including substantial increases in thalamocortical coherence [91]. Our VPC events were triggered by thalamic DBS, and initiated with increased low-frequency coherence along thalamocortical pathways, especially involving parietal cortex. More recently, evidence has also pointed to the importance of cortico-striatal [23,24,46,92] dynamics in absence seizures and, in rodent models, dopaminergic manipulation in the striatum can even modulate seizures [93]. Similarly, we found lower-frequency coherence changes along cortico-striatal pathways, and ongoing decreases in striatal firing rate were also characteristic of VPC.
The chief distinction between VPC and absence epilepsy is the presence of 3Hz spike-wave patterns in the EEG during absences. We observed increased oscillatory activity during VPC in the absence range (3-8Hz), though it lacked the characteristic spike-and-wave patterns attributed to the disorder in humans [43][44][45] and animal models [43][44][45][46]. One possibility for this discrepancy is that, while VPC seems to share mechanisms with absence, which can be triggered by thalamic DBS [88,89,94], it represents a pre-epileptic version of the syndrome that does not include spike-and-wave pathology. Previous work in cats implies that progression from behavioral arrest to 3Hz spike-and-wave absence events is possible [88] but difficult to replicate and may require stimulation of the thalamic reticular nucleus [26]. Our stimulation protocol involved a multi-microelectrode approach with relatively low-current (200μA) stimulations, which were based on parameters we have previously used to manipulate consciousness under anesthesia [4]. As other studies have shown, higher-current stimulations may be more likely to progress effects like VPC to epileptic activity [88]. In monkeys, absence responses with thalamic DBS are more common in drowsy animals [89] and after premotor cortical lesions [94,95], while our animals were intact, and recordings took place during the typical wakeful period. All of these results imply that VPC may be a sub-epileptic condition, useful for studies of consciousness. A further implication is that intralaminar thalamic DBS would more specifically model absence epilepsy if stimulation parameters were tailored to that effect, for example, using bilateral stimulations at higher current.
The mechanistic links between VPC and absence epilepsy are informative for clinical populations. Foremost, the similarity in symptoms between VPC and absence seizures, even without spike-and-wave patterns, could imply that these 3Hz rhythms are correlated with, but not necessarily causal towards, loss of consciousness during absences. In line with this, while rodent genetic models of absence epilepsy show spike-and-wave pathology [46], spike-andwave discharges have also been reported in wild-type rodents [96]. In humans, absence seizures can present with abnormal or disorganized spike-wave morphology [97], and their overall clarity and severity can vary both within and between patients [43]. Additionally, the frequency of spike-and-wave oscillations can be variable over the course of a seizure, typically starting faster and slowing over time [89,98,99]. Further, oscillations in animals are typically higher in frequency than those reported in humans, ranging up to 6Hz in macaques [89,100], or from 5-10Hz in rodents [46]. Thus, while spike-and-wave pathology is likely important for absence epilepsy, it may not be directly causal towards mechanisms of loss of consciousness. Rather, they may result as a consequence of the same abnormal thalamocortical dynamics which seem to trigger VPC in this study, and absence events in clinical populations [43]. If so, loss of consciousness is likely more closely linked to perturbed functional connectivity in CST networks.

Mechanisms of multi-microelectrode DBS
This study used multi-microelectrode DBS to manipulate consciousness in a frequency-dependent manner. Instead of larger DBS electrodes commonly used clinically, our stimulating electrodes were linear microelectrode arrays designed to conform to CL dimensions in macaques (Fig 1). Previously, we used this DBS method to increase consciousness in anesthetized macaques in a frequency-dependent manner that was CL-selective and maximally effective at 50Hz [4]. This result contrasts with another recent study achieving similar arousing effects using more traditional methodology, but at a much higher stimulation frequency (150Hz) and current (1.0-2.0mA vs our 200μA) [5]. These discrepancies could imply a fundamental difference in mechanism between the DBS methods. It is possible that the multi-microelectrode method more directly influences the targeted nucleus in ways that match the stimulation parameters.
This interpretation is supported by our current study, where thalamic DBS modulated VPC at both acute and longer-duration time scales. VPC was most strongly driven by 200Hz or 10Hz perturbations, both of which represent abnormal dynamics for CL in wakefulness, relative to 50Hz stimulation, which produced modulatory effects no stronger than no-stimulation controls. Further, VPC was most common when 200Hz stimulations followed 10Hz blocks. 10Hz activity is more similar to the dynamics of CL during sleep [4,29], while 200Hz represents, over longer time scales, non-physiological activity, or possibly the short inter-spike intervals within bursts generated by many other thalamic neurons in reduced conscious states. In comparison, 50Hz matches the dynamics of fast-firing CL neurons during wakefulness [4,28,29]. Clinically, absence seizures are associated with reduced arousal and are more common after sleep deprivation [101], and thalamic DBS more readily triggers absence-like responses in drowsy [89] or lightly anesthetized animals [102]. Similarly for VPC, 10Hz stimulations might reduce arousal and increase low-frequency oscillatory activity which propagates throughout the CST network. As our results show, such oscillations within CST areas predominate during VPC (Fig 12C). In addition, low-frequency coherence in the CST network often predates the onset of VPC (Fig 12A and 12B). Subsequent 200Hz stimulations, which might either mimic thalamic bursting or even lead to a shutdown of CL, capitalize on this vulnerability and thus induce VPC more often than other conditions. The fact that our results are frequency-specific, and these frequencies are functionally relevant for CL, implies that our method exerts effects by mimicking outputs of CL, and not through antidromic mechanisms.

Optimizing paradigms to produce VPC
VPC was associated with stimulation onset, implying that DBS directly influences CST dynamics to promote VPC. Further, VPC became more prevalent as stimulations were performed more frequently, even occurring without DBS application. These results imply that vulnerability to CST circuit disruptions linger and contribute to future disruptions in consciousness, not dissimilar to kindling effects in seizure models [87]. Our stimulations lasted for one minute, and were replicated three times per stimulation block, with multiple blocks performed each day. This paradigm could induce neural plasticity [103,104] promoting abnormal network dynamics that both exacerbate future responses to DBS and increase the likelihood for restingstate dynamics to spontaneously give way to VPC events. Similar mechanisms seem to contribute to epilepsy, as seizures commonly alter functional connectivity over the course of the disease, with the effect size dependent on seizure rate [105]. Interestingly for VPC, plasticity seemed to depend on exposure to the stimulation paradigm, as after experiments ceased, no abnormal events have been observed by experimenters or care staff during everyday behavior.
These results suggest that experimental parameters could be optimized to increase the rate of VPC beyond those we observed. Others have employed bilateral stimulation [3,106], which may produce stronger effects than our unilateral approach. Further, our experiments incorporated factors which decreased VPC likelihood (fixation task, no-stimulation blocks, 50Hz stimulations). Their inclusion in this study serve as controls that help link VPC to related states of decreased consciousness. For example, studies have shown that seizure rates tend to decrease with cognitive engagement in animal models [107], and VPC is similarly rare during fixation tasks. However, including the fixation task, no-stimulation blocks and 50Hz stimulations likely limited the overall occurrence of VPC well below what is possible with our stimulation method. Our results imply that the stimulation effects could increase with repeated batteries of 10Hz and 200Hz stimulations, which trended towards a synergistic effect. Such optimizations increasing VPC likelihood could further enable perceptual experiments in which sensory stimuli are presented during VPC or without perturbing consciousness. Increasing the VPC rate from DBS is likely to also increase the spontaneous VPC rate during non-stimulation transitions, which would further reduce potential confounds, such as artifacts caused by DBS.

VPC and neural correlates of consciousness in CST networks
This study builds on recent findings about CST circuitry during different conscious states. We previously showed that F � , a measure of neural integration, differentiates wakefulness, sleep, and anesthesia [20]. Here, F � further differentiates resting wakefulness, fixation tasks and stare controls from VPC. The time course of changes in F � map onto the time course of our ability to decode VPC from control stare events based on cortical and subcortical LFP power, as well as dynamic changes in CST network coherence. In this study, atypical broadband increases in coherence between cortex and striatum were linked to decreased consciousness at VPC onset, just as ongoing low-frequency intra-area coherence was linked to reduced consciousness during VPC events. Critically, this implies that different patterns of aberrant activity can contribute to reduced F � and consciousness in CST circuits.
Decreased consciousness also temporally aligned with changes in the integrated structure of the CST network (Fig 6), such that frontal regions, rather than parietal, became more affiliated with subcortical structures. This finding during VPC generalizes network configuration changes for reduced consciousness under two separate anesthetic agents and, as previously shown, during non-rapid-eye-movement sleep [20]. These findings challenge cortico-centric assumptions about consciousness. Affiliation of parietal and subcortical structures during wakefulness may represent parieto-striatal-thalamic interactions promoting sensory processing and perceptual representations. Shifts of subcortical affiliations to frontal cortex could represent breakdowns in the aforesaid more posterior associative networks and sensory pathways, resulting in decreases in consciousness.
The question remains how disruption of thalamic activity with DBS spreads across the CST network to perturb consciousness. The thalamus could disrupt cortico-striatal interactions by either directly influencing cortex, and consequently cortico-striatal projections, or via thalamo-striatal projections that directly modulate striatal responses [108]. Our results suggest the former has a greater influence. Thalamic DBS caused interference in thalamic outputs which translated to cortico-striatal, more so than thalamo-striatal, perturbances (Figs 11, 12 and 13). Abnormal cortico-striatal inputs likely alter basal ganglia output to the thalamus, further perturbing thalamo-cortical, cortico-cortical and cortico-striatal processes. Overall, our results highlight the influential role of the basal ganglia and higher-order thalamic structures, which control information processing and flow through CST networks in ways that support and manipulate consciousness [18,109]. This challenges common assumptions about the relevance of subcortical brain areas in consciousness research [17,110] which is of great importance to clinical practice.