γ-Aminobutyric Acid Receptor Type A Receptor Potentiation Reduces Firing of Neuronal Assemblies in a Computational Cortical Model

• Several major theories emphasize the importance of neuronal groups, sets of strongly connected neurons that fire in a time-locked fashion, in all aspects of brain function, particularly as a necessary substrate of consciousness

• With the use of computer modeling, γ-aminobutyric acid receptor type A receptor facilitation was shown to directly and indirectly inhibit the ability of neurons to form groups spontaneously

• A lack of group formation is consistent with some theories of anesthetic-induced loss of memory formation and consciousness

THE exact mechanisms by which general anesthetic agents induce unconsciousness, amnesia, and immobility remain unclear. In particular, loss of consciousness is poorly understood, reflecting our rudimentary understanding of the phenomena of consciousness. Considerable effort has been directed at identifying the cellular effects of anesthetics, with neurotransmitter-gated ion channels emerging as the leading candidate targets.1Anesthetics have widespread anatomic sites of action. Activity in the precuneus cortex and temporo–parietal–occipital junction appears particularly important in anesthesia-induced unconsciousness, whereas controversy exists as to the role of the thalamus.2What remains far from being understood is how anesthetic effects on individual cells scale to the whole-brain level to abolish consciousness.

This mystery is part of a larger problem: the underpinnings of consciousness. Prominent theories for the neural correlates of primary consciousness3,–,9share the idea that assemblies or coalitions of neurons, the formation of which is facilitated by extensive reentrant connections between cells, are integral to production of the conscious state. In contrast, during seizures or burst suppression, stereotypical firing patterns reduce the repertoire of responses with a corresponding loss of information.2,9It has even been proposed that the degree of consciousness attainable by a neuronal structure may be estimated by an analysis of the possible interconnections of all possible subassemblies within the larger structure.10Although this conjecture has been questioned in detail,9the general picture of consciousness arising from the interaction of dynamic patterns of neuronal activity at multiple scales seems inescapable. Evidence for the presence of neuronal assemblies at the macroscopic level comes from magnetoencephalographic data demonstrating enhanced synchrony between distinct brain regions during visual conscious perception6and coalesced spatial activation of voltage-sensitive fluorescent dyes in brain slices,11,12whereas at the microscopic level, analysis of multielectrode pyramidal neuron recordings has revealed the presence of elementary neuronal clusters.13 

A range of techniques has been brought to bear on the study of the mechanisms of anesthesia. The recording of field potentials and patch clamping has been used to give detailed temporal information about neuronal electrophysiology, but these modalities are unable to provide information about the behavior of larger neuronal groups. Electroencephalography and the noninvasive imaging modalities of functional magnetic resonance imaging and positron emission tomography can provide a global picture of neuronal activity but are limited by poor spatial or temporal resolution.12It is for these reasons that computational modeling of neuronal processes has become an invaluable means of attempting to bridge our knowledge gap between the behavior of individual neurons and the global behavior of networks of millions of neurons.14,–,19 

To better understand how anesthetics may influence the formation and firing of neuronal assemblies, we have examined a computer model of a cortical neuronal network with connections modulated by synaptic plasticity. We have focused on one class of anesthetic agent, the γ-aminobutyric acid receptor type A (GABA^A) potentiators, as exemplified by propofol, because of their widespread use. In our model, the response to excitatory and inhibitory neurotransmitters is modeled by enhancing channel conductances in response to single spiking events followed by decay over time in an exponential fashion. This model spontaneously forms neuronal groups that fire in a timing-dependent manner. We assume that the activity of these groups is representative of the type of group firing that has been implicated in conscious brain activity in the theories of consciousness to which we have referred. We hypothesize that the potentiation of GABA^Areceptors in the model will act directly and indirectly via  input drive inhibition to diminish the ability of a neuronal network to generate cell assemblies by reducing group sizes and their frequency of firing and that this reduction will lead directly to a diminution of consciousness and its associated observable behavioral manifestations.

Anesthesiologists modulate consciousness on a day-to-day basis. A greater understanding of general anesthesia might lead to improved monitors of the level of consciousness, helping to reduce dose-related side effects and the incidence of intraoperative awareness.

