Search tips
Search criteria 


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

Assessing Neuronal Interactions of Cell Assemblies during General Anesthesia


Understanding the way in which groups of cortical neurons change their individual and mutual firing activity during the induction of general anesthesia may improve the safe usage of many anesthetic agents. Assessing neuronal interactions within cell assemblies during anesthesia may be useful for understanding the neural mechanisms of general anesthesia. Here, a point process generalized linear model (PPGLM) was applied to infer the functional connectivity of neuronal ensembles during both baseline and anesthesia, in which neuronal firing rates and network connectivity might change dramatically. A hierarchical Bayesian modeling approach combined with a variational Bayes (VB) algorithm is used for statistical inference. The effectiveness of our approach is evaluated with synthetic spike train data drawn from small and medium-size networks (consisting of up to 200 neurons), which are simulated using biophysical voltage-gated conductance models. We further apply the analysis to experimental spike train data recorded from rats’ barrel cortex during both active behavior and isoflurane anesthesia conditions. Our results suggest that that neuronal interactions of both putative excitatory and inhibitory connections are reduced after the induction of isoflurane anesthesia.

I. Introduction

General anesthesia is a drug-induced state of reversible coma that is essential for performing most surgeries and many medical procedures [3]. In the United States, nearly 60,000 patients receive general anesthesia each day for surgical and medical procedures [1]. However, despite its ubiquity, the mechanisms by which anesthetic drugs induce the state of general anesthesia remains largely unknown. Developing a detailed understanding of these mechanisms may have significant functional implications for improved anesthesia practice, such as the design of safer drugs, or alternate methods of drug-delivery [3]. The study of general anesthesia from a systems neuroscience perspective is an approach that holds promise, using in vivo brain recordings to identify and assess alterations in neuronal circuitry during different conscious states. Previous investigations have used either non-invasive EEG or fMRI recordings to analyze the neural responses during general anesthesia (e.g., [7], [8], [13]).

Here, we focus on the problem of assessing functional connectivity of cell assemblies before or during general anesthesia using invasive spike train recordings. For this purpose, we apply a previously proposed point process generalized linear model (PPGLM) and variational Bayes (VB) algorithm [5], [6] to infer the neural interactions between population of neurons in a simulated network of size up to 100 neurons. The accuracy of our estimation results is validated on multivariate spike train data generated from a biophysical neuronal model under the baseline and anesthesia-like conditions. We further extend our analysis to experimental spike train data recorded from rats’ barrel cortex during both active behavior and isoflurane (inhalational) anesthesia conditions.

II. Methods

A. A Point Process Network Likelihood Model

A point process is a stochastic process with 0 and 1 observations. For the cth point process, let y1:Tc=(y1c,,yTc) denote the observed response variables during a (discretized) time interval [1, T], where ytc is an indicator variable that equals to 1 if there is a spike at time t and 0 otherwise. Therefore, multiple neural spike train data are completely characterized by a multivariate point process {y1:Tc}c=1C. To model the spike train point process data, we used the following logistic regression model with the logit link function. Specifically, let c be the index of target neuron, and let i = 1, … C be the indices of trigger neurons. The Bernoulli (binomial) logistic regression PPGLM is given by


where dim(βc) = d + 1 (where d = C × K) denotes total number of parameters in the augmented parameter vector βc={β0c,βi,kc}, and x(t) = {x0, xi, tk}. Here, x0 [equivalent] 1 and xi, tk denotes the spike count from neuron i at the kth time-lag history window. Since the spike count is nonnegative, therefore xi, tk ≥ 0. Alternatively, we can rewrite (1) as


which yields the probability of a spiking event at time t. It is seen from (2) that the spiking probability πt is a logistic sigmoid function of βcx(t); when the linear regressor βcx(t) = 0, πt = 0.5. Note that βcx (t) = 0 defines a (d + 1)-dimensional hyperplane that determines the decision favoring either πt > 0.5 or πt < 0.5.

