| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
Institute of Physiology, University of Bern, 3012 Bern, Switzerland
Submitted 26 March 2003; accepted in final form 7 May 2003
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIXDISCLOSURESACKNOWLEDGMENTSREFERENCES |
|---|
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIXDISCLOSURESACKNOWLEDGMENTSREFERENCES |
|---|
; McCormick et al. 1985; Rhodes 1999). Even the most detailed state-of-the-art model is unable to capture the entire phenomenology observed in the experiments. Such a richness calls for a model reduction that could provide a synthetic description of the response properties of the cells under particular conditions. We studied in vitro those features that are supposedly relevant when the cell is embedded in a large network of interconnected neurons, as it would be in in vivo conditions. The guidelines for selecting the relevant features were dictated by the theoretical framework developed in the last decade to study the dynamic properties of networks of integrate-and-fire neurons (see Gerstner and Kistler 2002 for a review). This approach was already successful in relating single neuron properties to several in vivo phenomena (Amit and Tsodyks 1991; Brunel 2000b) like the omnipresent spontaneous activity (Amit and Brunel 1997) or the persistent, selective delay activity observed in many areas of the cortex in behaving animals (Amit and Brunel 1997; Wang 2001; Yakovlev et al. 1998). In both cases the recorded spike activity is sustained throughout long intervals: low, spontaneous activity is always present, whereas elevated delay activity can last 
30s (Fuster 1995). During these intervals the statistics of the total synaptic input are likely to be stationary or quasi-stationary, that is, to vary on time scales that are much longer than the "reaction" time (Gerstner and Kistler 2002) of the assembly of neurons.
Despite the steadiness of the statistics, the total synaptic current results from a considerably irregular synaptic activity, and hence fluctuates all the time. Neurons on the presynaptic side emit spikes spontaneously at a frequency of a few spikes s–1 and the neocortical connectivity is rather high [approximately 104 synapses per neuron (Abeles 1991)]. As a consequence, during the interval between two successive spikes, every neuron integrates hundreds of excitatory and inhibitory postsynaptic potentials that arrive at random times. If the spike activities of the presynaptic neurons are statistically independent (see METHODS) then the resulting total synaptic conductance can be described as a random walk with a Gaussian distribution and finite time correlation length determined by the time development of unitary synaptic events [an Ornstein–Uhlenbeck process (Tuckwell 1988)]. Interestingly in vitro neocortical neurons produce in vivo–like activity when noisy current waveforms are injected (Destexhe et al. 2001; Mainen and Sejnowski 1995), and there is some preliminary indirect evidence, based on the analysis of the distribution of the membrane potential recorded intracellularly in vivo, that somatic currents could be modeled as an Ornstein–Uhlenbeck process (Destexhe and Paré 1999, 2000). The total synaptic current is also approximately Gaussian (Amit and Tsodyks 1992) (also see DISCUSSION and the APPENDIX) and, for fast synaptic currents (AMPA or GABAA), the time correlation length is short compared with other inherent time constants of the neuron (e.g., the membrane time constant). As a consequence the total input current is practically white noise and can be fully characterized by its average mI and SD sI.
In the present work we measured the response function of neocortical pyramidal cells to such a noisy current with stationary statistics. We first studied the time development of the neuronal response and characterized the functional role of the adaptation/facilitation components that determine the statistical properties of stationary responses. These components act on different time scales and, for some input currents, reach a steady regime in which they either disappear or are combined together to determine the stationary statistics of the train of spikes generated by the neuron. We then analyzed quantitatively the stationary responses by exploring systematically the whole parameter space {mI, sI} characterizing the input current and determining the resulting spike frequency f. The measured frequencies have been fitted by the theoretical response functions of two simple models of integrate-and-fire neurons with a single, effective component of adaptation/facilitation. These response functions have been computed analytically (Fusi and Mattia 1999; Ricciardi 1977) and have a relatively simple form. The value of this data modeling is multiple: besides providing a synthetic and efficient way of describing the whole data set, it allows reliable prediction of the output rate in response to currents not used in the experiment and, hence, to cover a wide class of experiments that study the response of the neuron when moving along some specific trajectory in the parameter space of the input current (see e.g., Chance et al. 2002). Moreover, the choice of the response function of model neurons instead of an arbitrary function gives a direct interpretation of the estimated parameters: they are the effective parameters of a model neuron that can re-create the measured response of pyramidal cells. The knowledge of these parameters also allows one to make quantitative predictions about the global dynamics of a network (in vivo) of this kind of cells. In addition, the response function to inputs with stationary statistics also gives important information about the reaction time of the network (Fourcaud and Brunel 2002; Gerstner and Kistler 2002; Mattia and DelGiudice 2002). If the parameters are tuned to reproduce the mean spike frequencies, the simulated neurons can also re-create the higher-order statistics of the interspike intervals expressing, for instance, the degree of irregularity of the spike train (e.g., the coefficient of variability of the interspike intervals). Finally, the statistics of the effective parameters across different cells provide a quantitative estimate of the heterogeneity of the functional properties of the cells.
|
METHODS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIXDISCLOSURESACKNOWLEDGMENTSREFERENCES |
|---|
Parasagittal slices of rat somatosensory cortex (300 µm thick) were prepared from 15- to 40-day-old female and male Wistar rats according to the institutional guidelines. The preparation was done in ice-cold extracellular solution using a Campden vibratome (752M; Campden Instruments, Loughborough, UK). Slices were incubated at 35°C for 25 min and afterward left at room temperature until being transferred to the recording chamber. The cells were visualized by infrared differential interference contrast videomicroscopy using a Newvicon camera (C2400, Hamamatsu City, Japan) and an infrared filter (RG9, Schott Mainz, Germany) mounted on an upright microscope (Axioscope FS, Zeiss, Germany).
We recorded in current-clamp whole cell configuration from the soma of layer 5 regular spiking (McCormick et al. 1985) pyramidal cells. Recordings and stimulations were made with an Axoclamp-2A amplifier (Axon Instruments, Burlingame, CA) in combination with Clampex 8 (Axon Instruments). The access resistance and the capacitance were compensated using the bridge balance and the capacitance neutralization after having established the whole cell configuration. The data were low-pass filtered at 2.5 kHz with sampling frequency twice the filter frequency. The temperature of the external solution was 31°C. Neurons were visually identified and some of them were filled with biocytin and then stained according to the avidin–biotin–peroxidase (ABC) procedure (Hsu et al. 1981).
Slices were continuously superfused with an artificial cerebrospinal fluid containing (in mM): 125 NaCl, 25 NaHCO3, 2.5 KCl, 1.25 NaH2PO4, 2 CaCl2, 1 MgCl2, 25 glucose, gassed with 95% O2-5% CO2.
The pipette solution for most of the analyzed cells contained (in mM): 110 K-gluconate, 30 KCl, 10 EGTA, 10 HEPES, 4 Mg-ATP, 0.3 Na2-GTP, 10 Na2-phosphocreatine, pH adjusted to 7.3 with KOH. This solution contains a relatively high concentration of EGTA, which allowed long stable recordings and made the cells produce more consistent responses (see Stimulation protocol and observables). We tried to identify undesired artifacts introduced by EGTA by studying the cell's response when two other pipette solutions were used. In what follows we will refer to the solution with high concentration of EGTA as the EGTA pipette solution. The other two pipette solutions were labeled KMeSO4 and KGluc and contained (in mM): KMeSO4: 135 K-methylsulfate, 20 KCl, 0.08 EGTA, 0.045 CaCl2, 10 HEPES, 4 Mg-ATP, 0.3 Na2-GTP, 10 Na2-phosphocreatine, pH adjusted to 7.3 with KOH. KGluc: 115 K-gluconate, 20 KCl, 10 HEPES, 4 Mg-ATP, 0.3 Na2-GTP, 10 Na2-phosphocreatine, pH adjusted to 7.3 with KOH. Pipette solution KGluc was always used with 10 mM biocytin. The measured osmolarity of all three pipette solutions was between 310 and 325 mOsm. Because pyramidal cells of the somatosensory cortex in vitro showed virtually no spontaneous activity, we did not systematically block synaptic input mediated by ligand-gated channels. Control experiments with blocked synaptic inputs by adding 50 µM D-APV, 10 µM CNQX, and 10 µM bicuculline to the extracellular solution did not change the spike frequency.
The input resistance of the neurons was calculated from the voltage transients in response to at least three different hyperpolarizing (600-ms duration, average of the last 300 ms) current pulses (amplitude, 0.05 nA). The membrane time constant 
m was estimated by injecting brief (0.5 ms) hyperpolarizing current pulses (–2.5 nA) into the soma. From the decaying averaged (n = 50) voltage transient after this current pulse, m was obtained from the slope of a straight line fitted through the tail portion of the semilogarithmic plot of the membrane voltage against time (Iansek and Redman 1973).
Model of in vivo–like input current
We assume that a large number of presynaptic neurons emit spikes at random times. On the postsynaptic side this heavy barrage is felt as a total synaptic current I that evolves as a random walk. If the synaptic currents are summed linearly and the different inputs are statistically independent (i.e., the emission of a spike by one presynaptic neuron does not affect or only slightly affects the initiation of an action potential in another presynaptic neuron) then the random walk can be replaced by a smoother version, the Ornstein–Uhlenbeck process (Tuckwell 1988), characterized by a Gaussian distribution (mean mI, SD sI) and by a time-correlation length I. If the synaptic evoked potentials sum linearly and there are Ne excitatory AMPA receptors, each activated at a mean rate fe and Ni inhibitory GABAA receptors (fi), then
 | (1) |
 | (2) |