Although this is not a detailed recreation of the cortical architecture, an attempt was made to incorporate the ratios and relative distances found in the mammalian cortex. The ratio of excitatory (E) to inhibitory (I) cells is 4:1 in the model.20To prevent boundary artifacts resulting from limited connectivity of neurons at the edge of the simulation field, the neurons are distributed randomly on a sphere of 8 mm radius (fig. 1). Each E neuron synapses on 150 neuronal compartments via  nonmyelinated local connections within a circle of 2 mm radius. Each I neuron has 160 synapses within a 1.5-mm radius via  nonmyelinated axonal collaterals. In addition, each E neuron has a randomly directed 10–15–mm myelinated axon that synapses on 50 neuronal compartments within 1.5 mm of the axonal termination.

The velocity of spike propagation is 0.15 m/s in the nonmyelinated local connections and 1 m/s in the myelinated axons.21 

The model comprises a network of 7,680 neurons in a single region, each with a soma and seven dendritic compartments (see Input section in Materials and Methods). There are 1,475,000 synapses. The spike timing dynamics of each compartment are simulated using Izhikevich's simple spiking neuron model,22in which neurons are modeled by the following equations:

where C  is the membrane capacitance, v  is the membrane potential (mV), v^r  is the resting potential, v^t  is the instantaneous threshold potential, u  is the recovery variable, c is the postspike membrane potential reset (mV), and I  represents the net dendritic and synaptic currents (pA). The parameters k, a , b , and d  are dimensionless and vary according to neuron type (table 1). Excitatory neuronal spiking patterns range from regular spiking to chattering type, with the regular spiking type predominant. Inhibitory neurons in the model exhibit a fast spiking response.

The neuronal equations are calculated at 1-ms time steps by a modification of the first order Euler method.23 

Each neuronal soma receives corticocortical synaptic input and current from dendritic compartments. In addition, each soma receives input drive from a random process that transiently depolarizes the membrane to a superthreshold level with a frequency of 1 Hz at baseline.

Synaptic transmission is modeled by increasing the conductance of the receptor channel in response to a spike event. The reversal potentials are 0 mV, 0 mV, −70 mV, and −90 mV for the α-amino-3-hydroxy-5-methylisoxazole-4-proprionic acid, N-  methyl-D-aspartate, GABA^A, and GABA^Breceptors, respectively,24,25and are incorporated into the following equation:

where g  denotes the conductance for each receptor type.26 

After synapse activation, the conductance at the channel decays in an exponential fashion over time according to the equation:

where τ denotes the decay time constant and is set at 5, 150, 6, and 150 ms for the α-amino-3-hydroxy-5-methylisoxazole-4-proprionic acid, N-  methyl-D-aspartate, GABA^A, and GABA^Breceptors, respectively, consistent with neurophysiological experimental data.27 

The ratio of α-amino-3-hydroxy-5-methylisoxazole-4-proprionic acid to N -methyl-D-aspartate synapses is set at 1 for the E neurons. Similarly, the ratio of GABA^Ato GABA^Bsynapses is set at 1 for all I neurons.

Currents coming from the upstream and downstream compartments determine the dendritic current for each compartment:

where G^down  and G^up  are conductances defined for each neuron type (see table 1) and v^up  and v^down  are the membrane potentials (mV) of individual upstream and downstream compartments, respectively.

At the start of simulation, excitatory synaptic weights c^E  are set within the range (0–1). These are then modified by spike-timing dependent plasticity (STDP; see Methods) and clipped to remain in the range (0–12). Inhibitory synaptic weights c^I  are assumed to be nonplastic in the model and are fixed at 4.

To model the short-term depression of synaptic weights that occurs with repeated firing of pyramidal cells and some interneurons, a phenomenological model matching in vitro  behavior was used,26in which the synaptic weight is modified by a per-synapse scalar factor x :

Thus, on each firing of an E neuron, the α-amino-3-hydroxy-5-methylisoxazole-4-proprionic acid and N -methyl-D-aspartate conductances of the postsynaptic cell are increased by xc^E . Similarly, on firing of an I neuron, conductances of GABA^Aand GABA^Bare increased by xc^I .

