Emergence of Physiological Oscillation Frequencies in a Computer Model of Neocortex

Introduction

Coordination of neocortical oscillations has been hypothesized to underlie the “binding” essential to cognitive function (Treisman, 1996; Engel and Singer, 2001; Singer, 2001; Wu et al., 2007; Fenton et al., 2010; Kelemen and Fenton, 2010). There is also evidence that neocortical oscillations contribute to working memory, and that alterations of brain oscillations are associated with diseases producing psychiatric disease (Spencer et al., 2003, 2004; Raghavachari et al., 2006; Uhlhaas and Singer, 2006; Uhlhaas et al., 2006; Düzel et al., 2010; Lisman, 2010). Recent experimental evidence suggests that human neocortex can generate theta locally (Raghavachari et al., 2006; Tsujimoto et al., 2006). At present, the mechanisms that generate neocortical oscillations in physiological frequency bands remain unknown, likely involving both intrinsic cortical properties and projections from other centers (Bazhenov et al., 2002). Understanding is complicated by the fact that neocortex is an area of great circuit complexity, which receives projections from multiple subcortical areas, including hippocampus, a potent generator of theta activity (Stewart and Fox, 1990; Hasselmo, 2006; Sirota et al., 2008; Cutsuridis and Hasselmo, 2010).

It is widely believed that neocortical circuits are arranged in columns and that these columns may be elemental information processing units (Mountcastle, 1997; Petersen and Sakmann, 2001). Although the columnar hypothesis remains controversial and lacks anatomical verification (Buxhoeveden and Casanova, 2002; Rockland and Ichinohe, 2004; Costa and Martin, 2010), it is at worst a convenient approximation. We therefore used a columnar organization in these studies. The columnar hypothesis suggests that individual columns can process and store information internally, and then send and receive information amongst themselves. Both intercolumnar and interlaminar connectivity provide a complex of circuitry with multiple levels of embedded and interlaced recurrent connections that may play complementary roles in information processing (Douglas et al., 1989; Nelson, 2002; Hill and Tononi, 2005; Sporns et al., 2005; Izhikevich and Edelman, 2008). These many levels of recurrence add to the difficulty in unraveling the dynamical contributions of layering and of cortical and regional organization.

The genesis of the electroencephalographic signal, first demonstrated over 80 years ago, remains little understood. We performed simulations of a group of nine neocortical columns from sensory cortex in order to determine the cellular, laminar and columnar relations that produce these characteristic oscillations. Our central hypothesis was that the characteristic spectral pattern would emerge from interrelations of the elements at multiple scales: cell types with differing intrinsic dynamics, connections via certain delays, complex topology involving interactions both within and between multiple layers, local intracolumnar and more remote intercolumnar connections. Due to the multiple nested and interlocking feedback loops, straightforward definition of a single oscillatory origin is impossible. In order to assess dynamical credit-assignment, we used minimal dynamical perturbation as a methodology for selectively altering individual ingredients in the dynamical system without pushing the system into a substantially different dynamical pattern. We minimally altered the neocortical topology by incrementally adding one hub (a cell with many more inputs/outputs than average) at a time to selected layers. This allowed us to demonstrate that layer 2/3 was critical for producing the specific spectral pattern characteristic of neocortex. Spectral power from our model’s population activity was found to be similar to spectral power from local field potential recordings from rat medial prefrontal cortex, both in profile and cross-frequency coupling.

Materials and Methods

Simulations

Simulations were run using NEURON (Hines and Carnevale, 2001; Carnevale and Hines, 2006). The cell model was an extension of an integrate-and-fire unit and was simulated in an event-driven fashion, since cell state variables were only calculated at input events. Complexity was added using rules with various dynamical features found in real neurons: adaptation, bursting, depolarization blockade, and voltage-sensitive NMDA conductance (Lytton and Stewart, 2005, 2006; Lytton and Omurtag, 2007; Lytton et al., 2008b). Briefly, each cell had a membrane voltage level state variable (V_m) and a resting membrane potential (V_RMP), which was the baseline value for V_m in the absence of external inputs. After synaptic input events, if V_m crossed spiking threshold (V_TH), the cell would emit an action potential and enter an absolute refractory period, lasting τ_refrac ms, during which it could not fire. Refractory periods were set to prevent a maximum firing frequency from being exceeded. After an action potential, an after-hyperpolarization voltage state variable was set (V_ahp) and then subtracted from V_m. Then V_ahp decayed exponentially with time-constant (τ_ahp) to 0. To simulate voltage blockade a cell could not fire if V_m surpassed V_block. Relative-refractory period was simulated after an action potential by increasing the threshold, V_TH, by W_RR. (V_block − V_TH) (to scale the increase in V_TH by a fraction of the voltage range where the cell could fire), where W_RR was a unitless weight parameter. V_TH then decayed exponentially to its baseline value with time-constant τ_RR.

Baseline wiring and number of cells per layer were based on published models and anatomical studies, which necessarily come from different species and different areas and are sometimes inferred rather than explicitly stated. We used data from cat visual cortex (Thomson et al., 2002; Binzegger et al., 2004), rat visual cortex (Song et al., 2005), and rat barrel cortex (Lefort et al., 2009) for intracolumnar connectivity information. We used mouse somatosensory cortex data for intercolumnar projections, which showed lateral inhibition as well as excitation between layers 2/3 and 5 of neighboring columns (Adesnik and Scanziani, 2010). In addition, we used data from previously published neocortical column models to ensure that our model was consistent with data used by others (Traub et al., 2005; Heinzle et al., 2007). Unfortunately, exact numbers are not available for all connections. Instead estimates are available from a number of techniques including light and electron microscopy, dual to quadruple intracellular recordings, laser scanning photostimulation, and others. Often these methods do not entirely agree, partially due to the difference in techniques but also due to the studies being done in different species and in different neocortical areas. Due to incomplete and conflicting information, a number of choices had to be made based on judgment. This approach of integrating multiple data sources and using best approximations to real parameters has generally been used in other large-scale neocortical models (Traub et al., 2005; Izhikevich and Edelman, 2008).

In our model, there were 13 cell types, arranged in 4 layers: 2/3, 4, 5, and 6. Interneurons were parameterized as fast-spiking (FS) or low-threshold spiking (LTS) interneurons. FS and LTS interneurons utilized somatic and dendritic GABA_A synapses, respectively. The following cell types were included: E2, I2, I2L, E4, I4, I4L, E5a, E5b, I5, I5L, E6, I6, and I6L. E (I) represent excitatory (inhibitory) cells, and the number following the cell type represents the layer, i.e., E2 represents pyramidal cells in layers 2/3. I6 represents FS interneurons in layer 6, and I6L represents LTS interneurons in layer 6. E5a and E5b are two subtypes of pyramidal neurons present in layer 5, with different connectivity patterns (Steriade, 2004; Brown and Hestrin, 2009; Groh et al., 2010). Parameters and counts for the cell types are listed in Table 1.

