A hidden Markov model reliably characterizes ketamine-induced spectral dynamics in macaque local field potentials and human electroencephalograms

Ketamine is an NMDA receptor antagonist commonly used to maintain general anesthesia. At anesthetic doses, ketamine causes high power gamma (25-50 Hz) oscillations alternating with slow-delta (0.1-4 Hz) oscillations. These dynamics are readily observed in local field potentials (LFPs) of non-human primates (NHPs) and electroencephalogram (EEG) recordings from human subjects. However, a detailed statistical analysis of these dynamics has not been reported. We characterize ketamine’s neural dynamics using a hidden Markov model (HMM). The HMM observations are sequences of spectral power in seven canonical frequency bands between 0 to 50 Hz, where power is averaged within each band and scaled between 0 and 1. We model the observations as realizations of multivariate beta probability distributions that depend on a discrete-valued latent state process whose state transitions obey Markov dynamics. Using an expectation-maximization algorithm, we fit this beta-HMM to LFP recordings from 2 NHPs, and separately, to EEG recordings from 9 human subjects who received anesthetic doses of ketamine. Our beta-HMM framework provides a useful tool for experimental data analysis. Together, the estimated beta-HMM parameters and optimal state trajectory revealed an alternating pattern of states characterized primarily by gamma and slow-delta activities. The mean duration of the gamma activity was 2.2s([1.7,2.8]s) and 1.2s([0.9,1.5]s) for the two NHPs, and 2.5s([1.7,3.6]s) for the human subjects. The mean duration of the slow-delta activity was 1.6s([1.2,2.0]s) and 1.0s([0.8,1.2]s) for the two NHPs, and 1.8s([1.3,2.4]s) for the human subjects. Our characterizations of the alternating gamma slow-delta activities revealed five sub-states that show regular sequential transitions. These quantitative insights can inform the development of rhythm-generating neuronal circuit models that give mechanistic insights into this phenomenon and how ketamine produces altered states of arousal.

neuronal circuit models that give mechanistic insights into this phenomenon and how ketamine produces altered states of arousal. 1 Author summary Monitoring brain activity during anesthesia can provide insights into the underlying mechanisms of how anesthetics elicit altered states of consciousness. Ketamine, a commonly used anesthetic, is known to cause short duration bursts of high frequency electrophysiological activity in the brain, but the neural mechanisms underlying this activity are not well understood. A key limitation in developing accurate models of the underlying mechanism is a lack of detailed knowledge of the dynamic structure and spectral properties of ketamine-induced oscillations. In this work, we address this limitation by developing a statistical framework to quantify ketamine-induced neural activity. Our framework is based on a hidden Markov model, which assumes that the neural activity switches among discrete states, each of which has its own distinctive probabilistic spectral representation. By estimating this versatile statistical model from electrophysiology data, we generated detailed descriptions of the dynamic properties and oscillatory signatures associated with ketamine-induced neurophysiological states in non-human primates and in human patients. Furthermore, we identified additional ketamine-induced states that have not yet been reported. In summary, our detailed quantitative descriptions of ketamine-induced spectra can aid further developments of neurophysiological mechanistic models of ketamine as well as biomarker discovery for clinical anesthesia monitoring.