Long-term synaptic plasticity of excitatory synapses is modified by STDP.28Differential synaptic strength changes are accumulated after presynaptic events and applied only once per second. If a spike from a presynaptic neuron i  arrives at a postsynaptic neuron j  before it fires, the accumulated synaptic change is modified in the direction of potentiation according to:

where tⁱ  is the time of arrival of the spike from neuron i  and t^j  is the time that neuron j  spikes. Similarly, if neuron j  spikes before the arrival of the spike from neuron I , accumulated synaptic change is modulated in the direction of depression according to:

Every 1,000 ms of simulation time, excitatory synaptic weights, and their derivatives are updated by:

The 0.01 term in equation x is used to potentiate synapses coming to silent neurons. Changes of synaptic weights in the model are slow and take many minutes of simulation time to develop.

For specificity, we took propofol as a representative example of those anesthetic agents that are thought to act primarily by facilitation of GABA^Achannels. The effects of propofol in the network are modeled in several ways.

The potentiation of GABA^Areceptor channels by propofol at a range of plasma concentrations (0.5–30 μM) covering a spectrum of clinical endpoints from subclinical effects to sedation and general anesthesia is simulated by increasing the conductance of the channels by between 100% and 800% and the decay times to between 12 ms and 48 ms, consistent with in vitro  data29,30(table 2). Propofol also acts indirectly by reducing input drive from other regions as a consequence of GABA^Areceptor potentiation. This was investigated by reducing the frequency of the random input by a percentage (range, 6.5–60%) equal to the percentage reductions in E cell firing rates under each of the GABA^Areceptor potentiation conditions. GABA^Areceptor potentiation and reduced input drive frequency were studied separately and in combination.

It is unclear to what degree propofol affects synaptic plasticity at clinical doses. At very high propofol concentrations, induction of long-term potentiation may be diminished.31,32We studied the potential effects of this by reducing the STDP potentiation term described by equation 9to between 0 and 75% of baseline. No effects on long-term depression have been demonstrated, and such effects are not modeled in this study.

On examination of plots displaying occurrence of action potentials in all neurons over time, the firing appears random and uncorrelated. However, certain groups of neurons fire together repeatedly with persistent spike timing patterns with millisecond precision. We term these neuronal collections polychronous groups (from Izhikevich 2006), meaning multiple timing. These groups form by the interaction between STDP in response to random spiking and axonal delays.33 

To find a polychronous group (fig. 2), we identify triplets of E neurons that have strong (more than 95% maximum synaptic weight) connections to a fourth E neuron at a given time point. The timing pattern of stimuli required to simultaneously elicit a response from the fourth neuron is determined by examining the axonal conduction delays. A recursive process is then used to find additional neurons that belong to the group. Pairs of neurons already in the group are identified that strongly connect to another neuron outside the group with spikes that arrive at the same time (±1 ms) based on their axonal conductances and the overall spike timing pattern. I neurons are also added to the group, but outgoing connections from these are not considered to be part of the recursive process. The process continues until no additional neurons can be added to the group. Thus, each group is a prediction of a possible stereotypical firing pattern based on the anatomical data.

Only polychronous groups that persist throughout the duration of the simulation period are used for firing frequency analysis. A group is determined to be activated when at least 50% of its component E neurons fire within ±1 ms of the previously determined spike-firing pattern.33Using each identified group as a template, the spike-firing data are scanned to determine the number of times each group is half-activated. The distribution of group sizes is compared using the Kolmogorov-Smirnov test with statistical significance defined as P < 0.05.

A simulated electroencephalographic signal was generated from the model using a mean of the currents in the vertically aligned dendritic compartments of the E neurons. Analysis of electroencephalogram frequencies was performed with spectral density plots generated by the fast Fourier transform.

The program simulating the neural network is purpose-written for parallel processing in CUDA (Compute Unified Device Architecture) extension to the C programming language and run on a personal computer based on an Intel Core i7 with dual NVIDIA GTX 580 graphic processing units (NVIDIA, Santa Clara, CA) with the Linux operating system (Red Hat, Raleigh, NC). Graphing functions and statistical analysis are performed by Octave.‡