Equation (1) essentially defines a spiking probability model for neuron c based on its own spiking history, and that of the other neurons in the ensemble. It has been shown that such a simple spiking model is powerful in inferring the functional connectivity among neuronal ensembles [4], [12], and in predicting single neuronal spikes from humans and primates based on collective population neuronal dynamics [15]. Here, exp(β0c) can be interpreted as the baseline firing probability of neuron c. Depending on the algebraic (positive or negative) sign of the coefficient βi,kc, exp(βi,kc) can be viewed as a “gain” factor (dimensionless, > 1 or < 1) that influences the current firing probability of neuron c from another neuron i at the previous kth time lag. A negative value of βi,kc will strengthen the inhibitory effect and move πt towards the negative side of the hyperplane; a positive value of βi,kc will enhance the excitatory effect, thereby moving πt towards the positive side of the hyperplance. Using this method, two neurons are said to be functionally connected if any of their pairwise connections is nonzero (or the statistical estimate is significantly nonzero).

Let θ = {β1, …, βC} be the ensemble parameter vector, where dim(θ) = C(1 + d). By assuming that the ensemble neuronal spike trains are mutually conditionally independent, the network log-likelihood of C-dimensional spike train data is given by [12]:


Note that the index c is uncoupled from each other in the network log-likelihood function, which implies that we can optimize the function L(βc) separately for individual spike train observations y1:Tc. For simplicity, from now on we will drop off the index c at notations ytc and βc when no confusion occurs.

B. Hierarchical Bayesian Modeling and VB Inference

It is known that maximum likelihood inference is prone to over-fitting by using a complex model. To address this issue, instead of maximizing the log-likelihood log p(y|x, β), we aim to maximize the marginal log-likelihood log p(y|x) or its lower bound L


where p(β|α) denotes the prior distribution of β, specified by the hyperparameter α. The variational distribution has a factorial form such that q(β, α) = q(β)q(α), which attempts to approximate the posterior p(β, α|y). This approximation leads to an analytical posterior form if the distributions are conjugate-exponential. The use of hyperparameters within the hierarchical Bayesian estimation framework employs a fully Bayesian inference procedure that makes the parameter estimate less sensitive to the fixed prior (unlike the empirical Bayesian approach) [14]. In our previous study [6], it was found that the VB approach produced excellent performance while dealing with sparse spiking data.

Specifically, we assume the following hierarchical Bayesian structure for the PPGLM [5], [6]:

p(β|α)~[mathematical script N](β|μ0,A1)[proportional, variant]exp(12(βμ0)[top top]A(βμ0)),p(α)=Πj=0dGamma(αj|a0,b0)

where Gamma(αj|a0,b0)=1Γ(a0)b0a0αja01eb0αj and A = diag {α} [equivalent] diag {α0, …, αd}. Here, we fix the hyperparameter as μ0 = 0 to favor a sparse solution.

Let ξ = {ξt} denote the data-dependent variational parameters (that are dependent on the input variables {xt}). In light of the variational approximation principle [10], one can derive a tight lower bound for the logistic regression likelihood, which will be used in the VB inference. Specifically, applying the VB inference yields the variational posteriors q(β|y) and q(α|y):

logq(β|y)=log[p with tilde](β,ξ)+Eq(α)[logp(β|α)]=log[mathematical script N](β|μT,ΣT),


which follow from updates from conjugate priors and posteriors for the exponential family (Gaussian and Gamma distributions). The derivations of {μT, ΣT} and {aT, bj, T} in equations (5) and (6) are given in [5]. The term [p with tilde](β, ξ) appearing in (5) denotes the variational likelihood bound for logistic regression:

logp(β,ξ)log[p with tilde](β,ξ)=t=1T(logσ(ξt)ξt2+ϕ(ξt)ξt2)[β[top top](ϕ(ξt)xtxt[top top])β]+β[top top](xtyt12xt),