Introduction
Ketamine is a phenylcyclidine derivative and one of the most commonly used anesthetics in clinical use world-wide [1][2][3][4]. It is classified as a WHO Essential Medicine [5], and is often the only anesthetic available in clinics of developing countries [4]. At low doses, ketamine is known to create a state of "dissociative anesthesia" characterized by altered sensory perception and analgesia [1,2]. At high doses it also leads to loss of consciousness, and therefore it is often used as an anesthetic agent during surgery in humans and animals [3,[6][7][8]. Under anesthetic dosages, ketamine produces distinct oscillatory signatures in the electroencephalogram (EEG) of healthy volunteers and patients [9][10][11]. In a recent retrospective study by Akeju et al. [10], the authors found that the frontal EEG from patients who received ketamine for the induction of general anesthesia showed distinct alternating periods of intermittent high power activity in the 27-40 Hz and 0.1-4 Hz frequency bands. Spontaneous gamma oscillations (i.e.  have also been demonstrated under high dose ketamine in mice [12], rats [13], cats [14], nonhuman primates (NHPs) [15][16][17], and sheep [18]. Recent neuroscience research in NHPs under high-dose ketamine revealed prominent gamma oscillations modulated by slow wave (0.3 Hz) activity [17], similar to the alternating periods of activity in gamma and slow-delta bands observed in humans [10]. Objective characterization of the intermittent band-limited spectral dynamics induced by ketamine would be useful both for monitoring patients during ketamine-induced unconsciousness in clinical settings, as well as for aiding computational and experimental neuroscience research aimed at building a mechanistic understanding of the phenomena. Therefore, the availability of an efficient analytic tool, which can objectively Hz). The power is averaged within each band and scaled between 0 and 1 to facilitate comparison among data sets that vary in signal amplitude but not in the spectral profiles they represent. To describe this multivariate observation sequence, we assume that the spectral dynamics are governed by a discrete-valued latent state process, where each state corresponds to a distinct spectral profile. We model these spectral profiles as realizations of state-specific and frequency-band specific beta distributions. From the estimated HMM, we are able to objectively characterize the dynamic evolution of neurophysiological activity induced by ketamine. By applying this tool to analyze neural data from both NHPs and humans, we obtain precise estimates of time-scales associated with neurophysiological phenomena observed in both species (e.g., alternating gamma slow-delta activities). Overall, our work provides methodological innovation to analyze switching spectral signatures as well as neuroscience inferences that result from applying this analytic tool to neural recordings from multiple NHPs and patients.
The paper is organized as follows. In the following section 3.2 we present our beta-HMM formulation and useful statistics based on the estimated beta-HMM. In section 3.3 we describe the electrophysiology datasets to which we apply our analysis framework. We present the results from our beta-HMM analyses in section 4. Finally, in section 5 we summarize the core contribution of our work, highlight the connections between our work and that of others, and suggest next steps.

Ethics statement
All NHP procedures reported here followed the guidelines of the National Institutes of Health and were approved by the Massachusetts Institute of Technology's Committee on Animal Care. The Partners Institutional Review Board (IRB) approved this human retrospective observational study. The IRB approved a waiver of consent for this study due to the anonymous nature of EEG data.

Observation model.
To facilitate the comparison of spectral observations across different data sets whose spectral power may vary in amplitude but represent comparable spectral profiles, we scale � S n ðωÞ to new vector y n 2 [0, 1] H whose h-th element is calculated as, where Q j ð � Sðo h ÞÞ denotes the function that operates on a real-valued sequence � Sðo h Þ ¼ f � S n ðo h Þg, n 2 [1, N] to yield the j-th quartile. The set of scaled power observations is denoted by Y = {y n }, n 2 [1, N]. The user-defined parameter λ h controls the mutual separation of points in the y nh -space relative to their corresponding map in the � S n ðo h Þ-space. We set this value to be The intuition behind the choice of parameters in the logistic function is that they result in scaling of the data such that the median is scaled to 0.5 and the data between the 1st and 3rd quartiles is scaled approximately linearly between [0.25, 0.75] (Sec A2 in S1 Appendix). Data below the 1st quartile and beyond the 3rd quartile are scaled non-linearly, such that very small values are scaled to be close to zero and very large values are scaled to be close to one. Compared to linear 0-1 scaling, the logistic scaling method is less sensitive to outliers which, in linear scaling, can lead to low variance observations and thus reduce the likelihood of detecting distinct spectral profiles. With logistic scaling, the data in � Sðo h Þ is scaled such that the standard deviations of the observations in each frequency band (y h ) are between [0.23, 0.31] for all data sets reported here.
In our model, we assume that for the n-th time-window, the observation y n is a manifestation of an underlying latent brain state, z n , at the same instant. We also assume that z n is a random process that transitions among K possible discrete states, each with its own characteristic observation probability distribution whose mean corresponds to a distinct spectral profile. We represent a state-specific and frequency band-specific probability density function (pdf) of the observations as a standard beta distribution: Here, p(y) denotes the pdf of a continuous-valued random variable y. In Eq (2), a hk > 0 and b hk > 0 are real-valued scalar parameters of the beta distribution that correspond to the scaled power in the h-th (h 2 [1, H]) frequency band when the state k 2 [1, K] is active. The normalization term Bða hk ; b hk Þ denotes the beta function [57].
Our assumption of statistical independence across the H distinct frequency bands conditioned on an underlying state, as implied in Eq (2), leads to the conditional joint pdf (jpdf) on y n , where ϕ k � {a hk , b hk } corresponds to the set of 2H beta distribution parameters associated with state k. The set of beta parameters across all K states is denoted as ϕ � {ϕ k } (Sec A3 in S1 Appendix). While the beta distribution can be multimodal when both a hk and b hk � 1, we require that each state-specific pdf is unimodal. A multimodal pdf associated with a latent state would indicate that the same state is associated with two different spectral profiles, which will confound our intended interpretation of the states. Therefore, to ensure the unimodality in the marginal pdf's we restrict the domain of the beta distribution parameters to avoid scenarios when both a hk � 1 and b hk � 1 for a given h and k. In summary, we chose to work with the beta distribution because it is a continuous probability distribution defined on the interval [0, 1]. Furthermore, it is part of the exponential family which leads to efficient maximum likelihood estimation (Sec A3 in S1 Appendix). The beta distribution is highly flexible as it allows for unimodal distributions where the modes can lie anywhere in [0, 1].

State transition model.
We assume the transition dynamics of the discrete-valued latent state process, z n , to be governed by first-order Markov transitions characterized by a constant K × K transition probability matrix, A = {A jk }, such that Here, Pr(z) denotes the probability mass function of a discrete-valued random variable z. The initial state probability at the first time-window is characterized by a vector π = {π k } such that Together, the observation model (Eq (3)), state transition matrix (Eq (4)), and initial state probabilities (Eq (5)) define our state-space model, which we refer to as the beta-HMM.

Statistical analysis using beta-HMM parameters.
The model parameters of the beta-HMM can be estimated from one or more independent EEG/LFP recording sessions (Sec A3 in S1 Appendix). The estimated beta-HMM provides reduced-order summaries of the scaled spectral dynamics. Between any two states, we can compare the spectral power in a given frequency band, as well as duration statistics derived from the transition matrix.
Summarizing scaled spectral power in a given frequency band. The set of beta distributions provides rich information about the scaled spectral observations from which the parameters are estimated. One helpful statistic to summarize each beta distribution is Pr(Y h > 0.5|Z = k) where (Y h |Z = k)*Beta(a hk , b hk ). The notation (Y h |Z = k)*Beta(a hk , b hk ) denotes a random variable (Y h |Z = k) that is distributed according to the beta distribution Beta(a hk , b hk ). A realization of (Y h |Z = k) would represent a scaled power observation in the h-th frequency band for a known latent state Z = k. Pr(Y h > 0.5|Z = k) describes the probability that an observation from the h-th frequency band and k-th state is greater than the median scaled power in the hth frequency band.
Comparing scaled spectral power in a given frequency band between any two beta-HMM states. One way to quantify how different any pair of beta distributions are from one another is to calculate the Kolmogorov-Smirnoff (KS)-distance from a large number of samples generated from the two distributions (we use 100,000 samples here). The KS-distance measures the maximum difference between two empirical CDFs, so values close to zero indicate that the underlying distributions are similar, and larger values indicate that the underlying distributions are dissimilar. Alternatively, we can leverage the analytic tractability of the beta distributions to analyze a difference random variable, Δ hjk = X hj − X hk , where X hj * Beta(a hj , b hj ) and X hk * Beta(a hk , b hk ). (In this section, we substitute (Y h |Z = k) with X hk to simplify the subsequent notation.) Here, the parameter sets {a hj , b hj } and {a hk , b hk } respectively correspond to a state j and another state k in the h-th frequency band. We use Pr(Δ hjk � 0) to denote the cumulative distribution function (cdf) value at 0 for Δ hjk 2 [−1, 1] such that where the pdf of the random variable Δ hjk can be expressed at a given realization δ as, where the joint pdf p(x hj , δ) is defined as, A general approach to calculating probability distributions of difference random variables can be found in the work by Cook and Nadarajah [58]. In this case, if Pr(Δ hjk � 0) = α, then there is probability of α that X hj � X hk , or equivalently a probability of 1 − α that X hj > X hk . In the ensuing discussion we adopt the following notational convention for clarity: For a withinsubject pairwise comparison of HMM states estimated for a single NHP (Sec 4.3), we add a superscript indicating the NHP (e.g. D MJ hjk ¼ X MJ hj À X MJ hk , for NHP MJ). When we compare states between two NHPs, we add superscripts indicating both NHPs (e.g. D MJ;LM hjk ¼ X MJ hj À X LM hk , when comparing an HMM state from NHP MJ to one from NHP LM). When we perform pairwise comparison of HMM states estimated from the human OR dataset (Sec 4.4), we add the superscript H, D H hjk . Comparing duration statistics between (i) any two beta-HMM states, and (ii) any two neurophysiological states each characterized by one or more beta-HMM states. We present several statistics to describe the dynamics of the model. If d k indicates a positive integer-valued random variable of the duration spent in state k after transitioning into it and before transitioning out to a different state, then the mean duration in that HMM state is given bŷ where Prðd k Þ � Prðz n:nþd k ¼ kjz n ¼ k; z nÀ 1 6 ¼ kÞ ¼ ð1 À A kk ÞA d k À 1 kk , which holds for any n > 1 [35].
In cases where a subset of beta-HMM states, k � � {1, � � �, K} corresponds to a neurophysiological phenomenon of interest (e.g. multiple HMM states corresponding to sub-states of gamma activity, or, multiple HMM states corresponding to the period between two consecutive gamma activities), we use the following Monte Carlo approach to estimate the consecutive time spent in this subset of beta-HMM states. Using the estimated transition matrixÂ and randomly sampled state z 1 , we generate a Markov sequence, z 1:N (N = 2000 in our study). For all z n 's taking values in k � , the corresponding state label is reassigned as z n = 1. The states that do not take values in k � are reassigned values as z n = 0. We then calculate the mean duration spent in the neurophysiological state,d k � , as the mean number of consecutive time points where z n = 1 holds. Similarly, for the same neurophysiological state k ? we also calculate a mean interval statistic as the average number of consecutive timepoints where z n = 0; this is an estimate of the time interval between the two consecutive occurrences of k ? . By repeating this procedure N MC times (= 4000 in our study) we calculate median and 95% confidence intervals of the state-specific mean duration and mean interval statistics, together referred to as the duration statistics in the ensuing discussion.

Experimental recordings
3.3.1 NHP LFP recordings under ketamine anesthesia. LFPs were recorded in multiple separate recording sessions from two rhesus macaques (Macaca mulatta) aged 14 years (NHP MJ, male, 13.0 kg) and 8 years (NHP LM, female, 6.6 kg). LFP recordings from 4 sessions in NHP MJ and 5 sessions in NHP LM were analyzed in this study. During each recording session, ketamine was administered as a single 20 mg/kg bolus intramuscular dose. This is a putative high dose for sedation and low dose for surgical anesthesia in NHPs [59][60][61]. Fifteen minutes prior to ketamine administration, glycopyrrolate (0.01 mg/kg) was delivered to reduce salivation and airway secretions. In both NHPs, LFPs were recorded from a 8 × 8 iridiumoxide contact microelectrode array ("Utah array", MultiPort: 1.0 mm shank length, 400 μm spacing, Blackrock Microsystems, Salt Lake City, UT) implanted in the frontal cortex (vlPFC). LFPs, recorded at 30 kHz, were low-pass filtered to 250 Hz and then downsampled to 1 kHz. LFPs were continuously recorded from 1-5 minutes prior to ketamine injection up to 18-20 minutes following ketamine injection. In the ensuing figures representing neural timeseries data, the time 0 denotes the time-point when the ketamine bolus was administered.

Human EEG recordings under ketamine anesthesia.
The data analyzed in this manuscript were acquired during an EEG study of ketamine-induced general anesthesia. We reviewed and selected 9 patients (Age: 51.8±11.7 years, Weight: 84.5±15.5 kg, mean ± standard deviation) from our database who were administered an intravenous bolus dose of ketamine as the sole hypnotic drug for the induction of general anesthesia. The ASA scores in all these 9 patients were less than or equal to 3, and none of the patients had a history of any neurological or psychiatric disorders. Prior to the ketamine bolus, patients were administered midazolam (n = 8; 1.81 ± 0.53 mg) and/or fentanyl (n = 7; 164.29 ± 80.18 mcg) for anxiolysis and to block the sympathetic response to laryngoscopy, respectively. To induce general anesthesia (GA) a bolus dose of ketamine (mean ± standard deviation; 182.22 ± 29.06 mg, 10 mg/ml) was administered, and intubation was carried out using succinylcholine, cisatracurium, or rocuronium for muscle relaxation. EEG data were acquired using Sedline monitor (Masimo Inc, USA). The standard Sedline Sedtrace electrode arrays were placed on the forehead that approximated the positions of Fp1, Fp2, F7, and F8, the ground electrode at Fpz, and the reference electrode 1 cm above Fpz. For our analysis, we use the data recorded from Fp2. Data were recorded with a pre-amplifier bandwidth of 0.5 to 92 Hz, a sampling rate of 250 Hz, with 16-bit, 29 nV resolution. Electrode impedance was less than 5 kO in each channel. We analyzed continuous EEG recordings starting from 2 minutes pre-ketamine bolus and up to 5.5 min-14 min post-ketamine bolus, and prior to the administration of any additional hypnotic drugs used to maintain GA.

Testing the beta-HMM analysis framework against simulated ground truth
We first validated the ability of our beta-HMM analysis framework to accurately retrieve a known Markov sequence and associated observation distributions by simulating spectrograms generated from pre-specifed Markov processes (see Table A2 in S1 Appendix for the simulation algorithm). We simulated spectral dynamics comparable to those caused by ketamine by using a spectrogram S({ω i }) calculated from one session of NHP experimental data (as described in Sec 3.2.1) as input to the algorithm. Additional inputs to this algorithm are userprescribed HMM parameters, the number of states K, π and A, as well as the duration M of the simulated data. This algorithm outputs a realization of the latent state path, z = {z m }, m 2 [1, M], corresponding simulated spectrogram, S sim ({ω i }), and Y. Salient features of the algorithm in Table A2 (S1 Appendix) are as follows. First, K distinct spectral clusters are created from the given NHP LFP spectrogram using an unsupervised clustering algorithm (different from beta-HMM) that does not respect the dependence across consecutive time-points. Then a Markov sequence is simulated using the known parameters, π and A. Finally, to generate the spectrogram with K underlying latent states with Markov state transitions, a spectrum is selected for a given instant by uniformly sampling from one of the K spectral clusters that corresponds to the realization of the latent state at the same instant. As described in Sec 3.2, Y is calculated from S sim ({ω i }) and used to estimate a K-state beta-HMM (Sec A3 in S1 Appendix).
For a given K, we generate 100 realizations of Y based on K spectral clusters. For each realization of Y, we estimate a K-state beta HMM to determine a set of estimated model parameterŝ Y ¼ fπ;Â;φg, and the estimated state path z � (calculated using the Viterbi algorithm). We assess goodness of fit for each realization of Y by comparing the estimated state path, z � , and model parameters,Ŷ, to the ground truth state path, z, and model parameters, Θ, respectively.
To assess the accuracy of the estimated path z � compared to the known path z, we calculated the fraction of time points for which the estimated and known path are identical: To assess the accuracy of the observation distribution estimation, we calculated the mean KS-distance between the HK ground truth and estimated beta distributions. KS-distances close to zero indicate that the estimated and ground truth distributions are similar; the maximum possible KS-distance is 1. We also calculated the absolute errors in the estimated initial state and state transition probabilities, denoted respectively by � A and � π , Note the the maximum possible value for Thus, the denominators scale the absolute errors to the range of [0, 1]. We report the median and 90% confidence intervals of each of these metrics calculated across 100 simulated spectrograms for each model order in Fig 2. Simulation analysis revealed that the model parameters of simulated spectrograms, generated from ketamine-induced neural data and characterized by known Markov transition dynamics, can be accurately estimated with the beta-HMM framework for low (K < 7) model orders. The performance of the beta-HMM analysis framework on a typical simulated spectrogram (Fig 2A) demonstrates that the estimated state trajectory (Fig 2C) matches well with the actual state trajectory (Fig 2B). Across 400 simulated datasets for models with 2 through 5 states (100 simulations each), the estimated state trajectory was estimated with high accuracy (PathAccuracy > 0.98 for all of 400 independent simulations), the beta distributions were estimated with high accuracy (mean KS-distance <8.78 × 10 −3 ), and the transition matrix was estimated with low relative error (� A < 0.01). For models with 2-5 states, the initial probability, which corresponds to a single data point, was also accurately estimated (� π < 3.09 × 10 −4 for all but three of the 400 independent simulations). Beyond 6 states, there was a significant drop-off in estimation accuracy. A potential reason for this occurrence is redundancy across some of the K > 6 spectral clusters derived from the original NHP spectrogram, which may not contain more than 6 clusters of distinct neural activity. Overall, this numerical experiment generated detailed insights on the level of inaccuracy in our beta-HMMs when estimated from spectral data derived from real neural activity induced by ketamine. The simulations revealed that the beta-HMM framework was able to estimate Markov state trajectories and model parameters with high accuracy for models with 2-5 states.

Demonstration of the beta-HMM analysis on a single observation sequence of NHP LFP data (L = 1 case)
Using the numerically tested beta-HMM analysis framework, we analyzed a single session high-dose ketamine NHP LFP recording (Sec 3.3.1). From qualitative analysis of the timeseries data and spectral estimates (S1 Fig), we identified two prominent neural signatures induced by ketamine: one is the gamma activity that has been previously described [10], and one is prominent slow-delta activity that also varies in power on a similar time scale to the gamma activity. The gamma and slow-delta activities tend to co-occur, resulting in four substates characterized by high or low gamma activity and high or low 0-4 Hz activity. We also sought to capture the transition from pre-to post-ketamine bolus neurophysiology, so we  (10)) in panel D, the KS-distance between the estimated and the actual state-specific and frequency band-specific beta distributions (averaged across all states and all frequency bands) in panel E, the relative error in transition matrix (Eq (11)) in panel F, and relative error in the initial state probabilities (Eq (12)) in panel G. For each specified number of states, 100 realizations of spectrograms were simulated and the beta-HMM analysis was performed on each.
The 5-state beta-HMM provided a quantitative characterization of the neurophysiological dynamics induced by ketamine (Fig 3A and 3B). The beta-HMM analysis objectively identified alternating gamma slow-delta dynamics in NHPs similar to those described by Akeju, et al. [10]. We found that the HMM identified alternating states which corresponded to periods of prominent gamma oscillations and prominent slow-delta oscillations (Fig 3C and 3D). From visual inspection of the the state-segmented spectrograms and further analysis of beta distributions (Fig 3E and 3G and S4 Fig), we identified relevant neurophysiological states with distinct spectral signatures. For example, in the state-segmented spectrogram for state 5 (Fig 3E), we can see that high spectral power in the 25-35 Hz frequency band (h = 6) is associated with a high probability that (Y 6 |Z = 5)>0.5 (Fig 3G). (See Sec A6 in S1 Appendix for a tutorial-styled explanation of how the estimated beta pdf's are used to analyze the spectral properties of the HMM states). Based on the beta distributions, we found that HMM states 4 and 5 both have high gamma power (25)(26)(27)(28)(29)(30)(31)(32)(33)(34)(35)
The transition matrix (Fig 3F) quantifies the alternating dynamics illustrated in Fig 3C and  3D. From the transition matrix, we can also determine the expected duration of each state, as denoted on the diagonal terms indicating self-transition probabilities (Eq 9). From the transition probabilities and optimal state trajectory, we infer there is a high probability of the state sequence 2 ! 3 ! 4 ! 5 ! 2. This confirms the alternating dynamics between time-localized high power activity in the gamma and low frequency bands induced by ketamine. The beta-HMM provides further evidence that the neural activity induced by ketamine is distinct from pre-ketamine neural activity. States 2, 4, and 5 do not occur prior to the ketamine bolus ( Fig  3A and 3B). Furthermore, HMM state 1 has sustained prevalence prior to the bolus and does not occur after a brief initial period (approximately 2 minutes) post-bolus. Therefore, HMM state 1 can be associated with a pre-ketamine neurophysiological state. This is further supported by the fact that the transition probability from any of the post-bolus states back to state 1 is less than 0.01.

Subject-specific beta-HMMs estimated from multiple recording sessions (L > 1 case) from 2 NHPs
A review of summary statistics from session-specific beta-HMMs estimated from multiple recording sessions of each NHP indicated the presence of subject-specific HMM states that have common features (both spectral and temporal) across sessions. This further suggested the utility of a subject-specific HMM for each NHP. Therefore, we fitted two 5-state beta-HMMs to multiple LFP recording sessions from 2 NHPs: 4 sessions for NHP MJ and 5 sessions for NHP LM (Sec 3.3.1). A key implementation feature in our beta-HMM analysis from multiple sessions (for a given group) is that we do not concatenate multiple sessions; rather we treat them as mutually independent in our EM algorithm implementation (as discussed in Sec A3 in S1 Appendix). Thus, the beta-HMM analysis respects the sequential nature of the observations within a session, while maintaining the distinction across multiple sessions. Since there were at least three days between each NHP session, we treated them as mutually independent and therefore utilized the EM algorithm based on Eq. (A4) in S1 Appendix to estimate the model parameters.
We present the optimal state trajectories and state-specific and frequency band-specific observation distributions estimated from 4 of 9 total sessions of 2 NHPs in Fig 4 (see S2 Fig for  all sessions). Although the amplitude of the spectral power tended to vary across sessions, our use of scaled power facilitated identification of common states with similar spectral profiles across multiple sessions (Fig 4A-4C and 4E-4G). The beta distributions provide detailed information about the spectral content of each HMM state, and we summarize the key results here. In both NHPs, we found that the beta distributions in each frequency band were unique across all states (KS-distance >0.04 for all pairwise comparisons; 140 comparisons = 10 pairs of states per frequency band per NHP × 7 frequency bands × 2 NHPs) (Fig 4D and 4H As demonstrated in Fig 3C and 3D, the multisession beta-HMM model captures the alternating dynamics induced by ketamine in both NHPs (Fig 4C and 4G). As noted in the single session analysis, there was some overlap between the gamma and slow-delta activity, resulting in We used the subject-specific HMMs to further analyze the duration statistics and spectral content across states, both within a subject as well as between the two subjects ( Fig 5). We found significant differences in the neurophysiological state dynamics between the two NHPs ( Fig 5A-5C). The mean duration of the gamma activity,  The transition matrices for NHP MJ and NHP LM with expected state duration (in seconds) indicated on the diagonal for the states that occur following the ketamine bolus. (D-F) Heatmaps indicating the probability of the event that scaled power in state j (of NHP MJ in panels D and F, and of NHP LM in panel E) is less than or equal to the scaled power in state k (of NHP MJ in panel D, and of NHP LM in panels E and F) for each of the 7 frequency bands. Note that in panels D and E, the upper triangle (1 < j < k < 5), omitted to avoid redundancy, is equal to 1 minus the lower triangle (1 < k < j < 5). The color indicates Pr(Δ � 0) which represents Pr(X hj − X hk � 0), where X hj is a random variable characterized by state j's beta pdf in frequency band h and X hk is a random variable characterized by state k's beta pdf in the same frequency band h.
https://doi.org/10.1371/journal.pcbi.1009280.g005 PLOS COMPUTATIONAL BIOLOGY (Fig 5B and 5C), it is also apparent that the expected duration in each HMM state in NHP MJ is greater than that of NHP LM, indicating faster state transition dynamics in the latter subject.
We used the ML estimated beta distributions to perform pairwise comparisons of the scaled power between any two states from either NHP (Fig 5D and 5E). By investigating Pr(Δ hjk � 0; a hj , b hj , a hk , b hk ) (per Sec 3.2.4), we were able to determine the probability that HMM state j had lower scaled power than HMM state k in the h-th frequency band. For this section, we simplify the notation by referring to Pr(Δ hjk � 0;a hj , b hj , a hk , b hk ) as Pr(Δ hjk � 0), and summarize the key findings. In both NHPs, we found that the scaled powers in the 0-1 Hz and 35-50 Hz bands during the pre-ketamine state (HMM state 1) were lower relative to every other state (0:00 < PrðD MJ hj;1 � 0Þ < 0:34 and 0:00 < PrðD LM hj;1 � 0Þ < 0:16 for h 2 {1, 7} and j 2 [2, K]). This provides further evidence that the oscillatory dynamics in these two frequency bands

Human-specific beta-HMM estimated using observations from multiple patient subjects (L > 1 case)
We investigated if the beta-HMM analysis framework could be applied to infer the common ketamine-induced neurophysiological state dynamics observed in the scalp EEG of multiple OR patients (Sec 3.3.2). We fit a single beta-HMM to independently collected EEG data from L = 9 human patients and generated a quantitative description of ketamine-induced neurophysiological dynamics for a typical patient (Fig 6). For this analysis, we chose a model order of K = 6 states, which allowed for classification of broadband high power, typical of noise artifacts in the OR, as an additional HMM state.
As with NHPs, we found that the beta distributions in each frequency band were unique across all states (KS-distance >0.04 for all pairwise comparisons; 105 comparisons = 15 pairs of states per frequency band × 7 frequency bands) (Fig 6D and S7 Fig) There were distinct neurophysiological differences between the human EEG and the NHP LFP (Sec 4.3). First, the gamma activity in humans (captured by states 4 and 5) was lower in frequency and extended into the beta (12-25 Hz, h = 5) frequency range (Pr(Y 5 > 0.5|Z = 4) = 0.81, Pr(Y 5 > 0.5|Z = 5) = 0.86). Second, in the human EEG, state 4 had both prominent gamma and slow-delta activity, indicating that the gamma and slow-delta activities overlapped rather than alternated. Third, there was high slow-delta power and, in some patients, high gamma power (S3 Fig) throughout the human EEG recordings, which resulted in HMM states that were present both before and after the ketamine bolus. Finally, soon after the time of the ketamine bolus, in several patients (Fig 6B and S3 Fig), there was broadband high power that, in the context of OR cases, is difficult to distinguish from noise. This activity is classified with the 6-th state, which also corresponds to periods of obvious noise artifacts (e.g. at t = 0 min in Fig 6A).
Quantitative analysis of the beta-HMM dynamics (Fig 7A)  cyclically through the HMM states induced by ketamine (Fig 7B)-the most probable PLOS COMPUTATIONAL BIOLOGY transition sequence, starting from state 1, is 1 ! 2 ! 3 ! 5 ! 4 ! 2. All states have a low probability of transitioning to state 6 (A j,6 < 0.012 for j 2 [1,5]), further indicating that this state captures spurious noise artifacts. Further quantitative analysis (Fig 7C) of human EEG beta-HMM parameters revealed higher scaled power in the 25-35 Hz range for HMM states 4 and 5 compared to states 1-3 (0:00 < PrðD H 6;j;k � 0Þ < 0:37 for j 2 {4, 5} and k 2 [1,3]). For state 5, the increased power extended into the 35-50 Hz range (0:01 < PrðD H 6;5;k � 0Þ < 0:24 for and k 2 [1,3]). States 2 (C) Heatmaps indicating probability of the event that scaled power in state j is less than or equal to the scaled power in state k for 1 � k � j � 6 and for each of the 5 frequency bands. Note that the upper triangle (1 < j < k < 6), omitted to avoid redundancy, is equal to 1 minus the lower triangle (1 < k < j < 5). The color indicates Pr (Δ � 0) which represents Pr(X hj − X hk � 0), where X hj is a random variable characterized by state j's beta pdf in frequency band h and X hk is a random variable characterized by state k's beta pdf in the same frequency band h. https://doi.org/10.1371/journal.pcbi.1009280.g007

Discussion
Our beta-HMM analysis framework provided parsimonious summaries of ketamine-induced slow-gamma alternating dynamics for NHP LFP and human EEG recordings following highdose ketamine administration. The input to the analysis framework was one or more sequences of instantaneous power in multiple canonical frequency bands. The power values were averaged across frequencies within each discrete band and scaled between 0 and 1. The scaling enabled identification of latent HMM states each with a distinct probability distribution on power values lying between a very low (corresponding to a 0) and a very high value (corresponding to a 1) in each of the chosen frequency bands, irrespective of the absolute value of spectral estimates. A key assumption incorporated in the HMM is that at any instant in time, the probability distribution from which the observation is sampled depends only on the discrete value of an underlying latent state process at the same instant. This state process itself is assumed to evolve according to a first-order Markovian transition dynamics. The key outputs of our analysis comprised the model parameters estimated using an EM algorithm and a corresponding state trajectory estimated using the Viterbi algorithm. The state trajectory provided an objective and efficient segmentation of LFP and EEG data. The estimated statespecific and frequency band-specific beta pdf's characterized the spectral representation of each HMM state, whereas the transition matrix characterized the underlying dynamics. Using the estimated observation pdf's, we objectively defined neurophysiological activities (e.g. gamma activity, slow-delta activity) in terms of underlying HMM states. Furthermore, by using the estimated state transition matrices, we calculated mean duration and mean interval statistics (as defined in Sec 3.2.4) corresponding to each of these neurophysiological activities.
Our analyses revealed an alternating pattern of states characterized primarily by gamma and slow-delta activities. The mean duration of the gamma activity was 2.2s ([1.7, 2.8]s) and 1.2s ([0.9, 1.5]s) for the two NHPs, and 2.5s ([1.7, 3.6]s) for the nine human subjects. The mean duration of the slow-delta activity was 1.6s ([1.2, 2.0]s) and 1.0s ([0.8, 1,2]s) for the two NHPs, and 1.8s ([1.3, 2.4]s) for the nine human subjects. Our characterizations of the alternating gamma slow-delta activity revealed five sub-states that show regular sequential transitions. Thus, our objective beta-HMM analysis framework enabled us to precisely characterize the dynamics of ketamine-induced gamma and slow-delta activities, and to compare the timescales of these neurophysiological states between subjects of the same species.
Our work advances the development and application of LFP or EEG data analysis tools. The beta-HMM analysis framework falls under the broad category of research aimed at developing time-frequency state space models to analyze neural dynamics [43,52,[62][63][64][65], and applying them in principled interpretation of neural time-series data for unconsciousness research [11,62]. A few key distinguishing features of the beta-HMM analysis framework are as follows. The assumption of state-specific and frequency-band specific beta distributions as marginal distributions of the instantaneous observation vectors allows for a flexible framework to characterize highly variable spectral profiles. This assumption, along with the assumption of Markovian state transitions, allowed us to directly characterize the discontinuous switching among multiple spectral profiles, as well as identify segments with similar spectral properties (Figs 3C, 3D; 4C, 4G; and 6C). Furthermore, our beta-HMM analysis framework accounts for the disjoint nature of data from separate sessions without requiring any concatenation of sessions or consequently any arbitrary choice of concatenation order. Potential extensions of this work on the statistical methods development front may include modeling the transition matrix as time-varying, which would account for longer time-scale variations (e.g. varying drug-levels). Furthermore, generalizing the current single-channel LFP or EEG data-based analysis framework to incorporate multi-channel LFP or EEG data would allow for neurophysiological characterization of ketamine-induced altered arousal states based on spatio-temporal statistical relationships in the observed data. Going beyond the context of ketamine-induced altered arousal states, the analysis framework itself can be relevant in general use-cases where LFP or EEG data demonstrate discontinuous, repeating transitions among time-limited, band-limited spectral signatures, such as those observed during sleep [43,66] or burst suppression in general anesthesia [67,68].
Our beta-HMM research can inform future neuroscience inquiries. While there are several hypotheses for the generation of gamma activity under ketamine [21], there is not yet a clear understanding of how the alternating gamma slow-delta activity is generated. Preliminary work indicates that there may be a cycle of cortical inhibition and dis-inhibition due to activity-dependent ketamine NMDAR inhibition of cortical interneurons and pyramidal neurons [22]. Our analyses of the NHP datasets using the beta-HMM led to the identification of neurophysiological activities that were primarily associated with the period following administration of a ketamine bolus. The summary statistics from each NHP (Fig 5) inform which frequency bands show significant activity, and how this activity varies over time. The duration and interval statistics of the alternating gamma and slow-delta activities that we estimated for NHP LFP can serve as quantitative constraints in the design of rhythm-generating neuronal microcircuit models that would mimic the neurophysiological dynamics caused by ketamine, similar to ones previously done for other anesthetics [23,24,26,69,70]. Furthermore, since ketamine is also used in treatment of depression [71,72], in pharmacologic models for schizophrenia [20,73,74], and in studies of altered states of consciousness [75], insights into ketamine's effects on brain state dynamics will provide neuroscientific insight beyond the field of anesthesia.
Our beta-HMM research can also inform future clinical inquiries. Ketamine is known to produce distinct oscillatory signatures in the electroencephalogram (EEG) of healthy volunteers and patients [9,10]. These oscillatory signatures are starkly different from those produced during propofol-induced unconsciousness [76,77]. These differences in oscillatory dynamics, particularly the presence of alternating gamma and slow-delta activities under ketamine, indicate that spectral markers of propofol-induced unconsciousness are not reliable markers for tracking ketamine-induced altered arousal states. However, in both the NHP LFP and human EEG recordings, we observed dramatic changes in the oscillatory state dynamics between the awake and post-ketamine states, suggesting that the alternating gamma slow-delta phenomena induced by ketamine is a marker of altered arousal. The existence of consistent neurophysiological activity observed in EEG following ketamine bolus administration in multiple patients indicates promise in precise clinical monitoring of these states. Towards this goal, as demonstrated in this work, the beta-HMM analysis framework can provide precise, objective characterization of the neurophysiological dynamics associated with ketamineinduced altered states of arousal. In summary, our work can be regarded as a part of the ongoing neuroscience research efforts to investigate alternating dynamical states of the brain and how these states transition from one to the other [78]. Alternating brain states have been previously identified in anesthesia-induced altered states of consciousness across different imaging modalities, species, analytic tools and anesthetic drugs [11,79,80]. The experimental and analytic framework developed here can be used to investigate differences in brain state dynamics between awake and altered states of consciousness. This investigation could be relevant in future studies looking into the theories of consciousness (e.g. [81]) as well as to improve the understanding of mechanism of actions of these drugs in therapeutic applications (e.g. [13,14]).
There are two primary limitations in both the data set that we analysed here and in our assumptions. The first limitation is a need for prospective studies in humans and animals with carefully recorded simultaneous neural, behavioral response and physiological activities. Future experimental studies with larger cohorts can address the current study's limitation due to low number of subjects in both the NHP LFP and human EEG analyses. The consistency of the alternating gamma-slow spectral activity from 9 patients and 2 NHPs during ketamineinduced altered arousal states indicate that transition dynamics between the gamma and slowdelta neurophysiological states can be a candidate biomarker of ketamine-induced unconsciousness. However, to reliably establish these EEG or LFP correlates of unconsciousness, experimental studies in healthy volunteers and animal subjects with simultaneous monitoring of behavioral response, physiological parameters and neural activity will be essential. The second limitation is a need for extension to real-time tracking of neurophysiological states during ketamine-induced general anesthesia. Our current approach utilizes the entire data set to calculate appropriate scaling factors for transforming spectral power to real numbers lying between 0 and 1, and thus (in the present formulation) cannot be executed in real time. A larger study may allow for the development of a generalized scaling approach. This would facilitate objective identification of ketamine neurophysiological states in real-time. While a common limitation of EEG analysis in real-time OR settings is the presence of motion and other noise artifacts in recordings, the beta-HMM framework is well poised to handle these artifacts, as observations that fall far outside the distribution of the signal can be assigned to a distinct HMM state.
In conclusion, we have developed a detailed analysis framework for ketamine-induced neurophysiological phenomena. Our work provides a methodological innovation to analyze switching spectral dynamics in single-channel electrophysiological data. The scaling of the frequency band-wise power to [0, 1] interval followed by its analyses in the beta-HMM framework facilitates estimation across multiple recording sessions, as well as comparisons between sessions. To our knowledge, this is the first application of an HMM with beta observation distributions for characterizing neural data. Our work also provides insights into the neural dynamics due to ketamine anesthesia. By applying our methodological innovation to LFP data from NHP subjects and EEG data from human patients collected during ketamine anesthesia, we have identified distinct neurophysiological states by their spectral signatures and their duration statistics. Our quantitative findings will inform future neurophysiological models as well as clinical biomarker search. Also, the generalizability of the beta-HMM framework indicates utility beyond what has been reported here. Future work will investigate how the spectral dynamics revealed by our analysis contribute to ketamine-induced altered states of arousal.
Supporting information S1 Appendix. In the appendix, we provide additional detail in the calculation of band-wise power (Sec A1), scaling of the beta-HMM observations (Sec A2), estimation of model parameters (Sec A3), and EM algorithm (Sec A4). We also provide the algorithm for simulating a spectrogram with known Markov model parameters (Sec A5, Table A2), and a tutorial on interpreting beta distributions (Sec A6). In Table A1, we provide a glossary of mathematical symbols. (PDF) S1 Fig. LFP time series and corresponding spectrograms. A typical 20 second epoch (following a high-dose ketamine bolus) of LFP from a single electrode of a multi-electrode array located in the frontal cortex (vlPFC) of NHP MJ is presented in panel A. Panel B shows the multitaper spectrogram for the same electrode, and panel C shows the average multitaper spectrogram for the vlPFC electrode array. Time 0 corresponds to 400 seconds after the ketamine bolus was administered. We fit our beta-HMM to the reduced-order representation of the corresponding spectrogram (derived using Eqs. (1) and (A1) in S1 Appendix), presented in panel C. The resulting optimal segmentation is represented by the colored vertical bars overlaid on the neural time-series in panel A. State 2 is shown in blue, state 3 in red, state 4 in yellow, and state 5 in green. (State 1, which corresponds to the time before the ketamine bolus, is not present in this epoch. The state-specific and frequency-band specific beta pdfs corresponding to a 5-state beta-HMM were estimated from L = 1 session of LFP recording from NHP MJ. For each state (viewed row-wise) and a frequency band (viewed column-wise), the corresponding subplot presents (1) the empirical pdf plotted as a histogram based on the observations that correspond to the optimal segmentation of the data sequence, and (2) the continuous beta pdf with parameters estimated by the EM algorithm. Each color distinguishes a state as in The state-specific and frequency-band specific beta pdfs corresponding to a 5-state beta-HMM were estimated from L = 4 sessions of LFP recording from NHP MJ. For each state (viewed row-wise) and a frequency band (viewed column-wise), the corresponding subplot presents (1) the empirical pdf plotted as a histogram based on the observations that correspond to the optimal segmentation (Viterbi algorithm) across the L = 4 sessions, and (2) the continuous beta pdf with parameters estimated by the EM algorithm. Each color distinguishes a state as in The state-specific and frequency-band specific beta pdfs corresponding to a 5-state beta-HMM were estimated from L = 5 sessions of LFP recording from NHP LM. For each state (viewed row-wise) and a frequency band (viewed column-wise), the corresponding subplot presents (1) the empirical pdf plotted as a histogram based on the observations that correspond to the optimal segmentation (Viterbi algorithm) across the L = 5 sessions, and (2) the continuous beta pdf with parameters estimated by the EM algorithm. Each color distinguishes a state as in Fig 4. (TIF) S7 Fig. State-and frequency-specific beta pdfs for 9 human sessions. The state-specific and frequency-band specific beta pdfs corresponding to a 6-state beta-HMM were estimated from L = 9 sessions of EEG recording from 9 human OR patients. For each state (viewed row-wise) and a frequency band (viewed column-wise), the corresponding subplot presents (1) the empirical pdf plotted as a histogram based on the observations that correspond to the optimal segmentation (Viterbi algorithm) across the L = 9 sessions, and (2) the continuous beta pdf with parameters estimated by the EM algorithm. Each color distinguishes a state as in Fig 6. (TIF) S8 Fig. Patient variability in state-wise mean scaled spectral power. In each state plot, one dot represents the mean scaled power in a specific frequency band for one patient, where the mean is calculated from the observations assigned to that state via the Viterbi algorithm. The mean scaled power across frequencies for each patient are connected with a line. The shaded region indicates the 95% confidence interval for the mean scaled power across 10000 samples of the mean of the corresponding beta distribution, where the mean is calculated from 200 independent samples of the beta distribution. Note that while some observed means exceed the 95% confidence interval of the expected mean of the corresponding estimated beta distribution, the overall trends are consistent across patients. (TIF)

S9 Fig. Exploration of how frequency resolution affects the outcome of the beta-HMM.
Panels A and B present analysis from NHP MJ sessions 1-4, panels C and D present analysis from NHP LM sessions 1-5, and panels E and F present analysis from the 9 human EEG sessions. In panels A, C, and E, a portion of the spectrogram is displayed on the top, and two state trajectories are displayed below, estimated from models with H = 51 frequency bands (i.e. the highest possible frequency resolution equal to that of the estimated spectrogram) and H = 7 � frequency bands (where 7 � frequency bands indicate the canonical frequency bands described in Sec 3.2.2). Panels B, D, and F present the key duration statistics presented in Sec 4.3 and 4.4 across varying frequency-band resolution. With the exception of the models corresponding to 7 � , the canonical frequency bands, all other models utilized evenly spaced frequency bands. For each model, 4000 Markov sequences of length N = 2000 were simulated using the estimated transition matrix. The mean duration and interval corresponding to the gamma and slow-delta activities were calculated for each realization of these sequences. Median and 95% confidence bounds across simulated sequences for models with equally spaced frequency bands are indicated in grey and, for the model utilizing canonical frequency bands, in blue. The mean duration and interval calculated from the estimated state trajectory (output of the Viterbi algorithm which uses the maximum likelihood beta-HMM parameters and the observations as input) are indicated by a cross (×) symbol. Note that the 95% confidence intervals of the durations across the varying frequency resolutions almost always overlap. This indicates that the specification of the frequency bands does not have a critical effect on the inferences related to the model dynamics. Thus, in the main manuscript, we chose to present the model estimated from canonical frequency bands that are commonly used to describe neural oscillations. (TIF)