ampa = gaba = I). 
e (i) is the average peak postsynaptic current evoked by the arrival of a single excitatory (inhibitory) spike. The rise time to this peak current is zero, and then the current decays exponentially: I(t) = e,ie–t/1. (See the APPENDIX for more details on how these constants are related to the mean synaptic conductances.) Slow NMDA components can be added to mI only, given that they leave sI almost unaffected (Brunel and Wang 2001). The total synaptic current I(t) can be generated with a single equation
 | (3) |
). Stimulation protocol and observables
The noisy input current was generated as an Ornstein–Uhlenbeck stochastic process by iterating the following expression
 | (4) |
t is a unitary Gauss distributed random variable, updated at every time step. The process was generated and injected at a rate of 5 kHz (
t = 0.2 ms) and the correlation length I was usually 1 ms (for 10 cells we also used I = 5 ms). The resulting current I(t) has a stationary Gauss distribution with mean mI = µII and variance 
(Cox and Miller 1965).
The space {mI, sI} was systematically explored as follows: data points were collected at fixed sI (ranging from 0 to 500 pA), stepwise increasing mI from a subthreshold value up to nonstationary frequencies. This protocol was used to determine the threshold mean current (the rheobase current). Then, the whole space {mI, sI} was discretized and then explored in random order to prevent correlations between time and one of the two parameters mI, sI. In both protocols, the duration of the stimulation depended on the cell response: it was 10 s long, or shorter if 
150 spikes were collected. For those cells used to characterize the response over time (see RESULTS), a stimulus duration of 10 s was always used. The first transient part of the neuronal response (2 s if the stimulation time was longer than 4 s; 0.5 s otherwise) was discarded when estimating the mean spike frequency and the coefficient of variability (see following text). The intervals between successive stimulations were 50–60 s (Fig. 1).
|
The mean spike frequency (response) was estimated as the ratio between the total number of action potentials Nsp and the stimulus duration T. The confidence intervals (68%) of the experimentally measured frequencies were given by = (f+ + f–)/2, with (see, e.g., Meyer 1965)
 | (5) |