Initialization of the model takes approximately 2 min, and each second of simulation time takes 1 s of computing time. Each simulation is run for 6 h of simulation time, divided into three epochs as follows: The network is first run for 3 h of simulation time to allow stabilization of synaptic weights under the STDP rule (fig. 3). Data from this “run-in” epoch time are not included in subsequent analyses. The model is then run for a 2-h “baseline” epoch. Network connection data from this baseline epoch are analyzed for the presence of polychronous groups. Only groups that are present throughout are included in the firing frequency analysis. An additional hour of simulation time (“induction and maintenance” or “test” epoch) is performed with either control conditions or varying degrees of GABA^Areceptor potentiation and input frequency reduction or abolition of long-term potentiation. The frequency of occurrence of the polychronous groups identified from the 2-h baseline time period is determined for each simulated state during the “test” epoch. Group activation frequency is reported as the total number of groups activating per second averaged over the entire “test” epoch.

Initially, the cell activity in the model approximates the 1 Hz random excitatory input. As the run-in epoch progresses, synaptic weights are strengthened or depressed, depending on the relative timing of the random input spikes at connected neuronal compartments. Within the first 90 min of simulation time, the synaptic weights distribute in a bimodal manner. Most of the synaptic weights are weak, eliciting an excitatory postsynaptic potential of less than 0.1 mV. At steady state, just less than 20% of synaptic weights are within 5% of maximal strength, with a mean excitatory postsynaptic potential of approximately 5 mV. Once the synaptic weights stabilize, E cells fire at a mean of 4.2 Hz, and I cells fire at a mean of 20 Hz (fig. 4A). The proportionally higher firing rate of I cells relative to their reduced numbers leads to an approximate balance of excitation and inhibition in the model.

During the following 2 h of baseline simulation, the mean firing rates of cells and the frequency of the predominant rhythms remain stable. The simulated electroencephalogram shows the spontaneous development of low theta (4 Hz) and gamma (42 Hz) rhythmicity (fig. 4B). The appearance of these rhythms coincides with the development of a bimodal distribution of synaptic weights.

GABA^Afacilitation and input frequency reduction independently reduce the firing rate of cells in the network in a dose-dependent manner (figs. 5and 6). The combination of these two effects is approximately multiplicative (fig. 7).

Increasing degrees of GABA^Areceptor potentiation and input frequency reduction diminish the spectrographic amplitude of the electroencephalographic theta and gamma rhythms. The maximal reduction in spectral power by GABA^Areceptor potentiation alone is 75%, by input frequency reduction alone is 68%, and in combination is 94%. There is a modest reduction of the dominant gamma frequency from 42 to 39 Hz at the higher range of simulated propofol effects.

If the model is run at baseline without synaptic plasticity, no polychronous groups form. In the presence of STDP, analysis of the neuronal connection data identified more than 8 million polychronous groups. Many of these groups are transient, appearing and disappearing in response to ongoing synaptic plasticity. However, approximately 168,000 groups persist throughout the entire 2-h baseline simulation. These groups contain as many as 54 neurons, and their activity spans from 25 to 90 ms. The neurons connect via  a mixture of local nonmyelinated axonal collaterals and distant axonal synapses. There is considerable overlap of the groups, with each neuron having membership in a number of different collections.

A scan of the spike train data from the “test” epoch shows that the mean activation frequency of polychronous groups under control conditions is approximately 12 times per hour, with approximately 30% of the groups activating at least once per hour. Groups that are active at the start of the epoch are not necessarily the same groups that are active toward the end, reflecting ongoing dynamic changes in the strength of connections between the neurons. However, the total number of groups active in the first 100 s of the epoch was very close to the number active in the last 100 s, demonstrating that under baseline conditions, the network has achieved a form of equilibrium during this time period.

