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.
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.
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
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.
The coefficient of variation, Cv = σISI/μISI, where μISI and σISI denote the mean and standard deviation of the ISI distribution.
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.
(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.
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
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
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
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 = Vsoma −ri[Iinj + ∑ gi(Ei− Vsoma)], 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.
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.
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).
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).
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).
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.)
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).
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.
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.
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 (PEISI/σISI) 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 PEISI/σ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 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.