where σ(·) is a logistic sigmoid function. The variational likelihood bound at the right-hand side of (7) has a quadratic form in terms of β, and therefore it can be approximated by a Gaussian likelihood. We can also further derive the variational lower bound of marginal log-likelihood L. (see [5]). The VB inference alternately updates (5) and (6) to monotonically increase L. The criterion for algorithmic convergence is set until the consecutive change of L is sufficiently small. Upon completing the VB inference, the confidence bounds of the estimates can be derived from the posterior mean and the posterior variance.

C. Measuring Goodness-of-fit of the Model

The goodness-of-fit of the point process models estimated from all algorithms is evaluated based on the Kolmogorov-Smirnov (KS) test [6]. The KS statistics are used to measure the maximum deviance between the empirical cumulative distribution function (cdf) of the time-rescaled data and the theoretical uniform cdf. If the curve falls within the 95% of the confidence intervals (CIs) of the KS plot, it suggests that the model has a sufficiently accurate characterization of the spike train data. In simulation studies, in addition to the KS test, we also compute the mis-identification error rate, which is the sum of the the false positive (FP) and false negative (FN) rates.

III. Data And Results

A. Simulations

In simulating a neuronal network with excitatory pyramidal cells and inhibitory interneurons, we used the single-compartment voltage-gated conductance models of the Hodgkin-Huxley type [9]. Specifically, the membrane potential of each cell is governed by a system of nonlinear differential equations of the form [11]:


where Iion, Isyn and Iapp denote, respectively, the ionic, synaptic and external applied currents specific to each cell, and cm is the membrane capacitance. Ionic currents are induced by the flow of charged particles across the cell membrane and are described by the equation


where m=m(V)mτm(V) and h=h(V)hτh(V). The parameters gx and Ex are constant, while p and q are non-negative integers. The term gxmphq is the ionic conductance and determines the behavior of the current as the voltage deviates from equilibrium. The most common ionic currents are related to the flow of sodium and potassium, which are responsible for generating action potentials (i.e., spikes). In addition, connectivity between neurons is established through the synaptic currents Isyn. Two such currents are considered herein: the excitatory AMPA current (IAMPA) and the inhibitory GABAA current (IGABA). Both are described by an equation of the form


where s is an activation variable that depends on the voltage of the presynaptic cell υpre, and Es is a constant. Since the simulated network is of limited size, each cell receives a current Iapp that mimics the exogenous excitatory background drive that would be present in vivo. It is assumed that each pyramidal cell excites some subset of interneurons, and each interneuron makes reciprocal inhibitory synapses onto a subset of pyramidal cells (see Fig. 1 right panel and Fig. 2b and 2c for the synaptic connectivity maps).

Fig. 1
The 5-s snapshot of spike rasters of 14 cells and their synaptic connectivity map. Cells #1–10 are pyramidal cells, and Cells #11–14 are interneurons. Red/blue color at (i, j)-entry implies the presence of an excitatory/inhibitory synaptic ...
Fig. 2
The average firing rates of selected 25 cells (a) and the synaptic connectivity maps that were used to generate spike trains during baseline (b) and general anesthesia (c). Cells #1–20 are pyramidal cells, and Cells #21–25 are interneurons. ...

Two networks were constructed as follows. The first network was simulated at only a baseline state, while the second network was simulated at both baseline and anesthesia states. To simulate the anesthesia-like state we made parametric changes that are consistent with the molecular targets of well-known drugs. For instance, the general anesthetic drug propofol acts by potentiating the GABAA synaptic current [2]. In our simulation, we accounted for this effect by making a threefold increase in the synaptic conductance g in the decay dynamics of spre) in equation (10) [7], [11].

Setup-1: The first simulation is a small network of 14 cells, with 10 pyramidal neurons (i.e., regular-spiking, or RS cells) and 4 interneurons (i.e., fast-spiking, or FS cells). The simulation length is 2 minutes, with sampling rate 1 kHz. The mean±SD firing rates of the pyramidal neurons and interneurons are 3.86±0.18 Hz and 20.31±0.41 Hz, respectively. The spike rasters and the synaptic connectivity map are shown in Fig. 1.

