Cerebral Cortex Advance Access originally published online on June 1, 2006
Cerebral Cortex 2007 17(4):894-908; doi:10.1093/cercor/bhk044
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Dynamical Basis of Irregular Spiking in NMDA-Driven Prefrontal Cortex Neurons
Centre for Theoretical and Computational Neuroscience, University of Plymouth, Plymouth, UK
Address correspondence to Dr Daniel Durstewitz, Centre for Theoretical and Computational Neuroscience, University of Plymouth, Portland Square, A220, Plymouth, PL4 8AA, UK. Email: daniel.durstewitz{at}plymouth.ac.uk.
 |
Abstract |
|---|
 Top Abstract IntroductionMaterials and MethodsResultsDiscussionReferences |
|---|
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.
Key Words: chaos • neural coding • computation • bursting • slice electrophysiology • biophysical model
|
Introduction |
|---|
|
TopAbstractIntroductionMaterials and MethodsResultsDiscussionReferences |
|---|
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 |
|---|
|
TopAbstractIntroductionMaterials and MethodsResultsDiscussionReferences |
|---|
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 x40 water immersion objective. Thick-walled borosilicate patch pipettes (2–6 M
) were filled with (in mM) K-gluconate (115), 4-(2-hydroxyethyl)-1-piperazineethanesulfonic acid (10), MgCl2 (2), KCl (20), MgATP (2), Na2-ATP (2), and Guanosine 5'-triphosphate (0.3) (pH = 7.3, 285 ± 5 mOsm). Ag/AgCl wires and pellets (Warner, Hamden, CT) were used for the patch and reference electrodes, respectively. Whole-cell current-clamp recordings were performed with an Axoclamp 2B amplifier (Axon Instr., Foster City, CA) and acquired with Labview-based software, written by Lee Campbell (Salk Institute, La Jolla, CA), at sampling rates of 2–5 kHz. Pipette series resistance was <20 M, and was compensated. Electrode potentials were adjusted to zero before recording, without correcting the liquid junction potential.
NMDA (7–8.5 µM; Sigma-Aldrich, Munich, Germany; a concentration which is within a physiological range, see Melendez et al., 2005) and the selective D1 agonist SKF 38393 (2–5 µM; Tocris, Ellisville, MO) were bath applied within the ACSF. In some experiments, 12–20 µM DNQX (Sigma-Aldrich) and 24–30 µM bicuculine (Sigma-Aldrich) were added to block 
-Amino-3-hydroxy-5-methylisoxazole-4-propionic acid (AMPA) and 
-aminobutyric acid A (GABAA) currents, respectively. In the cells probed with current steps for control (see Results), additionally all NMDA currents were blocked by 2-Amino-5-phosphonopentanoic acid (APV) (50 µM, Sigma-Aldrich). All drugs were mixed into ACSF from frozen stock solution aliquots. All experiments were conducted at 34–36 °C in ACSF at a perfusion rate of 2–3 ml/min.
Analysis Methods
Recordings from single cells consisted of one to several samples of 30–300 s length, separated by brief nonrecording periods. To judge stationarity both within and between samples of the same cell, the original voltages traces and the instantaneous firing rate (=1/interspike-interval [ISI]) as a function of time were first inspected by eye. Traces were then divided into consecutive segments of 15 or 30 s (depending on total recording duration, and down to 3 s for the NMDA-omission controls [see Results]) for which ISI distributions together with their mean and standard deviation were plotted (for all segments overlaid within the same graph). All ISI distributions were compared through Kolmogorov–Smirnov tests, probing for significant differences between any 2 distributions. For none of the samples judged as stationary, the overall corrected (for number of statistic comparisons) P value fell below 10%, indicating that all ISI distributions belonged to the same process. Note that this stationarity requirement implies stability of the spiking statistics across the entire length of recording periods used for analysis, unconfounded by potential drifts caused by intracellular washout etc.
For quantifying bimodality in the membrane potential distributions, spikes were first cut out within the range [tsp – 2 ms, tsp + 4 ms], for all tsp which denote the times of crossing a –20-mV threshold from below. Next, the sum of 2 Gaussian-like exponentials was fitted to the distribution, which had the general form
 |