Particular care was taken to ensure that the response of the cell was consistent throughout the whole recording session. Usually the cells showed a transient phase at the beginning, followed by a long time interval (10–90 min) during which the response was consistent (when the same current was injected, the spike frequency did not change within the statistical error of Eq. 5), and by a final unstable phase. Cells with a stable period shorter than 40 min were not included in the analysis. Response consistency depended on the pipette solution. The cells were classified into three groups depending on the level of consistency: 1) consistent cells: less than a third of the repeated recordings were out of range; 2) poorly consistent cells: more than a third of the repeated recordings were out of range; and 3) clearly inconsistent cells: errors as for poorly consistent cells plus there were clear inconsistencies at low frequencies (e.g., the frequency decreased with sI for some points and increased for others). The number of consistent or poorly consistent cells was (out of total number): 41/44 for EGTA pipette solution, 14/19 for KGluc pipette solution, 6/12 for KMeSO4 pipette solution.
Model of the integrate-and-fire neurons
The subthreshold dynamics of the only dynamic variable, the membrane voltage V, obey
 |
and initiates an action potential. After the spike emission, Vr is reset to some value Vr from which the neuron starts again integrating the input current after a refractory period r. For the leakage L(V) we studied two models: 1) the classical leaky integrate-and-fire neuron (LIF), which has a leakage proportional to the membrane depolarization L(V) = –VC/ ( is the membrane time constant); and 2) the constant leakage integrate-and-fire neuron with a floor (CLIFF), which has a constant leakage L(V) = –
. For this last neuron, to have the same behavior as for the LIF neuron, it is necessary to impose a rigid lower boundary for the depolarization, which we chose to be the cell resting potential (see Response functions of the model neurons). Adaptation/facilitation components
The stationary response results from the combination of several dynamic processes occurring on different time scales. Some of them are likely to be transient and to disappear in a few hundreds of milliseconds (see RESULTS). Others are long-lasting and might reach a steady state on time scales of the order of seconds. When the effects of long-lasting processes merge together to produce a neuronal response with stationary statistics, we model them by introducing an additional current I
proportional to the mean output spike frequency f. This current is meant to imitate any combination of adaptation and facilitation processes that determine the neuronal response in stationary condition. The linear dependency on f was suggested by the discrepancies observed between the best-fit theoretical response without adaptation/facilitation and the data. Moreover it is supported by the experiments that focus on specific components of frequency-dependent modifications of the input current. For example the calcium-dependent potassium current, which is responsible for fast adaptation, is usually modeled as a negative current proportional to the intracellular calcium concentration [Ca2+]i, which in turn is proportional to the spike frequency f. A possible implementation in terms of a detailed spike driven dynamics is as follows (see also Ermentrout 1998; Fuhrmann et al. 2002; Liu and Wang 2001 for related models): calcium (or whatever ion species is responsible for the phenomenon; see Powers et al. 1999; Sanchez-Vives et al. 2000) concentration is increased on every spike emission and then decays exponentially to its resting value
 | (6) |
