Dynamical Basis of Irregular Spiking in NMDA-Driven Prefrontal Cortex Neurons

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.

chaos, neural coding, computation, bursting, slice electrophysiology, biophysical model

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), NaHCO₃ (26), glucose (10), KCl (3.5), NaH₂PO₄ (1.25), CaCl₂ (0.5), MgCl₂ (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 CaCl₂ increased to 2 mM and MgCl₂ 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% O₂/5% CO₂ 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), MgCl₂ (2), KCl (20), MgATP (2), Na₂-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 (GABA_A) 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 [t_sp − 2 ms, t_sp + 4 ms], for all t_sp 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

where h_V denotes the V_m distribution and μ₁ ≤ μ₂. So whereas one exponential (with μ₂ as mean) is a pure Gaussian, the other exponential (with μ₁ 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 V_m 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, μ₁ = μ₂ 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 V_m distribution was not detected as such by this algorithm and had to be enforced “by hand.” Based on these fits, a parameter d_V quantifying the degree of bimodality within the membrane potential distribution was defined as d_V = |μ₁ − μ₂|/σ₂.

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 = V_m distribution, black line = function fit, see Materials and Methods). (C) Index for the amount of bimodality within the membrane potential distribution (d_V), various indices for the amount of spiking irregularity (C_v, C_vL, L_v, H_ISI), and average duration of bursts (〈T_bs〉 in s).

The average duration of bursts (〈T_bs〉) 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 V_m distributions, the average burst duration was set to 0.

Four different measures for the irregularity of the spike trains were obtained:d_V, 〈T_bs〉, 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, C_v = σ_ISI/μ_ISI, where μ_ISI and σ_ISI denote the mean and standard deviation of the ISI distribution.

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

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

4. H_ISI = − ∑nP(ISI ∈ B_n) log₂P(ISI ∈ B_n) with bins B₀ = [0, 2] ms and B_n = [(1 − α) c_n, (1 + α) c_n] for n > 0. The bin centers (c_n) were determined as c_n = c_n−1 (1 + α)/(1 − α) (with α = 0.1). This procedure ensures that bin size increases proportionally to the bin center (i.e., |B_n|/c_n = const. = 2α) because the ISI standard deviation can be expected to scale with the mean (but other choices for H_ISI yielded similar results). H_ISI 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 H_ISI 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)

where ISI_n is the m-dimensional ISI vector under consideration with components (ISI_n−m+1 ⋯ ISI_n) 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 ≥Δk_max future values) in m-dimensional state space:

(2)

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)