Setup-2: The second simulation is a medium network of 200 cells (180 RS cells and 20 FS cells). In order to imitate a realistic recording condition where only a small portion of population neurons can be accessed, we randomly selected 25 cells (i.e., 1/8 of the population) that consist of 20 RS cells and 5 FS cells. A total of 5-min data were generated, with sampling rate 1 kHz. Specifically, the averaged firing rates during baseline (anesthesia) are 3.3 (2.1) Hz and 15.5 (3.9) Hz for the RS and FS cells, respectively. Summary of neurons’ firing rates and synaptic connectivity maps from two states are shown in Fig. 2.

In modeling the spike train data in these simulation studies, the spikes were binned with 1-ms resolution. We selected six firing history temporal windows that consist of the spike counts in the past 1–3, 4–8, 9–15, 16–25, 26–40, 41–60 ms. As a general rule, we set the hyperprior parameters to a0 = 10−3, b0 = 10−3 for FS cells, and a0 = 10−2, b0 = 10−4 for RS cells. Because neuronal spiking activity was assumed to be conditionally independent, individual neurons were fit with separate PPGLMs, followed by a KS test. To create a functional map of inferred neuronal interactions, two cells are said to be interacting if there is a nonzero (in the statistically significant sense) spiking dependent coefficient (at any time lag) between a trigger cell and the target cell. To determine the “functional” connection being excitatory or inhibitory, we counted the majority of the nonzero coefficients at all six time lags—if the majority of the coefficients are greater than 0, or there are more positive connections than negative connections, then the trigger cell is concluded to have an excitatory connection to the the target cell. A similar rule holds for the inhibitory connection.

In the experiment setup-1, The estimation accuracy of the inferred connection map is 94.9%. Checking the estimation errors showed that most errors were induced by FP: a few weak connections between some pyramidal cell pairs were mistakenly identified.

In the experiment setup-2, we first investigated the baseline condition. We have tested the impact of the data length to the estimation accuracy. Using 1-min, 2-min, and 5-min recordings, it was found that the mis-identification (FP+FN) error rates were 16.5% (99/600), 14.7% (88/600), 12.5% (75/600), respectively (the cell self-interactions were excluded). Hence, increasing the length of the recordings improved the estimation accuracy, but there was still a fundamental bottleneck because of the overall sparsity of the spiking data and the limit of the statistical model. It was also found that most of (either FP or FN) errors occurred in the RS→FS connections. By examining the results, we suspected that it was probably due to the imbalance of the firing rates between the RS and FS cells, since a low-spiking RS cell would likely be mistakenly estimated with an inhibitory effect on the FS cell (which causes FP). The FN error might be either because of the lack of sufficient spiking data, or because of the insufficient detectability of the model (in terms of the window size and length) or the inaccurate assumption of statistical model. Next, the simulated anesthesia-like condition are investigated. Using the full 5-min recordings (Anesthesia-1), the obtained mis-identification error rate was 18.9% (113/600). Therefore, the reduction of ensemble firing rates increased the sparsity level of the spiking data and consequently caused an increase of the FP/FN error (even the KS plot still falls within the 95% CIs). To further investigate the impact of the firing rate on the estimation accuracy, we fixed the network connectivity and data length, and further reduced the cell firing rates (~20% reduction at the population level; Anesthesia-2). In this case, we observed an increased mis-identification error rate (23.3%, 140/600). See Table I for a summary of the results.

Summary of Setup-2 results based on 5-min simulated data.

B. Experimental Data