Ca >> 1/f), then [Ca2+]i 
Caf and the current turns out to be proportional to the spike count in some temporal window, I = –
Caf. Any other model that, in stationary conditions, produces a negative current proportional to the mean spike rate would be equivalent. Theoretical response functions
If the input current I is Gauss distributed and delta-correlated, then the equation for V can be written as
 |
is a unitary Gauss-distributed variable that is updated every time step. For simplicity we focus on the CLIFF neuron, the same considerations apply to the LIF neuron. µ and 
2 are the instantaneous mean and variance of the change in V(t) per unit time and characterize the statistics of V on short time scales (see, e.g., Cox and Miller 1965). For a delta-correlated input current the variance of V grows linearly with time on short time scales and is proportional to 2. For such an input current, the response function can be computed analytically for both models (Fusi and Mattia 1999; Ricciardi 1977) and the expressions for the average firing rate f = 
(mI, sI) is reported in Table 1. The CV of the interspike intervals can also be computed analytically. The expression can be found in Brunel (2000a) for the LIF neuron and in Fusi and Mattia (1999) and Salinas and Sejnowski (2002) for the CLIFF neuron.
|
The equations of Table 1 express the mean spike frequency f as a function of the mean µ and variance 2 in unit time of the input current (erf = error function). These parameters are related to mI, sI and the time correlation length I as indicated in the bottom row: these expressions are based on the assumption that the input current has only fast synaptic components and hence the time correlation length I is much shorter than the typical interspike interval. Indeed when the current of Eq. 3 is injected, the variance of V is (Cox and Miller 1965)
 |
t >> scales as 
, as expected from a delta-correlated process. For a more detailed analysis of the effects of time correlated currents on the response function see Brunel and Sergi (1998) and Fourcaud and Brunel (2002). These effects can be reabsorbed in an effective spike-emission threshold , which depends on I.
The effects of adaptation/facilitation on the stationary response are introduced as an extra current (mI 
mI – f) proportional to the mean spike frequency. For each combination of mI, sI, the mean output frequency determines the amount of negative current that should be added to mI and that modifies the output frequency. The process is iterated until the equation for the frequency becomes self-consistent and the frequency that appears in the negative current is the same as the output frequency
 |
The asymptotic stationary response function and the response function with = 0 for the two model neurons are plotted in Fig. 5 and discussed in Response functions of the model neurons.
|
|
The theoretical f – mI curves have been fitted to the stationary data. The fit was achieved through a Monte Carlo minimization (see e.g., Press et al. 1992) of the distance between the measured mean spike frequency 
and the theoretical spike frequency 
predicted by the model
 | (7) |