where the η_k are independent Gaussian random numbers (〈η_k〉 = 0, 〈η_k²〉 = 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 (I_Na, I_NaP, I_DR, I_KS, I_HVA, I_C), 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 (〈g_AMPAmax〉 = 0.02 μS, 〈g_NMDAmax〉 = 〈g_AMPAmax〉/50), or GABA_A-like inhibitory (〈g_GABAAmax〉 = 0.006μS) connections, respectively (see Durstewitz and Seamans 2002). The NMDA current is given by

(4a)

with

(4b)

where g_NMDA(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(V_m) 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 g_NMDA(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 U_SE = 0.3, recovery from depression τ_rec = 0.8 s, facilitation τ_fac = 0.8 s; for details, see Tsodyks and Markram 1997), IN→PC (U_SE = 0.25, τ_rec = 0.6 s, τ_fac = 1 s) and IN→IN (U_SE = 0.12, τ_rec = 0.6 s, τ_fac = 1.2 s) connections weak short-term depression (Reyes and others 1998), and PC→IN (U_SE = 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(V_m) and reversal potential as the synaptic NMDA conductance (eq. 4), but the product of g_NMDAmax and the time-dependent term g_NMDA(t) modeling the rise and decay of the conductance after receptor activation was replaced by a constant, g_NMDAcons (with a mean of 0.095 mS/cm² 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 I_NaP, reducing I_KS everywhere by 60%, reducing dendritic I_HVA by 15%, and reducing g_leak 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 I_NaP, I_KS, I_HVA, and I_C maximum conductances (all 20%), size of PCs and INs (with all morphological dimensions changed by the same factor with 20% standard deviation), g_NMDAcons (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⁻⁸ and absolute tolerances of 10⁻¹⁰. All simulation code is available from the following URL: www.pion.ac.uk/∼daniel.

For the single PC model, 2-dimensional (V_soma, dV_dend/dt) state space projections were obtained analytically. For given V_soma, all somatic gating variables and Ca²⁺ and K⁺ ionic concentrations were solved for their steady-state values. By setting dV_soma/dt = 0, the dendritic membrane potential can be obtained as V_dend = V_soma −r_i[I_inj + ∑ g_i(E_i− V_soma)], where r_i is the internal (axial) resistance, the g_i are the total conductances (including gating), and E_i are the corresponding reversal potentials. Using V_dend, steady-state values for dendritic gating variables and ionic concentrations could be computed and the curve dV_dend/dt = F(V_soma) (for which dV_soma/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 Ca²⁺ current I_HVA, the dendritic [Ca²⁺]_i concentration which controls the Ca²⁺-dependent activation of the K⁺ current I_C, and the voltage-dependent inactivation (h) of the dendritic persistent Na⁺ current I_NaP (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 V_slow: v_fixed = v_∞(V_slow), h_fixed = h_∞(a₁V_slow + b₁), and [Ca²⁺]_fixed = [Ca²⁺]_∞(a₂V_slow + b₂), where the steady-state value of [Ca²⁺]_i is a function of voltage via Ca²⁺ influx through I_HVA. Note that while v is set to its steady-state value for V_slow directly, the other 2 steady-state values are computed for linear functions of V_slow. The parameters a_i and b_i 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 V_slow (v) = v_∞⁻¹ (v) versus V_slow([Ca²⁺]_i) = [Ca²⁺]_∞⁻¹ ([Ca²⁺]_i), and V_slow(v) versus V_slow(h) = h_∞⁻¹ (h) (where ⁻¹ 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, V_slow, cutting through the 3-dimensional parameter space (v_fixed, [Ca²⁺]_fixed, h_fixed) 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 (C_v), the L_v (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 C_v even when spiking is otherwise highly regular, a local version of the C_v (C_vL) that likewise minimizes contributions of long ISIs between bursts, and the entropy (H_ISI) of the ISI-distribution with bin sizes proportional to the interval centers (for details, see Materials and Methods). As shown in Figure 1C, C_vL, L_v, and H_ISI all attain a maximum for the irregular spiking pattern, while the C_v is closest to 1 (the C_v 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 (d_V) 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 L_v and H_ISI attain a maximum (while the standard C_v is closest to 1, not shown) for a d_V range indicative of V_m distributions that just start to separate into 2 discernable peaks. Figure 2C gives the dependence of C_vL on d_V for all stationary recording sets as dots, to better visualize the variance in the data (clearly unimodal distributions always have d_V = 0, clustering along the ordinate). Analyses of variance comparing the 3 groups d_V < 1, 1 ≤ d_V < 4, and d_V ≥ 4 yielded highly significant main effects (all F_2,105 > 53, all P < 10⁻¹⁵) 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 d_V, revealing that brief-duration bursts just start to appear for the range of d_V values for which the irregularity is highest (〈T_bs〉 ∼ 164 ms for 1 ≤ d_V < 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 L_v as a measure of spike train irregularity attains a maximum for small to intermediate d_V values, indicating V_m distributions at the border of bimodality. (B) The entropy of the ISI distribution (H_ISI) is also maximal for the same range of d_V values. (C) Local C_vL for all stationary samples (n = 108, from 62 cells) as a function of d_V. (D) The duration of bursts is highly but not perfectly linearly correlated with d_V. For d_V 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 (ISI_n+1, ISI_n) 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 ISI_n- and ISI_n+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 ISI_n+1 = F(ISI_n), 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) ISI_n+1 as a function of preceding ISI_n: 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 (PE_ISI, 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 PE_ISI/σ_ISI 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 GABA_A 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 d_V values within an intermediate range, as indicative of V_m distributions at the edge of bimodality. Again, all post hoc comparisons for the d_V groups as defined above were significant for all 4 irregularity indices, except for H_ISI, where the intermediate and the high d_V group were statistically en par (possibly partly due to the lack of cells with d_V > 7 in this sample), although there was, once again, a tendency for cells with intermediate d_Vs to carry the highest entropies (not shown).

Figure 4.

Same as Figure 1 for PFC slices after blockade of AMPA (DNQX: 12 μM) and GABA_A (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 (d_V), various indices for the amount of spiking irregularity (C_v, C_vL, L_v, H_ISI), and average duration of bursts (〈T_bs〉 in s).

Figure 5.

Summary graphs for PFC slices after blockade of AMPA and GABA_A synapses. Both the L_v (A) and the local C_vL (B) still peak for d_V 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 (C_v) versus average ISI for all steady-state traces from 8 recorded PFC neurons. Except for one trace, the C_vs are all very low, indicating highly regular single-spiking activity. The one trace with a higher C_v (∼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 L_v and C_vL, 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 GABA_A-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-(Mg²⁺)-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 GABA_A 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/GABA_A 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 (d_V), various indices for the amount of spiking irregularity (C_v, C_vL, L_v, H_ISI), and average duration of bursts (〈T_bs〉 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).

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 I_HVA). 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 (g_NMDAcons 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 dV_dend/dt as a function of V_soma (black line, labeled dV_soma/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 dV_dend/dt curve crosses the dV_dend/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 g_NMDAcons. (A) Nonlinear NMDA current: For most values of g_NMDAcons, the cell exhibits only a few different ISIs and hence fires periodically with different periods, but for some g_NMDAcons range (∼0.103–0.114 mS/cm²), 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 [Ca²⁺]_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 dV_soma/dt = 0 in Fig. 9B bottom), assuming that all gating variables and ionic concentrations are in their steady states, and that dV_soma/dt = 0 (i.e., V_soma = const.). Consequently, all points where in addition dV_dend/dt = 0 mark the fixed points of the whole system. Within these 2-dimensional spaces spanned by V_soma and dV_dend/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 g_NMDAcons (not shown) and, above some critical g_NMDAcons, 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 g_NMDAcons, 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 V_m distribution and the merging of its 2 peaks as evidenced in Figures 1A, 4A, and 7A (middle). With further increases of g_NMDAcons, 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 V_m 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 V_m distribution with merging peaks (intermediate d_V values). It is also interesting to note that a hyperpolarizing current injection would lower the dV_soma/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 Mg²⁺ 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 Mg²⁺-dependent voltage-gate s (eq. 4b) to a constant value (we set s_fixed = 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 g_NMDAcons axis). Figure 10B shows that after removal of the NMDA nonlinearity, the model neuron passes from silence directly to regular spiking with increasing g_NMDAcons 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 Mg²⁺-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 Ca²⁺ current I_HVA, the dendritic [Ca²⁺]_i concentration (steering Ca²⁺-dependent activation of the K⁺ current I_C), and the inactivation gate (h) of the dendritic persistent Na⁺ current I_NaP to fixed values parameterized through a single voltage variable denoted V_slow (the parameterization was done such that the resulting 1-dimensional cut through the 3-dimensional parameter space (v, [Ca²⁺]_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 V_slow (higher values corresponding to stronger inactivation of I_HVA and I_NaP and stronger activation of I_C), 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 V_slow) 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 V_slow, as well as the lower and upper edge of the limit cycle (dotted) and the (V_slow(ν), V_soma) 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 V_slow 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, I_HVA and I_NaP slowly inactivate, whereas the Ca²⁺ concentration slowly rises and enhances I_C. Due to this increasing loss of depolarizing force, the spiking process breaks down at some point and the (V_slow(v), V_soma) trajectory (gray line in Fig. 11A), projected onto the (V_slow, V_soma) 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 (V_slow > −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 I_HVA and I_NaP and slow drop in [Ca²⁺]_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 V_slow(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 V_m 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 Ca²⁺ 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 (PE_ISI/σ_ISI) was significantly lower for the transition regime (1 ≤ d_V < 4) than for the steady spiking (d_V < 1) or bursting (d_V ≥ 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 PE_ISI/σ_ISI-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 GABA_A 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 (Ca²⁺ 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 Mg²⁺ block) that inverts the I_NMDA 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 L_vs 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 H_ISI), 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.