µ2. So whereas one exponential (with µ2 as mean) is a pure Gaussian, the other exponential (with µ1 as mean) was allowed different slopes right versus left from the mean, to account for the fact that unimodal and left-hand parts of bimodal Vm distributions were often strongly skewed (see e.g., Fig. 1B bottom). Fits with 6 free parameters were performed in 2 different ways. For the first type of fit, k = 2 was fixed. For the second type of fit, µ1 = µ2 was enforced, but k was allowed to vary in the range [0, 3]. Hence, the second type of fit was unimodal but with the same number of free parameters as the bimodal fit. Fitting was performed by constrained optimization using MatLab's (The MathWorks, Inc., Natick, MA) lsqcurvefit-routine. In general, whichever of these 2 fits produced a smaller (absolute) error was accepted. Only in <10% of cases, a clearly uni- or bimodal Vm distribution was not detected as such by this algorithm and had to be enforced "by hand." Based on these fits, a parameter dV quantifying the degree of bimodality within the membrane potential distribution was defined as dV = |µ1 – µ2|/
2.
|
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:
- 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 = –
n P(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.
Tbs, and all measures of irregularity were first computed separately for all samples and then averaged across all samples within a stationary set.
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) |
ISIn) for n = 1 and U is the neighborhood of size of this vector (|·| denotes its cardinality; for each point was iteratively adjusted such that at least 5 neighbors could be found). Temporal neighbors < 10 ISIs away from the ISI vector under consideration were excluded, as these may simply pick up temporal correlations not due to reoccurring deterministic structure. A local prediction error is then calculated as squared discrepancy between predicted and actual future values and averaged across all N points (with kmax future values) in m-dimensional state space:
 | (2) |
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) |
k are independent Gaussian random numbers (k = 0, 
= 1) and f is a monotonic but possibly nonlinear function (for details, see Schreiber and Schmitz 2000). Briefly, surrogates that obey this null hypothesis are created by "rescaling" original values to a Gaussian distribution (empirically inverting f), performing a fast Fourier transform (after reducing discontinuities between end points to minimize periodicity artifacts), randomizing the phases (but retaining the amplitudes), transforming back into the time domain, and rescaling to the original distribution (because of the rescaling, some "polishing" may be necessary and distribution and power spectrum may have to be adjusted iteratively in several steps; see Schreiber and Schmitz 2000). The resulting data sets have the same power spectrum (autocorrelation function) and distribution as the original ISI series but destroy any nonlinear dependencies and hence, predictability that would result specifically from nonlinear deterministic processes. For performing the predictions and creating the surrogate data, the TISEAN 2.1 software package (Hegger and others 1999; Schreiber and Schmitz 2000) was used, where the prediction routine "zeroth" was adapted to accommodate multiple temporally discontinuous samples. Computational Methods
Networks of 2-compartment pyramidal cells (PCs), equipped with 6 different ionic conductances (INa, INaP, IDR, IKS, IHVA, IC), and single-compartment interneurons (IN), equipped with just fast Na+ and DR K+ channels, were constructed. Gating variables for all ionic conductances followed Hodgkin-Huxley–like channels kinetics. The basic model cells and synapses were taken from Durstewitz and Seamans (2002), and details regarding parameters and equations can be found there and in Durstewitz and others (2000). Neurons were connected through AMPA- and NMDA-like excitatory (gAMPAmax = 0.02 µS, gNMDAmax = gAMPAmax/50), or GABAA-like inhibitory (gGABAAmax = 0.006µS) connections, respectively (see Durstewitz and Seamans 2002). The NMDA current is given by
 | (4a) |
 | (4b) |