To ensure that group firing is not appearing by chance in randomly spiking neurons, we determined the frequency of group firing in a surrogate data group by inverting the times of the spike data. This procedure preserves many of the relevant statistics of the spiking data.33Although group activation is still detected in the surrogate data, the rates are much lower (data not shown), demonstrating that group activation is a statistically significant event.

The potentiation of GABA^Areceptors as described in Methods and table 2acts directly to reduce the mean activation frequency of polychronous groups (fig. 5). Similarly, GABA^Areceptor potentiation acting indirectly via  reduced input drive substantially reduces polychronous group activation (fig. 6). The direct and indirect effects combine multiplicatively to almost completely prevent group firing at levels consistent with high plasma propofol concentrations (fig. 7).

Reduction of group activation occurs immediately in the “test” epoch on the introduction of the direct or indirect effects of GABA^Areceptor potentiation. At lower levels of GABA^Areceptor potentiation and input drive reduction, group activation frequency is stable throughout the epoch. Under conditions consistent with high plasma propofol concentrations, group activation frequency increases by approximately 50% during the epoch, reflecting ongoing synaptic plasticity that acts to partially normalize the perturbed group firing function of the network.

Of the groups that did activate under conditions of GABA^Areceptor potentiation or input drive reduction, their distribution is skewed toward those containing smaller numbers of neurons (fig. 8) (P < 0.02). In addition, there are many fewer occasions in which multiple groups fire concurrently to form larger assemblies (fig. 9).

Polychronous group activation is relatively insensitive to inhibition of synaptic plasticity via  reduction of long-term potentiation. At STDP potentiation levels as low as 5% of baseline, dominant rhythms are maintained and group firing rates reduce by no more than 7%. In the absence of potentiation consistent with very high plasma propofol concentrations, the ability of the network to maintain and activate polychronous groups rapidly declines.

This study used a computational model of cortical neurons incorporating important neurophysiological details, including synaptic plasticity, excitatory and inhibitory neurotransmitters, and accurate phenomenological descriptions of spiking activity. The action of GABA^A-acting anesthetics such as propofol was simulated by increasing the conductance and decay times of GABA^Achannels and by reducing the input drive frequency. At baseline, the model exhibited firing of neuronal groups in precise spatiotemporal patterns known as polychronous groups. We have shown that the simulated effects of propofol, separately and in combination, substantially disrupted the ability of the cortical network to display polychronous group activity. This was reflected by reductions in group activation frequency, size, and the formation of larger collections of groups. These results do not imply that existing groups are necessarily destroyed by inhibition of group activity; their activity may resume when GABA^Ainhibition levels normalize.

It has been hypothesized that neuronal groups constitute working memories33or building blocks of perception.13These groups may assemble7or undergo selection,6giving rise to repertoires of such complexity that they reflect the richness and diversity of conscious experience.

On a much smaller scale in this model, a far greater number of groups arose than the number of neurons. In this manner, the network is able to represent many more states than one would expect on the basis of its limited neuronal number. The input to the network is random, so the groups that form are also random, representing “memories” of arbitrary external experiences.

The reduction of group activation frequency for combined GABA^Areceptor potentiation and input drive reduction at levels consistent with sedation and general anesthesia is about 70% and 90%, respectively. Although there are no data to indicate whether these levels of inhibition of group formation and activity could explain the abolition of new memory formation and consciousness, it certainly seems reasonable that the observed almost 10-fold reduction in group activity in the model would reduce dramatically the information-carrying capacity of the network and, in some theoretical frameworks,6,7could account for loss of consciousness. Although removal of long-term potentiation alone is sufficient to prevent polychronous group activation, it may occur only at very high plasma concentrations of propofol31,32and was not required in the model to account for the observed lack of group activation.

Currently available techniques to study anesthesia-induced unconsciousness in humans are hampered by a lack of temporal resolution (functional magnetic resonance imaging) or spatial resolution (positron emission tomography and electroencephalography). Nevertheless, important insights consistent with our findings have been provided by these modalities. Functional magnetic resonance imaging studies have demonstrated that propofol diminishes total integration of the brain34and may reduce the number of states that can be discriminated.35Studies of multichannel electroencephalographic recordings36,37using analyses derived from graph theory found that propofol reduces network connections and diminishes the entropy of the connection matrix while preserving global network organization.

