Abstract

Slow N-Methyl-D-aspartic acid (NMDA) synaptic currents are assumed to strongly contribute to the persistently elevated firing rates observed in prefrontal cortex (PFC) during working memory. During persistent activity, spiking of many neurons is highly irregular. Here we report that highly irregular firing can be induced through a combination of NMDA- and dopamine D1 receptor agonists applied to adult PFC neurons in vitro. The highest interspike-interval (ISI) variability occurred in a transition regime where the subthreshold membrane potential distribution shifts from mono- to bimodality, while neurons with clearly mono- or bimodal distributions fired much more regularly. Predictability within irregular ISI series was significantly higher than expected from a noise-driven linear process, indicating that it might best be described through complex (potentially chaotic) nonlinear deterministic processes. Accordingly, the phenomena observed in vitro could be reproduced in purely deterministic biophysical model neurons. High spiking irregularity in these models emerged within a chaotic, close-to-bifurcation regime characterized by a shift of the membrane potential distribution from mono- to bimodality and by similar ISI return maps as observed in vitro. The nonlinearity of NMDA conductances was crucial for inducing this regime. NMDA-induced irregular dynamics may have important implications for computational processes during working memory and neural coding.

Introduction

Networks in prefrontal cortex (PFC) actively maintain goal-related information in the absence of external cues, a capacity termed working memory, through persistently elevated neural firing rates (Fuster 1973; Funahashi and others 1989; Goldman-Rakic 1995; Rainer and others 1998). Active networks in PFC are characterized by both strong recurrent N-Methyl-D-aspartic acid (NMDA) synaptic input and dopaminergic modulation through D1 receptors (Sawaguchi and others 1990; Sawaguchi and Goldman-Rakic 1994; Dudkin and others 1997; Watanabe and others 1997; Aura and Riekkinen 1999; Lewis and O'Donnell 2000; Durstewitz and Seamans 2002; Seamans and others 2003; Phillips and others 2004). These ingredients could strongly support the persistent firing of prefrontal neurons through the almost constant synaptic drive provided by slowly decaying NMDA currents (Wang 1999; Compte and others 2000; Durstewitz and others 2000; Brunel and Wang 2001; Koulakov 2001; Durstewitz and Seamans 2002), which are further amplified by D1 receptor activation (Zheng and others 1999; Seamans and others 2001; Wang and O'Donnell 2001; Chen and others 2004; Tseng and O'Donnell 2004).

PFC spiking activity is not just uniformly enhanced during working memory tasks but—like in many other cortical areas (Dean 1981; Softky and Koch 1993; Holt and others 1996)—exhibits a high amount of irregularity that significantly increases throughout the memory periods (Compte, Constantinidis, and others 2003), that is, the periods associated with a rise in NMDA input. This high spiking irregularity might be a signature of other important dynamical processes that characterize working memory activity, depending on whether it is mainly due to inherently stochastic (Shadlen and Newsome 1994, 1998; Brunel and Wang 2003) or deterministic factors. Deterministic chaos (Strogatz 1994) at the single neuron or network level could give rise to highly irregular firing even in the absence of any noise or random factors (Hansel and Sompolinsky 1996; van Vreeswijk and Sompolinsky 1996, 1998). If PFC networks would reside in a close to or slightly chaotic regime with complex deterministic temporal structure, this may strongly facilitate the processing of temporal sequences (Bertschinger and Natschläger, 2004; Jaeger and Haas, 2004). The necessity to integrate temporal sequence information is a crucial feature of any working memory task, and part of the PFC capacity (Ninokura and others, 2003, 2004). Thus, the sort of dynamics prefrontal neurons exhibit, and its biophysical basis, will profoundly affect information processing during working memory.

Because NMDA conductances and D1 receptors are assumed to play such a vital role in shaping PFC attractor dynamics and activity states (Wang 1999; Durstewitz and others 2000; Lewis and O'Donnell 2000; Brunel and Wang 2001; Koulakov 2001; Durstewitz and Seamans 2002; Seamans and others 2003; Fellous and Sejnowski 2003), we were interested in dynamical properties of the spiking of PFC neurons induced by NMDA and D1 receptor activation. PFC slices bathed in a combination of NMDA and D1 receptor agonists have been reported previously to exhibit spontaneous activity with interesting spatiotemporal signatures (Wang and O'Donnell 2001; Beggs and Plenz 2003; Tseng and O'Donnell 2005). Here we used the same sort of preparation and observed that high spiking irregularity occurred within a certain regime characterized by a transition in the membrane potential distribution. Spike trains within this regime exhibited short-term predictability that was incompatible with a purely linear process driven by random inputs. This suggested a potential mechanism for the observed irregularity that was further analyzed through computational modeling. Modeling results confirmed that high irregularity could be due to chaotic dynamics within a regime characterized by a transition in the membrane potential distribution. Adding nonlinear NMDA conductances to the dendrites of single model neurons was sufficient to induce such a regime. Based on these simulation results, it is therefore suggested that NMDA conductances may lift PFC neurons into a chaotical regime with a characteristic signature in the membrane potential distribution. Parts of the data have been presented previously in abstract form (Gabriel and others 2004).

Materials and Methods

In Vitro Methods

Long-Evans rats (P > 55) were deeply anesthetized with isoflurane and decapitated. Brains were rapidly removed into ice-cold artificial cerebrospinal fluid (ACSF) containing (in mM) NaCl (125), NaHCO3 (26), glucose (10), KCl (3.5), NaH2PO4 (1.25), CaCl2 (0.5), MgCl2 (3); pH 7.45, osmolarity 300 ± 5 mOsm. Coronal slices (450 μm) containing the medial PFC (mPFC) were cut on a Vibratome and transferred into warm (±37 °C) ACSF solution with CaCl2 increased to 2 mM and MgCl2 decreased to 1 mM. After 20 min, ACSF was allowed to cool down to room temperature for at least 40 min before recording. All solutions were constantly oxygenated with 95% O2/5% CO2 throughout the experiments.

Pyramidal neurons located in deep layers of the medial PFC were identified under visual guidance using infrared-differential interference contrast video microscopy with a ×40 water immersion objective. Thick-walled borosilicate patch pipettes (2–6 MΩ) were filled with (in mM) K-gluconate (115), 4-(2-hydroxyethyl)-1-piperazineethanesulfonic acid (10), MgCl2 (2), KCl (20), MgATP (2), Na2-ATP (2), and Guanosine 5′-triphosphate (0.3) (pH = 7.3, 285 ± 5 mOsm). Ag/AgCl wires and pellets (Warner, Hamden, CT) were used for the patch and reference electrodes, respectively. Whole-cell current-clamp recordings were performed with an Axoclamp 2B amplifier (Axon Instr., Foster City, CA) and acquired with Labview-based software, written by Lee Campbell (Salk Institute, La Jolla, CA), at sampling rates of 2–5 kHz. Pipette series resistance was <20 MΩ, and was compensated. Electrode potentials were adjusted to zero before recording, without correcting the liquid junction potential.

NMDA (7–8.5 μM; Sigma-Aldrich, Munich, Germany; a concentration which is within a physiological range, see Melendez et al., 2005) and the selective D1 agonist SKF 38393 (2–5 μM; Tocris, Ellisville, MO) were bath applied within the ACSF. In some experiments, 12–20 μM DNQX (Sigma-Aldrich) and 24–30 μM bicuculine (Sigma-Aldrich) were added to block α-Amino-3-hydroxy-5-methylisoxazole-4-propionic acid (AMPA) and γ-aminobutyric acid A (GABAA) currents, respectively. In the cells probed with current steps for control (see Results), additionally all NMDA currents were blocked by 2-Amino-5-phosphonopentanoic acid (APV) (50 μM, Sigma-Aldrich). All drugs were mixed into ACSF from frozen stock solution aliquots. All experiments were conducted at 34–36 °C in ACSF at a perfusion rate of 2–3 ml/min.

Analysis Methods

Recordings from single cells consisted of one to several samples of 30–300 s length, separated by brief nonrecording periods. To judge stationarity both within and between samples of the same cell, the original voltages traces and the instantaneous firing rate (=1/interspike-interval [ISI]) as a function of time were first inspected by eye. Traces were then divided into consecutive segments of 15 or 30 s (depending on total recording duration, and down to 3 s for the NMDA-omission controls [see Results]) for which ISI distributions together with their mean and standard deviation were plotted (for all segments overlaid within the same graph). All ISI distributions were compared through Kolmogorov–Smirnov tests, probing for significant differences between any 2 distributions. For none of the samples judged as stationary, the overall corrected (for number of statistic comparisons) P value fell below 10%, indicating that all ISI distributions belonged to the same process. Note that this stationarity requirement implies stability of the spiking statistics across the entire length of recording periods used for analysis, unconfounded by potential drifts caused by intracellular washout etc.

For quantifying bimodality in the membrane potential distributions, spikes were first cut out within the range [tsp − 2 ms, tsp + 4 ms], for all tsp which denote the times of crossing a −20-mV threshold from below. Next, the sum of 2 Gaussian-like exponentials was fitted to the distribution, which had the general form 

graphic
where hV denotes the Vm distribution and μ1 ≤ μ2. So whereas one exponential (with μ2 as mean) is a pure Gaussian, the other exponential (with μ1 as mean) was allowed different slopes right versus left from the mean, to account for the fact that unimodal and left-hand parts of bimodal Vm distributions were often strongly skewed (see e.g., Fig. 1B bottom). Fits with 6 free parameters were performed in 2 different ways. For the first type of fit, k = 2 was fixed. For the second type of fit, μ1 = μ2 was enforced, but k was allowed to vary in the range [0, 3]. Hence, the second type of fit was unimodal but with the same number of free parameters as the bimodal fit. Fitting was performed by constrained optimization using MatLab's (The MathWorks, Inc., Natick, MA) lsqcurvefit-routine. In general, whichever of these 2 fits produced a smaller (absolute) error was accepted. Only in <10% of cases, a clearly uni- or bimodal Vm distribution was not detected as such by this algorithm and had to be enforced “by hand.” Based on these fits, a parameter dV quantifying the degree of bimodality within the membrane potential distribution was defined as dV = |μ1 − μ2|/σ2.

Figure 1.

Highest spiking irregularity within PFC slices occurs in a regime where the membrane potential distribution is at the edge from mono- to bimodality. (A) Example traces from 3 different prefrontal neurons with long-duration bursts (top), highly irregular spiking (middle), and continuously regular spiking (bottom). Pharmacological conditions were similar in all 3 cases (NMDA: 8 μM; SKF 38393: 3 μM (top) or 4 μM [middle and bottom]). (B) Membrane potential distributions of these neurons indicating a shift from a clearly bimodal (top), through an intermediate (middle), to a clearly unimodal (bottom) distribution (gray = Vm distribution, black line = function fit, see Materials and Methods). (C) Index for the amount of bimodality within the membrane potential distribution (dV), various indices for the amount of spiking irregularity (Cv, CvL, Lv, HISI), and average duration of bursts (〈Tbs〉 in s).

Figure 1.

Highest spiking irregularity within PFC slices occurs in a regime where the membrane potential distribution is at the edge from mono- to bimodality. (A) Example traces from 3 different prefrontal neurons with long-duration bursts (top), highly irregular spiking (middle), and continuously regular spiking (bottom). Pharmacological conditions were similar in all 3 cases (NMDA: 8 μM; SKF 38393: 3 μM (top) or 4 μM [middle and bottom]). (B) Membrane potential distributions of these neurons indicating a shift from a clearly bimodal (top), through an intermediate (middle), to a clearly unimodal (bottom) distribution (gray = Vm distribution, black line = function fit, see Materials and Methods). (C) Index for the amount of bimodality within the membrane potential distribution (dV), various indices for the amount of spiking irregularity (Cv, CvL, Lv, HISI), and average duration of bursts (〈Tbs〉 in s).

The average duration of bursts (〈Tbs〉) was determined with a 2-threshold method (Anderson and others 2000), after cutting spikes within 35 ms windows and linear interpolation between end points. Durations briefer than 10 ms were excluded. For unimodal Vm distributions, the average burst duration was set to 0.

Four different measures for the irregularity of the spike trains were obtained:dV, 〈Tbs〉, and all measures of irregularity were first computed separately for all samples and then averaged across all samples within a stationary set.

  1. The coefficient of variation, Cv = σISIISI, where μISI and σISI denote the mean and standard deviation of the ISI distribution.

  2. A local measure of the Cv (denoted CvL), which measures the Cv within temporal windows containing on average 4 spikes and then averages across all these local measurements. CvL reduces the impact of long ISIs during burst-like activity which could lead to a high Cv despite otherwise completely regular activity.

  3. forumla (Shinomoto and others 2003). The Lv is another measure introduced to quantify the true amount of irregularity, minimizing the contribution of burst-like but otherwise very regular activity.

  4. HISI = − ∑nP(ISI ∈ Bn) log2P(ISI ∈ Bn) with bins B0 = [0, 2] ms and Bn = [(1 − α) cn, (1 + α) cn] for n > 0. The bin centers (cn) were determined as cn = cn−1 (1 + α)/(1 − α) (with α = 0.1). This procedure ensures that bin size increases proportionally to the bin center (i.e., |Bn|/cn = const. = 2α) because the ISI standard deviation can be expected to scale with the mean (but other choices for HISI yielded similar results). HISI gives the entropy (the amount of information) of the ISI distribution and is not obscured by burst-like but regular activity. Regular activity with only a narrow range of different ISIs will always yield relatively low HISI values.

For determining nonlinear deterministic structure in the data, a nonlinear prediction algorithm (Kantz and Schreiber 2004) was applied to original and surrogate ISI data. Briefly, the method is as follows: First, a delay embedding of the ISI series is performed by combining m temporally consecutive measurements (separated by Δn time steps in the ISI sequence) into m-dimensional vectors, yielding a trajectory in an m-dimensional state space (for details, see Kantz and Schreiber 2004; we used Δn = 1 and m = 3, but also checked m = 6; for the present purpose the choice of m is not so crucial). Next, for each point in this m-dimensional space, neighboring points are sought within some neighborhood ε. (One might think of this procedure simply as looking up similar m-dimensional ISI patterns in the data base.) These neighbors are used to form a prediction of Δk (≥1) time steps (ISIs) ahead into the future according to the formula 

(1)
graphic
where ISIn is the m-dimensional ISI vector under consideration with components (ISIn−m+1 ⋯ ISIn) for Δn = 1 and Uε is the neighborhood of size ε of this vector (|·| denotes its cardinality; for each point ε was iteratively adjusted such that at least 5 neighbors could be found). Temporal neighbors < 10 ISIs away from the ISI vector under consideration were excluded, as these may simply pick up temporal correlations not due to reoccurring deterministic structure. A local prediction error is then calculated as squared discrepancy between predicted and actual future values and averaged across all N points (with ≥Δkmax future values) in m-dimensional state space: 
(2)
graphic

To check for the possibility that the predictability within ISIs could be due to a simple autoregressive linear process with random inputs (Autoregressive Moving Average), phase-randomized surrogate data were created (Schreiber and Schmitz 1996, 2000). Formally, the null hypothesis tested against states that the time series can be described by any process of the form 

(3)
graphic
where the ηk are independent Gaussian random numbers (〈ηk〉 = 0, 〈ηk2〉 = 1) and f is a monotonic but possibly nonlinear function (for details, see Schreiber and Schmitz 2000). Briefly, surrogates that obey this null hypothesis are created by “rescaling” original values to a Gaussian distribution (empirically inverting f), performing a fast Fourier transform (after reducing discontinuities between end points to minimize periodicity artifacts), randomizing the phases (but retaining the amplitudes), transforming back into the time domain, and rescaling to the original distribution (because of the rescaling, some “polishing” may be necessary and distribution and power spectrum may have to be adjusted iteratively in several steps; see Schreiber and Schmitz 2000). The resulting data sets have the same power spectrum (autocorrelation function) and distribution as the original ISI series but destroy any nonlinear dependencies and hence, predictability that would result specifically from nonlinear deterministic processes. For performing the predictions and creating the surrogate data, the TISEAN 2.1 software package (Hegger and others 1999; Schreiber and Schmitz 2000) was used, where the prediction routine “zeroth” was adapted to accommodate multiple temporally discontinuous samples.

Computational Methods

Networks of 2-compartment pyramidal cells (PCs), equipped with 6 different ionic conductances (INa, INaP, IDR, IKS, IHVA, IC), and single-compartment interneurons (IN), equipped with just fast Na+ and DR K+ channels, were constructed. Gating variables for all ionic conductances followed Hodgkin-Huxley–like channels kinetics. The basic model cells and synapses were taken from Durstewitz and Seamans (2002), and details regarding parameters and equations can be found there and in Durstewitz and others (2000). Neurons were connected through AMPA- and NMDA-like excitatory (〈gAMPAmax〉 = 0.02 μS, 〈gNMDAmax〉 = 〈gAMPAmax〉/50), or GABAA-like inhibitory (〈gGABAAmax〉 = 0.006μS) connections, respectively (see Durstewitz and Seamans 2002). The NMDA current is given by 

(4a)
graphic
with 
(4b)
graphic
where gNMDA(t) is a sum of double exponentials as in Durstewitz and others (2000) that implements the time course of all recurrent NMDA inputs to the cell (parameters as in Durstewitz and Seamans 2002) and s(Vm) implements the very fast nonlinear voltage dependence (made instantaneous here, without any impact on results). All synaptic conductances were multiplied by a short-term plasticity factor (absorbed into gNMDA(t) in eq. 4a) that developed according to the formalism given in Tsodyks and Markram (1997; see also eq. 6 in Durstewitz 2003). PC→PC connections exhibited short-term depression (Markram and others 1998; utilization USE = 0.3, recovery from depression τrec = 0.8 s, facilitation τfac = 0.8 s; for details, see Tsodyks and Markram 1997), IN→PC (USE = 0.25, τrec = 0.6 s, τfac = 1 s) and IN→IN (USE = 0.12, τrec = 0.6 s, τfac = 1.2 s) connections weak short-term depression (Reyes and others 1998), and PC→IN (USE = 0.05, τrec = 0.15 s, τfac = 1.2 s) connections short-term facilitation (Markram and others 1998).

The standard network consisted of 100 PCs and 25 INs, with a random 20% connectivity between different neurons. For the basic observations reported here, precise connectivity and network size did not seem of major relevance. Constant (i.e., time independent) NMDA activation was added to the dendritic compartment of all PCs and to INs. This conductance had the same voltage dependence s(Vm) and reversal potential as the synaptic NMDA conductance (eq. 4), but the product of gNMDAmax and the time-dependent term gNMDA(t) modeling the rise and decay of the conductance after receptor activation was replaced by a constant, gNMDAcons (with a mean of 0.095 mS/cm2 for the network simulations reported here). The inclusion of this conductance was meant to mimic the continuous application of NMDA to the bath solution in the in vitro experiments. D1 application was mimicked by doubling the inactivation time constant of INaP, reducing IKS everywhere by 60%, reducing dendritic IHVA by 15%, and reducing gleak of IN by 25% (Yang and Seamans 1996; Gorelova and Yang 2000; Gorelova and others 2002), with regards to the “standard” parameters given in Durstewitz and Seamans (2002). We stress, however, that such D1-like manipulations are not important for the basic phenomena reported here (similar results were obtained without them).

Many parameters followed a Gaussian distribution within the network, with the following standard deviations given as percentage of reference (mean) values in brackets: Somatic and dendritic INaP, IKS, IHVA, and IC maximum conductances (all 20%), size of PCs and INs (with all morphological dimensions changed by the same factor with 20% standard deviation), gNMDAcons (20%), and all synaptic conductances (50%). Network simulations were run with NEURON (Hines and Carnevale 1997).

Single neuron simulations were also performed with NEURON and were often double-checked with MatLab-based code, using a model cell picked from the network. Single cell results were further confirmed by drawing at random 10 cells from the network distribution and in model cells without D1-like changes. NEURON was used with its implicit variable time-step solver CVODE (Cohen and Hindmarsh 1996), and for the MatLab simulations, the implicit multi-step solver ode15s (Shampine and Reichelt 1997) was taken. Simulation results were probed with relative tolerances of down to 10−8 and absolute tolerances of 10−10. All simulation code is available from the following URL: www.pion.ac.uk/∼daniel.

For the single PC model, 2-dimensional (Vsoma, dVdend/dt) state space projections were obtained analytically. For given Vsoma, all somatic gating variables and Ca2+ and K+ ionic concentrations were solved for their steady-state values. By setting dVsoma/dt = 0, the dendritic membrane potential can be obtained as Vdend = Vsomari[Iinj + ∑ gi(EiVsoma)], where ri is the internal (axial) resistance, the gi are the total conductances (including gating), and Ei are the corresponding reversal potentials. Using Vdend, steady-state values for dendritic gating variables and ionic concentrations could be computed and the curve dVdend/dt = F(Vsoma) (for which dVsoma/dt = 0) could be derived, as used in Figure 9B. Wherever this function becomes 0, the whole neuron is in a fixed point. Exact solutions for the fixed points were obtained numerically in the vicinity of the graphically identified points (Fig. 9B) using MatLab's fzero function. The Jacobian matrix of the complete PC system of 26 differential equations was then evaluated at these fixed points, that is, eigenvalues were computed, to determine their stability (for an introduction, see Strogatz 1994).

To gain insight into the bursting process of single NMDA-equipped model neurons, a reduced model was constructed by singling out 3 major slow variables (time constants > 100 ms) and replacing them by constants. These were the voltage-dependent inactivation gate (v) of the high-voltage–activated dendritic Ca2+ current IHVA, the dendritic [Ca2+]i concentration which controls the Ca2+-dependent activation of the K+ current IC, and the voltage-dependent inactivation (h) of the dendritic persistent Na+ current INaP (see Durstewitz and others 2000; Durstewitz and Seamans 2002). The values of all 3 variables were fixed to their steady-state values according to a single given voltage parameter Vslow: vfixed = v(Vslow), hfixed = h(a1Vslow + b1), and [Ca2+]fixed = [Ca2+](a2Vslow + b2), where the steady-state value of [Ca2+]i is a function of voltage via Ca2+ influx through IHVA. Note that while v is set to its steady-state value for Vslow directly, the other 2 steady-state values are computed for linear functions of Vslow. The parameters ai and bi of these functions were determined from simulations of the corresponding full 26-variable model by least-squared-mean-error fits to the covariation plots of Vslow (v) = v−1 (v) versus Vslow([Ca2+]i) = [Ca2+]−1 ([Ca2+]i), and Vslow(v) versus Vslow(h) = h−1 (h) (where −1 denotes the inverse function). This way, bifurcation graphs could be obtained (with fixed points and their stability determined similarly as described above) as functions of just a single parameter, Vslow, cutting through the 3-dimensional parameter space (vfixed, [Ca2+]fixed, hfixed) in an optimal manner with reference to the full model behavior.

Results

Irregular Activity in NMDA-Driven PFC Neurons

To examine the impact of NMDA- and D1 receptor activation on the spiking dynamics of PFC neurons, adult rat (d > 55) PFC slices were treated with a cocktail of a D1 agonist (SKF 38393, 2–5 μM) and NMDA (7–8.5 μM), added to the ACSF (see also Beggs and Plenz 2003; Tseng and O‘Donnell 2005). Somatic current-clamp patch recordings from PCs within this preparation revealed a variety of different spontaneous activity patterns, ranging from long-duration bursts occurring at regular (or sometimes irregular) intervals, highly irregular spiking, to quite regular continuous spiking (Fig. 1A; despite their sometimes long duration, we will usually refer to events as in the top graph of Fig. 1A as “bursts,” to stress that we do not imply that they represent the same phenomenon as seen during sleep cycles in vivo, see Discussion). These different sorts of patterns could be observed for different neurons even for exactly the same concentrations of agonists (although different concentrations may shift the balance, see Tseng and O’Donnell 2005) and hence, presumably reflect largely differences in biophysical and morphological properties of the recorded cells, and/or differences in total (including recurrent) NMDA input.

Different spiking patterns were associated with different membrane potential distributions (Fig. 1B). For burst-like activity patterns, it was clearly bimodal, whereas for continuous regular spiking in most cases, this distribution was clearly unimodal. The highly irregular spiking patterns, however, in most cases had a membrane potential distribution that ranged somewhere in between, where 2 peaks just seemed to separate (Fig. 1B, middle) and continuous spiking was intermingled with brief-duration (<200 ms) “mini-upstates” (Fig. 1A, middle). To quantify the amount of spiking irregularity, we obtained a number of different indices: The standard coefficient of variation (Cv), the Lv (Shinomoto et al., 2003) which captures more of the true (ISI-to-ISI) irregularity and reduces the impact of bursting which can produce a high Cv even when spiking is otherwise highly regular, a local version of the Cv (CvL) that likewise minimizes contributions of long ISIs between bursts, and the entropy (HISI) of the ISI-distribution with bin sizes proportional to the interval centers (for details, see Materials and Methods). As shown in Figure 1C, CvL, Lv, and HISI all attain a maximum for the irregular spiking pattern, while the Cv is closest to 1 (the Cv of a Poisson process) for this pattern.

Figure 2 summarizes the results for n = 108 stationary recording sets (from a total of 62 cells). The membrane potential distributions for all cells were characterized through a single parameter (dV) that quantifies the degree to which the distribution splits into 2 underlying components with different means. In Figure 1C, this parameter increases from bottom to top, in accordance with the tendency of the distribution toward bimodality. Figure 2A,B shows as histograms that the Lv and HISI attain a maximum (while the standard Cv is closest to 1, not shown) for a dV range indicative of Vm distributions that just start to separate into 2 discernable peaks. Figure 2C gives the dependence of CvL on dV for all stationary recording sets as dots, to better visualize the variance in the data (clearly unimodal distributions always have dV = 0, clustering along the ordinate). Analyses of variance comparing the 3 groups dV < 1, 1 ≤ dV < 4, and dV ≥ 4 yielded highly significant main effects (all F2,105 > 53, all P < 10−15) and significant post hoc effects for all group comparisons for all 4 irregularity indices. Figure 2D, finally, gives the average duration of bursts as a function of dV, revealing that brief-duration bursts just start to appear for the range of dV values for which the irregularity is highest (〈Tbs〉 ∼ 164 ms for 1 ≤ dV < 4). In summary, the highest spiking variability occurs in a transition regime where the membrane potential moves from uni- to bimodality and where continuous spiking and doublets start to intermingle with brief bursts.

Figure 2.

Summary graphs: Highest spiking irregularity occurs close to the emergence of bimodality. (A) The Lv as a measure of spike train irregularity attains a maximum for small to intermediate dV values, indicating Vm distributions at the border of bimodality. (B) The entropy of the ISI distribution (HISI) is also maximal for the same range of dV values. (C) Local CvL for all stationary samples (n = 108, from 62 cells) as a function of dV. (D) The duration of bursts is highly but not perfectly linearly correlated with dV. For dV values within the transition regime (<4), brief-duration bursts just start to appear. Error bars = standard error of the mean.

Figure 2.

Summary graphs: Highest spiking irregularity occurs close to the emergence of bimodality. (A) The Lv as a measure of spike train irregularity attains a maximum for small to intermediate dV values, indicating Vm distributions at the border of bimodality. (B) The entropy of the ISI distribution (HISI) is also maximal for the same range of dV values. (C) Local CvL for all stationary samples (n = 108, from 62 cells) as a function of dV. (D) The duration of bursts is highly but not perfectly linearly correlated with dV. For dV values within the transition regime (<4), brief-duration bursts just start to appear. Error bars = standard error of the mean.

Next, we asked whether the irregular spiking behavior in the prefrontal slices bears any signs characteristic of a chaotic process. Figure 3A shows the first-return maps of the ISI time series for the 3 cells illustrated in Figure 1. The first-return map plots the length of each ISI as a function of the length of the preceding ISI. For the bursting (Fig. 1A top) and the regular continuous spiking cell (Fig. 1A bottom), the (ISIn+1, ISIn) pairs are clustered within 3 disconnected or 1, respectively, regions of the return map (Fig. 3A, top and bottom), indicating a considerable degree of periodicity (albeit, expectedly, the absolute jitter increases for large ISIs between bursts). In contrast, for the highly irregular behavior (Fig. 1A, middle), the points are scattered all along the ISIn- and ISIn+1 axes within a single continuous band (Fig. 3A, middle), indicating much less periodicity. On the other hand, the points are not just randomly distributed across the whole map but are still confined within a quite narrow band. Such characteristics, quite aperiodic yet not completely random behavior in some cases while more periodic in others, and a noninvertible return map ISIn+1 = F(ISIn), are typical for chaotic processes. In fact, it will be shown later that chaotic model neurons in a noise-free network exhibited the same sort of structure (Fig. 8A). Yet we caution that these observations are still compatible with a much simpler (nonchaotic) deterministic process involving noise (see Discussion).

Figure 3.

ISI first-return maps and nonlinear predictability within ISI time series for the traces from Figure 1. (A) ISIn+1 as a function of preceding ISIn: For the bursting cell (top), clustering tends to occur within 3 disconnected regions of the map, indicating successions of brief ISIs occasionally interrupted by long ISIs. For the continuously regular spiking cell (bottom), ISIs cluster within one region on the diagonal indicating that the cell fires at a fairly constant rate. For the irregular cell (middle), ISIs smear out along the 2 axes within a continuous band, indicating the absence of clearly periodic components. (B) Prediction error (PEISI, normalized to the standard deviation, σISI, of the data) as a function of prediction time step for the original ISI series (black solid), the mean of the 99 surrogates (gray solid), and the 90% interval around this mean (gray dashed) (The expectancy value of PEISIISI for a stationary renewal process will range between 1 and sqrt(2) depending on neighborhood size.).

Figure 3.

ISI first-return maps and nonlinear predictability within ISI time series for the traces from Figure 1. (A) ISIn+1 as a function of preceding ISIn: For the bursting cell (top), clustering tends to occur within 3 disconnected regions of the map, indicating successions of brief ISIs occasionally interrupted by long ISIs. For the continuously regular spiking cell (bottom), ISIs cluster within one region on the diagonal indicating that the cell fires at a fairly constant rate. For the irregular cell (middle), ISIs smear out along the 2 axes within a continuous band, indicating the absence of clearly periodic components. (B) Prediction error (PEISI, normalized to the standard deviation, σISI, of the data) as a function of prediction time step for the original ISI series (black solid), the mean of the 99 surrogates (gray solid), and the 90% interval around this mean (gray dashed) (The expectancy value of PEISIISI for a stationary renewal process will range between 1 and sqrt(2) depending on neighborhood size.).

The hallmark of deterministic processes is predictability. If the irregular spiking is due to nonlinear deterministic processes, subsequent ISIs should be predictable to some degree on the basis of other observations with a similar history of preceding intervals. Because predictability in chaotic systems exponentially vanishes, the prediction horizon may be limited, but predictability should still be better than chance at short timescales. More precisely, we tested the hypothesis that predictability within the irregular spike trains is higher than expected from a linear autoregressive process driven by Gaussian random inputs and a potentially nonlinear but monotonic instantaneous measurement (or response) function. For this, we constructed surrogate data (using the TISEAN package, Hegger and others 1999; Schreiber and Schmitz 2000) that preserve the distribution and autocorrelation (power spectrum) of the original ISI sequences. The results for the 3 cells from Figure 1 are shown in Figure 3B. Figure 3B (black lines) shows the uncertainty about a future ISI as a function of the prediction step (the number of ISIs extrapolated into the future). The gray solid line gives the average of 99 phase-randomized surrogate data, and the gray dashed lines mark the 90% interval around this mean assuming a Gaussian distribution of the data (cutting off 5% of the distribution below and above, as we are interested in 1-tailed probabilities). As can be seen, although predictability quickly vanishes for the irregular and the bursting cell, it is significantly higher in the original than in the surrogates at least for the first 2–3 time steps. This holds also true for a nonparametric comparison, as 99/100 of the data sets (99 surrogates + original series) had a larger prediction error than the original ISI series for the first 2 prediction steps (for the irregular cell, also the third step was significant). We obtained similar results for a number of cells for which the quality and amount of data (>500 ISIs) were sufficient to perform this kind of analysis: According to nonparametric statistics, 21/26 stationary recording sets (from 19 cells) yielded a significantly better prediction than the surrogates for the first time step, 20/26 for the first and second, 13/26 up to the third, 10/26 up to the fourth, and still 2/26 were consistently significant even up to the tenth time step (parametric comparisons based on the normal distribution assumption gave very similar results). This provides further evidence that the temporal variation in the spike trains is at least partly due to complex nonlinear deterministic processes, rather than to purely probabilistic factors or linear processes with random inputs.

Next, we checked the contribution of AMPA and GABAA synaptic currents to the observed spiking irregularity (a total of 27 cells providing n = 52 stationary sets). Although the density of functional synaptic inputs is very much reduced in the PFC slices, it was still surprising that after blockade of these currents with DNQX (12 μM) and bicuculine (24 μM) essentially the same range of phenomena could be observed as before (Fig. 4). The summary statistics (Fig. 5) also gave similar distributions as for the nonblocked condition, with maximum irregularity for dV values within an intermediate range, as indicative of Vm distributions at the edge of bimodality. Again, all post hoc comparisons for the dV groups as defined above were significant for all 4 irregularity indices, except for HISI, where the intermediate and the high dV group were statistically en par (possibly partly due to the lack of cells with dV > 7 in this sample), although there was, once again, a tendency for cells with intermediate dVs to carry the highest entropies (not shown).

Figure 4.

Same as Figure 1 for PFC slices after blockade of AMPA (DNQX: 12 μM) and GABAA (bicuculine: 24 μM) synapses. (A) Example traces of 3 different PFC neurons with bursting activity (top), highly irregular spiking (middle), and continuously regular spiking (bottom). (B) Membrane potential distributions and function fit for these neurons. (C) Index for the amount of bimodality within the membrane potential distribution (dV), various indices for the amount of spiking irregularity (Cv, CvL, Lv, HISI), and average duration of bursts (〈Tbs〉 in s).

Figure 4.

Same as Figure 1 for PFC slices after blockade of AMPA (DNQX: 12 μM) and GABAA (bicuculine: 24 μM) synapses. (A) Example traces of 3 different PFC neurons with bursting activity (top), highly irregular spiking (middle), and continuously regular spiking (bottom). (B) Membrane potential distributions and function fit for these neurons. (C) Index for the amount of bimodality within the membrane potential distribution (dV), various indices for the amount of spiking irregularity (Cv, CvL, Lv, HISI), and average duration of bursts (〈Tbs〉 in s).

Figure 5.

Summary graphs for PFC slices after blockade of AMPA and GABAA synapses. Both the Lv (A) and the local CvL (B) still peak for dV values within an intermediate range (cf., Fig. 2). Error bars = standard error of the mean.

Figure 5.

Summary graphs for PFC slices after blockade of AMPA and GABAA synapses. Both the Lv (A) and the local CvL (B) still peak for dV values within an intermediate range (cf., Fig. 2). Error bars = standard error of the mean.

Finally, as a control, we omitted application of NMDA to the bath but kept the D1 agonist (SKF 38393, 3 μM), blocked all synaptic input (DNQX 12–20 μM, bicuculine 24–30 μM, APV 50 μM), and probed cells with a series of depolarizing and a few hyperpolarizing current steps (20–30 s). Without NMDA in the bath, all cells (n = 9/9) displayed only single-spiking trains after the adaptation transient at the onset of a step (Fig. 6), that is, none of them exhibited sustained bursting or an irregular transition regime as in Figures 1 and 4 (in line with Tseng and O'Donnell, 2005, who also reported that burst-like patterns within a similar sort of PFC preparation crucially depended on NMDA).

Figure 6.

Single cell responses in the absence of NMDA. (A) Example of the steady-state responses of a PFC pyramidal neuron in vitro to long depolarizing current steps when all synaptic currents were blocked (DNQX 20 μM, bicuculine 30 μM, APV 50 μM) and NMDA was omitted from the bath solution, but D1 receptors were stimulated (SKF 38393: 3 μM). (B) Average steady-state ISIs as a function of injected current for the same cell (first step from A omitted since only one ISI was observed). Error bars give the standard deviations that are small compared with the mean, reflecting the regular single-spiking process. (C) Coefficient of variation (Cv) versus average ISI for all steady-state traces from 8 recorded PFC neurons. Except for one trace, the Cvs are all very low, indicating highly regular single-spiking activity. The one trace with a higher Cv (∼0.6) comes from a cell with “stuttering” single-spike responses, that is, regular spiking with occasional gaps but no bursting (the same cell had a very low Lv and CvL, both <0.11).

Figure 6.

Single cell responses in the absence of NMDA. (A) Example of the steady-state responses of a PFC pyramidal neuron in vitro to long depolarizing current steps when all synaptic currents were blocked (DNQX 20 μM, bicuculine 30 μM, APV 50 μM) and NMDA was omitted from the bath solution, but D1 receptors were stimulated (SKF 38393: 3 μM). (B) Average steady-state ISIs as a function of injected current for the same cell (first step from A omitted since only one ISI was observed). Error bars give the standard deviations that are small compared with the mean, reflecting the regular single-spiking process. (C) Coefficient of variation (Cv) versus average ISI for all steady-state traces from 8 recorded PFC neurons. Except for one trace, the Cvs are all very low, indicating highly regular single-spiking activity. The one trace with a higher Cv (∼0.6) comes from a cell with “stuttering” single-spike responses, that is, regular spiking with occasional gaps but no bursting (the same cell had a very low Lv and CvL, both <0.11).

Computational Model

We constructed networks of 100 2-compartment conductance-based PCs and 25 single-compartment IN, randomly connected through AMPA-, NMDA-, and GABAA-like synapses, respectively (the precise connectivity seemed not crucial for the present observations). All neurons received an additional source of constant NMDA receptor activation (time independent, but voltage-(Mg2+)-gated NMDA conductances), supposed to mimic the constant application of NMDA to the bath solution. D1-like changes in cell parameters (Yang and Seamans 1996; Gorelova and others 2002) with regard to a standard reference model (Durstewitz and others 2000; Durstewitz and Seamans 2002) were also included but were not crucial for the principal phenomena to occur. Model neurons varied with regard to their morphological dimensions and strength of several voltage-gated conductances, as well as the strength and number of synaptic inputs, and the constant NMDA activation (reflecting natural variation among neurons). Importantly, there was no source of noise or randomness in the running model network.

The whole range of phenomena as observed in vitro could be reproduced in the network, from long-duration bursts, to highly irregular spiking, to regular continuous spiking (Fig. 7). Note that the high irregularity in the model cells must be due to deterministic chaos because there was no noise source in the simulations. Like in the in vitro preparations, the highest irregularity occurred in a regime where a bimodal distribution of membrane potential just started to form (compare Fig. 7 middle with top and bottom). Figure 8 furthermore illustrates that for these model cells similar ISI return maps and evidence for nonlinear predictability are obtained as in the real neurons (compare with Fig. 3). Finally, removing AMPA and GABAA synapses from the network, highly irregular firing and long-duration bursts could still be observed, confirming the presence of those phenomena in the in vitro preparations with AMPA/GABAA currents blocked (not shown). (Interestingly, the remaining NMDA connections were also sufficient to synchronize bursts between the bursting neurons.)

Figure 7.

Same as Figure 1 for a purely deterministic biophysical network model of the PFC slice. (A) Example traces of 3 different pyramidal neurons (with randomly drawn parameters, see Materials and Methods) from the same simulated network with bursting activity (top), highly irregular spiking (middle), and continuously regular spiking (bottom). (B) Membrane potential distributions and function fits for these neurons. (C) Index for the amount of bimodality within the membrane potential distribution (dV), various indices for the amount of spiking irregularity (Cv, CvL, Lv, HISI), and average duration of bursts (〈Tbs〉 in s).

Figure 7.

Same as Figure 1 for a purely deterministic biophysical network model of the PFC slice. (A) Example traces of 3 different pyramidal neurons (with randomly drawn parameters, see Materials and Methods) from the same simulated network with bursting activity (top), highly irregular spiking (middle), and continuously regular spiking (bottom). (B) Membrane potential distributions and function fits for these neurons. (C) Index for the amount of bimodality within the membrane potential distribution (dV), various indices for the amount of spiking irregularity (Cv, CvL, Lv, HISI), and average duration of bursts (〈Tbs〉 in s).

Figure 8.

Same as Figure 3 for the chaotic biophysical network model. (A) ISI return maps for the model neurons from Figure 7 are similar to the ones shown for real prefrontal neurons in Figure 3A. (B) Prediction error as a function of prediction time step behaves similarly for the 3 different types of model neuron as for the real neurons (Fig. 3B), indicating a significant nonlinear contribution to predictability for the irregular (middle) and bursting (top) cell but less so for the continuously regular cell (bottom).

Figure 8.

Same as Figure 3 for the chaotic biophysical network model. (A) ISI return maps for the model neurons from Figure 7 are similar to the ones shown for real prefrontal neurons in Figure 3A. (B) Prediction error as a function of prediction time step behaves similarly for the 3 different types of model neuron as for the real neurons (Fig. 3B), indicating a significant nonlinear contribution to predictability for the irregular (middle) and bursting (top) cell but less so for the continuously regular cell (bottom).

Next, we examined single model PCs in isolation, trying to get closer to the potential mechanism of the irregular chaotic dynamics. Adding a time-independent NMDA conductance to the dendritic compartment of a single model cell was sufficient to induce irregular behavior with a phenomenology quite closely resembling that observed in vitro (Fig. 9A). Varying the strength of this NMDA conductance, the model cell traversed different regimes, including bursting and regular spiking in addition to the irregular behavior. Such a transition from bursting to regular spiking through chaos was also observed in 10/10 other model cells with randomly drawn biophysical and morphological parameters (see Materials and Methods), as well as in model cells without D1-like changes. To many of these, a small negative biasing current (≥−0.1 nA) had to be added to reveal these transitions (the amplitude required seemed to inversely correlate with the strength of the dendritic IHVA). While in the network every cell was irregular to some degree, however, chaotic firing in single model cells was observed only within a rather restricted regime of constant NMDA activation (Fig. 10A).

Figure 9.

Neurodynamical basis of spiking irregularity at the single neuron level. (A) A single deterministic model cell (same cell as in Fig. 7A middle) exhibits similar spiking regimes as the full network for different levels of constant NMDA activation (gNMDAcons values as in B), with one regime of chaotic irregularity with frequent membrane potential transitions and mini-upstates (middle). (B) Two-dimensional state space projection of the biophysical model cell, showing dVdend/dt as a function of Vsoma (black line, labeled dVsoma/dt = 0 in bottom graph), assuming that all other model variables (gates and ion concentrations) are in their steady states. Hence, the neuron is in a fixed point wherever the dVdend/dt curve crosses the dVdend/dt = 0 axis, but none of them is stable to small perturbations. Gray lines indicate the trajectories corresponding to the time graphs in (A). The chaotic regime (middle), characterized by a trajectory that densely fills a whole region of state space, appears around the bifurcation point where 2 of the 3 fixed points coalesce and vanish.

Figure 9.

Neurodynamical basis of spiking irregularity at the single neuron level. (A) A single deterministic model cell (same cell as in Fig. 7A middle) exhibits similar spiking regimes as the full network for different levels of constant NMDA activation (gNMDAcons values as in B), with one regime of chaotic irregularity with frequent membrane potential transitions and mini-upstates (middle). (B) Two-dimensional state space projection of the biophysical model cell, showing dVdend/dt as a function of Vsoma (black line, labeled dVsoma/dt = 0 in bottom graph), assuming that all other model variables (gates and ion concentrations) are in their steady states. Hence, the neuron is in a fixed point wherever the dVdend/dt curve crosses the dVdend/dt = 0 axis, but none of them is stable to small perturbations. Gray lines indicate the trajectories corresponding to the time graphs in (A). The chaotic regime (middle), characterized by a trajectory that densely fills a whole region of state space, appears around the bifurcation point where 2 of the 3 fixed points coalesce and vanish.

Figure 10.

Bifurcation graphs of a single model cell with constant NMDA activation showing all ISIs for given gNMDAcons. (A) Nonlinear NMDA current: For most values of gNMDAcons, the cell exhibits only a few different ISIs and hence fires periodically with different periods, but for some gNMDAcons range (∼0.103–0.114 mS/cm2), aperiodic (irregular) behavior appears where the ISI can assume all possible values within a large range. (B) Linear NMDA current: When the nonlinearity of the NMDA current is removed, the same model cell just responds with regular single-spike trains.

Figure 10.

Bifurcation graphs of a single model cell with constant NMDA activation showing all ISIs for given gNMDAcons. (A) Nonlinear NMDA current: For most values of gNMDAcons, the cell exhibits only a few different ISIs and hence fires periodically with different periods, but for some gNMDAcons range (∼0.103–0.114 mS/cm2), aperiodic (irregular) behavior appears where the ISI can assume all possible values within a large range. (B) Linear NMDA current: When the nonlinearity of the NMDA current is removed, the same model cell just responds with regular single-spike trains.

We next examined the fixed points of the single pyramidal neuron with constant NMDA activation, that is, the points where neither the membrane potentials nor any of the channel gating variables or [Ca2+]i and [K+]o change anymore. These points can be determined by setting all temporal derivatives of the differential equation system describing the model neuron to 0. The fixed points are graphically illustrated in Figure 9B by plotting the temporal derivative of the dendritic membrane potential (i.e., the total dendritic current divided by capacitance) as a function of the somatic membrane potential (curve labeled dVsoma/dt = 0 in Fig. 9B bottom), assuming that all gating variables and ionic concentrations are in their steady states, and that dVsoma/dt = 0 (i.e., Vsoma = const.). Consequently, all points where in addition dVdend/dt = 0 mark the fixed points of the whole system. Within these 2-dimensional spaces spanned by Vsoma and dVdend/dt, we also plotted the temporal evolution of the system (the trajectory) for the respective parameter configurations (gray lines), which correspond to the time graphs in Figure 9A.

We now examined what happens at the transition to chaos. Below the regime where chaos sets in (Fig. 9B top), the system has 3 fixed points whose stability can be determined from linear stability analysis of the full 26-dimensional PC differential equation system. All 3 fixed points are unstable, that is, any small perturbation leads the system away from these fixed points—the system never precisely settles into one of them and would stay there only if placed exactly at one of these points. The lowest fixed point (around −56 mV) is, however, a stable equilibrium for low gNMDAcons (not shown) and, above some critical gNMDAcons, changes to an unstable spiral point (through a Hopf bifurcation), that is, trajectories are repelled from this point in an outward spiral-like fashion (and in the limit cycle around it). This causes trajectories approaching the neighborhood of this point being either propelled far to its left (Fig. 9B top) or sliding fast along its right side, giving rise to the subthreshold membrane potential bimodality.

Chaos sets in as, with increasing gNMDAcons, the unstable spiral point moves to the right into the path of the approaching trajectory where it coalesces with the second fixed point in a saddle node bifurcation (Fig. 9B middle). Within the vicinity of this bifurcation, the cycling trajectory starts to densely fill up the whole space between approximately −58 and approximately −46 mV which accounts for the broad Vm distribution and the merging of its 2 peaks as evidenced in Figures 1A, 4A, and 7A (middle). With further increases of gNMDAcons, as the bifurcation region finally disappears, chaos vanishes (Fig. 9B bottom). Instead, without any fixed points left within the range below −40 mV, trajectories pass right through and spikes start shooting off from a rather uniform membrane potential level, creating the narrow unimodal Vm distribution witnessed in Figures 1A, 4A, and 7A (bottom).

To summarize, the single cell model equipped with dendritic NMDA conductances can explain the variety of spiking patterns observed in the PFC slices (Fig. 1A), in particular the occurrence of a highly irregular (chaotic) regime, and why this is accompanied by a broad Vm distribution with merging peaks (intermediate dV values). It is also interesting to note that a hyperpolarizing current injection would lower the dVsoma/dt nullcline, hence moving the model cell from a regular spiking to a bursting regime. In agreement herewith, we observed that hyperpolarizing current steps in vitro could change a continuously spiking mode into the bursting mode (not shown).

Is the nonlinearity of the NMDA current, imposed by the fast voltage dependence of its Mg2+ block, crucial for its ability to induce bursting and chaotic spiking? In the model, it is possible to specifically remove the voltage dependence by setting the Mg2+-dependent voltage-gate s (eq. 4b) to a constant value (we set sfixed = s(−51 mV), approximately the value where the saddle node bifurcation happens in Figure 9B, but the exact value does not matter because changing it just corresponds to a rescaling of the gNMDAcons axis). Figure 10B shows that after removal of the NMDA nonlinearity, the model neuron passes from silence directly to regular spiking with increasing gNMDAcons values (the same tendency was observed for all 10 other randomly drawn neurons, of which only one still started off with spike-doublets before reverting to regular single-spiking, and for different amounts of biasing current, from −0.1 to 0 nA) Moreover, strongly slowing down the Mg2+-dependent gating by equipping the s gate with a very slow time constant (5 s) had essentially the same effect, that is, model neurons displayed only regular spiking, neither bursting nor chaos. This is expected because strongly slowing down the voltage-dependent gating effectively linearizes the NMDA current–voltage relationship on short timescales.

Finally, to understand more precisely the nature of NMDA-induced bursting, we singled out the 3 major slow variables of the model and studied its behavior fixing those variables as parameters (see e.g., Rinzel and Ermentrout 1998, for this approach to the study of bursting neurons). Specifically, we set the inactivation (v) of the dendritic high-voltage–activated Ca2+ current IHVA, the dendritic [Ca2+]i concentration (steering Ca2+-dependent activation of the K+ current IC), and the inactivation gate (h) of the dendritic persistent Na+ current INaP to fixed values parameterized through a single voltage variable denoted Vslow (the parameterization was done such that the resulting 1-dimensional cut through the 3-dimensional parameter space (v, [Ca2+]i, h) optimally aligned (in a least-squared-error sense) with the free fluctuation of these variables in the full model; see Materials and Methods for details). Fixing these 3 variables, the model settles into either silence or regular single spiking. Figure 11A shows a bifurcation graph of this reduced 23-variable model. All fixed points of the system are plotted as a function of Vslow (higher values corresponding to stronger inactivation of IHVA and INaP and stronger activation of IC), with stable fixed points indicated by continuous and unstable ones by broken black lines. The lower line of dots which terminates on the middle unstable branch of fixed points demarks the lower edge of a limit cycle (stable periodic orbit), that is the lowest somatic membrane potential reached during each cycle of sustained spiking in the stable regular spiking regime, and the upper line of dots demarks its upper edge.

Figure 11.

Bifurcation diagram of a reduced single cell model where 3 major slow variables (parameterized by Vslow) have been separated from the remaining dynamics (see text). The diagrams show all the fixed points (stable ones as continuous black lines, unstable ones as dashed lines) as a function of Vslow, as well as the lower and upper edge of the limit cycle (dotted) and the (Vslow(ν), Vsoma) trajectory (gray solid line) of the full model projected onto this plane. (A) Nonlinear NMDA current: With the NMDA nonlinearity in place, a broad hysteresis region appears (from approximately −58.7 to −47 mV) where the limit cycle and the lower stable fixed point coexist, and the trajectory cycles a few times (the burst) before crawling back below the middle unstable arch toward the next bursting cycle. (B) Linear NMDA current: The hysteresis region is dramatically reduced, and the trajectory leaves the limit cycle region already during the first spike, giving rise to a simple single-spike train.

Figure 11.

Bifurcation diagram of a reduced single cell model where 3 major slow variables (parameterized by Vslow) have been separated from the remaining dynamics (see text). The diagrams show all the fixed points (stable ones as continuous black lines, unstable ones as dashed lines) as a function of Vslow, as well as the lower and upper edge of the limit cycle (dotted) and the (Vslow(ν), Vsoma) trajectory (gray solid line) of the full model projected onto this plane. (A) Nonlinear NMDA current: With the NMDA nonlinearity in place, a broad hysteresis region appears (from approximately −58.7 to −47 mV) where the limit cycle and the lower stable fixed point coexist, and the trajectory cycles a few times (the burst) before crawling back below the middle unstable arch toward the next bursting cycle. (B) Linear NMDA current: The hysteresis region is dramatically reduced, and the trajectory leaves the limit cycle region already during the first spike, giving rise to a simple single-spike train.

The crucial point about this graph (Fig. 11A) is the existence of a broad hysteresis region, that is, the coexistence of the limit cycle with a stable resting point over the Vslow range approximately −58.7 to −47 mV. The bursting process in the full model may now be explained as follows: During spiking within a burst, IHVA and INaP slowly inactivate, whereas the Ca2+ concentration slowly rises and enhances IC. Due to this increasing loss of depolarizing force, the spiking process breaks down at some point and the (Vslow(v), Vsoma) trajectory (gray line in Fig. 11A), projected onto the (Vslow, Vsoma) graph of the reduced model, passes over to the right into the region where the lower branch of fixed points constitutes the only stable solution (Vslow > −47 mV). Due to the hysteresis region, the trajectory then crawls back within the basin of the lower branch of stable fixed points (due to reactivation of IHVA and INaP and slow drop in [Ca2+]i), which accounts for the interburst interval, until it pops onto the limit cycle again as the stable resting solution vanishes in a saddle node bifurcation where the stable and unstable branches meet (at approximately −58.7 mV). The remarkable consequence of such a (square-wave-type) bursting process is that chaotic solutions between the bursting and regular single-spiking regimes have been shown to be a generic property of such systems under rather general assumptions (Terman 1992).

Figure 11B depicts the situation with a linear NMDA current (s fixed): The hysteresis region is dramatically diminished, in fact so much that Vslow(v) passes over to the monostable resting region during the upstroke of the first spike and bursting is no longer supported. Hence, the cell engages only in regular single spiking after removal of the NMDA nonlinearity.

Discussion

We have shown here that highly irregular spiking in NMDA-driven prefrontal neurons occurs within a transition regime, where the activity exhibits a mixture of various spiking modes, and the membrane potential distribution hovers at the edge from mono- to bimodality. Both in cells that exhibit no sign of bimodality and in cells that quite clearly do so, firing is more regular. Within the irregular regime, ISIs still maintained a degree of predictability that is incompatible with a purely linear or random process, therefore indicating the presence of nonlinear deterministic structure. Moreover, the various spiking patterns observed in vitro, together with their signatures in the ISI return maps and membrane potential distributions, could be reproduced in purely deterministic model neurons. A detailed analysis of single model cells furthermore revealed a crucial role for the nonlinearity of the NMDA current in establishing bursting and chaos and provided insights into the nature of these processes and their association with particular Vm distributions. The combination of these various observations therefore suggests that the empirically observed irregularity could be due to deterministic chaos induced by NMDA currents.

Neural Chaos and Mechanisms of Spiking Irregularity

A variety of neural models, both single cells and networks, exhibits chaos within certain parameter regimes (Chay and Rinzel 1985; Fan and Holden 1992; Hansel and Sompolinsky 1996; van Vreeswijk and Sompolinsky 1996, 1998; Bernasconi and others 1999; Clay 2003; Compte, Sanchez-Vives, and others 2003). In many of the single neuron systems, chaos occurs at the transition from bursting to single-spiking modes (Chay and Rinzel 1985; Holden and Fan 1992; Rinzel and Ermentrout 1998; Shuai and others 2003; Shilnikov and others 2005), as observed here. Chaos is indeed an inevitable consequence of a certain class of bursting processes (Terman 1992). Bursting leading into chaos may be established through a variety of different biophysical means in the nervous system, including for instance mechanisms depending on persistent Na+ or dendritic Ca2+ currents (Magee and Carruth 1999; Shuai and others 2003; Traub and others 2003). In rat PFC with a paucity of intrinsic bursters (Yang and others 1996), recurrent NMDA inputs may provide an interesting network-regulated mechanism for pushing many neurons into a chaotic regime which they would not attain based on their intrinsic properties alone.

For the present study, an interesting alternative interpretation to deterministic chaos is that cells within the highly irregular domain are at the edge from true mono- to true bistability, with noise driving them back and forth between 2 “barely” stable states. Such a regime, approximating a line attractor configuration, would strongly facilitate noise-induced fluctuations (Durstewitz 2003), thereby causing high irregularity. During clear monostability, in contrast, there would be no transitions between 2 states, and during clear bistability transitions would be rare as noise is dampened and has a harder time to drive the cell from one state into the other. Strictly, this scenario is not ruled out by the present time series analysis results, as they still allow for the possibility that the time series were generated by nonlinear yet completely periodic processes, irregularized by external noise.

Although this question cannot yet be decided based on the present experimental data, for three reasons we feel this alternative scenario at least cannot be the sole explanation for the observed irregularity. First, as noted above, such a configuration close to the edge of true bistability strongly facilitates noise-induced fluctuations, approaching a random walk situation. If this is indeed the cause of the high irregularity, much of the deterministic structure would be destroyed, and predictability should be lowest within this regime. In contrast, the relative prediction error (PEISIISI) was significantly lower for the transition regime (1 ≤ dV < 4) than for the steady spiking (dV < 1) or bursting (dV ≥ 4) regimes, as revealed by Tukey post-hoc tests. This supports the idea that much of the large ISI variance in the transition regime may be due to cell-intrinsic deterministic processes, while the smaller fluctuations in the other two regimes may be caused mainly by synaptic inputs or noise superimposed on otherwise periodic processes like those in Fig. 9 top and bottom (this interpretation is furthermore backed by the slower leveling off of the PEISIISI-curve for the transition compared to the other two regimes in Fig. 3B). Second, removing one major noise source, namely, synaptic inputs mediated by AMPA and GABAA currents, the highly irregular regime remained still in place (Figs 4 and 5). Third, we showed computationally that NMDA currents were capable of inducing chaotic behavior with the same signatures as observed experimentally and, moreover, that they may induce bursting through a mechanism that implies chaos at the transition from bursting to regular spiking regimes (Terman 1992).

Notwithstanding, even if the chaotic interpretation were basically correct, external noise will most likely add to it, and NMDA-induced chaos and random transitions between ‘quasi-stable’ states (cf. Fig. 11A; see Compte, Sanchez-Vives, and others 2003) caused by noise may work hand in hand to generate the empirically observed irregularity. A similar high irregularity at the transition from bursting or down/up-state oscillatory to regular spiking modes is also apparent in a chaotic network model by Compte, Sanchez-Vives, and others (2003), and in recent cortical cell culture studies where transitions may partly be regulated through cholinergic mechanisms (Tateno and others 2005; Fujisawa and others 2006).

A number of other mechanisms have been proposed to explain the high spiking irregularity observed in vivo, most of them not in conflict with the data presented here. Many of these are based on the idea that neurons within a high conductance regime are endowed with a brief effective membrane time constant (Softky and Koch 1993; Softky 1994; Rudolph and Destexhe 2003; Leger and others 2005), rendering them highly sensitive to correlated input (Softky 1994; König and others 1996; Stevens and Zador 1998; Salinas and Sejnowski 2000), whereas other authors favored a probabilistic interpretation (Shadlen and Newsome 1994, 1998; Brunel and Wang 2003). Although correlated activity and inherently stochastic factors (like channel noise or probabilistic synaptic release; Bekkers and others 1990; White and others 1998) are likely to contribute to spiking irregularity and the large membrane potential fluctuations observed in vivo, this has often not been demonstrated in a self-consistent account that proves that the conditions giving rise to high irregularity also constitute an attracting state of the recurrent network dynamics (but see Brunel 2000). Chaos would provide such a self-consistent account (van Vreeswijk and Sompolinsky 1996, 1998).

Other Evidence for Deterministic Temporal Structure in Neural Activity

To provide strong evidence for chaos in neurophysiological preparations has in general proven difficult (Guevara and others 1981; Holden and others 1982; Freeman 1994, 2003; Hayashi and Ishizuka 1995; Schindler and others 1997; Lovejoy and others 2001; Fujisawa and others 2004; see especially Canavier and others 2004). Other evidence for deterministic structure, however, comes from studies demonstrating precise temporal spiking patterns within neural activity in vitro, in cultures, and in vivo (Abeles and others 1993; Riehle and others 1997; Tsodyks and others 1999; Grün and others 2002; Beggs and Plenz 2004; Fellous and others 2004; Ikegaya and others 2004). Spiking patterns of single or among multiple neurons may reoccur with high temporal precision aligned to sensory stimulation or significant task events (Bair and Koch 1996; Rieke and others 1996; Riehle and others 1997; Grün and others 2002; Shmiel and others 2005), and may form ‘local attractor states’ of the neural dynamics (Tsodyks and others 1999; Tiesinga and others 2002; Fellous and others 2004). Spatio-temporal activity patterns resembling those evoked by stimulation also occur spontaneously and allow to predict spiking responses of single neurons with high accuracy (Tsodyks and others 1999). Reoccurring temporal patterns imply above-chance predictability within spike trains and hence support our present findings and their potential significance for in vivo situations. Recently, reoccurring and hierarchically organized spatiotemporal patterns, susceptible to NMDA antagonists, within large-scale population signals (Ca2+ imaging, multi-electrode arrays) and synaptic input sequences have also been described in neocortical slices and cultures (Beggs and Plenz 2004; Ikegaya and others 2004).

Relation between NMDA-Induced Spiking Patterns and Activity In Vivo

NMDA synaptic inputs are believed to strongly shape recurrent dynamics in PFC networks during working memory (Wang 1999; Compte and others 2000; Durstewitz and others 2000; Durstewitz and Seamans 2002). NMDA receptors reach a particularly high density in PFC (Scherzer and others 1998), both delay period activity and working memory behavior are disrupted by prefrontal NMDA blockade (Dudkin and others 1997; Aura and Riekkinen 1999), and the maintenance of upstates in PFC cultures and in vivo crucially depends on NMDA conductances (Seamans and others 2003), as do bursting or up-state-like phenomena in vitro (Tseng & O'Donnell, 2005; see also Golomb et al., 2006 for a recent computational study). NMDA channels do not only exhibit a slow decay but are also equipped with a fast nonlinear voltage dependence (due to their Mg2+ block) that inverts the INMDA current–voltage relation in the subthreshold regime and could give rise to dendritic spikes (Schiller and others 2000). These properties seem crucial for the current's propensity to induce bursting and chaos. Given the variance of biophysical and morphological parameters of real neurons, at least some of them should reside in the critical, irregular transition regime described here during recurrent network activity in vivo. Our network simulations further suggested that such NMDA-induced chaotic irregularity present in some cells could be conveyed, to some degree, through synaptic connections to other neurons (in fact, NMDA-like synaptic connections seemed sufficient for this).

But is there also more direct evidence pointing to the occurrence of such a transition regime in vivo? One signature of this regime is the broad membrane potential distribution with 2 merging peaks. Such a distribution may well characterize activity during upstates in vivo, which are maintained through strong recurrent NMDA inputs (Seamans and others 2003) and are likely to constitute the baseline state in awake behaving animals (Steriade and others 2001; Timofeev and others 2001; Durstewitz and Seamans 2005): Within such up-states, it appears that two membrane potential levels can be further discerned (e.g. Figs. 1 & 9a in Azouz and Gray, 1999), but they are not nearly as distinct as the transitions between down- and up-states that occur during sleep or anesthesia (Stern et al., 1997; Steriade et al., 2001; Timofeev et al., 2001). An in vivo study by Léger et al., (2005; see traces in their Figs. 1 & 2) shows this quite clearly. The authors divided within-up-state activity further into two different membrane potential levels which, however, do not seem well separated, as is the case for the transition regime described here. A tendency towards burstiness with close-to or slightly bimodal ISI distributions (Compte et al., 2003a), as well as Lvs fairly above one (Shinomoto et al., 2005), as indicative of the irregular regime, have also been documented in a large number of neurons recorded during working memory delays and other tasks in vivo.

Thus, to summarize, there are some indications that the NMDA-dependent close-to-bimodal transition regime described here may also contribute to the irregular firing of PFC neurons recorded in vivo. To further test these ideas, one should 1) examine in more detail membrane potential distributions in awake animals using measures as introduced here and 2) perform in vivo experiments where spike train irregularity is studied before and after partial blockade of NMDA receptors.

Potential Functional Implications

We have suggested that an NMDA-induced chaotic transition regime may account for some of the spiking irregularity observed in vivo, but could it also serve a computational purpose? Many roles for chaos in neural computations have been suggested, including, for instance, its utilization in memory search or pattern recognition processes (Kushibe and others 1996; De Feo 2003; Nara 2003). In terms of neural coding, one advantage of settling into chaotic irregular spiking may be that on the one hand, the information content is quite high (indicated by the HISI), yet on the other hand, the series of ISIs still exhibits short-term predictability that could be exploited for coding (Davies and others 2006, provided some evidence that irregularity may increase with task complexity). Another potential advantage in the context of working memory (in the sense of ‘non-Markovian’ temporal sequence tasks where the correct choice depends on the recent history of inputs) is that complex system dynamics at the edge of chaos may strongly enhance computation on time series (Bertschinger and Natschläger, 2004). The complexity in dynamics allows for well (linearly) separable encoding of stimulus patterns, while close to a bifurcation at the edge of chaos separability will be retained across quite long periods due to long effective time constants within this regime. Only a weakly chaotic or complex limit-cycle (but in any case deterministic) dynamics, as backed by our nonlinearity test, would be required. However, whether the irregular regime observed here could play such an instrumental role rather than just being a side-effect of NMDA-induced bursting, still remains to be shown.

All slice physiology and much of the computational work were done at the Institute for Cognitive Neuroscience, Ruhr-University Bochum, Germany, funded by a research grant from the Deutsche Forschungsgemeinschaft to DD (Du 354/2-3, 2-4). Many thanks to Emili Balaguer-Ballester, Roman Borisyuk, Marc-Thorsten Hütt, Sven Kröner, Jeremy Seamans, and Thomas Wennekers for feedback and discussions on this manuscript, and to Maxim Volgushev for kindly providing us with his lab space for some of the recordings. Conflict of Interest: None declared.

References

Abeles
M
Bergman
H
Margalit
E
Vaadia
E
Spatiotemporal firing patterns in the frontal cortex of behaving monkeys
J Neurophysiol
 , 
1993
, vol. 
70
 (pg. 
1629
-
1638
)
Anderson
J
Lampl
I
Reichova
I
Carandini
M
Ferster
D
Stimulus dependence of two-state fluctuations of membrane potential in cat visual cortex
Nat Neurosci
 , 
2000
, vol. 
3
 (pg. 
617
-
621
)
Aura
J
Jr
Riekkinen P
Blockade of NMDA receptors located at the dorsomedial prefrontal cortex impairs spatial working memory in rats
Neuroreport
 , 
1999
, vol. 
10
 (pg. 
243
-
248
)
Azouz
R
Gray
CM
Cellular mechanisms contributing to response variability of cortical neurons in vivo
J Neurosci
 , 
1999
, vol. 
19
 (pg. 
2209
-
2223
)
Bair
W
Koch
C
Temporal precision of spike trains in extrastriate cortex of the behaving macaque monkey
Neural Comput
 , 
1996
, vol. 
8
 (pg. 
1185
-
1202
)
Beggs
JM
Plenz
D
Neuronal avalanches in neocortical circuits
J Neurosci
 , 
2003
, vol. 
23
 (pg. 
11167
-
11177
)
Beggs
JM
Plenz
D
Neuronal avalanches are diverse and precise activity patterns that are stable for many hours in cortical slice cultures
J Neurosci
 , 
2004
, vol. 
24
 (pg. 
5216
-
5229
)
Bekkers
JM
Richerson
GB
Stevens
CF
Origin of variability in quantal size in cultured hippocampal neurons and hippocampal slices
Proc Natl Acad Sci USA
 , 
1990
, vol. 
87
 (pg. 
5359
-
5362
)
Bernasconi
C
Schindler
K
Stoop
R
Douglas
R
Complex response to periodic inhibition in simple and detailed neuronal models
Neural Comput
 , 
1999
, vol. 
11
 (pg. 
67
-
74
)
Bertschinger
N
Natschlager
T
Real-time computation at the edge of chaos in recurrent neural networks
Neural Comput
 , 
2004
, vol. 
16
 (pg. 
1413
-
1436
)
Brunel
N
Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons
J Comput Neurosci
 , 
2000
, vol. 
8
 (pg. 
183
-
208
)
Brunel
N
Wang
XJ
Effects of neuromodulation in a cortical network model of object working memory dominated by recurrent inhibition
J Comput Neurosci
 , 
2001
, vol. 
11
 (pg. 
63
-
85
)
Brunel
N
Wang
XJ
What determines the frequency of fast network oscillations with irregular neural discharges? I. Synaptic dynamics and excitation-inhibition balance
J Neurophysiol
 , 
2003
, vol. 
90
 (pg. 
415
-
430
)
Canavier
CC
Perla
SR
Shepard
PD
Scaling of prediction error does not confirm chaotic dynamics underlying irregular firing using interspike intervals from midbrain dopamine neurons
Neuroscience
 , 
2004
, vol. 
129
 (pg. 
491
-
502
)
Chay
TR
Rinzel
J
Bursting, beating, and chaos in an excitable membrane model
Biophys J
 , 
1985
, vol. 
47
 (pg. 
357
-
366
)
Chen
G
Greengard
P
Yan
Z
Potentiation of NMDA receptor currents by dopamine D1 receptors in prefrontal cortex
Proc Natl Acad Sci USA
 , 
2004
, vol. 
101
 (pg. 
2596
-
2600
)
Clay
JR
A novel mechanism for irregular firing of a neuron in response to periodic stimulation: irregularity in the absence of noise
J Comput Neurosci
 , 
2003
, vol. 
15
 (pg. 
43
-
51
)
Cohen
S
Hindmarsh
A
CVODE, a stiff/nonstiff ODE solver in C
Comput Phys
 , 
1996
, vol. 
10
 (pg. 
138
-
143
)
Compte
A
Brunel
N
Goldman-Rakic
PS
Wang
XJ
Synaptic mechanisms and network dynamics underlying spatial working memory in a cortical network model
Cereb Cortex
 , 
2000
, vol. 
10
 (pg. 
910
-
923
)
Compte
A
Constantinidis
C
Tegner
J
Raghavachari
S
Chafee
MV
Goldman-Rakic
PS
Wang
XJ
Temporally irregular mnemonic persistent activity in prefrontal neurons of monkeys during a delayed response task
J Neurophysiol
 , 
2003
, vol. 
90
 (pg. 
3441
-
3454
)
Compte
A
Sanchez-Vives
MV
McCormick
DA
Wang
XJ
Cellular and network mechanisms of slow oscillatory activity (<1 Hz) and wave propagations in a cortical network model
J Neurophysiol
 , 
2003
, vol. 
89
 (pg. 
2707
-
2725
)
Davies
RM
Gerstein
GL
Baker
SN
Measurement of time-dependent changes in the irregularity of neural spiking
J Neurophysiol
  
in press
Dean
AF
The variability of discharge of simple cells in the cat striate cortex
Exp Brain Res
 , 
1981
, vol. 
44
 (pg. 
437
-
440
)
De Feo
O
Self-emergence of chaos in the identification of irregular periodic behavior
Chaos
 , 
2003
, vol. 
13
 (pg. 
1205
-
1215
)
Dudkin
KN
Kruchinin
VK
Chueva
IV
Effect of NMDA on the activity of cortical glutaminergic structures in delayed visual differentiation in monkeys
Neurosci Behav Physiol
 , 
1997
, vol. 
27
 (pg. 
153
-
158
)
Durstewitz
D
Self-organizing neural integrator predicts interval times through climbing activity
J Neurosci
 , 
2003
, vol. 
23
 (pg. 
5342
-
5353
)
Durstewitz
D
Seamans
JK
The computational role of dopamine D1 receptors in working memory
Neural Netw
 , 
2002
, vol. 
15
 (pg. 
561
-
572
)
Durstewitz
D
Seamans
JK
Beyond bistability: Biophysics and temporal dynamics of working memory
Neuroscience
 , 
2006
, vol. 
139
 (pg. 
119
-
133
)
Durstewitz
D
Seamans
JK
Sejnowski
TJ
Dopamine-mediated stabilization of delay-period activity in a network model of prefrontal cortex
J Neurophysiol
 , 
2000
, vol. 
83
 (pg. 
1733
-
1750
)
Fan
Y
Holden
AV
From simple to complex oscillatory behaviour via intermittent chaos in the Rose-Hindmarsh model for neuronal activity
Chaos Solitons Fractals
 , 
1992
, vol. 
2
 (pg. 
349
-
369
)
Fellous
JM
Sejnwoski
TJ
Regulation of persistent activity by background inhibition in an in vitro model of a cortical microcircuit
Cereb Cortex
 , 
2003
, vol. 
13
 (pg. 
1232
-
1241
)
Fellous
JM
Tiesinga
PH
Thomas
PJ
Sejnowski
TJ
Discovering spike patterns in neuronal responses
J Neurosci
 , 
2004
, vol. 
24
 (pg. 
2989
-
3001
)
Freeman
WJ
Neural networks and chaos
J Theor Biol
 , 
1994
, vol. 
171
 (pg. 
13
-
18
)
Freeman
WJ
Evidence from human scalp electroencephalograms of global chaotic itinerancy
Chaos
 , 
2003
, vol. 
13
 (pg. 
1067
-
1077
)
Fujisawa
S
Matsuki
N
Ikegaya
Y
Single neurons can induce phase transitions of cortical recurrent networks with multiple internal states
Cereb Cortex
 , 
2006
, vol. 
16
 (pg. 
639
-
654
)
Fujisawa
S
Yamada
MK
Nishiyama
N
Matsuki
N
Ikegaya
Y
BDNF boosts spike fidelity in chaotic neural oscillations
Biophys J
 , 
2004
, vol. 
86
 (pg. 
1820
-
1828
)
Funahashi
S
Bruce
CJ
Goldman-Rakic
PS
Mnemonic coding of visual space in the monkey's dorsolateral prefrontal cortex
J Neurophysiol
 , 
1989
, vol. 
61
 (pg. 
331
-
349
)
Fuster
JM
Unit activity in prefrontal cortex during delayed-response performance: neuronal correlates of transient memory
J Neurophysiol
 , 
1973
, vol. 
36
 (pg. 
61
-
78
)
Gabriel
T
Seamans
JK
Durstewitz
D
The basis of irregular spontaneous activity in a prefrontal cortex slice preparation. Program No. 58.2, 2004 Abstract Viewer and Itinerary Planner
2004
Washington, DC
Society for Neuroscience
Goldman-Rakic
PS
Cellular basis of working memory
Neuron
 , 
1995
, vol. 
14
 (pg. 
477
-
485
)
Golomb
D
Shedmi
A
Curtu
R
Ermentrout
GB
Persistent synchronized bursting in cortical tissues with low magnesium concentration: a modeling study
J Neurophysiol
 , 
2006
, vol. 
95
 (pg. 
1049
-
1067
)
Gorelova
N
Seamans
JK
Yang
CR
Mechanisms of dopamine activation of fast-spiking interneurons that exert inhibition in rat prefrontal cortex
J Neurophysiol
 , 
2002
, vol. 
88
 (pg. 
3150
-
3166
)
Gorelova
NA
Yang
CR
Dopamine D1/D5 receptor activation modulates a persistent sodium current in rat prefrontal cortical neurons in vitro
J Neurophysiol
 , 
2000
, vol. 
84
 (pg. 
75
-
87
)
Grün
S
Diesmann
M
Aertsen
A
Unitary events in multiple single-neuron spiking activity: I. Detection and significance
Neural Comput
 , 
2002
, vol. 
14
 (pg. 
43
-
80
)
Guevara
MR
Glass
L
Shrier
A
Phase locking, period-doubling bifurcations, and irregular dynamics in periodically stimulated cardiac cells
Science
 , 
1981
, vol. 
214
 (pg. 
1350
-
1353
)
Hansel
D
Sompolinsky
H
Chaos and synchrony in a model of a hypercolumn in visual cortex
J Comput Neurosci
 , 
1996
, vol. 
3
 (pg. 
7
-
34
)
Hayashi
H
Ishizuka
S
Chaotic responses of the hippocampal CA3 region to a mossy fiber stimulation in vitro
Brain Res
 , 
1995
, vol. 
686
 (pg. 
194
-
206
)
Hegger
R
Kantz
H
Schreiber
T
Practical implementation of nonlinear time series methods: the TISEAN package
Chaos
 , 
1999
, vol. 
9
 (pg. 
413
-
435
)
Hines
ML
Carnevale
NT
The NEURON simulation environment
Neural Comput
 , 
1997
, vol. 
9
 (pg. 
1179
-
1209
)
Holden
AV
Fan
Y
From simple to complex oscillatory behaviour via intermittent chaos in the Rose-Hindmarsh model for neuronal activity
Chaos Solitons Fractals
 , 
1992
, vol. 
2
 (pg. 
349
-
369
)
Holden
AV
Winlow
W
Haydon
PG
The induction of periodic and chaotic activity in a molluscan neurone
Biol Cybern
 , 
1982
, vol. 
43
 (pg. 
169
-
173
)
Holt
GR
Softky
WR
Koch
C
Douglas
RJ
Comparison of discharge variability in vitro and in vivo in cat visual cortex neurons
J Neurophysiol
 , 
1996
, vol. 
75
 (pg. 
1806
-
1814
)
Ikegaya
Y
Aaron
G
Cossart
R
Aronov
D
Lampl
I
Ferster
D
Yuste
R
Synfire chains and cortical songs: temporal modules of cortical activity
Science
 , 
2004
, vol. 
304
 (pg. 
559
-
564
)
Jaeger
H
Haas
H
Harnessing nonlinearity: predicting chaotic systems and saving energy in wireless communication
Science
 , 
2004
, vol. 
304
 (pg. 
78
-
80
)
Kantz
H
Schreiber
T
Nonlinear time series analysis
2004
2nd ed.
Cambridge, UK
Cambridge University Press
König
P
Engel
AK
Singer
W
Integrator or coincidence detector? The role of the cortical neuron revisited
Trends Neurosci
 , 
1996
, vol. 
19
 (pg. 
130
-
137
)
Koulakov
AA
Properties of synaptic transmission and the global stability of delayed activity states
Network
 , 
2001
, vol. 
12
 (pg. 
47
-
74
)
Kushibe
M
Liu
Y
Ohtsubo
J
Associative memory with spatiotemporal chaos control
Phys Rev E
 , 
1996
, vol. 
53
 (pg. 
4502
-
4508
)
Leger
JF
Stern
EA
Aertsen
A
Heck
D
Synaptic integration in rat frontal cortex shaped by network activity
J Neurophysiol
 , 
2005
, vol. 
93
 (pg. 
281
-
293
)
Lewis
BL
O‘Donnell
P
Ventral tegmental area afferents to the prefrontal cortex maintain membrane potential ’up' states in pyramidal neurons via D(1) dopamine receptors
Cereb Cortex
 , 
2000
, vol. 
10
 (pg. 
1168
-
1175
)
Lovejoy
LP
Shepard
PD
Canavier
CC
Apamin-induced irregular firing in vitro and irregular single-spike firing observed in vivo in dopamine neurons is chaotic
Neuroscience
 , 
2001
, vol. 
104
 (pg. 
829
-
840
)
Magee
JC
Carruth
M
Dendritic voltage-gated ion channels regulate the action potential firing mode of hippocampal CA1 pyramidal neurons
J Neurophysiol
 , 
1999
, vol. 
82
 (pg. 
1895
-
1901
)
Markram
H
Wang
Y
Tsodyks
M
Differential signaling via the same axon of neocortical pyramidal neurons
Proc Natl Acad Sci USA
 , 
1998
, vol. 
95
 (pg. 
5323
-
5328
)
Melendez
RI
Vuthiganon
J
Kalivas
PW
Regulation of extracellular glutamate in the prefrontal cortex: focus on the cystine glutamate exchanger and group I metabotropic glutamate receptors
J Pharmacol Exp Therapeutics
 , 
2005
, vol. 
314
 (pg. 
139
-
147
)
Nara
S
Can potentially useful dynamics to solve complex problems emerge from constrained chaos and/or chaotic itinerancy?
Chaos
 , 
2003
, vol. 
13
 (pg. 
1110
-
1121
)
Ninokura
Y
Mushiake
H
Tanji
J
Representation of the temporal order of visual objects in the primate lateral prefrontal cortex
J Neurophysiol
 , 
2003
, vol. 
89
 (pg. 
2868
-
2873
)
Ninokura
Y
Mushiake
H
Tanji
J
Integration of temporal order and object information in the monkey lateral prefrontal cortex
J Neurophysiol
 , 
2004
, vol. 
91
 (pg. 
555
-
560
)
Phillips
AG
Ahn
S
Floresco
SB
Magnitude of dopamine release in medial prefrontal cortex predicts accuracy of memory on a delayed response task
J Neurosci
 , 
2004
, vol. 
24
 (pg. 
547
-
553
)
Rainer
G
Asaad
WF
Miller
EK
Selective representation of relevant information by neurons in the primate prefrontal cortex
Nature
 , 
1998
, vol. 
393
 (pg. 
577
-
579
)
Reyes
A
Lujan
R
Rozov
A
Burnashev
N
Somogyi
P
Sakmann
B
Target-cell-specific facilitation and depression in neocortical circuits
Nat Neurosci
 , 
1998
, vol. 
1
 (pg. 
279
-
285
)
Riehle
F
Grun
S
Diesmann
M
Aertsen
A
Spike synchronization and rate modulation differentially involved in motor cortical function
Science
 , 
1997
, vol. 
278
 (pg. 
1950
-
1953
)
Rieke
A
Warland
D
vanSteveninck
RRD
Bialek
W
Spikes: exploring the neural code
 , 
1996
Cambridge, MA
MIT Press
Rinzel
J
Ermentrout
B
Koch
C
Segev
I
Analysis of neural excitability and oscillations
Methods in neuronal modeling
 , 
1998
Cambridge, MA
MIT Press
(pg. 
251
-
292
)
Rudolph
M
Destexhe
A
The discharge variability of neocortical neurons during high-conductance states
Neuroscience
 , 
2003
, vol. 
119
 (pg. 
855
-
873
)
Salinas
E
Sejnowski
TJ
Impact of correlated synaptic input on output firing rate and variability in simple neuronal models
J Neurosci
 , 
2000
, vol. 
20
 (pg. 
6193
-
6209
)
Sawaguchi
T
Goldman-Rakic
PS
The role of D1-dopamine receptor in working memory: local injections of dopamine antagonists into the prefrontal cortex of rhesus monkeys performing an oculomotor delayed-response task
J Neurophysiol
 , 
1994
, vol. 
71
 (pg. 
515
-
528
)
Sawaguchi
T
Matsumura
M
Kubota
K
Effects of dopamine antagonists on neuronal activity related to a delayed response task in monkey prefrontal cortex
J Neurophysiol
 , 
1990
, vol. 
63
 (pg. 
1401
-
1412
)
Scherzer
CR
Landwehrmeyer
GB
Kerner
JA
Counihan
TJ
Kosinski
CM
Standaert
DG
Daggett
LP
Velicelebi
G
Penney
JB
Young
AB
Abstract expression of N-methyl-D-aspartate receptor subunit mRNAs in the human brain: hippocampus and cortex
J Comp Neurol
 , 
1998
, vol. 
390
 (pg. 
75
-
90
)
Schiller
J
Major
G
Koester
HJ
Schiller
Y
NMDA spikes in basal dendrites of cortical pyramidal neurons
Nature
 , 
2000
, vol. 
404
 (pg. 
285
-
289
)
Schindler
KA
Bernasconi
CA
Stoop
R
Goodman
PH
Douglas
RJ
Chaotic spike patterns evoked by periodic inhibition of rat cortical neurons
Z Naturforsch Sect A-A J Phys Sci
 , 
1997
, vol. 
52
 (pg. 
509
-
512
)
Schreiber
T
Schmitz
A
Improved surrogate data for nonlinearity tests
Phys Rev Lett
 , 
1996
, vol. 
77
 (pg. 
635
-
638
)
Schreiber
T
Schmitz
A
Surrogate time series
Physica D
 , 
2000
, vol. 
142
 pg. 
346
 
Seamans
JK
Durstewitz
D
Christie
BR
Stevens
CF
Sejnowski
TJ
Dopamine D1/D5 receptor modulation of excitatory synaptic inputs to layer V prefrontal cortex neurons
Proc Natl Acad Sci USA
 , 
2001
, vol. 
98
 (pg. 
301
-
306
)
Seamans
JK
Nogueira
L
Lavin
A
Synaptic basis of persistent activity in prefrontal cortex in vivo and in organotypic cultures
Cereb Cortex
 , 
2003
, vol. 
13
 (pg. 
1242
-
1250
)
Shadlen
MN
Newsome
WT
Noise, neural codes and cortical organization
Curr Opin Neurobiol
 , 
1994
, vol. 
4
 (pg. 
569
-
579
)
Shadlen
MN
Newsome
WT
The variable discharge of cortical neurons: implications for connectivity, computation, and information coding
J Neurosci
 , 
1998
, vol. 
18
 (pg. 
3870
-
3896
)
Shampine
LF
Reichelt
MW
The MATLAB ODE suite
SIAM J Sci Comput
 , 
1997
, vol. 
18
 (pg. 
1
-
22
)
Shilnikov
A
Calabrese
RL
Cymbalyuk
G
Mechanism of bistability: tonic spiking and bursting in a neuron model
Phys Rev E Stat Nonlin Soft Matter Phys
 , 
2005
, vol. 
71
 pg. 
056214
 
Shinomoto
S
Shima
K
Tanji
J
Differences in spiking patterns among cortical neurons
Neural Comput
 , 
2003
, vol. 
15
 (pg. 
2823
-
2842
)
Shinomoto
S
Miyazaki
Y
Tamura
H
Fujita
I
Regional and laminar differences in in vivo firing patterns of primate cortical neurons
J Neurophysiol
 , 
2005
, vol. 
94
 (pg. 
567
-
575
)
Shmiel
T
Drori
R
Shmiel
O
Ben Shaul
Y
Nadasdy
Z
Shemesh
M
Teicher
M
Abeles
M
Neurons of the cerebral cortex exhibit precise interspike timing in correspondence to behavior
Proc Natl Acad Sci USA
 , 
2005
, vol. 
102
 (pg. 
18655
-
18657
)
Shuai
J
Bikson
M
Hahn
PJ
Lian
J
Durand
DM
Ionic mechanisms underlying spontaneous CA1 neuronal firing in Ca2+-free solution
Biophys J
 , 
2003
, vol. 
84
 (pg. 
2099
-
2111
)
Softky
W
Sub-millisecond coincidence detection in active dendritic trees
Neuroscience
 , 
1994
, vol. 
58
 (pg. 
13
-
41
)
Softky
WR
Koch
C
The highly irregular firing of cortical cells is inconsistent with temporal integration of random EPSPs
J Neurosci
 , 
1993
, vol. 
13
 (pg. 
334
-
350
)
Steriade
M
Timofeev
I
Grenier
F
Natural waking and sleep states: a view from inside neocortical neurons
J Neurophysiol
 , 
2001
, vol. 
85
 (pg. 
1969
-
1985
)
Stern
EA
Kincaid
AE
Wilson
CJ
Spontaneous subthreshold membrane potential fluctuations and action potential variability of rat corticostriatal and striatal neurons in vivo
J Neurophysiol
 , 
1997
, vol. 
77
 (pg. 
1697
-
1715
)
Stevens
CF
Zador
AM
Input synchrony and the irregular firing of cortical neurons
Nat Neurosci
 , 
1998
, vol. 
1
 (pg. 
210
-
217
)
Strogatz
SH
Nonlinear dynamic and chaos
 , 
1994
Reading, MA
Addison-Wesley
Tateno
T
Jimbo
Y
Robinson
HPC
Spatio-temporal cholinergic modulation in cultured networks of rat cortical neurons: spontaneous activity
Neurosci
 , 
2005
, vol. 
134
 (pg. 
425
-
437
)
Terman
D
The transition from bursting to continuous spiking in excitable membrane models
J Nonlinear Sci
 , 
1992
, vol. 
2
 (pg. 
135
-
182
)
Tiesinga
PHE
Fellous
JM
Sejnowski
TJ
Attractor reliability reveals deterministic structure in neuronal spike trains
Neural Comp
 , 
2002
, vol. 
14
 (pg. 
1629
-
1650
)
Timofeev
I
Grenier
F
Steriade
M
Disfacilitation and active inhibition in the neocortex during the natural sleep-wake cycle: an intracellular study
Proc Natl Acad Sci USA
 , 
2001
, vol. 
98
 (pg. 
1924
-
1929
)
Traub
RD
Buhl
EH
Gloveli
T
Whittington
MA
Fast rhythmic bursting can be induced in layer 2/3 cortical neurons by enhancing persistent Na+ conductance or by blocking BK channels
J Neurophysiol
 , 
2003
, vol. 
89
 (pg. 
909
-
921
)
Tseng
KY
O'Donnell
P
Dopamine-glutamate interactions controlling prefrontal cortical pyramidal cell excitability involve multiple signaling mechanisms
J Neurosci
 , 
2004
, vol. 
24
 (pg. 
5131
-
5139
)
Tseng
KY
O'Donnell
P
Post-pubertal emergence of prefrontal cortical up states induced by D1-NMDA co-activation
Cereb Cortex
 , 
2005
, vol. 
15
 (pg. 
49
-
57
)
Tsodyks
M
Kenet
T
Grinvald
A
Arieli
A
Linking spontaneous activity of single cortical neurons and the underlying functional architecture
Science
 , 
1999
, vol. 
286
 (pg. 
1943
-
1946
)
Tsodyks
MV
Markram
H
The neural code between neocortical pyramidal neurons depends on neurotransmitter release probability
Proc Natl Acad Sci USA
 , 
1997
, vol. 
94
 (pg. 
719
-
723
)
van Vreeswijk
C
Sompolinsky
H
Chaos in neuronal networks with balanced excitatory and inhibitory activity
Science
 , 
1996
, vol. 
274
 (pg. 
1724
-
1726
)
van Vreeswijk
C
Sompolinsky
H
Chaotic balanced state in a model of cortical circuits
Neural Comput
 , 
1998
, vol. 
10
 (pg. 
1321
-
1371
)
Wang
XJ
Synaptic basis of cortical persistent activity: the importance of NMDA receptors to working memory
J Neurosci
 , 
1999
, vol. 
19
 (pg. 
9587
-
9603
)
Wang
J
O'Donnell
P
D1 dopamine receptors potentiate NMDA-mediated excitability increase in layer V prefrontal cortical pyramidal neurons
Cereb Cortex
 , 
2001
, vol. 
11
 (pg. 
452
-
462
)
Watanabe
M
Kodama
T
Hikosaka
K
Increase of extracellular dopamine in primate prefrontal cortex during a working memory task
J Neurophysiol
 , 
1997
, vol. 
78
 (pg. 
2795
-
2798
)
White
JA
Klink
R
Alonso
A
Kay
AR
Noise from voltage-gated ion channels may influence neuronal dynamics in the entorhinal cortex
J Neurophysiol
 , 
1998
, vol. 
80
 (pg. 
262
-
269
)
Yang
CR
Seamans
JK
Dopamine D1 receptor actions in layers V-VI rat prefrontal cortex neurons in vitro: modulation of dendritic-somatic signal integration
J Neurosci
 , 
1996
, vol. 
16
 (pg. 
1922
-
1935
)
Yang
CR
Seamans
JK
Gorelova
N
Electrophysiological and morphological properties of layers V-VI principal pyramidal cells in rat prefrontal cortex in vitro
J Neurosci
 , 
1996
, vol. 
16
 (pg. 
1904
-
1921
)
Zheng
P
Zhang
X-X
Bunney
BS
Shi
W-X
Opposite modulation of cortical N-methyl-D-aspartate receptor-mediated responses by low and high concentrations of dopamine
Neuroscience
 , 
1999
, vol. 
91
 (pg. 
527
-
535
)