that implements the time course of all recurrent NMDA inputs to the cell (parameters as in Durstewitz and Seamans 2002) and s(Vm) implements the very fast nonlinear voltage dependence (made instantaneous here, without any impact on results). All synaptic conductances were multiplied by a short-term plasticity factor (absorbed into gNMDA(t) in eq. 4a) that developed according to the formalism given in Tsodyks and Markram (1997; see also eq. 6 in Durstewitz 2003). PC
PC connections exhibited short-term depression (Markram and others 1998; utilization USE = 0.3, recovery from depression 
rec = 0.8 s, facilitation fac = 0.8 s; for details, see Tsodyks and Markram 1997), INPC (USE = 0.25, rec = 0.6 s, fac = 1 s) and ININ (USE = 0.12, rec = 0.6 s, fac = 1.2 s) connections weak short-term depression (Reyes and others 1998), and PCIN (USE = 0.05, rec = 0.15 s, fac = 1.2 s) connections short-term facilitation (Markram and others 1998).
The standard network consisted of 100 PCs and 25 INs, with a random 20% connectivity between different neurons. For the basic observations reported here, precise connectivity and network size did not seem of major relevance. Constant (i.e., time independent) NMDA activation was added to the dendritic compartment of all PCs and to INs. This conductance had the same voltage dependence s(Vm) and reversal potential as the synaptic NMDA conductance (eq. 4), but the product of gNMDAmax and the time-dependent term gNMDA(t) modeling the rise and decay of the conductance after receptor activation was replaced by a constant, gNMDAcons (with a mean of 0.095 mS/cm2 for the network simulations reported here). The inclusion of this conductance was meant to mimic the continuous application of NMDA to the bath solution in the in vitro experiments. D1 application was mimicked by doubling the inactivation time constant of INaP, reducing IKS everywhere by 60%, reducing dendritic IHVA by 15%, and reducing gleak of IN by 25% (Yang and Seamans 1996; Gorelova and Yang 2000; Gorelova and others 2002), with regards to the "standard" parameters given in Durstewitz and Seamans (2002). We stress, however, that such D1-like manipulations are not important for the basic phenomena reported here (similar results were obtained without them).
Many parameters followed a Gaussian distribution within the network, with the following standard deviations given as percentage of reference (mean) values in brackets: Somatic and dendritic INaP, IKS, IHVA, and IC maximum conductances (all 20%), size of PCs and INs (with all morphological dimensions changed by the same factor with 20% standard deviation), gNMDAcons (20%), and all synaptic conductances (50%). Network simulations were run with NEURON (Hines and Carnevale 1997).
Single neuron simulations were also performed with NEURON and were often double-checked with MatLab-based code, using a model cell picked from the network. Single cell results were further confirmed by drawing at random 10 cells from the network distribution and in model cells without D1-like changes. NEURON was used with its implicit variable time-step solver CVODE (Cohen and Hindmarsh 1996), and for the MatLab simulations, the implicit multi-step solver ode15s (Shampine and Reichelt 1997) was taken. Simulation results were probed with relative tolerances of down to 10–8 and absolute tolerances of 10–10. All simulation code is available from the following URL: www.pion.ac.uk/
daniel.
For the single PC model, 2-dimensional (Vsoma, dVdend/dt) state space projections were obtained analytically. For given Vsoma, all somatic gating variables and Ca2+ and K+ ionic concentrations were solved for their steady-state values. By setting dVsoma/dt = 0, the dendritic membrane potential can be obtained as Vdend = 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
(v) versus Vslow([Ca2+]i) = [Ca2+] ([Ca2+]i), and Vslow(v) versus Vslow(h) = h (h) (where –1 denotes the inverse function). This way, bifurcation graphs could be obtained (with fixed points and their stability determined similarly as described above) as functions of just a single parameter, Vslow, cutting through the 3-dimensional parameter space (vfixed, [Ca2+]fixed, hfixed) in an optimal manner with reference to the full model behavior. |
Results |
|---|
|
TopAbstractIntroductionMaterials and MethodsResultsDiscussionReferences |
|---|
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).
|
Computational Model
We constructed networks of 100 2-compartment conductance-based PCs and 25 single-compartment IN, randomly connected through AMPA-, NMDA-, and GABAA-like synapses, respectively (the precise connectivity seemed not crucial for the present observations). All neurons received an additional source of constant NMDA receptor activation (time independent, but voltage-(Mg2+)-gated NMDA conductances), supposed to mimic the constant application of NMDA to the bath solution. D1-like changes in cell parameters (Yang and Seamans 1996; Gorelova and others 2002) with regard to a standard reference model (Durstewitz and others 2000; Durstewitz and Seamans 2002) were also included but were not crucial for the principal phenomena to occur. Model neurons varied with regard to their morphological dimensions and strength of several voltage-gated conductances, as well as the strength and number of synaptic inputs, and the constant NMDA activation (reflecting natural variation among neurons). Importantly, there was no source of noise or randomness in the running model network.
The whole range of phenomena as observed in vitro could be reproduced in the network, from long-duration bursts, to highly irregular spiking, to regular continuous spiking (Fig. 7). Note that the high irregularity in the model cells must be due to deterministic chaos because there was no noise source in the simulations. Like in the in vitro preparations, the highest irregularity occurred in a regime where a bimodal distribution of membrane potential just started to form (compare Fig. 7 middle with top and bottom). Figure 8 furthermore illustrates that for these model cells similar ISI return maps and evidence for nonlinear predictability are obtained as in the real neurons (compare with Fig. 3). Finally, removing AMPA and GABAA synapses from the network, highly irregular firing and long-duration bursts could still be observed, confirming the presence of those phenomena in the in vitro preparations with AMPA/GABAA currents blocked (not shown). (Interestingly, the remaining NMDA connections were also sufficient to synchronize bursts between the bursting neurons.)
|
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.
|
), Vsoma) trajectory (gray solid line) of the full model projected onto this plane. (A) Nonlinear NMDA current: With the NMDA nonlinearity in place, a broad hysteresis region appears (from approximately –58.7 to –47 mV) where the limit cycle and the lower stable fixed point coexist, and the trajectory cycles a few times (the burst) before crawling back below the middle unstable arch toward the next bursting cycle. (B) Linear NMDA current: The hysteresis region is dramatically reduced, and the trajectory leaves the limit cycle region already during the first spike, giving rise to a simple single-spike train.The crucial point about this graph (Fig. 11A) is the existence of a broad hysteresis region, that is, the coexistence of the limit cycle with a stable resting point over the Vslow range approximately –58.7 to –47 mV. The bursting process in the full model may now be explained as follows: During spiking within a burst, IHVA and INaP slowly inactivate, whereas the Ca2+ concentration slowly rises and enhances IC. Due to this increasing loss of depolarizing force, the spiking process breaks down at some point and the (Vslow(v), Vsoma) trajectory (gray line in Fig. 11A), projected onto the (Vslow, Vsoma) graph of the reduced model, passes over to the right into the region where the lower branch of fixed points constitutes the only stable solution (Vslow > –47 mV). Due to the hysteresis region, the trajectory then crawls back within the basin of the lower branch of stable fixed points (due to reactivation of IHVA and INaP and slow drop in [Ca2+]i), which accounts for the interburst interval, until it pops onto the limit cycle again as the stable resting solution vanishes in a saddle node bifurcation where the stable and unstable branches meet (at approximately –58.7 mV). The remarkable consequence of such a (square-wave-type) bursting process is that chaotic solutions between the bursting and regular single-spiking regimes have been shown to be a generic property of such systems under rather general assumptions (Terman 1992
). Figure 11B depicts the situation with a linear NMDA current (s fixed): The hysteresis region is dramatically diminished, in fact so much that Vslow(v) passes over to the monostable resting region during the upstroke of the first spike and bursting is no longer supported. Hence, the cell engages only in regular single spiking after removal of the NMDA nonlinearity.
|
Discussion |
|---|
|
TopAbstractIntroductionMaterials and MethodsResultsDiscussionReferences |
|---|
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