Recent animal data from multiple cortical implanted electrodes have shown that desflurane suppresses long latency flash responses, suggesting a deintegration of cortical processing.38Voltage-sensitive dye studies of cortical slices exposed to intravenous and volatile anesthetics suggest a prolongation of neuronal activation after stimulation, perhaps reducing the brain's ability to respond rapidly to incoming sensory information by forming new neuronal collections.11,12 

Although the random input into the model has a frequency of 1 Hz, the simulated electroencephalogram displayed prominent low-frequency activity (4–6 Hz) and high-frequency activity (42 Hz) that developed within 100 min of simulation time. These rhythms are not built into the network but develop in response to synaptic potentiation and depression under the STDP rule. Similar emergent frequencies have been identified in an anatomically detailed model of the neocortex containing modified integrate-and-fire neurons.39 

The GABA^Areceptor potentiation and input frequency reduction in our model reduced the amplitude of these emergent rhythms and was accompanied by a modest slowing of the gamma rhythm. Ching et al.  19recently reported modeling of propofol-induced alpha rhythm in thalamocortical neurons, arising from a slowing of the gamma rhythm that predominates at baseline. That the same magnitude of gamma slowing did not occur in our model perhaps emphasizes the importance of the thalamus in this phenomenon and may reflect differences in the magnitude of the neuronal input current between the two models. Certainly, regional differences exist in the electroencephalographic response to propofol anesthesia,36and loss or slowing of gamma may not be a universal finding.40 

We have not modeled all of the potential cellular effects of propofol. Tonic GABA^Acurrents are likely enhanced, particularly at the higher end of the dose spectrum.41Propofol may also cause inhibition of the glutamate release machinery.42It is likely that the cumulative effect of these potential mechanisms will lead to additional impairment of cell firing and group formation in the model.

In this study, we have not examined the effects of thalamic activity or modeled area-specific specializations of the cortex. As discussed, many electroencephalographic and imaging studies34,37,43,44have shown significant regional variations in the effects of anesthesia. To fully model the phenomenon of anesthesia-induced unconsciousness, it will be important to include such anatomic detail. Nevertheless, what the model shows is a way that the known direct and indirect GABA^Areceptor effects of propofol and similar agents can lead to a disruption of the ability of local cortical circuits to form polychronous groups. Additional work will help to elucidate how these local phenomena translate to the whole brain scale.

The approach that we use to identify polychronous groups requires knowledge of axonal delays and synaptic weights for all neurons in the network. For this reason, in vivo  recording of neuronal membrane potentials would not lend itself to such an analysis. In addition, the number of concurrent neuron recordings required to have a reasonable chance of detecting a polychronous group is beyond the capabilities of currently available multielectrode arrays.21 

The modeling approach described in this article could be easily adapted to study other agents that act predominantly on GABA^Achannels, in particular etomidate and the benzodiazepines. Although volatile anesthetics may have similar effects on neuronal group formation to those seen here for the predominantly GABA^Apotentiating anesthetics, modeling is complicated by a more diverse range of targets, including ion channels not explicitly modeled in our current formulation. In subsequent modeling, we plan to study a larger scale model of the thalamocortical system, closer to that described by Izhikevich and Edelman,26that recreates anatomic details, including global thalamocortical white matter anatomy, thalamic nuclei, cortical microcircuitry, and a diverse range of neuronal subtypes.

We have studied a cortical network with synaptic plasticity that spontaneously self-organizes into groups that repeatedly fire in precise spatiotemporal patterns. The direct and indirect phasic GABA^Aeffects of propofol substantially disrupt the ability of the network to form these neuronal assemblies. These findings are consistent with theories of memory and consciousness that emphasize the importance of neuronal group formation.

The authors thank Hugh Hemmings, M.D., Ph.D. (Professor in Anesthesiology, Department of Anesthesiology, Weill Cornell Medical College, New York, New York), and Kane Pryor, M.D. (Assistant Professor, Department of Anesthesiology, Weill Cornell Medical College), for reviewing the manuscript.