The spike train data were recorded from the barrel cortices of two rats during active behavior (maze running) and during the administration of isoflurane anesthesia. Experimental protocol details are referred to [16]. Recordings from two rats at two days were used here, each with about 25–30 min recording time during run or anesthesia. For consistency, we used 20-min recordings for each rat for both run and induction of anesthesia (excluding the first 5 min data) conditions. The depth of anesthesia was assessed by breathing rate and the hindpaw withdrawal reflex. When the depth of anesthesia was reached, the delivery of isoflurane was stopped. The same cells were tracked through the active behavior and anesthesia. Classification of cell types (RS vs. FS) were determined by the firing rate and the peak-to-trough width [16]. Snapshots of spike train recordings from one rat are shown in Fig. 3. Summary statistics of the experimental data are shown in Table II.

Fig. 3
Snapshots of 1-min raster plots of 13 RS cells recorded from barrel, cortex (Rat 1) during maze running (left) and induction of anesthesia (right).
Summary statistics of real spike train data.

The spike trains were binned with 2-ms resolution. Considering the large variability of firing rates and inter-spike intervals (ISIs) across all neurons, to account for low firing rates we used five firing history windows that consist of the spike counts in the past 1–100, 101–200, 201–500, 501–1000 ms. No extra sensory or behavior covariate was used for fitting the spike train data during active behavior or induction of anesthesia. The hyperprior parameters a0 = 10−3, b0 = 10−4 are chosen for PPGLM fitting (but the results are insensitive to these hyperprior values). Due to imbalance of firing rates and high irregularity of ISIs, none of PPGLM fit fell within the 95% CIs of the KS plots. This might also be due to the insufficiency of the covariates used in our PPGLM, given the potentially non-stationary nature of the spike dynamics (e.g., large chunks of burst suppression periods) within the 20-min timescale. Based on the similar methodology in simulations, we inferred the functional connectivity of the cell assemblies. Results are shown in Fig. 4. As seen, the neuronal interactions (mostly RS-RS connections) among cell assemblies reduced from the active behavior to the induction of anesthesia. Specifically in Rat 2, the dominant inhibitory effect of one FS cell (#6 in Fig. 4c) was suppressed by the isoflurane, which further induced some excitatory or inhibitory interactions among other cells (Fig. 4d). Note that based on our experimental data, the inferred functional connectivity of neuronal ensembles during active behavior was rather low in both animals (Fig. 4a,c). This might be due to the fact that many recorded neurons were indeed physically far apart (resulting in a decreasing chance of connectivity) because of the placement of electrodes at different layers of the barrel cortex.

Fig. 4
Inferred functional connectivity maps for two rats during active, behavior (a,c) and induction of anesthesia (b,d).

IV. Discussion and Future Work

In this paper, we have used a PPGLM to assess the neuronal interactions of cell assemblies during baseline and anesthesia. By first testing on the synthetic data generated by biophysical neuronal models, our statistical model shows reasonably good estimation results on the neuronal interactions among cell assemblies. In the simulation study, it is observed that the reduction of neuronal firing rates pose a challenge in accurately estimating the network functional connectivity. Therefore, improving the detection accuracy and robustness in the presence of sparse spiking data requires more research effort in future investigation. In the real data analysis, it is observed that the interactions of RS-RS and FS-RS cells decreased after the induction of isoflurane anesthesia. As a future goal, we will apply the methods developed here to more experimental spike train data recorded from animals under different anesthetic drugs. Our approach suggests a way to characterize the effects of anesthestic drugs on single neurons and neuronal ensembles.


Supported by NIH Grants DP1-OD003646.

Contributor Information

Zhe Chen, Department of Brain and Cognitive Sciences, MIT, Cambridge, MA 02139, USA. Neuroscience Statistics Research Lab, MGH/HMS, Boston, MA 02114.

Sujith Vijayan, Department of Brain and Cognitive Sciences, MIT, Cambridge, MA 02139, USA.

ShiNung Ching, Department of Brain and Cognitive Sciences, MIT, Cambridge, MA 02139, USA. Neuroscience Statistics Research Lab, MGH/HMS, Boston, MA 02114.

Greg Hale, Department of Brain and Cognitive Sciences, MIT, Cambridge, MA 02139, USA.

Francisco J. Flores, Department of Brain and Cognitive Sciences, MIT, Cambridge, MA 02139, USA. Neuroscience Statistics Research Lab, MGH/HMS, Boston, MA 02114.

Matthew A. Wilson, Department of Brain and Cognitive Sciences, MIT, Cambridge, MA 02139, USA.

Emery N. Brown, Department of Brain and Cognitive Sciences, MIT, Cambridge, MA 02139, USA. Neuroscience Statistics Research Lab, MGH/HMS, Boston, MA 02114.


1. Sentinel event alert: preventing, and managing the impact of, anesthesia awareness. Oakbrook Terrace, IL: The Joint Commission; 2004.
2. Bai D, Pennefather PS, MacDonald JF, Orser BA. The general anesthetic propofol slows deactivation and desensitization of GABAA receptors. J. Neurosci. 1999;vol. 19:10635–10646. [PubMed]
3. Brown EN, Lydic R, Schiff ND. General anesthesia, sleep, and coma. New Engl. J. Med. 2010;vol. 363:2638–2650. [PMC free article] [PubMed]
4. Chornoboy ES, Schramm LP, Karr AF. Maximum likelihood identification of neural point process systems. Biol. Cybern. 1988;vol. 59:265–275. [PubMed]
5. Chen Z, Kloosterman F, Wilson MA, Brown EN. Proc. IEEE ICASSP’10. Dallas, TX: 2010. Variational Bayesian inference for point process generalized linear models in neural spike trains analysis; pp. 2086–2089.
6. Chen Z, Putrino DF, Ghosh S, Barbieri R, Brown EN. Statistical inference for assessing neuronal interactions and functional connectivity with sparse spiking data. IEEE Trans. Neural Syst. Rehab. Engr. 2011;vol. 19(no. 2):121–135. [PMC free article] [PubMed]
7. Ching S, Cimenser A, Purdon PL, Brown EN, Kopell N. Thalamocortical model for a propofol-induced α-rhythm associated with loss of consciousness. Proc. Natl. Acad. Sci. USA. 2010;vol. 107:22665–22670. [PubMed]
8. Cimenser A, Purdon PL, Pierce ET, et al. Tracking brain states under general anesthesia by using global coherence analysis. Proc. Natl. Acad. Sci. USA. 2011 in press. [PubMed]
9. Cronin J. Mathematical Aspects of Hodgkin-Huxley Neural Theory. Cambridge Univ. Press; 1987.
10. Jaakkola TS, Jordan MI. Bayesian parameter estimation via variational methods. Statist. Comput. 2000;vol. 10:25–37.
11. McCarthy MM, Brown EN, Kopell N. Potential network mechanisms mediating electroencephalographic beta rhythm changes during propofol-induced paradoxical excitation. J. Neurosci. 2008;vol. 28:13488–13504. [PMC free article] [PubMed]
12. Okatan M, Wilson MA, Brown EN. Analyzing functional connectivity using a network likelihood model of ensemble neural spiking activity. Neural Computat. 2005;vol. 17:1927–1961. [PubMed]
13. Purdon PL, Pierce ET, Bonmassar G, et al. Simultaneous electroencephalography and functional magnetic resonance imaging of general anesthesia. Ann. NY Acad. Sci. 2009;vol. 1157:61–70. [PMC free article] [PubMed]
14. Stevenson IH, et al. Bayesian inference of functional connectivity and network structure from spikes. IEEE Trans. Neural Syst. Rehab. Engr. 2009;vol. 17:203–213. [PMC free article] [PubMed]
15. Truccolo W, Hochberg LR, Donoghue JP. Collective dynamics in human and monkey sensorimotor cortex: predicting single neuron spikes. Nat. Neurosci. 2010;vol. 13:105–111. [PMC free article] [PubMed]
16. Vijayan S, Hale GJ, Moore CI, Brown EN, Wilson MA. Activity in the barrel cortex during active behavior and sleep. J. Neurophysiol. 2010;vol. 103:2073–2084. [PubMed]