Intracolumnar synapses utilized our previously developed just-in-time synapses, optimized for large networks supporting high-frequency synaptic events (Lytton et al., 2008a). Intercolumnar synapses used NEURON’s NetCon synapses (Lytton, 1996), included to permit future parallelization of single columns per node via the message passing interface (MPI). Each cell had a voltage state variable associated with a synapse type, V_syn, one for each of excitatory AMPA, NMDA, and two inhibitory GABA_A s, which simulated GABA_A at soma (fast time-constant) and GABA_A at dendrite (slower time-constant). Synaptic inputs were simulated by step-wise changes in V_syn and then added to the cell’s V_m level. To allow for dependence on V_m, synaptic inputs changed V_syn by δ_V = w_syn.(1 − V_m/E_syn), where w_syn is the synaptic weight and E_syn is the reversal potential, relative to V_RMP. E_syn took the following values (in mV): AMPA 65, NMDA 90, GABA_A −15. w_syn was positive for excitatory synapses and negative for inhibitory synapses. NMDA synapses also had an additional voltage-dependent scaling factor based on Jahr and Stevens(1990a,b). For all synapses, after synaptic input events, V_syn decayed exponentially toward 0 with time-constant τ_syn. τ_syn took the following values (in milliseconds): AMPA 20, NMDA 300, somatic GABA_A 10, dendritic GABA_A 20. Synaptic weights were constant between a given set of populations. Dendritic synapses (AMPA, NMDA, dendritic GABA_A) utilized delays chosen from a uniform distribution ranging from 3 to 5 ms, while somatic synapses (somatic GABA_A) had delays ranging from 1.8 to 2.2 ms.

The network consisted of nine cortical columns (from sensory neocortex) with each column connected to all other columns. We refer to simulations with interconnected columns as multicolumnar simulations. Each of the nine cortical columns had 470 cells, for a total of 4230 cells. The number of synapses for a particular cell was chosen from a normal distribution using average connection densities between different classes of cells. Intra- and intercolumnar connectivity information is shown in Tables 2 and 3, respectively. Each cell’s spatial location is determined by its column and cell type. Together, these determine the cell’s probability of connecting to other cells.

Biological connection densities across and within cell type groupings are believed to cover the wide range of 1–50% in different areas (Thomson et al., 2002; Song et al., 2005; Orman et al., 2008). A typical number of synapses for a baseline simulation was 831,967 out of 17,888,670 (full number of possible connections: n² − n for a fully connected network without self-connections), for an overall density of 4.65%, consistent with densities found in the literature. Within a column, the average number of synapses was 8390, 5767, 7449, and 3810 for E → E (excitatory → excitatory), E → I (excitatory → inhibitory), I → E, and I → I synapses, respectively. Between columns, the average number of synapses was 1974 and 1331, for E → E and E → I connections, respectively.

Descriptions of the dynamics emerging from structural layouts of neural tissue must be seen in the context of the vertical wiring layout between different layers and the horizontal layout across areas and columns. Although both wiring and component specifications remain incompletely understood, there is growing knowledge on which we base our model, depicted in Figure 1 (Douglas et al., 1989; Thomson and Bannister, 2003; Binzegger et al., 2004; Douglas and Martin, 2004, 2007a,b; Markram et al., 2004; Traub et al., 2005; Thomson and Lamy, 2007; Izhikevich and Edelman, 2008; Weiler et al., 2008; Brown and Hestrin, 2009; Lefort et al., 2009; Adesnik and Scanziani, 2010; Anderson et al., 2010; Groh et al., 2010; Neymotin et al., 2010).

Figure 1A shows the cell types in the various layers, with connectivity depicted by red arrows for excitatory and blue for inhibitory connections. Intergroup connectivity density within the column is moderate (43%; 72 out of 169 possible connections for the 13 cell groups). Total E → E connections were 88%; E → I connections 25%; I → E connections 45%; I → I connections 34%. All populations have self-connections (note that these are not connections from a cell to itself but from cells in a group to different cells in the same group).

Figure 1B shows a graph-theoretic visualization of the intercolumnar group-to-group wiring between neocortical columns. There are again 169 possible connections, but the connection probability is 7% (12/169), far lower than within column. Although all projection pathways are excitatory, most project onto inhibitory cells and thereby produce feed-forward inhibition, largely laterally to the same layer across columns. An additional major lateral inhibitory projection goes from E5b → I2L. E5, the main output layer, is also the source of intercolumnar feed-forward excitatory connections, including E5b → E2 and connections from both E cell groups in layer 5.

Figure 1C shows a graph-theoretic schematic of intracolumnar connectivity. Each population is arranged in a circular layout in order to minimize the number of edge crossings (Six and Tollis, 1999a,b; Gansner and North, 2000). This has the effect of placing the nodes that are more densely connected closer to each other. The size of the circles shows the number of cells in that population and the width of the arrows indicates the strength of a projection (consisting of the product of the connection density and the average synaptic weights). The red circles are larger than the blue circles, indicating that the majority of neurons are excitatory. E2, the circle representing layer 2/3 pyramidal cells, is the largest, foreshadowing its large role in determining network dynamics. Additionally, the E2 → E2 connecting edge is very thick, indicating the strong recurrent connectivity within E2. E2 also strongly feeds into the I2 inhibitory population and gets strong feedback inhibition from there. E2 activity is thereby modulated through strong self-coupling via both excitation and inhibition.

The thick edge from E2 → E5a reflects the strength of this major intracolumnar, inter-layer excitatory pathway. E5a also has relatively strong recurrence, and feeds strongly onto E5b, reinforcing output layer 5’s activity. In most of the populations, inhibitory feedback loops are prominent, but with the E → I projection generally weaker than the I → E back-projection. We suspect that the stronger I → E might reflect the dual role of inhibitory cells as part of both feedback and feed-forward inhibition.

The graph visualization of network structure suggested that E2, I2, and E5 would have particularly large influence on the dynamics of the system. We tested this in the simulations.

Development of Hubbed Networks

To increase the efficacy of a chosen subpopulation, we added hubs to that subpopulation. Hubs are individual cells with a larger number of inputs and outputs than the average cell of the given type. Wiring changes were performed by first randomly selecting 1–10 cells from a group to be hubs. To preserve the overall number of synapses in the network, we randomly moved existing synapses from other cells and set the new hubs as their target or as their source. We avoided unconnected cells, holes, in the network by preventing wiring changes that resulted in a cell having no inputs, or no outputs. We set the hub cells to have 3× the post-rewiring average convergence and divergence of the non-hub cells of the same type. For excitatory cell hubs, the number of inputs and outputs to and from interneurons was not altered. This allowed these hubs to control excitation in the network, without changing inhibition. For inhibitory cell hubs, we increased the number of excitatory inputs from and inhibitory outputs to excitatory cells, allowing this type of hub to more effectively inhibit.

Data Analysis