= {r, Vr, C, }, plus m or in case of the LIF or CLIFF neuron, respectively. k represents the confidence intervals defined in Stimulation protocol and observables. Because only 5 of the 6 parameters of the model neurons are independent, we set the threshold arbitrarily at = 20 mV. In fact, both response functions (see Table 1) are invariant under the scaling 
, Vr Vr, C C/, with > 0.
The minimum 
of Eq. 7 with respect to the parameters set follows approximately a 
2-distribution with N – M degrees of freedom, where N is the number of experimental points, and M = 5 is the number of free parameters. The fit was accepted whenever the probability 
was greater than 0.1 for consistent cells, and greater than 0.01 for poorly consistent cells.
In some cases we also report the absolute discrepancy d, defined as the average (across all points) absolute difference between the measured and the theoretical frequencies of the best-fit curves. This number is not correlated with the goodness-of-fit test and gives a useful indication of the error made in considering a theoretical response function which does not strictly pass the least-squares test.
|
RESULTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIXDISCLOSURESACKNOWLEDGMENTSREFERENCES |
|---|
We measured the cell mean spike frequency in response to noisy currents with stationary statistics (see Fig. 1 and METHODS for details about the protocol). We identified at least three components of adaptation/facilitation, already known in the literature, which determine the features of the stationary response. What follows in this section is not meant to be an extensive and systematic analysis of all the factors that determine the stationary response. The goal is rather to clearly define what we mean by stationary response, to understand when it is possible to have it, and to summarize euristically all the observable mechanisms that affect the transient on different time scales and that might contribute to determine the stationary response. Such stationary response was usually reached in 1–2 s and then sustained until the end of the stimulation. For strong enough injected currents the cells were unable to sustain the elevated activity imposed by the stimulation and no stationary response was possible. The three components are illustrated in Fig. 2. Their features depend on the pipette solution and are summarized in Table 2. They are as follows.
|
|
): fast and invariably present for all the different preparations, it manifests itself for high enough spike frequencies: when the neuron is injected with a constant current, the second and successive interspike intervals are longer than the first one (see Fig. 2, a3–a5). This component, known in the literature to be attributed to calcium-dependent potassium current (see e.g., McCormick et al. 1985; Sah 1996), reaches a steady state in a few spikes and it is usually modeled as negative spike frequency–dependent current, which clearly affects the stationary response. This adaptation component was present for all the cells and all the preparations and in particular it was observed also in the strong presence of EGTA, a calcium chelator. Buffering of calcium concentration transients by EGTA is probably not fast enough to disrupt fast initial adaptation (Smith et al. 1984). The spike-frequency dependency was different depending on the pipette solution. The same effect was produced for lower currents in the case of EGTA pipette solution when compared with the case of KGluc pipette solution: 1.2 ± 0.4 nA (KGluc pipette solution), 0.6 ± 0.2 nA (EGTA pipette solution), corresponding to 36 ± 10 Hz (EGTA pipette solution) and 21 ± 8 Hz (KGluc pipette solution). For KMeSO4 pipette solution: 0.8 ± 0.2 nA and 24 ± 7 Hz. ) and its effects on the stationary response compete with those of the late adaptation component (described below) when the spike frequency is high enough. For this reason it is difficult to isolate its effects and to characterize its dependency on the preparation. However, it was clear that early facilitation was mostly prominent when KMeSO4 pipette solution was used and this would be compatible with the known fact (Zhang et al. 1994) that methylsulfate is the least disruptive to intracellular structures of calcium homeostasis. ; Powers et al. 1999; Sanchez-Vives et al. 2000; Sawczuk et al. 1997; Schwindt et al. 1989) and it is hypothesized to be attributed to slow inactivation of sodium channels (Fleidervish et al. 1996; Powers et al. 1999; Sawczuk et al. 1997) and to an outward Na+-dependent K+ current (Sanchez-Vives et al. 2000; Schwindt et al. 1989). The spike frequency decay was always accompanied by a progressive drift of the maximal upstroke velocity of the membrane potential, indicating that the spike shape was continuously and slowly broadening. This phenomenon is best displayed by using longer stimuli that elicit high initial spike rates. An example is shown In Fig. 3B), where 60-s-long stimuli were used.
|
Many of the recorded cells (19/44) showed an initial doublet of spikes, resembling the beginning of a burst. This doublet sometimes masks the initial adaptation component and makes it difficult to analyze. Notice that all the recorded cells were selected not to be inherently bursting. The doublet was observed for all the preparations for strong enough stimulation (usually stronger than the one that reveals fast spike adaptation) and it is shown in Fig. 2, a4–a5. The threshold current depended on the pipette solution: 1.2 nA for the 4/13 in EGTA pipette solution, about 1.8 nA for 13/19 cells in KGluc pipette solution and for 2/12 cells in KMeSO4 pipette solution. Because of its transitory nature, the mechanism responsible for this doublet probably does not affect the stationary response.
Maximal stationary response
The late adaptation component makes the cell unable to sustain a response with stationary statistics above a critical output spike frequency, which depended on the cell and on the pipette solution. Although the observed modifications in the spike mean frequency and in the spike shape were usually quite substantial at the end of long stimulations, all cells could recover during the interstimulus intervals, indicating that the phenomenon is reversible.
For most cells we did not stimulate for very long intervals (60 s) as for the cell shown in Fig. 3. Hence we do not have enough data to determine the maximal stationary response for all cells. This analysis would have required the injection of several currents, strong enough to produce nonstationary responses, and it would have limited even further the number of stationary data points. Nonstationary responses are not included in the following analysis and are not modeled.
Usually we restricted the input currents to the range in which the slow frequency decay was below 1 Hz/s. This criterion restricted our analysis of the response function to a limited range of frequencies that depended on the pipette solution. The maximal frequency of this range could be determined for those cells for which the stimulation was strong enough to reveal a frequency decay above 1 Hz/s and it was (see Table 2): 44 ± 12 Hz (n = 11/13) for EGTA pipette solution and 56 ± 9 Hz (n = 5/19) for KGluc pipette solution. In KMeSO4 pipette solution the cells did not display maximal stationary response with the exception of only one out of 12, in which case the maximal stationary frequency was 64 Hz.
Experimentally measured response functions
In Fig. 4 we show the response function as measured in the experiment. The measurements are shown in the form of f – mI curves that represent the output frequency f (the response) as a function of the mean current mI at constant SD sI. The CV of the interspike interval is also reported in the same form. The qualitative behavior was the same for all cells: for constant currents (sI = 0), there is obviously no activity for average inputs that are below the rheobase current (i.e., the minimal constant current that makes the neuron fire), whereas the response curve is approximately linear for supra-rheobase currents. For noisy inputs (sI > 0) there are essentially two regimes. 1) A sub-rheobase, fluctuation-dominated regime in which the mean current mI is below the rheobase current and hence not sufficient to drive the membrane voltage across the threshold for emitting a spike. In this case the action potentials are sporadically initiated by fluctuations and hence the spike activity is very irregular and the CV is high (Fig. 4C). 2) A supra-rheobase, drift-dominated regime, in which the membrane potential fluctuates around a rising ramp that leads to the emission threshold at a more regular pace (Fig. 4B). The introduction of noise permits sub-rheobase activity, which in turn smooths out the response curves at the rheobase. The curves, convex for weak input currents, sometimes bend at high frequencies (see e.g., Fig. 6, B and C) and hence change curvature.
|
Response functions of the model neurons
The data points have been fitted by the response functions of two simple model neurons sharing a similar qualitative behavior (see METHODS). The two models differ in the form of the leakage L(V): for the first, classical leaky integrate-and-fire (LIF) neuron the leakage is proportional to V [L(V) = –VC/m, where m is the membrane time constant], whereas for the second model—the CLIFF neuron—the leakage is constant [L(V) = –] and the input current is integrated linearly. To obtain two qualitatively similar response functions for the two models, it is essential to limit the membrane potential from below in the case of the CLIFF neuron (Fusi and Mattia 1999; Salinas and Sejnowski 2002). Indeed, when both neurons are injected with a sub-rheobase noisy current, the variance of the fluctuations of the membrane voltage tends to increase linearly with time. However, for the LIF neuron the leakage compensates for this tendency, and, depending on the parameters of the current, can lead to a stationary sub-rheobase regime in which V fluctuates around an asymptotic value. For the CLIFF neuron this is not possible because the current is integrated linearly and all the fluctuations accumulate: for negative net currents (mean synaptic current minus leakage) the mean voltage decreases and the fluctuations increase boundlessly. In such a sub-rheobase regime the average spike frequency of the CLIFF neuron is always zero (Fusi and Mattia 1999; Gerstein and Mandelbrot 1964). The rigid barrier that limits the membrane voltage from below permits the neuron to fluctuate steadily around a value near the lower barrier, waiting for a fluctuation that drives V across the threshold. In this case a sub-rheobase regime with nonzero frequency is possible.
In both neurons the adaptation/facilitation components which determine the asymptotic stationary response were modeled by assuming that the effective mean current driving the cell is reduced by a term proportional to the cell's own spike frequency (mI mI – f; see METHODS). There are several simple and detailed biophysical models that correspond to this phenomenological model: the simple models based on calcium dependent potassium channels (see e.g., Wang 1998) and the more realistic Hodgkin–Huxley type model in Traub and Miles (1991) are but two examples (see also Ermentrout 1998). All these models exploit the fact that the adaptation/facilitation processes are slow compared with the average interspike intervals. Note that these mechanisms do not account for the nonstationarity observed in Fig. 3B, given that our model neuron can always have a stationary response, no matter what the intensity of stimulation.
The response functions corresponding to the two models are plotted in Fig. 5 in the same way as the measured curves are shown in Fig. 4. The main difference between the two model neurons is in the dependency of the spike frequency on sI, the fluctuations amplitude of the noisy current: the CLIFF neuron is more sensitive to large sI and this fact manifests itself in a progressively increasing distance between the f – mI curves. This behavior can be exposed by plotting the spike frequency as a function of sI at constant mI (insets in Fig. 5): the curvature has a different sign for the two neurons. Moreover the response function of the LIF neuron without adaptation/facilitation components is highly nonlinear around the rheobase current for low levels of noise and, for sI = 0, the derivative of the f – mI curve is infinite. These nonlinearities are attenuated by the frequency-dependent feedback (Ermentrout 1998).
The CV of the interspike intervals behave in the same way for the two model neurons: in the sub-rheobase fluctuations dominated regime the CV is close to 1 for a wide range of input currents. This corresponds to the maximal irregularity and the train of spikes has Poisson statistics. As the current drives the neuron toward a drift-dominated regime, the CV approaches 0 with a speed that depends on sI (i.e., the fluctuations in the input current).
Fitting the theoretical response functions to the data
Both models could faithfully reproduce most of the cells stationary responses, at least for those neurons that had consistent enough responses to allow for a quantitative analysis (see METHODS). For each cell we fitted the theoretical response functions to the data points by searching the space of the 5 independent parameters characterizing each model: the capacitance C, the refractory periodr, the reset potential Vr, the adaptation/facilitation constant , and the time membrane constant m for the LIF neuron and C, r, Vr, , and the leakage for the CLIFF neuron (see METHODS). We adopted a Monte Carlo technique that minimizes the mean square distance between the data points and the predicted values, each point being weighted by the inverse of the confidence interval. The 2 test passed either with one model or with the other (or with both) with a probability P > 0.1 for 27 consistent cells out of 37. For all these cells, the EGTA pipette solution was used. For other pipette solutions the consistency was usually poor and the 2 test passed with a probability P > 0.01 for 16 cells out of 24.1 The responses could anyway be well approximated by the theoretical functions even when the models did not pass the test: the average discrepancy between the measured and the model frequencies were <2.5 Hz; <1.5 Hz if frequencies below 50 Hz were considered. That is, the average difference between a single data point and its theoretical match was below 3 Hz even when the test did not pass. Examples of fitted response functions are shown in Fig. 6. The LIF model was best suited to reproduce f – mI curves that were almost equally spaced (Fig. 6B), whereas the CLIFF model could capture better those cells that were less sensitive to low levels of noise as in Fig. 6A.
Adaptation/facilitation components of the model dynamics turned out to be essential to fit the response of the models to the data. For the LIF neuron it is necessary to linearize the response function at the rheobase. For the CLIFF neuron, which already has an almost linear f – mI curve at the rheobase, one might wonder if the fitting would have been possible without adaptation/facilitation components. Increasing the membrane capacity C as well as increasing reduces the slope of the f – mI curves. However, in contrast to adaptation/facilitation, an increased capacity also reduces the effect of fluctuations on the spike frequency and hence diminishes the sensitivity to sI. Thus for the CLIFF neuron, adaptation/facilitation are the only mechanisms that can change the slope of the f – mI curves almost without affecting the distance between f – mI curves that correspond to different values of sI (La Camera et al. 2002).
For high frequencies (>40–50 Hz) there is a preliminary evidence of a phenomenon that is not captured by our model neurons: the f – mI curves corresponding to different sI have a slight tendency to diverge, as if the sensitivity to sI would increase for strong enough stimulation (see, e.g., Fig. 6C). This is not captured by our models of integrate-and-fire neurons for which asymptotically the f – mI curves always tend to merge. This divergence effect is a source of discrepancy that usually does not compromise the 2 test, but that clearly indicates the activation of an extra process that is not modeled. The increased sensitivity to sI is usually accompanied by an increase in the CV of the interspike intervals (see the previous section) and it might be attributable to the activation of calcium spikes in the distal dendrite (Larkum et al. 1999). The study of this phenomenon will require further investigation.
Some cells (n = 6, all of them in EGTA pipette solution) could be fitted only by searching a region of limit parameter values, letting the reset potential be very close to the threshold. This was the only way to capture sets of f – mI curves that do not tend to converge to a common asymptotic value for elevated currents (Fig. 6D). Although the response function could be described by the theoretical response function, we do not consider the model as appropriate for describing these cells. The theoretical response functions are computed under the assumption that the current is delta-correlated, and this is not a good approximation when the reset potential is very close to the threshold. Hence we consider these cells belong to a different class that the simple integrate-and-fire models are probably inappropriate to describe.
We recorded the response of 10 cells to noisy currents with a longer time correlation length (I = 5 ms). The f – mI curves were qualitatively the same as in the case of shorter I and they could be fitted by the theoretical response functions of Table 1 (not shown), even if they are valid only for delta-correlated currents. A full account of the analysis of these cells will be reported elsewhere.
Higher-order statistics of the interspike intervals
The main goal of this investigation was to study the neuronal response in terms of mean spike rates. However, it is also interesting to consider the higher-order statistics of the distribution of the interspike intervals. A detailed and systematic analysis would require more data to estimate the higher-order moments of the distribution, although our preliminary analysis of the variance of the interspike intervals already gives interesting indications. In Fig. 7 we plotted, for two cells, the CV of the interspike intervals as a function of the spike frequency. Each curve is obtained by keeping sI fixed and then sweeping along mI. The parameters of the model are tuned to fit only the mean firing frequencies. Then, with the best-fit parameters we computed the CV by means of simulations and compared it to the measured one.
|
The cell in the top row of Fig. 7 represents the typical case, whereas the one below is our best fit. The agreement between the predictions and the data points is in general reasonably good for most of the points and for most of the cells, even if it would probably not pass a 2 test. The largest discrepancies are observed at high firing frequencies, for elevated values of the variance of the input current. There are at least two reasons for these discrepancies. First, these points might not be completely stationary, given that the firing frequency approaches the maximal allowed value for the cell. As expected, this produces a higher CV than predicted by the theoretical models for which the nonstationarity is not modeled. Second, at high frequencies and high sI we often observed doublets or bursts of spikes (see the bottom left inset in Fig. 7). This is also not captured by the model and increases the measured CV.
Moreover there is often a large discrepancy around the rheobase current, when the amount of noise in the current is small. Near the threshold for spike emission, the response sometimes consisted of sequences of tonic, regular firing interrupted by periods of subthreshold activity (see the bottom right inset in Fig. 7). This behavior, not captured by a simple integrate-and-fire model, artificially increases the variability of a train of spikes that otherwise would be rather regular.
Effective parameters
The parameters corresponding to the best fit must be considered as effective parameters, that is, as those parameters that provide the best description of the stationary response function. Other observables might not be captured by the same parameters (also see DISCUSSION). Fitting the f – mI curves for three different values of sI was already restricting the model parameters to a small region of the parameter space (we had 3 to 5 f – mI curves for each cell). The range in which single parameters can vary and the 2 test still passes depends on the sensitivity of the spike frequency on the different parameters. A rough estimate indicates that the least determined parameters are r (up to about ±70% error) and Vr (about ±25% error), whereas the other parameters can vary at most in intervals of the order of ±20% (, CLIFF neuron), ±5% (C, , , LIF neuron; C, , CLIFF neuron). The refractory period is clearly the effective parameter to which the response function is least sensitive because it affects mostly the very high-frequency region (f 
1/r), usually beyond the observed range of output rates.
The effective parameters for the two model neurons are reported in Table 3 for all the cells that were consistent or poorly consistent and that did pass the 2 test (the number of cells is reported in column N). For the analysis of effective parameters we mainly focused on the cells with EGTA pipette solution because they were stable enough to give consistent responses for a large number of data points. The parameters are reported as an average across all cells with EGTA pipette solution and separately for the other two pipette solutions. Interestingly, the differences between different pipette solutions were usually comparable to the differences across cells (reported as SD).
|
Two passive parameters characterizing the cell (m, C) could be measured directly (see METHODS), and their values were in general different from the effective parameters. The cross-correlations were negligible for the CLIFF neuron and weak for the LIF neuron: the correlation between the directly measured capacitance and the effective capacitance of the LIF neuron was 0.59 (0.03 for the CLIFF neuron), and the correlation between the membrane time constants was 0.30. Because real neurons are not integrate-and-fire neurons, such a mismatch is not very surprising. Indeed, the effective parameters provide a good fit to the mean spike frequencies, while the passive membrane parameters are estimated from very different observables (the subthreshold response to short pulses), in different experimental conditions. No clear patterns between effective or directly measured parameters and the model that could fit the data emerged.
|
DISCUSSION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIXDISCLOSURESACKNOWLEDGMENTSREFERENCES |
|---|
; Amit and Tsodyks 1991; Chance et al. 2002; Fusi and Mattia 1999) or when a network of interacting neurons has to respond fast to a time varying input (Rudolph and Destexhe 2001; Silberberg, Bethge, Markram, Tsodyks, and Pawelzik, unpublished observations, 2002).
The agreement between the theoretical and the measured response function is remarkably good for a wide range of statistics of the input currents, on the whole {mI, sI} space and for all the sustainable spike frequencies of the analyzed cells. Hence, our results cover a variety of physiological scenarios that can be studied by measuring the neural response when moving along specific trajectories in the {mI, sI} space. For instance balanced synaptic input, as studied e.g., in Chance et al. 2002, would correspond to increasing the variability of the current sI, while keeping the mean current mI fixed, and the outcome could be predicted by studying the theoretical response function of integrate-and-fire neurons.
Moreover the knowledge of the response functions to in vivo–like currents allows one to study many dynamic properties of networks of interacting neurons and to relate single neuron properties to the activity observed in in vivo experiments. For instance, imposing the stability of a network state, like spontaneous activity (Amit and Brunel 1997), can constrain the possible architecture of the network, thus providing an indirect measurement of the synaptic efficacies in vivo. This will require a more extensive analysis of different types of cells, in different layers and different areas. The analysis of the dynamic states of a network of neurons will also require an additional study of the dependency of the parameters mI, sI of the current on the mean output rate of the presynaptic neurons. This dependency can be quite complicated, involving for instance the short-term dynamics of single synaptic responses (Tsodyks et al. 1998), the effects of a particular morphological structure (Rhodes 1999), or calcium dynamics (Larkum et al. 1999). So far the analysis of the in vivo phenomena based on single neuron response functions has been successfully performed with a simple linear relation between both m and 
and the spike frequency of the presynaptic neurons (Amit and Brunel 1997;