Simulation spike trains were organized into multiunit activity (MUA) vectors, defined for a cell population within a column, as the number of spikes in the population over a time interval (bin). Bin sizes were set to 5 ms (200 Hz sampling rate). We evaluated bin sizes ranging from 0.1 to 5 ms and found similar results. We also generated simulated local field potential (sLFP) for easier visual comparison with experimental data. The sLFP is a summed and filtered result of synaptic activity that produces dipoles. With our limited, “point-neuron” model, sLFPs were approximated by summing excitatory-inhibitory synaptic input differences across excitatory cells and adding a small component of spike activity. The following expression was used for the synaptic contribution from each excitatory cell to the sLFP at each time step: V_cell = V_AMPA + V_NMDA − V_GABAAsoma−V_{GABAA dend}. Organization of simulations and analysis of output data was done with Neural Query System (Lytton, 2006; Lytton and Stewart, 2007). Spectral analyses used mean-subtracted vectors and were performed using Matlab’s implementation of the multitaper spectral power estimation method and FFT (pmtm() and FFT() functions; The MathWorks, Natick, MA, USA 2009b). The resulting power spectra were smoothed using a moving average filter for easier visualization and comparison of multiple power spectra. Similar peaks were visible in the unsmoothed spectra as well.

To measure regularity of individual spike trains, we used the local variation of interspike intervals method, defined as:

(Shinomoto et al., 2005), where T_i represents the i^th interspike interval, and n is the number of intervals. L_v is 0 for interspike intervals formed from regular spike trains with a fixed interval between spikes and L_v is 1 for irregular spike trains with interspike intervals that are independently exponentially distributed, such as those generated by a Poisson process. Higher values tend to suggest clustered or burst firing.

Electrophysiology

All procedures received IACUC approval and were in accordance with institutional and NIH guidelines. Recordings were part of an experiment on the effect of hippocampal lesions using injections of ibotenic acid. The single recording used for this paper was from a control animal that had a saline injection. Briefly, timed-pregnant Long-Evans rats were obtained at 13–14 days of gestation from Charles River Laboratories (Wilmington, MA). On postnatal day 7 (P7), pups weight 14–18 g were anesthetized by hypothermia. A puncture hole was made bilaterally on the skull for injection sites (relative to bregma AP: 3.0 mm, ML: ±3.5 mm), and 0.3 μl of saline was infused bilaterally into the ventral hippocampal formation (relative to the skull surface DV: 5.0 mm) at 0.15 μl/min. The pups were warmed and then returned to their mothers. Since only saline was injected, this animal was considered a control animal which had no significant brain lesion. On P21, animals were weaned. At P60, the rats were implanted and recorded for local field potentials (LFP, procedure adapted from Olypher et al., 2006). The LFP recordings used in this paper were made bilaterally from the medial prefrontal cortex AP: +3.0 mm, ML: ±1.0 mm, DV: 4.0 mm (only left side recording used). Electrodes were referred to a reference electrode implanted in the cerebellar white matter. The rat was recorded in its homecage for 30 min. The signal was digitized (24-bits, 48 kHz; Axona Ltd, St. Albans, U.K.) and transmitted galvanically along a counter-balanced cable. The digital signals were received by a set of dedicated digital signal processors, bandpass filtered (1–500 Hz), amplified digitally, and then down-sampled (2000 Hz). Comparisons with simulation were done using normalized spectral data since simulation spectra were based on MUAs of arbitrary amplitude.

Results

This study involved over 1200 20 s simulations, using five different random wirings and five different input streams, as well as many briefer simulations used during tuning. Results shown were consistent across the variations in input streams and random wirings. Each simulation had nine columns, connected in some cases and not in others. Simulations were run on Linux on a 2.7 GHz quad-core Intel XEON CPU. A 20-s simulation ran in 2–4 min, depending on the simulation type. The full model is available on ModelDB (https://senselab.med.yale.edu/ModelDB/ShowModel.asp?model=138379).

Emergent Rhythms

To maintain activity throughout the simulation, each synapse was independently provided with background inputs: sub-threshold stimulation from external, independent, Poisson-random inputs at rates of 240–360, 40–60, and 100–150 Hz, for AMPA, NMDA, and GABA_A synapses respectively. The frequency spectra of the interspike intervals of the background inputs were relatively flat, indicating white-noise spectra. Interspike interval local variation (ISI L_v, see Materials and Methods) of each input was near 1, further verifying that intervals between inputs were exponentially distributed (Shinomoto et al., 2005). In addition to external inputs, synapses received inputs from other cells in the model. Initially, with the cells of the network not connected, the weights of the background inputs were adjusted to provide low levels of activity. The cells were then connected up and the networks hand-tuned through adjustment of E → E, E → I, I → E, and I → I gains till physiologically realistic firing rates were obtained: E cells firing 0.5–3 Hz and I cells 1.5–10 Hz (Table 4). In the resulting network, all excitatory postsynaptic potentials remained sub-threshold, hence cell firing would only occur through integration of several inputs.

A network with strong E → E connectivity is at risk of runaway excitation that will produce epileptiform activity, either immediately self-limiting (interictal spike), or continuing (seizure). In order to avoid these pathologies, we did further E/I balancing so that there were no multilayer population spikes seen in the raster or simulated local field potentials (sLFPs), as shown in Figure 2. In addition, inhibitory weights onto E cells were set so as to prevent oversynchronization and stereotyped firing. In this way, we avoided the positive feedback latch-up seen with increased synaptic strength and increased excitatory firing rates (Lytton and Omurtag, 2007).

The network produced a complex repertoire of activity patterns, including periods of low and high synchrony between neurons. In Figure 2, this alteration can be appreciated both in the raster, with its diffuse vertical stripes alternating with periods of asynchrony, and in the amplitude alterations of the sLFP. Additionally, individual cells fired with irregular patterns (Figure 2C). Note that the E2 pyramidal cell, like other excitatory cells, fired at a significantly lower rate than the inhibitory cells. We quantified the irregularity of single cell firing by using the local variation (L_v) measure used physiologically (Shinomoto et al., 2005). The background spike driving, being Poisson, is very tightly clustered around L_v = 1. Activity of the cells had a lower mean (L_v = 0.9; differs p < 0.001 by Student’s T-test) indicating more regular firing on average, and a much broader range of values (SD 0.32 compared to 0.01). This indicates that some cells were showing regular firing (L_v ≪ 1) while others were fairly bursty (L_v ≫ 1), despite there being no intrinsic bursting mechanism for the cells in this network. Increasing the strength of E → E coupling led to increased regularity of firing, with L_v decreasing from 0.9 to 0.7 (50% E → E weight increase) and to 0.3 (75% increase). This greater single cell firing regularity was associated with oversynchrony and fairly regular recurrence of population spikes. A somewhat more subtle form of synchronization pathology was the occurrence of synchronization across columns. A slight correlation between the sLFPs from different columns can be see in Figure 2B. We tuned intercolumnar excitatory connection strength to bring baseline correlations down to ∼0.09 (Pearson correlation of sLFPs). Note that there was no correlation when columns were unconnected, since each was driven with independent noise.

We tuned the networks using only two criteria, realistic firing rates and absence of pathological activity. When connecting the cells, all networks showed emergent activity with MUA frequency peaks close to physiologically relevant frequency bands (Figure 3). We used MUA rather than sLFP for spectral analysis since MUA is a representation of population output while sLFP is a representation of input to a population. MUA was formed by binning spikes with bin size 5 ms (200 Hz sampling rate); further decreasing bin size incrementally down to 0.1 ms did not change the spectra. Simulations with disconnected cells showed a fairly flat spectrum inherited from the white-noise driving (Figure 3A black traces). Strong theta/alpha (our model is not species specific, and one man’s alpha is another rat’s theta), emerged in the excitatory cells of the connected network with a prominent peak in E cell MUA near 7 Hz (Figure 3A, left, red). I cell MUA of connected cells showed a broad high gamma peak between 40 and 80 Hz (Figure 3A, right, red), consistent with evidence that interneurons generate gamma oscillations physiologically (Lytton and Sejnowski, 1991; Traub et al., 1996; Wang and Buzsaki, 1996; Tamás et al., 2000; Bartos et al., 2007; Tukker et al., 2007; Cardin et al., 2009; Sohal et al., 2009). This gamma peak was the result of selective attenuation of low-frequency oscillations due to synchrony within the I cell populations produced by the faster-firing I cells. The I2 population reflected the 7 Hz E2 peak when looked at in isolation (not shown), but this peak was lower than that of the E2 population. Both the E and I peaks demonstrated that individual subnetworks filtered out some frequency bands while augmenting or resonating to others. Note that all the emergent frequencies were considerably higher than the firing frequencies of the constituent cells. This demonstrated that these frequencies were produced by coordination within and across populations. Different cells fired on different cycles of the dominant frequencies, consistent with physiological data (Fries et al., 2007).

In Figure 3A we showed the emergence of specific frequencies when cells were connected within the column. In Figure 3B we show the alteration in spectrum when the columns were connected to one another. The theta peak in E cell MUA near ∼7 Hz was dampened (Figure 3B left, blue). Dampening was also seen in the gamma band out to about 35 Hz. By contrast, I cell MUA was augmented when columns were connected consistent with the strong intercolumnar excitatory projections onto inhibitory cells (Figure 3B right, blue). This then provides the feed-forward inhibition that reduced the power in the E cell spectrum.

Rewiring to Determine Dynamical Causality

Dynamical system evaluation offers a chicken and egg problem: given the existence of multiple loops, and loops within loops, there is no clear way to elucidate who drives whom. Those who drive are driven in turn. If we remove a subpopulation from a network (an ablation experiment) in an effort to determine whether this subpopulation was the major driver, we are left with a radically different network, that will show very different dynamics.

Our approach to this conundrum was to use a minimal dynamical perturbation: gently alter the network and determine how dynamics gradually changed. We slowly increased the efficacy of particular subpopulations through the addition of intracolumnar hubs. In a directed network (directed graph) such as this one, hubs are nodes that both receive projections from, and in turn project to, a large number of other nodes (as in an airline hub). In this case, the hubs were individual cells with 3× the post-rewiring average of convergence and divergence of non-hub cells of the same type. Connections remained individually sub-threshold, but large convergence means that these cells fired more. To preserve the overall number of synapses in the network (connection density), we randomly moved existing synapses from other cells and set the new hubs as their target or as their source. We avoided unconnected cells by preventing changes that would result in a cell having 0 inputs or 0 outputs.

Using individual (unconnected) columns, we increased the number of hubs in each excitatory population independently and calculated the E MUA spectral power from the remaining subset of cells. Power greatly increased above baseline with addition of hubs to E2 but not with addition of hubs to other populations (Figure 4). This increase was evident both in the E MUA power and the I MUA power (not shown), across frequencies. E2 cell parameters did not differ from those of other E cell types, so the change was due to the position of the E2 population in the overall network. The increase was particularly prominent in the theta/alpha and low beta bands (note different y-axis scales). Detectable increases in oscillatory power took place with as few as 5 of the 150 E2 cells being wired as hubs.

Local Structural Changes Produced Dynamical Changes

In the 10-E2-hubs columns (multi-hub columns), individual E2 hubs had average excitatory convergence of 95 and divergence of 91, approximately 3× the 31 average for the non-hub cells in the network. E2 hub cells were activated by more cells so showed average firing rates of 2.68, compared to 1.01 Hz for non-hub E2 cells. The E2 hubs also increased the average firing rates of other populations (in Hz): overall E2 1.124; I2 2.988; I2L 4.83; E4 0.804; E5a 0.914; E5b 0.659; I5 1.46; I5L 3.48 (compare to Table 4). The increased rate of non-hub cells was due to the excitatory positive feedback resulting from enhanced bidirectional connectivity between the hubs and their targets/sources. As a result, when a hub cell fired, a much larger number of cells in the network had their probability of firing increased after synaptic delay. Since the hub received a higher number of feedback projections from its targets than the average cell, it was also much more likely to fire within a short interval of its targets firing, potentially receiving disynaptic activation from itself. This positive feedback resulted in enhanced synchrony among hubs and between hubs and targets.

Individual columns of unconnected networks showed population spikes (Figure 5), where population spikes were defined as the co-firing of ≥35% of all E cells during a 20-ms interval (Neymotin et al., 2008). As seen from the raster plot, these extend across pyramidal cells of Layers 2/3 and 5. The population spikes are also apparent in the sLFP (Figure 5B). Individual columns had an average of 1.6 population spikes per minute. The number of spikes was negligible for columns without E2 hubs.

We then connected the columns as before (Figure 6). Therefore, the E2 hubs remained as strictly intracolumnar hubs and did not have any increased probability of projecting laterally to neighboring columns. However, E2 is a major lateral-projection source, so some of these E2 hubs would project from the column. As before, connecting the columns produced feed-forward inhibition. This reduced the overall hyperexcitability that was seen with the unconnected columns. Hence E2 firing rates were augmented less in comparison to the isolated columns. Rates of the other excitatory cells groups were also relatively depressed. Due to increased feed-forward inhibition from lateral connections, hub-induced population spikes were also reduced to only 0.25/min.

I2 Hubs Sculpt Out Additional High Gamma

We hypothesized that the I2 interneurons played a prominent role in determining gamma activation in the network. We tested the effects of increases in I2 hub density in the individual unconnected columns. I2 hubs were created with respect to excitatory cells: convergence from and divergence to excitatory cells was increased 3× relative to other I2 cells. As with the E2 case, the I2 cells had greatest influence over the spectrum compared to I cells in other layers. However, where the E2 hubs had created increases in spectral strength, I2 hubs caused prominent reductions in theta and beta power (Figure 7). These reductions were consistent with the expected dampening effect of inhibitory interneurons. Interestingly, dampening was far less in the gamma band, where maximum reduction was <3%. Hence the presence of the I2 hubs produced an overall shift of frequencies toward gamma through selective dampening, sculpting the response by selective frequency removal. This selective dampening is consistent with the effects of low-frequency attenuation of I MUA oscillations seen when connecting the cells (Figure 3A right).

Additional dynamical effects were also noted in this 10-I2-hub network. I2 hub rates were 2.93 ± 0.04 Hz, compared to 1.24 ± 0.02 Hz for non-hub I2 cells. E2 firing rates decreased slightly to 0.70 Hz, down from 0.76 Hz without I2 hubs. Note that the dampening effects became complex since I hubs reduced activation of E cells which then fed back to produce reduced I cell activation, mostly in the non-hub cells. At the same time, inhibition from hubs to other I cells caused further suppression of interneuronal activity.

Frequency Homeostasis

We looked at the effect of driving on the network by increasing the strength of the random inputs onto the E2 subpopulation (Figure 8). This would correspond to activation due to the back-projections from higher cortical areas, postulated to accompany increased attentional focus (Bollimunta et al., 2008). Firing rates increased substantially with this manipulation: nearly doubling to 1.1 Hz for the E2 cells, more than doubling for I2 cells to 10.19 Hz; and increasing 18, 9, 16% for E5a (0.8 Hz), E5b (0.76 Hz), and I5 (3.31 Hz) populations (compare to Table 4). Small increases were also seen in most other cell subpopulations. Remarkably, these large changes in individual cell firing rates did not produce substantial changes in the E cell spectral profile. Spectral power increased across all frequencies for the E cell population with no appreciable shift in peak, demonstrating that the network compensated in order to maintain spectral consistency. This compensation, a form of homeostasis, appeared to involve feedback from the I cell population which did show substantial shifts, with greatly increased strength in the theta and high gamma bands.

We also assessed the effects of changes in synaptic delays (Figure 9). We ran 75 20-s simulations increasing and decreasing mean values of synaptic delays by 50% relative to baseline: dendritic delays centered at 2, 4 (baseline simulation), 6 ms; somatic delays at 1, 2 (baseline), 3 ms. Increasing the delays to 150% reduced MUA spectral power and produced a significant shift in the theta peak down to 6.1 Hz, which could be partially explained by sharpening of the peak at the left side. The power reduction was explained by noting a slight but significant (p < 0.001) increase in the firing rate of I cells from 4.29 to 4.38 Hz. The added I cell activity may also be responsible for the shift in gamma toward higher frequencies. Firing rates of E cells across the three sets of simulations did not change significantly. Similarly, excitatory MUA spectral power was increased in the 50% delay simulations, with some decrease in I cell firing rate. The relative alterations in population frequency are small compared to the 50% change in synaptic delays, further demonstrating frequency homeostasis.

Coordination Across Columns

Shared spectrum across columns has been hypothesized to allow for various forms of phase locking that permits intercolumnar and interareal communication and coordination (von der Malsburg and Schneider, 1986; Gray and Singer, 1989; Engel et al., 1991; von der Malsburg, 1995). These experiments are largely done using LFP recordings. Therefore we looked at pairwise correlations between our simulated LFPs across connected columns of multiple multicolumnar networks of each simulation type (Figure 10). As expected, all network types showed substantial increases in correlations compared to the random correlations seen in the unconnected network. There was a significant further increase in correlations in the 10-E2-hub models, with an average of 0.117 compared to 0.089 (p < 0.001). This was due to the added excitation spreading first within, and then across, columns, as a result of the hub activity.

Comparison to Experimental Data

We compared the spectral power of each of our multicolumn simulations against the pattern of spectral power taken from a 30-min duration LFP recording from left medial prefrontal cortex from an awake rat (Figure 11). Although this was a characteristic spectrum, analysis of spectra from 1-min sections showed considerable variation over this period, with substantial shifts in both amplitude and position of the major frequency bands. Best matches to the spectrum shown were found by minimizing mean-squared error across different scalings of the arbitrarily-scaled simulation power spectra. The best overall match came from an E2_inc simulation with a mean-squared error of 1.4·10⁻⁴. The power spectra shown lines up with the experimental at the low-frequency peak and at the high-frequency peak in gamma near 60 Hz. Best fit error for the other simulations were 1.6·10⁻⁴ for the E2_hub simulation, 2.3·10⁻⁴ for a baseline simulation, and 2.5·10⁻⁴ for a I2_hub simulation. Other exemplars of our simulation variants also showed similarity to the experimental spectrum over the full range of frequencies, many showing a small gamma peak similar to that of the experimental spectrum. The distribution of errors for each simulation type is shown in Figure 11B. Overall, the E2_hub and E2_inc simulations matched the best, with errors of 2.76·10⁻⁴ ± 4·10⁻⁶ and 2.86·10⁻⁴ ± 4·10⁻⁶ respectively (mean ± SEM). The added E2 activity levels in these simulations supports the role of E2 in generating physiological oscillatory activity. The baseline, and I2_hub simulations had higher errors of 3.71·10⁻⁴ ± 4·10⁻⁶ and 3.99·10⁻⁴ ± 4·10⁻⁶, respectively. As noted, biological recordings also show considerable variability among LFP power spectra measured at different times or in different areas, with theta and gamma peaks always recognizable despite shifts and changes in amplitude.

Experimentally, theta and gamma show cross-frequency coupling, with activations coincident in time (Palva et al., 2005; Canolty et al., 2006). We were able to capture a similar pattern of relationship in multiple simulations (Figure 12). To quantify the coupling of power changes between pairs of frequencies across time, we used power-fluctuation correlation matrices. This involved calculating pairwise Pearson correlation between the MUA spectrogram’s time-series of power from each pair of frequencies, f₁ and f₂. The values were then stored in a matrix and plotted in a two-dimensional heat-map with the color at coordinates f₁, f₂ in the matrix representing the correlation between frequencies f₁ and f₂. High off-diagonal values indicated that the power of one frequency fluctuated in tandem with another. In Figure 12, the frequencies that show coupling, and those that do not, are similar in both experiment and simulation. Our simulations suggest that this pattern would be expected to emerge from the relation of E2 generation of theta with I2 sculpting of gamma. E2 hubs produced augmented theta directly (Figure 4) and at the same time activated I cells that fired in the gamma range (Figure 8). The I cells, I2 in particular, then dampened activity but preferentially allowed gamma expression (Figure 7).

Discussion

We have demonstrated that a surprisingly realistic pattern of firing frequencies and frequency relationships (Figures 11 and 12) emerges from a relatively simple multicolumnar network: highly simplified cells subdivided into 13 cell populations organized using anatomically and physiologically determined intra- and intercolumnar projection densities (Figure 1). Function follows form insofar as layer 2/3, containing the largest subpopulations with the strongest projections, can be identified as a primary driver (for E cells) or sculptor (for I cells) of overall spectral patterns (Figures 4 and 7). Feed-forward inhibition between columns in the full multicolumn network reduces activity, but does not change the form of the spectrum (Figure 3). The network possesses internal homeostatic mechanisms: the spectral form remains intact when stressed by forcing the network with added activation (Figure 8) and altered synaptic delays (Figure 9). Finally, the network produces a set of intercolumnar correlation relationships that could provide the basis for the neocortical phase relations putatively used in cognitive binding (Figure 10).

We make a number of testable predictions from the model:

1 Layers of high recurrent connectivity and population density, E2/3 in our current networks, will play the role of a control layer. Such a layer will largely determine network dynamics, drive other layers, and set the spectral profile. We suggest that different Brodmann areas may have different dominant layers: some primarily supragranular as here, others infragranular. This would be testable by combinations of unit and field recording in different neocortical layers.

2 Hub cells will be present in these control layers. This increases the control layer’s synchronization of network dynamics. Experimentally, hub existence and prevalence would be testable with paired/multiple simultaneous single cell recording in slice. Further simulation would allow us to test hub influence using causal measures such as Granger causality and transfer entropy (Neymotin et al., 2011).

3 Frequency homeostasis: spectral profile is resilient to changes in network activation level as well as to variation in delay. Activation of cortex testable with modulating inputs indirectly in vivo using attention modulation behavioral experimental paradigms (Bollimunta et al., 2008). Delay alteration might be testable in vitro using temperature changes.

4 Fast-spiking interneurons enhance the proportion of gamma power by selectively dampening low-frequency activity. Testable using optogenetics to selectively alter firing rates of different interneuron classes (Cardin et al., 2009).

5 Excitatory populations produce low-frequency oscillations in vivo and drive interneuron production of gamma. The feedback between these populations allows for intrinsic network dynamics to create the commonly observed phase relationship of gamma and theta oscillations. Testable using optogenetics to selectively alter firing rates of different cell classes (Cardin et al., 2009).

Graphical Rendering of Neuronal Circuits and Dynamics

Neuronal circuits are extremely complex, and not well represented by simply categorizing as feed-forward or recurrent as done by artificial neural network theory (Lytton, 2002). Understanding of emergent dynamics will require multiple depictions of various aspects of cell typology, of wiring, and of dynamics at different scales. Graphical depiction techniques may be drawn from neuroanatomy, from graph theory, and from other domains. These are multiscale models, so that different depiction scales are necessary in order to make predictions at the various levels of organization (Sejnowski et al., 1988). In our explorations of circuitry, we found that a graph-theoretic schematic provided the best clues to the major dynamical features (Figures 1B,C), but necessarily obscured the underlying cortical structure that we depicted in Figure 1A. In our exploration of dynamics, we found that rasters, MUAs, sLFPs and correlation plots were all useful for seeing different scales and different aspects of the dynamics.

The graph-theoretic schematic that we used for anatomy provided visual indicators of population proportions (nodes) and of the strength and type of connections (edges). Network features that stood out in this depiction played an important role in shaping dynamics. Important aspects of neocortical architecture that were revealed included the centrality of layer 2/3 pyramidal cells, the descending excitatory pathway from layer 2/3 to 5, and negative feedback patterns featuring stronger I → E compared to E → I projections. Additional schematics (not shown) suggested a dominant lateral inhibitory effect between columns. The graph visualizations led us to our initial hypotheses and simulation experiments which focused on tests for the role of E2 and I2.

Similar graphical techniques may be particularly beneficial when detailed simulations are not feasible, and can help play a linkage role in the development of multiscale models. Working backward, it is also valuable to use these visualization techniques on experimental (or simulation) dynamical data to infer functional connectivity, as is routinely done for functional magnetic resonance imaging data. Assessing graphs of functional connectivity at different time-periods may also allow for measuring similarity of correlation structure across time in an intuitive manner.

Spectrum as an Emergent Property

We set up our network from what is known from neurophysiology and neuroanatomy. We then had to tune the network by adjusting weights in order to set approximate known firing rates of the cells in the different cell subpopulations. This tuning is typically required for simulations, in part because the numbers of cells and convergence of inputs are so much less than those found in the brain (Izhikevich and Edelman, 2008). Part of this tuning involves the provision of continuous ongoing random synaptic inputs (white-noise). These white-noise inputs are meant to very imperfectly represent the incoming activations and excitations from other parts of the brain, which may not have a white-noise spectral profile under more realistic conditions (Kerr et al., 2011).

Tuning our network required adjustment of local properties utilizing only global parameters: the strength of external excitation, and four tuning factors for weight of connection between excitatory and inhibitory populations: E → E, E → I, I → E and I → I connections. Cell properties themselves remained unchanged. The local property that was the main tuning target was firing rates for all cell types, keeping them within the bounds of propriety. Some further tuning was then done for a global property: weights were reduced in order to avoid any appearance of impropriety at the network level: hypersynchrony and epileptiform activity. After tuning based on these two criteria, physiologically realistic power spectra and cross-frequency coupling appeared as an emergent property through the cooperativity of the individual neuronal elements in the system. The general spectral form revealed was similar to that of rat brain in terms of both direct comparison (Figure 11) and frequency interrelationships (Figure 12). These results suggest that our model may have captured aspects of neocortical dynamics.

In order to evaluate how readily a comparable broad spectrum arose from other models, we surveyed several neocortical and thalamocortical models from the literature. Many papers, including those describing our prior models, do not show power spectra (Lytton and Omurtag, 2007; Lytton et al., 2008b; Neymotin et al., 2011). Generally this is because most models are focused on some particular phenomenon that does not require exploration of this facet of model behavior. Some others, that do concern themselves with frequencies, model the sleep state, or a slice, and produce spectra not similar to the whole-brain awake state. As an example, a neocortical slice simulation by Traub et al. models kainate-induced gamma oscillations and shows a single peak at 30 Hz (Figure 2 of Traub et al., 2005). However, this same paper has a model with a much broader spectrum (Figure 11A of Traub et al., 2005) that is comparable to in vivo activity. That spectrum has additional high-frequency peaks due to modeling of epileptiform activity. Another complex model is that of Izhikevich and Edelman(Izhikevich and Edelman, 2008). Although they do not show spectra, they noted that areas show local fast (40–50 Hz) traveling gamma rhythms, and note that other areas (of this very large whole-brain model) are dominated by 20 Hz beta. However, they state that the overall model was dominated by 1–3 Hz delta rhythms with intermixed alpha (Figure 5 of Izhikevich and Edelman, 2008), rhythms that, as they point out, do not typically co-occur at high power in vivo.

Most other neocortical models are far simpler than the two just described. Despite its simplicity, a model of mu (10 Hz motor) rhythm by Jones et al. shows a fairly realistic awake spectrum (Figure 7 of Jones et al., 2009). The absence here of high gamma frequencies is justified, since the experimental data being modeled also lacks activity in that band. Other simplified multineuronal models generally show tight frequency peaks. A generalized cortical network by Rulkov et al. shows a tight peak at 45–50 Hz (Figure 13 of Rulkov et al., 2004), as does a generalized network model by Brunel and Wang (Figures 1 and 6 of Brunel and Wang, 2003). A model by Karameh et al. has a more complex power spectrum, but still shows a relatively simple spectral structure with three separate peaks at 8, 15, and 24 Hz (Figure 11 of Karameh et al., 2006). By contrast, simplified “lumped” models, based on differential equations representing brain regions rather than individual brain cells, may readily demonstrate far more complex spectral fidelity (e.g., Figure 9 of von Albada et al., 2010).

From this brief survey, we conclude that simple models produce simple spectra, generally narrowly focused on a single frequency. Broad spectra are likely to emerge from complex multilayer networks. However, we note that different simulations produced different spectral patterns that sometimes involved physiologically disjoint frequencies or abrupt frequency cut-offs.

Assessing Causality in a Complex Dynamical System

In highly interconnected, high-dimensional dynamical systems, such as neuronal networks, bidirectional and nested loops are omnipresent, making it extremely difficult to ascertain which elements are causally responsible for each other’s activity pattern (Douglas and Martin, 2007b). Traditional methods for determining causality in non-linear dynamics, such as phase plane analysis, are not useful when applied directly to such high-dimensional networks, requiring an initial dimensional reduction which is of itself a major research project and which necessarily eliminates much of the dynamics. Other methods that have been developed to assess causality include transfer entropy (Schreiber, 2000; Gourevitch and Eggermont, 2007; Neymotin et al., 2011; Quinn et al., 2011), Granger causality (Roebroeck et al., 2005; Chen et al., 2006; Seth, 2008), Bayesian dynamical modeling (Baldi et al., 1998; Rajapakse and Zhou, 2007; Eldawlatly et al., 2008), and combinations of modeling and optimization (Makarov et al., 2005). Although these methods are useful in many scenarios, a unique benefit of simulation is that it allows the experimenter to have full control of the system. This allows the use of a more direct method to ascertain dynamic causality, by performing perturbations to a system and observing the effects (Paus, 2005; Chen et al., 2008; Daunizeau et al., 2009). However, a caveat is in order. If the perturbation is too drastic, then one is likely to enter a new dynamical regime rather than accurately probing the dynamics of the current regime. An example of an overly drastic intervention would be any perturbation producing a parameter change that crosses a bifurcation boundary, since this will create a new dynamical pattern that was not there before.

One of the most extreme, yet most commonly used, such perturbations is the ablation experiment. In this type of experiment, one proceeds around the system removing one component at a time until some oscillation or activity pattern of interest, pattern A, disappears. Let’s say that removal of component Q abolishes pattern A. One then declares victory: component Q is the producer, cause, or primary driver of pattern A. However, what if the remaining network, network minus Q, has moved into an entirely different pattern, pattern B? Should we now define component Q as driver of pattern A, or as suppressor of pattern B? An anecdote from the era of vacuum tubes bears repeating. Removing a particular vacuum tube from the radio “motherboard” would transform the white-noise static of the untuned radio to a high amplitude, narrow-band signal (a howl). That evidence did not justify calling the particular tube the “static-generator unit,” nor the “howl-suppressor unit.”

To evaluate dynamics, we avoided drastic network changes by using only minimal dynamical perturbation as a methodology. We chose hub formation as a subtle perturbation that could be applied gradually to the system, starting at a low level which would have no detectable effect. Gradual addition of hubs could then be done without changing the number of elements or changing the total number of connections, thereby maintaining the basic structural attributes of the system. With addition of hubs, the influence of the hubbed subpopulation gradually increased as activity spread more effectively. For example, addition of hubs in the E2 subpopulation had no observable effect until a threshold was reached at about five hubs (3.3% of population). Above this point, the addition of additional hubs progressively augmented power in all bands tested (Figure 4). Note that though the number of hubs were changed identically in all subpopulations, the proportion of hubs at any given number was lowest in the case of E2 due to the larger size of this subpopulation. In our inhibitory subpopulation experiments, the addition of even one hub produced alterations in power, suggesting that we should have set our hub’s convergence and divergence level smaller in order to get less effect with one hub. Despite this, we were able to demonstrate a clear distinction between hub addition in the I2 subpopulation and hub addition in the other inhibitory subpopulations (Figure 7).

Experimental Evidence for Hubs in the Brain

Graph-theoretic hubs have been implicated in the spread of epidemics, synchronization, and enhanced communication in complex networks across a wide variety of fields and spatiotemporal scales (Strogatz, 2001; Wang and Chen, 2003). Examination of the connectivity distribution of a nematode’s neurons revealed that it obeyed a power-law distribution with rarely occurring, but highly connected hub cells present (Barabási and Albert, 1999; Strogatz, 2001). The role of hubs in synchronizing neuronal network activity was shown in a series of experimental and simulation studies of dentate gyrus (Dyhrfjeld-Johnsen et al., 2007; Morgan and Soltesz, 2008). In these studies, it was shown that granule cell hubs played a contributing role in epileptogenesis even when only a few hub cells were added. This is similar to the higher frequency of population spikes observed in our model when several E2 hub cells were present. These studies also showed that efficacious hubs had increases in both incoming and outgoing connections: hubs with only increased incoming connections produced much less hyperexcitability, while hubs with only increased outgoing connections produced no hyperexcitability. The existence of these hubs was explained by noting that some few cells had both large dendritic and large axonal arborizations.

Recent experimental work demonstrated that in neocortical networks, layer 5 pyramidal neurons showed non-random connectivity suggestive of hubs: a low probability of high-connection strength between pairs of neurons (Song et al., 2005). This study found a large asymmetry in the strength of synapses: 17% of the synaptic connections were shown to contribute half of the total synaptic strength, suggesting that some neurons may be significantly more effective in causing cells to fire than others. In addition, there was a correlation between the strength of incoming and outgoing connections. From these observations, the authors suggested the existence of pyramidal hub neurons in the neocortex.

Interneuron hubs have been demonstrated in developing hippocampus (Bonifazi et al., 2009). Functional and effective connectivity analysis was performed and revealed a power-law distribution. The hubs were found to have much larger axonal arborizations than the average cell in the network, accounting for large divergence. Dynamically, the hubs were able to synchronize network activity, either when stimulated or spontaneously, producing population spikes. Similarly, we observed the synchronizing effects of hubs in our simulations.

At the larger spatial scale of networks of brain regions, hubs have been found to play an important role in functional integration (Hagmann et al., 2008). Analysis of networks of brain areas with diffusion-weighted MRI also revealed the presence of hubs (Iturria-Medina et al., 2008). Other recent work in fMRI-based functional networks has demonstrated the presence of hubs in regions of unimodal or multimodal association cortex, suggesting these types of hubs can integrate different forms of information effectively (Achard and Bullmore, 2007).

The Centrality of Layer 2/3 in this Network

In our simulations, the E2 subpopulation, the excitatory cells in layer 2/3, appeared to provide a primary oscillatory drive to the network, particularly in the theta/alpha band. In addition, E2 appeared to strongly drive the I2 cells to produce gamma, sculpted out of the network through differential damping. E2 hubs were able to synchronize the intracolumnar activity via the mechanism of positive feedback between convergent and divergent cells. Activity spread first within, and then between columns. Our simulations comparing unconnected and connected columns demonstrated that columns may communicate by modulating the excitation within their neighboring columns via feed-forward activation of both E and I cells.

Layer 2/3 is implicated in the generation of theta oscillations in entorhinal cortex (Ylinen et al., 1995; Sirota and Buzsaki, 2005), and motor cortex (Weiler et al., 2008). Additionally, layer 2/3 has been suggested to be a major driver and determinant of activity in a number of studies of sensory cortex (Fransen and Lansner, 1998; Lundqvist et al., 2010). By contrast, a classical physiological study, Silva et al., does not support this model (Silva et al., 1991). This in vitro study evaluated smaller and smaller fragments of sensorimotor neocortical slice and found that only layer 5 could produce activity in isolation. Clearly, this is a form of the type of ablation experiment criticized above. Additionally, slicing will injure cells, and injured cells may be prone to hyperexcitability.

However, rather than describe this as an either/or situation, we would note that there are considerable anatomical differences across cortical regions, as originally described by Brodmann. These feature different thicknesses of layers, numbers of cells in each subpopulation, and apparently different connectivity densities. These anatomical differences are likely to be associated with dynamical differences, including alterations in spectral profile. Additionally, there is physiological evidence to suggest that the direction of driving differs from one area to another. In the awake, behaving animal study of Bollimunta et al. (2008), Granger causality inferred infragranular (layer 5/6) to supragranular (layer 2/3) driving in sensory areas (V2, V4), but supragranular to infragranular driving in area IT (inferotemporal cortex). This suggested the possibility of dynamical differences between lower-level sensory areas (V2, V4) and higher cortical association areas (IT). Based on this distinction, we would propose our model as a model for higher centers. We will then need to produce variant models for primary sensory cortex, and for motor cortex. As with most areas of biology, there will be no single cortical model (hence, alas, no direct model falsification) but rather complex multiscale combinations of model variants.

Network Stabilization: Frequency Homeostasis

Homeostasis, the maintenance of stability under varying environmental conditions, is a feature of any biological system. In the brain, one form of homeostasis entails maintenance of baseline activity patterns across varying levels of inputs and neuromodulators, avoiding seizures despite maintenance of the high levels of correlation needed for intracortical communication. We call this frequency homeostasis, the maintenance of the correct spectrum of frequencies and frequency interrelationships despite alternation of drive. This has been demonstrated physiologically (Bollimunta et al., 2008; Rajagovindan and Ding, 2011). In our simulations, we saw that adding power to the primary driving center, layer 2/3, a major perturbation that produced increased frequencies in most cell groups, caused the alpha/theta peak to grow but not shift substantially (Figure 8). This pattern of added driving might be produced in vivo by attentional focus provided through strong feedback from higher cortical centers to a lower center, an influence that arrives primarily in supragranular layers (Quilodran et al., 2008). We hypothesize that this homeostatic effect in the simulation would be provided by a balance between I cell domination of the high-end of the spectrum and E cell domination of the low-end of the spectrum. The increased high-band inhibition counteracted the tendency of the existing peak to start spilling over into higher frequency bands, preserving the overall shape of the spectral power profile. An additional homeostatic mechanism was shown with the preservation of population frequencies despite alterations in the delays between neurons (Figure 9).

The beginning of a failure of homeostasis is shown in Figures 5 and 6 – the occurrence of population spikes, harbinger of seizure activity. Notably, the strong feed-forward inhibition between columns served to reduce the number of population spikes considerably, partly suppressing population spikes and restoring something of a normal appearance to the sLFPs (Figure 6). Additional homeostatic processes remain to be explored and explained. For example, we expect that the spectral relationship shown in Figure 12 will also be resistant to increased activation. It will also be important to check the effects of activation by way of other pathways, most notably via sensory activation through layer 4.

Cognitive Implications

Brains receive information from different sensory modalities as well as internal sources and integrate them into coherent representations (Treisman, 1996). One hypothesis for how the brain solves this problem is via the synchronized activation of distributed neuronal ensembles (von der Malsburg and Schneider, 1986; Engel and Singer, 2001; Singer, 2001; Wu et al., 2007). Additionally, synchronization across brain areas has been suggested as a prerequisite for optimal information transfer between them (Buehlmann and Deco, 2010). Oscillations have been suggested as a mechanism for fine temporal tuning of neural firing times which could provide coding and learning via spike timing (Canolty et al., 2006; Fries et al., 2007; Clopath et al., 2010). Specifically, gamma and theta oscillations have been suggested as partners in a multiplexing mechanism via subsets of cells that fire on a particular gamma subcycle superimposed on a theta cycle (Lisman and Idiart, 1995). Human neocortical theta has also been implicated in working memory tasks (Raghavachari et al., 2001; Jensen and Tesche, 2002). Co-occurrence of theta and gamma, demonstrated in the present study in both experiment and simulation (Figure 12), would provide a substrate for these types of mechanisms.

While fairly high levels of synchronization are required for cognition, pathological increase in synchronization may lead to seizures (Osterhage et al., 2007; Schindler et al., 2007; Lytton, 2008; Lehnertz et al., 2009; Neymotin et al., 2010). This suggests a tradeoff: synchronization between components is required for information processing, but too much leads to pathology. Frequency homeostasis might then be important to maintain conditions required for a cognition that relies on tuned brain oscillations (von der Malsburg and Schneider, 1986; Gray and Singer, 1989; von der Malsburg, 1995; Kopell et al., 2000; Spencer et al., 2003; Buzsaki and Draguhn, 2004; Uhlhaas et al., 2006).

In our model, E2 cells had the strongest recurrent connectivity and had excitatory projections laterally via E5, as well as strong lateral inhibitory projections. We therefore hypothesize that layer 2/3, together with layer 5, will be primarily responsible for the intercolumnar coordination we demonstrate in Figure 10. When we added hubs to E2, we produced population spikes in the isolated columns, an epileptiform indicator (Figure 5). These hubs were strongly suppressed by lateral inhibition in the connected network but were still apparent (Figure 6). The enhanced synchrony due to hubs also increased correlation across columns, possibly allowing for enhanced information transfer, but equally possibly associated with excessive coordination that would actually reduce information transfer (Olypher et al., 2006).

Our simulations modeled both the generation and the coupling of oscillations in the theta and gamma bands, and also suggest a possible mechanism for this coupling. Theta oscillations emerged primarily from E cell activity. This activity then provided drive to the interneurons, which then produced gamma oscillations. Feedback of this activity onto E cells then sends a gamma timing signal into the E population. Interestingly our columnar sLFPs showed fluctuations of synchrony. We hypothesize that states of low intercolumnar synchrony may be characterized by intracolumnar computation, and states of high intercolumnar synchrony may be periods when the columns communicate. We speculate that synchrony fluctuation may allow for isolating periods of local and global coordination.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

Supported by: DARPA grant N66001-10-C-2008, National Institute of Mental Health (R01MH084038); National Institute of Neurological Disorders and Stroke (R42NS064474). The authors would like to thank the reviewers for their insightful comments and for help in improving the manuscript, Mingzhou Ding, Gordon M. G. Shepherd, Ben Suter, Kimberle Jacobs for assistance with the neocortical model, Ahmet Omurtag for discussions on data analysis, Larry Eberle for computer support and administration, and Cliff Kerr for helpful comments on the manuscript.

References