Abstract
Wang, XiaoJing. Calcium coding and adaptive temporal computation in cortical pyramidal neurons. J. Neurophysiol. 79: 1549–1566, 1998. In this work, we present a quantitative theory of temporal spikefrequency adaptation in cortical pyramidal cells. Our model pyramidal neuron has twocompartments (a “soma” and a “dendrite”) with a voltagegated Ca^{2+} conductance (g _{Ca}) and a Ca^{2+}dependent K^{+} conductance (g _{AHP}) located at the dendrite or at both compartments. Its frequencycurrent relations are comparable with data from cortical pyramidal cells, and the properties of spikeevoked intracellular [Ca^{2+}] transients are matched with recent dendritic [Ca^{2+}] imaging measurements. Spikefrequency adaptation in response to a current pulse is characterized by an adaptation time constant τ_{adap} and percentage adaptation of spike frequency F _{adap} [% (peak − steady state)/peak]. We show how τ_{adap} and F _{adap} can be derived in terms of the biophysical parameters of the neural membrane and [Ca^{2+}] dynamics. Two simple, experimentally testable, relations between τ_{adap} and F _{adap} are predicted. The dependence of τ_{adap} and F _{adap} on current pulse intensity, electrotonic coupling between the two compartments, g _{AHP} as well the [Ca^{2+}] decay time constant τ_{Ca}, is assessed quantitatively. In addition, we demonstrate that the intracellular [Ca^{2+}] signal can encode the instantaneous neuronal firing rate and that the conductancebased model can be reduced to a simple calciummodel of neuronal activity that faithfully predicts the neuronal firing output even when the input varies relatively rapidly in time (tens to hundreds of milliseconds). Extensive simulations have been carried out for the model neuron with random excitatory synaptic inputs mimicked by a Poisson process. Our findings include 1) the instantaneous firing frequency (averaged over trials) shows strong adaptation similar to the case with current pulses; 2) when the g _{AHP} is blocked, the dendritic g _{Ca} could produce a hysteresis phenomenon where the neuron is driven to switch randomly between a quiescent state and a repetitive firing state. The firing pattern is very irregular with a large coefficient of variation of the interspike intervals (ISI CV > 1). The ISI distribution shows a long tail but is not bimodal. 3) By contrast, in an intrinsically bursting regime (with different parameter values), the model neuron displays a random temporal mixture of single action potentials and brief bursts of spikes. Its ISI distribution is often bimodal and its power spectrum has a peak. 4) The spikeadapting current I _{AHP}, as delayed inhibition through intracellular Ca^{2+} accumulation, generates a “forward masking” effect, where a masking input dramatically reduces or completely suppresses the neuronal response to a subsequent test input. When two inputs are presented repetitively in time, this mechanism greatly enhances the ratio of the responses to the stronger and weaker inputs, fulfilling a cellular form of lateral inhibition in time. 5) The [Ca^{2+}]dependent I _{AHP} provides a mechanism by which the neuron unceasingly adapts to the stochastic synaptic inputs, even in the stationary state following the input onset. This creates strong negative correlations between output ISIs in a frequencydependent manner, while the Poisson input is totally uncorrelated in time. Possible functional implications of these results are discussed.
INTRODUCTION
Cortical neurons display a large repertoire of voltage and calciumgated potassium ion channels with kinetic time constants ranging from milliseconds to seconds (Llinás 1988; Rudy 1988; Storm 1990). The diversity and richness of K^{+} conductances indicate that they likely contribute to neuronal inputoutput computation in ways more complex than sculpturing the waveform of action potentials or regulating the overall membrane excitability. For example, slow K^{+} currents, in interplay with Ca^{2+} and/or Na^{+} currents, can generate rhythmic firing patterns intrinsic to single neurons (Llinás 1988; Wang and Rinzel 1995). Or a slowly inactivating K^{+} current can integrate synaptic inputs in a temporalhistory–dependent manner (Storm 1988; Turrigiano et al. 1996; Wang 1993). Moreover, K^{+} channels at dendritic sites are capable of modifying cable properties and may regulate synaptic transmission (Hoffman et al. 1997) and prevent input saturation (Bernander et al. 1994; Wilson 1995).
Spikefrequency adaptation that depends on a Ca^{2+}gated K^{+} conductance is a conspicuous neuronal firing characteristic exhibited by a majority of (“regular spiking”) pyramidal neurons in neocortex and hippocampus (Avoli et al. 1994; Connors et al. 1982; Foehring et al. 1991; Gustafsson and Wigström 1981; Lanthorn et al. 1984; Lorenzon and Foehring 1992; Mason and Larkman 1990; McCormick et al. 1985). In response to a constant current pulse, the firing frequency of an adapting neuron is initially high then decreases to a lower steadystate plateau level within hundreds of milliseconds. This phenomenon has been studied intensively in in vitro slice experiments (as is the case for all aforecited references). Recently, Ahmed et al. (1993; B. Ahmed, C. Anderson, R. J. Douglas; K.A.C. Martin, unpublished results) observed and quantified spikefrequency adaptation of in vivo cortical neurons with intracellular recordings from the primary visual cortex of the anesthetized cat. They found that when subjected to a injected current pulse, the adaptation time course of cortical cells can be fitted empirically by an exponential time course (Ahmed et al. 1993; unpublished results), i.e., the instantaneous firing rate f(t) = f _{ss} + (f _{0} − f _{ss}) exp(−t/τ_{adap}), where f _{0} is the initial firing rate, f _{ss} is the steadystate firing rate, and τ_{adap} is an adaptation time constant. Thus this time course is characterized by two quantities: τ_{adap} and the percentage adaptation of firing frequency F _{adap} = (f _{0} − f _{ss})/f _{0}. Ahmed et al. (1993; unpublished results) found that τ_{adap} = 10–50 ms and F _{adap} = 50–70% with a significant difference between superficial and deep layer neurons. They also performed computer simulations that reproduced many of their observations.
In this modeling work, we present a quantitative study of spikefrequency adaptation temporal dynamics, which, in particular, yields analytical expressions for τ_{adap} and F _{adap} in terms of the cellular biophysical parameters. We also explore possible implications of this phenomenon in the realtime inputoutput computation of cortical neurons. Compared with an early quantitative work modeling spikefrequency adaptation in motoneurons by Baldsissera and Gustafsson (1974), the present study benefitted from a number of recent experimental findings and quantitative data about the cellular mechanisms underlying the spikefrequency adaptation phenomenon. First, it is well known that the spikefrequency adaptation is produced mainly by a voltageindependent, Ca^{2+}dependent K^{+} current, although other K^{+} currents (such as the M current) also are involved to a lesser degree (Madison and Nicoll 1984; Madison et al. 1987; McCormick and Williamson 1989). This current is associated with the slow after hyperpolarization (AHP) after a burst of spikes, hence is called the AHP current (I _{AHP}) (Hotson and Prince 1980; Lancaster and Adams 1986; Schwindt et al. 1988). Second, it has been demonstrated by photolytic manipulation of Ca^{2+} that the intrinsic gating of I _{AHP} is rapid; its slow activation is thus attributable to the kinetics of the cytoplasmic calcium concentration [Ca^{2+}] (Lancaster and Zucker 1994). Third, spikeevoked [Ca^{2+}] transients now can be measured by fluorescence imaging techniques (see Yuste and Tank 1996 for a review). Recently, to overcome the problem that a [Ca^{2+}] indicator dye like Fura2 is also a [Ca^{2+}] buffer, Helmchen et al. (1996) used increasingly low concentrations of Fura2 and, by extrapolation to zero dye concentration, obtained measurements of putatively intrinsic [Ca^{2+}] transient signals. Their estimated spikeevoked [Ca^{2+}] transient from dendrites of cortical pyramidal neurons are larger and faster than previously reported. In the model pyramidal neuron of the present paper, the calcium dynamics (spikeevoked influx and decay) is constrained by accurate measurements of Helmchen et al. (1996).
Can calcium signaling perform interesting sensory computation in cortical neurons? In previous computational models, spikefrequency adaptation has been incorporated as a gain control mechanism for neuronal excitability (Barkai and Hasselmo 1994; Douglas et al. 1995). However, the adaptation temporal dynamics, i.e., its role in momenttomoment neural computation in response to timevarying inputs, has not been emphasized. Through spikefrequency adaptation, [Ca^{2+}] dynamics produces a “forward masking” phenomenon: the neuronal response to a stimulus may be masked due to another stimulus that precedes it in time as was demonstrated experimentally in the cricket auditory neurons (Sobel and Tank 1994). Results reported here suggest that such an effect also exists in cortical pyramidal cells and may be used to selectively respond to temporal input patterns. In particular, when two or several competing inputs are presented, the neuronal output is sharply “tuned” to the strongest input, and the responses to weaker inputs are greatly suppressed. Furthermore, if the input consists of temporally uncorrelated excitatory postsynaptic potentials (EPSPs) (a Poisson process), spikefrequency adaptation leads to strong anticorrelation between the consecutive interspike intervals of the output spike train. These results indicate that the adaptation mechanism is operative even in the stationary state after the input onset, and suggest a direct means to assess its efficacy from extracellularly recorded spike trains.
METHODS
Model
The neuron model has two compartments, representing the dendrite and the soma plus axonal initial segment, respectively (Pinsky and Rinzel 1994). Many of our results can be obtained with a single compartment. However, we used a twocompartment model for three reasons. 1) Ca^{2+} imaging measurements from pyramidal cells show spikeevoked [Ca^{2+}] transient that is much larger at the dendrite than at the soma (Jaffe et al. 1994; Schiller et al. 1995; Svoboda et al. 1997; Yuste et al. 1994), and the [Ca^{2+}] influx is produced primarily by Ca^{2+} entry through voltagegated channels (Miyakawa et al. 1992). In the present model, we focus mainly on spikefrequency adaptation that is caused by a dendritic [Ca^{2+}]dependent I _{AHP}. 2) We wanted to see whether our theoretical analysis can be carried out even with two (or more) compartments. When I _{AHP} is present both at soma and dendrite, we show that there are two “calcium modes” and the spikefrequency adaptation time course should be described as a sum of two exponentials. And 3) with an appropriate choice of parameters, the neuron model displays burst firing patterns that require weak electrotonic interactions between the two compartments.
The dendritic compartment has a highthreshold calcium current (L type), I
_{Ca}, and a calciumdependent potassium current I
_{AHP}. The somatic compartment contains spike generating currents (I
_{Na} and I
_{K}) and possibly also I
_{Ca} and I
_{AHP}. The somatic and dendritic membrane potentials V
_{s} and V
_{d} obey the following currentbalance equations
The voltagedependent currents are described by the HodgkinHuxley formalism (Hodgkin and Huxley 1952). Thus a gating variable x satisfies a firstorder kinetics
The highthreshold calcium current I
_{Ca} = g
_{Ca}
m
_{∞}(V − V
_{Ca}), where m is replaced by its steadystate m
_{∞}(V) = 1/{1 + exp[−(V + 20)]/9} (Kay and Wong 1987). The voltageindependent, calciumactivated potassium current I
_{AHP} = g
_{AHP}[[Ca^{2+}]/([Ca^{2+}] + K
_{D})](V − V
_{K}), with K
_{D} = 30 μM. The intracellular calcium concentration [Ca^{2+}] is assumed to be governed by a leakyintegrator (Helmchen et al. 1996; Tank et al. 1995; Traub 1982)
The parameter values for the spikegenerating currents and the I _{AHP} were chosen so that the model displayed the initial as well as the adapted fI curves similar to the data of McCormick et al. (1985) and Mason and Larkman (1990): g _{L} = 0.1, g _{Na} = 45,g _{K} = 18, dendritic g _{Ca} = 1, and g _{AHP} = 5 (in mS/cm^{2}); V _{L} = −65, V _{Na} = +55, V _{K} = −80, V _{Ca} = +120 (in mV). The somatic g _{Ca} and g _{AHP} are zero except for Fig. 7. The resting state (with I = 0 and g _{syn} = 0) is at V _{s} = −64.8 and V _{d} = −64 mV.
The model was simulated on a Silicon Graphics Workstation, using a fourthorder RungeKutta method, with time step dt = 0.02–0.05 ms.
Statistical analysis
With Poisson synaptic inputs, each random output train of spike times is converted into a sequence of interspike intervals {t
_{1}, t
_{2}, t
_{3}, . . . , t
_{N}}. The interspike interval (ISI) histogram is tabulated, with a mean μ and a standard deviation σ given by
The ISI variability is quantified by the coefficient of variation (CV) (Softky and Koch 1993; Tuckwell 1988)
The temporal correlations are further assessed by the power spectrum of the spike train (with a time bin of 1 ms), using the subroutine Spctrm.c from Numerical Recipes (Press et al. 1989), modified by Yinghui Liu.
RESULTS
Time course of spikefrequency adaptation
In response to a depolarizing current pulse, the model neuron initially fires at a high frequency, then adapts to a lower steadystate frequency (Fig. 1 A). Spikefrequency adaptation is accompanied by a gradual increase of the fast spike AHP (from −53 to −57 mV, see Fig. 1 A, inset). This firing pattern is in parallel with the time course of Ca^{2+} accumulation, at a rate of ≃200 nM/spike [comparable with the [Ca^{2+}] imaging measurements from proximal apical dendrites of cortical layer V pyramidal cells (Helmchen et al. 1996)]. The I _{AHP} increases with [Ca^{2+}], hence the cell is gradually hyperpolarized and the firing frequency is decreased in time. In the steady state, an equilibrium is reached in the [Ca^{2+}] dynamics, when the spikeevoked [Ca^{2+}] influx rate is balanced with the [Ca^{2+}] decay rate. After the current pulse, there is a longlasting AHP that mirrors the Ca^{2+} (hence the I _{AHP}) decay (Fig. 1 A).
Frequencycurrent fI curves are shown in Fig. 1 B (left), for the initial first, third, and fifth interspike intervals, as well as the steady state. At the onset of repetitive firing (rheobase I ≃ 0.5), the firing frequency starts at zero, through a homoclinic bifurcation of the saddlenode type (see also Crook et al. 1997). It is noticeable that the initial fI curves are quite nonlinear, but the steadystate fI relation is very close to linear, similar to regular spiking pyramidal neurons (cf. Figs. 89 in Mason and Larkman 1990). Intuitively, the adapting AHPcurrent provides a delayed negative feedback to the cell. It is larger at higher firing frequencies, thus the difference between the initial f _{0} and the final f _{ss} increases with the current intensity. As a result, the steadystate inputoutput relation is linearized by the I _{AHP}. We also computed the mean dendritic [Ca^{2+}] as I was varied. Its steadystate plateau level depends linearly on f _{ss} (Fig. 1 B, right), with a slope ≃ 17 nM/Hz (compared with the measured 16 nM/Hz from dendrites of layer V pyramidal cells) (Helmchen et al. 1996). In this sense, the dendritic Ca^{2+} level encodes the neuronal firing activity (Helmchen et al. 1996; Johnston 1996).
For the spike train in Fig. 1, the instantaneous firing rate f(t) is defined as the reciprocal of ISIs (Fig. 2 A, top). Its time course can be well fitted by a monoexponential curve, f(t) = A + B exp(−t/τ_{adap}), where A = f _{ss}, and B = f _{0} − f _{ss}. The empirical best fit is given by f(t) = 116 + 156 exp(−t/33). Thus τ_{adap} = 33 ms, and the percentage adaptation F _{adap} = (f _{0} − f _{ss})/f _{0} = B/(A + B) = 57%. The [Ca^{2+}] also follows an exponential time course with the same time constant τ_{adap}, the steadystate plateau is [Ca^{2+}]_{ss} = 1.74 μM. Note that τ_{adap} is much shorter than the decay time constant τ_{Ca} = 80 ms. We now show how this adaptation time course, given a constant current pulse, can be predicted quantitatively from the biophysics of the membrane dynamics Eqs. 14 . We shall see further that this description leads to a calciumcoding model of neuronal output, even when the input varies temporally.
The fastslow variable dissection method (Baer et al. 1995; Ermentrout 1994; Guckenheimer et al. 1997; Rinzel 1987; Wang and Rinzel 1995) is based on the observation that Ca^{2+} evolves much more slowly than the membrane potential and the other channel gating variables of the model and that this slow [Ca^{2+}](t) determines the adaptation time course. Thus we can first analyze the fast subsystem and the slow [Ca^{2+}] dynamics separately then put them back together. This is done in three steps (see for details). In step 1, the fast electrical subsystem is analyzed while considering the slow [Ca^{2+}] as if it was a static parameter rather than a dynamical variable. This allows us to determine the functional [Ca^{2+}] dependence of the firing frequency f([Ca^{2+}]) as well as other voltagedependent quantities like 〈I
_{Ca}〉, where 〈x〉 denotes an average of x over a typical ISI. The functional [Ca^{2+}] dependence of f([Ca^{2+}]) and of 〈I
_{Ca}〉 was found to be approximately linear (Figs. 2, B and C)
In step 2, we turn to the [Ca^{2+}] dynamics Eq. 4
. Because it does not vary a lot within a sufficiently short ISI, I
_{Ca} is substituted by 〈I
_{Ca}〉, which is a function of [Ca^{2+}]. The resulting equation now only depends on [Ca^{2+}]
Finally, in step 3, we insert the time course for [Ca^{2+}] into the expression f = f
_{0} − G
_{f}[Ca^{2+}], which yields the adaptation process for the firing frequency as function of time
Note that the percentage adaptation
Calcium model of neuronal activity
We have described above how to reduce the biophysical membrane model of the neuron to a calciummodel, for a given applied current I
Equations 14 and 15 completely describe the calcium model of neuronal activity. For this model to be useful, it should be able to predict the output firing rate f(t), even when the input current I(t) varies temporally in an arbitrary fashion. In Helmchen et al. (1996), the dendritic [Ca^{2+}] signal was shown to encode firing frequency of a cortical layer V pyramidal neuron, when the input changed at the time scale of seconds. Intuitively, one expects that the calciummodel of neuronal discharges to be valid only when the temporal change of I(t) is slower than both the individual ISIs and the [Ca^{2+}] dynamics. Because of the adapting I _{AHP}, the effective [Ca^{2+}] time constant is given by τ_{adap}, typically much shorter than τ_{Ca}. Hence, one can expect that the [Ca^{2+}] code could remain effective even for input changes much more rapidly than in seconds. An example is shown in Fig. 4. Driven by a chaotic input current I(t) (Fig. 4, bottom), the membrane potential fires spikes irregularly (Fig. 4, top). The instantaneous firing rate (the reciprocal of the ISI) and [Ca^{2+}] as function of time are shown in Fig. 4, middle (in blue), superimposed with the predictions from the calciummodel (in red). This example suggests that the calcium model can indeed predict the instantaneous firing rate, even for relatively rapid time changes (within tens to hundreds of milliseconds) of the input current.
Dependence on the electrotonic coupling g_{c}
We next consider how the spikefrequency adaptation properties depend on the various biophysical parameters of the model. First, the spikefrequency adaptation is influenced strongly by the electrotonic coupling g _{c}, which controls the twoway current flow between the somatic and dendritic compartments. An example is illustrated in Fig. 5 A, with two different values of g _{c} and a same current pulse to the soma. With a larger g _{c}, there is greater current loss to the dendrite, hence the initial firing frequency is lower (Fig. 5 A, left). On the other hand, the dendritic membrane potential repolarizes more rapidly after the somatic spike, the dendritic spike width is narrower, and the [Ca^{2+}] influx per spike is reduced (Fig. 5 A, right). This leads to a slower [Ca^{2+}] accumulation, larger τ_{adap}, and smaller percentage adaptation F _{adap} (Fig. 5 B). As in the case of varying the applied current I (Fig. 3 B), the two quantities τ_{adap} and F _{adap} change in an antagonistic manner (faster adaptation time course is correlated with a greater degree of adaptation).
On the other hand, burst firing of spikes can be observed when the electrotonic coupling g _{c} between the two compartments is sufficiently small, and the dendritic membrane area is significantly larger than the somatic one so that the coupling is asymmetric and the somatodendrite influence is weak (Mainen and Sejnowski 1996; Pinsky and Rinzel 1994; Rhodes and Gray 1994; Traub 1982; Traub et al. 1994). Similar to the intrinsically bursting pyramidal cells in layer V neocortex (Kim and Connors 1993; Mason and Larkman 1990; McCormick et al. 1985; Nishimura et al. 1996), with moderate current intensity the model neuron fires doubles of spikes repetitively at low frequencies (∼4 Hz); and (more typically) current injection of higher intensities elicits an initial burst of spikes followed by a train of single spikes (Fig. 5 C). This firing pattern can be viewed as an extremely strong form of spikefrequency adaptation, produced by the same set of ion channels as the adaptation phenomenon (if these ion channels are located at dendritic sites sufficiently isolated from the soma).
Relations between τ_{adap} and F_{adap}
In addition to the neuronal electrotonic structure, spikefrequency adaptation depends also on the channel conductances g
_{Ca} and g
_{AHP}, as well as the [Ca^{2+}] kinetic parameters α and τ_{Ca}. These dependences were explored within the framework of our calciummodel. First, the initial firing rate f
_{0} should not depend on any of these parameters. Second, because 〈I
_{Ca}〉 = 〈I
_{Ca}〉_{0} + G
_{cc}[Ca^{2+}], 〈I
_{Ca}〉_{0} and G
_{cc} should be proportional to g
_{Ca}. Third, both gain parameters G
_{f} and G
_{cc} must be proportional to g
_{AHP}. In summary, we can write
By contrast, if τ_{Ca} is held constant and αg _{Ca} g _{AHP} is increased (Fig. 6 C), adaptation is faster (τ_{adap} is smaller) and the percentage adaptation F _{adap} is larger. The plot of F _{adap} versus τ_{adap} shows a linear relation with a negative slope (Fig. 6 D). In Fig. 6 D, we also plotted the F _{adap} versus τ_{adap} simulation data, obtained when the input current intensity I is varied (Fig. 3 B) and when the electrotonic coupling g _{c} is varied (Fig. 5 B). Quite surprisingly, in all cases, the F _{adap}τ_{adap} curve is linear with approximately the same slope ≃ 1/65, which is close to, but differs from, the reciprocal of the [Ca^{2+}] decay time constant 1/τ_{Ca} = 1/80.
To gain some insights about this general relation between F
_{adap} and τ_{adap}, let us suppose that the ratio between the initial value and the gain parameter is roughly the same for f and 〈I
_{Ca}〉, i.e., f
_{0}/G
_{f} ≃ (−〈I
_{Ca}〉_{0})/G
_{cc}. With this assumption, we then can write −〈I
_{Ca}〉_{0} G
_{f}/f
_{0} ≃ G
_{cc}. Inserting this relation into Eq. 13
, we have
Two calcium modes
We have developed our calcium model of neuronal activity Eqs. 14
and
15
with the ion currents I
_{Ca} and I
_{AHP} localized at the dendrite. In real cortical pyramidal neurons, of course, these channels are distributed widely on the dendritic trees as well as the soma (Jaffe et al. 1994; Johnston et al. 1996; Yuste et al. 1994). It is of interest to investigate whether our approach can be generalized to such cases where the [Ca^{2+}] signaling and [Ca^{2+}]dependent spikefrequency adaptation are distributed across multiple compartments. For our twocompartment model, suppose that the distributions of g
_{Ca} and g
_{AHP} are uniform at the soma and dendrite, g
_{Ca} = 1 and g
_{AHP} = 5 (in mS/cm^{2}) for both compartments. We then have two equations like Eq. 4
for somatic and dendritic [Ca^{2+}], respectively
The theoretical analysis for the adaptation time course can be generalized in this case, but only approximately. As is detailed in , our liner approach yields good estimates for the two adaptation time constants τ_{adap,1} and τ_{adap,2}. However, the steady state plateaus and the actual time courses cannot be predicted accurately unless nonlinear interactions between the two [Ca^{2+}] modes are taken into account.
Adaptation to stochastic Poisson input
So far, spikefrequency adaptation of cortical pyramidal neurons has been investigated with current pulse stimulation. However, in in vivo conditions, pyramidal cells are driven by synaptic inputs that are stochastic and vary unceasingly in time. As a first step toward addressing the question of spikefrequency adaptation to natural stimuli, we stimulated the response of our model neuron to random synaptic inputs generated by a Poisson process with rate λ. As illustrated in Fig. 8
A, when the Poisson input is turned on, the cell initially fires rapidly (at 180 Hz), then the firing rate is reduced in parallel with the accumulation of [Ca^{2+}]. The instantaneous firing frequency calculated by averaging over trials shows a time course similar to that with a constant current pulse. This time course can be predicted using the same analysis as before, except that now we use trialaveraged firing rate f and calcium current I
_{Ca}. As shown in Fig. 8
B, both f and 〈I
_{Ca}〉 are approximately linear with [Ca^{2+}] when the latter is considered as a parameter
Driven by random EPSPs, the neuronal firing activity displays fluctuations. We calculated the coefficient of variation of the ISIs (cf. methods). The CV displays a time course that mirrors the mean instantaneous firing rate (Fig. 8 A): its initial value is low but not zero (∼0.1), increases as the firing rate decreases, and plateaus at a value close to 0.5. The ISI variability is expected to be larger with lower firing frequencies when the cell acts more like a coincidence detector than a temporal integrator of excitatory inputs (Perkel et al. 1967; Softky and Koch 1993; Stern 1965).
Because the firing is now random, the membrane potential shows two noteworthy features. First the initial response appears as a burst, even though the neuron shows typical adaptation time course and no burst with a current pulse (Fig. 1). Second, in the steady state after the transient burst, fast doublets can be observed typically after a relatively long silent time. This phenomenon can be understood in terms of temporal correlations caused by the spikeadapting current I _{AHP}. Statistically, if by chance the cell happens to fire rapidly in a short time window, significant [Ca^{2+}] influx will be produced, the enhanced I _{AHP} will hyperpolarize the cell, and the firing rate subsequently will be reduced. Conversely, during a time period of relative low firing, I _{AHP} decreases as a result of the [Ca^{2+}] decay, and the probability of firing increases. In short, the adaptation generates negative statistical temporal correlations between ISIs. This is shown in Fig. 8 C by the ISI return map (computed in the steady state) where the (n + 1)th ISI is plotted against the nth ISI. Given ISI_{n}, the average ISI_{n+1} is calculated (blue curve) and yields a linear relation with a slope κ = −0.3. This slope is the same as the coefficient of correlation between successive ISIs (the CC, see Methods). It is about a third of the maximum negative correlation (κ = −1) possible.
The steadystate inputoutput relation, i.e., the mean firing rate f versus the Poisson input rate λ, is linear (Fig. 9 A, red curve). We also computed the ISI CV and CC as function of λ; they are plotted against the mean output frequency f (Fig. 9, B and C). One observes that the CV changes slightly with f, ranging 0.3–0.5 (Fig. 9 B, red curve). The CC shows a much stronger frequency dependence with a large negative peak at f ≃ 20 Hz (Fig. 9 C, red curve). This strong frequency dependence can be explained in terms of two time constants: τ_{Ca} = 80 ms and τ_{adap} ≃ 15 ms (for the sake of argument, we assume that τ_{adap} remains about the same as λ is varied). The negative ISI correlation is produced by the spikeadapting current I _{AHP}, and the effect is large only in a neuronal firing rate range such that the mean ISI (1/f) is between τ_{adap} and τ_{Ca}. This is because when the firing rate f is very low (f < 1/τ_{Ca} ∼ 10 Hz), [Ca^{2+}] decays back to baseline between spikes so that I _{AHP} cannot accumulate and κ ≃ 0. On the other hand, at very high firing rates (f > 1/τ_{adap} ∼ 70 Hz), no significant adaptation dynamics takes place during two consecutive ISIs, and their correlation should again be small.
This result demonstrates that, when the input is stochastic, the [Ca^{2+}]dependent current I _{AHP} strongly sculptures the firing patterns by producing negative correlations between ISIs in a frequencydependent manner, even in the stationary state where the trialaveraged firing rate is constant.
Stochastic burst firing
We confirmed that the ISI correlation (CC) is essentially zero when the spikefrequency adaptation is absent by blocking either g _{AHP} or g _{Ca} (Fig. 9 C, blue and green curves). Interestingly, with g _{AHP} = 0, the ISI variability (CV) is >1 at low frequencies, whereas it is always <1 if g _{Ca} = 0 (Fig. 9 B). Note that the behavior with g _{AHP} = 0 should be the same as the initial unadapted firing state when g _{AHP} is not blocked (but I _{AHP} = 0). Therefore, the fλ curves with or without g _{AHP} blockade in Fig. 9 A can be considered as the initial and the plateau inputoutput relations, respectively, under control conditions.
The firing is highly irregular (with CV > 1) when the g _{AHP} is blocked because of a hysteresis phenomenon that is produced by a dendritic calcium plateau potential (cf. Marder et al. 1996). Indeed, when the neuron model is excited by a constant input current (within an appropriate range), there is a bistability where both the resting state and a state of repetitive spiking are possible (data not shown). Because Poisson synaptic inputs are not constant in time, the model neuron is driven to switch randomly between the two states (Fig. 10 A). There is no significant ISI correlation in the absence of I _{AHP} (Fig. 10 B). The switching is not periodic but stochastic, and the ISI distribution is not bimodal (Fig. 10 C). The large CV is related to the presence of a long tail in the ISI distribution, with large ISIs corresponding to the time epochs when the cells is in the silent state (Fig. 10 C). Note that hysteresis between two states has been suggested previously to be a cause of large ISI variability, although in these cases the ISI distribution is bimodal (Stern et al. 1997; Wilbur and Rinzel 1983). Such dynamics with many long time scales also is reflected in the power spectrum of the spike train, which displays large amplitudes at low frequencies (Fig. 10 D, blue curve). By contrast, in an adapting state with g _{AHP} ≠ 0, powers at lower frequencies are strongly suppressed (Fig. 10 D, red curve). Such highpass filtering behavior is a common and important signature of temporally adapting dynamics.
Although random switching between two states gives rise to an apparently bursty pattern of neuronal firing (Fig. 10), we emphasize that it is very different from an intrinsically bursting behavior that does not require random switching by fluctuating inputs. We simulated the bursting state of Fig. 5 C with Poisson inputs (Fig. 11). In this case, the membrane potential displayed a temporal mixture of single spikes and bursts of spikes (Fig. 11 A). Unlike the bistable dynamics, bursting state showed a large negative correlation between ISIs due to the g _{AHP} (Fig. 11 B). Furthermore, the ISI distribution can be bimodal: a broad peak at typical ISI values between single spikes, whereas a second sharp peak atISI ≃ 3–5 ms corresponds to the short ISIs within a burst (Fig. 11 C). The bimodality is most evident at low mean firing rates when the two maxima are well separated (Fig. 11 C, blue and red curves). The bimodality also is reflected by the presence of a distinct peak in the power spectrum of spike trains (Fig. 11 D). At higher mean firing rates, the broad peak moves closer to, and eventually overlaps with, the sharp peak of the ISI distribution (Fig. 11 C, black curve), and the hump in the power spectrum disappears (Fig. 11 D).
To summarize, stochastic switching between a resting state and a firing state is characterized by a long tail in the ISI distribution and a high power at low frequencies. By contrast, an intrinsically bursting neuron is characterized by a random temporal mixture of bursts and single spikes, a strong negative ISI correlation, a bimodal ISI distribution and a peak in the power spectrum that depends on the level of the excitatory drive.
Forward masking
We investigated possible computational implications of the spikefrequency adaptation, in particular a “forward masking” effect suggested by the experiments of Sobel and Tank (1994) on the cricket Omega auditory neurons. Namely, when two or several inputs are presented sequentially in time, neuronal response to the first input would activate an I _{AHP} with a delay, thereby inhibiting responses to subsequent inputs. To see if such an effect also occurred in our pyramidal neuron model, we first presented to the cell a “masking” input pulse, then a test input pulse after a time interval of τ ms (Fig. 12). Indeed, the response to the test input was reduced dramatically by the masking input (Fig. 12, A and B). The peak response to the test input was decreased in proportion with the amplitude of the masking input, i.e., there was a linear relationship between the two with a negative slope (Fig. 12 C). The masking effect was produced by the residual [Ca^{2+}], which was accumulated by the neuronal firing during the masking input and which did not have enough time to return to its baseline between the two pulses as long as τ was not too large. Therefore, a significant amount of residual I _{AHP} inhibited the neuronal firing at the onset of the test pulse. Interestingly, there is often a time lag between the test input onset and neuronal response, which allows [Ca^{2+}] to decay further and leads to a gradually increasing time course of the neuronal response (Fig. 12 B). The requirement for the relative timing of the two inputs can be quantified by plotting the ratio of the two peak responses f _{2}/f _{1} versus τ. The recovery follows a decay time course like ∼exp(− τ/τ_{Ca}) (Fig. 12 D). This is expected because it should be the same as the decay time course of [Ca^{2+}] during the time τ, which, in the absence of spiketriggered influx, satisfies d[Ca^{2+}]/dt = −[Ca^{2+}]/τ_{Ca}.
When two or several inputs of different amplitudes are presented to the neuron model, one expects that the forward masking mechanism can produce a “selective attention” effect where responses to all inputs but the strongest one are suppressed (Pollack 1988). This is illustrated with two periodic trains of pulses with different amplitudes (Fig. 13 A). As is shown in Fig. 13 B, the inputoutput relation for the stimulus 2 is approximately linear; the presence of stimulus 1 shifts this response curve to higher input rates, but does not significantly change its slope (cf. Fig. 1C in Sobel and Tank 1994). Note that because the inputs were repetitive and [Ca^{2+}] accumulated between pulses, firing responses to both stimuli were affected by the adaptation dynamics: each was smaller than what it would have been in the absence of the other input. However, the effect was differential and the reduction was much smaller to the response to stronger input than that to the weaker one. This was shown by plotting the ratio of responses f _{2}/f _{1} versus that of inputs λ_{2}/λ_{1}, which yields a very nonlinear curve (Fig. 13 C). For example, if the stimulus 1 input rate is twice that of stimulus 2, the response to stimulus 1 was 25 times that of response 2 (Fig. 13 C). Therefore, in the presence of two or more inputs, the [Ca^{2+}]dependent adaptation process can selectively suppress the neuronal responses to weaker inputs so that the response to the strongest input “pops up” in time. This represents a cellular selective attention mechanism operating in the temporal domain.
DISCUSSION
The main findings of this work are twofold and are summarized in the following text.
Quantitative theory of spikefrequency adaptation and calcium coding
Cortical neurons and networks display many forms of adaptation to sensory inputs. Here we focused on the spikefrequency adaptation and developed a method to predict its time course (in response to a current pulse) characterized by two quantities: the adaptation time constant τ_{adap} and the percentage adaptation F _{adap}. This procedure led to a reduction of the twocompartment conductancebased model to a simple calcium model of neuronal firing activity. It was shown that this calcium model can predict the instantaneous neuronal output as the input changes rapidly in time in an arbitrary manner partly because the effective time constant τ_{adap} for the [Ca^{2+}] dynamics is much shorter than τ_{Ca} due to the adaptation feedback. This result confirms and strengthens the suggestion (Helmchen et al. 1996) that the intracellular Ca^{2+} is capable of encoding the temporal output computation of cortical pyramidal neurons. Of course, this description is limited to a rate code of neural information, the knowledge about the precise spike timing within a typical ISI being lost. Thus our calcium model is a firingrate model of neuronal activity that is derived from a conductancebased model and that incorporates the spikefrequency adaptation, a main firing characteristic of regular spiking cortical pyramidal neurons. Such a calciumbased model can be generalized to a network of interconnected pyramidal neurons.
Our analysis permitted us to derive τ_{adap} and F _{adap} in terms of the biophysical parameters of the membrane and intracellular Ca^{2+} dynamics. With plausible approximations, two relations between τ_{adap} and F _{adap} were obtained. The first relation is given by Eq. 18 : F _{adap} ≃ αG _{cc}τ_{adap}, where αG _{cc} is a “negative feedback gain parameter” and is proportional to g _{Ca} and g _{AHP} but independent of all [Ca^{2+}] extrusion and buffering processes. The second relation is given by Eq. 19 : F _{adap} ≃ 1 − τ_{adap}/τ_{Ca}, which predicts a negative proportionality between the two quantities as long as τ_{Ca} is maintained fixed. We emphasize that these relations are not exact but only approximate (see Fig. 6). They are not expected to hold if the linear theory is not valid and the adaptation time course is not exponential, for example, when adaptation is so strong that the firing is blocked completely in the steady state (Guckenheimer et al. 1997; Madison and Nicoll 1984). Moreover, it is important to assess how these predictions depend on details of the model. We performed some limited computer simulations to address this tissue, using a different functional form for I _{AHP} = g _{AHP} ([Ca^{2+}]/([Ca^{2+}] +K _{D})^{q} × (V − E _{K}), with q = 2 (Koch et al. 1995) or q = 3 (Lytton and Sejnowski 1993; Zhang et al. 1995) instead of q = 1. In this case, the reduced calcium equation was no longer linear. Numerically, we found that the second relation between τ_{adap} and F _{adap} still held roughly with the slope comparable with, but larger than, 1/τ_{Ca}. Another issue concerned other ionic currents such as the [Ca^{2+}]independent I _{M} (Madison and Nicoll 1984) or/and the nonspecific cation I _{H} (Lorenzon and Foehring 1992), which in some pyramidal neurons contribute to spikefrequency adaptation although to a much lesser extent than I _{AHP}. In principle, their effects could be analyzed in the same way as we did for I _{AHP}; the predicted relations between τ_{adap} and F _{adap} would need to be modified accordingly.
To the first order of approximation, however, we feel that the two relations are useful predictions that can be tested experimentally. The second relation can be assessed in pyramidal neurons, e.g., by measuring τ_{adap} and F _{adap} while a current pulse is applied with different intensities or when the g _{AHP} is blocked pharmacologically at various drug concentrations. These manipulations preserve τ_{Ca}, and our theory predicts that F _{adap} and τ_{adap} vary in such a way that they satisfy F _{adap} ≃ 1 − τ_{adap}/τ_{Ca}. Such measurements should provide an independent (if only crude) estimate for the [Ca^{2+}] decay time τ_{Ca}. Furthermore, because g _{AHP} is known to be modulated by transmitters such as acetylcholine, norepinephrine, serotonin, and histamine (McCormick and Williamson 1989; Nicoll 1988; Pedarizani and Storm 1993), one can ask if neuromodulation of the g _{AHP} could strongly regulate the spikefrequency adaptation dynamics while preserving the simple relation between τ_{adap} and F _{adap}.
In addition, these relations could be used to assess how the biophysical mechanism for spikefrequency adaptation differs from cell to cell in neocortical circuits. Ahmed et al. (1993, 1997) have quantitatively described spikefrequency adaptation in in vivo pyramidal neurons of the cat primary visual cortex. It was found that there was a conspicuous difference between superficial and deep layer neurons: superficial layer cells adapt faster and more strongly (τ_{adap} =11.5 ± 1.3 ms, and F _{adap} = 67 ± 3%) than deep layer cells(τ_{adap} = 51.4 ± 6.4 ms and F _{adap} = 51 ± 5%). According to our theoretical results, this finding cannot be explained primarily by a difference in the [Ca^{2+}] extrusion/buffering processes between superficial and deep layer neurons. If that was the case, the relation F _{adap} ≃ αG _{cc}τ_{adap} would predict that a smaller τ_{adap} should be correlated with a smaller F _{adap}, contrary to the observations.
Cellular selective attention and decorrelation
We investigated in detail the effect of spikefrequency adaptation to stochastic synaptic inputs of the Poisson type. In particular, we showed that the [Ca^{2+}]gated I _{AHP} can produce a forward masking effect similar to that shown experimentally in the Omega auditory neurons of the cricket (Sobel and Tank 1994). Essentially, the I _{AHP} provides a delayed negative feedback that is intrinsic to single pyramidal neurons and hence can be viewed as a cellular mechanism of lateral inhibition in time. Thus when two (or more) competing inputs converge onto a pyramidal neuron, the I _{AHP} can differentially suppress the responses to all but the strongest input (Fig. 13).
Because awake states are associated with activation of the brain stem cholinergic system and acetylcholine is a potent inhibitor of the I _{AHP} (and I _{M}) (McCormick and Williamsson 1989; Nicoll 1988), it is important to know how strong the spikefrequency adaptation effect is in pyramidal neurons during awake behavioral conditions. This question can be addressed, if the degree of neuronal adaptation can be assessed directly from extracellularly recorded spike trains. We found that the spikefrequency adaptation mechanism was manifested by a large negative coefficient of correlations of the ISIs (CC) (Figs. 8 and 9). Because the CC is readily computable from spike trains, it may subserve a probe for assessing the strength of adaptating ion currents (especially the I _{AHP}), under different in vivo conditions of the intact brain. Note that the CC is computed in the stationary state after transients. Therefore, in contrast to constant current pulses, with stochastic synaptic inputs the I _{AHP} is operative all the time and is responsible for generating negative temporal correlations between output ISIs in a frequencydependent manner, whereas the Poisson input is totally uncorrelated in time. This result suggests that, if the inputs actually are correlated strongly in time, such positive correlations can be reduced (inputs are decorrelated) by a subtraction mechanism through spikefrequency adaptation. Computational implications of this observation will be explored in a separate work.
Acknowledgments
I thank G. Turrigiano for comments on the manuscript.
This work was supported by National Institute of Mental Health Grant MH5371701), the Alfred P. Sloan Foundation, and the W. M. Keck Foundation.
Footnotes

Address reprint requests to X.J. Wang.
 Copyright © 1998 the American Physiological Society
APPENDIX A: ISI RETURN MAP AND COEFFICIENT OF CORRELATION
For a neuronal spike train converted into a sequence of ISIs {t
_{1}, t
_{2}, t
_{3}, . . ., t
_{N}} (with a mean μ), the coefficient of correlation between consecutive ISIs (CC) is defined in Eq. 7
of methods. Here we show that the CC is the same as the slope of the conditional average 〈t
_{n+1}‖t
_{n}〉 versus t
_{n}. Let p(t
_{n+1}, t
_{n}) be the joint probability for t
_{n+1} and t
_{n}. Then by definition
APPENDIX B: DERIVATION OF SPIKEFREQUENCY ADAPTATION TIME COURSE
The method is outlined in results. Here, we describe the procedure step by step. In step 1, we focus on the fast membrane dynamics Eqs. 13
, where [Ca^{2+}] is considered as a parameter rather than a variable. For a constant [Ca^{2+}], the I
_{AHP} = g
_{AHP}[[Ca^{2+}]/([Ca^{2+}] + K
_{D})](V − E
_{K}) is a “passive” outward current. Given I, the larger is [Ca^{2+}], so is I
_{AHP}, and the less the neuron should fire. We simulated the membrane Eqs. 13
with different [Ca^{2+}] values. For each [Ca^{2+}] value, the firing rate f was computed as well as the average 〈I
_{Ca}〉 over an ISI; they are plotted in Fig. 2, B and C, respectively. These curves are nonlinear, in fact f = 0 for [Ca^{2+}] > 2.2 μM. However, if we restrict ourselves to the range of [Ca^{2+}] between 0 and [Ca^{2+}]_{ss} = 1.74 μM of Fig. 2
A, the two curves in Fig. 2, B and C, can be approximated well by straight lines. This yields
In step 2, we turn our attention to the slowly evolving dynamics of [Ca^{2+}], Eq. 4. Because [Ca^{2+}] varies slowly, we substitute I
_{Ca} by its average 〈I
_{Ca}〉 which is a function of [Ca^{2+}] itself (Eq. EB1
). Therefore, we have
APPENDIX C: ADAPTATION TIME CONSTANTS BY TWO COUPLED CALCIUM MODES
With I
_{Ca} and I
_{AHP} in both somatic and dendritic compartments, we attempted to generalize our fastslow variable analysis to predict the biphasic spikefrequency adaptation time course. In step 1, we hold [Ca^{2+}]_{s} and [Ca^{2+}]_{d} as parameters. We make the simple assumption that the two adapting currents act independently and additively, so we numerically compute the dependence of the firing rate f and the calcium currents I
_{Ca,s} and I
_{Ca,d} on each of the [Ca^{2+}]_{s} and [Ca^{2+}]_{d} separately. For instance, if we consider the effect of the somatic I_{AHP} on spikefrequency adaptation, we let [Ca^{2+}]_{d} = 0 while [Ca^{2+}]_{s} is varied, we have (Fig. B1A)
The combined effect of both I
_{AHP,s} and I
_{AHP,d} is expressed as
The solution of Eq. EC4
has the form
The fact that c
_{d1} and c
_{d2} have different signs implies that the time course [Ca^{2+}]_{d}(t) is not monotonic and exhibits a maximum. Indeed, we can find the time t
_{max} when the maximum occurs by letting d[Ca^{2+}]_{d}/dt = 0, which yields
Finally, in step 3, by inserting the solutions for [Ca^{2+}]_{s}(t) and [Ca^{2+}]_{d}(t) into the expression of the firing rate (Eq. EC3
), we obtain
Compared with the empirical fits given by Eq. 21, we see that the linear theory predicts reasonably well the two adaptation time constants but not the steadystate values and the weighting factors for the two modes (the coefficients in Eqs. C5 and C9 ). We found that the linear analysis is not accurate, because there are significant interactions between the two modes which would introduce crossproduct terms like [Ca^{2+}]_{s} × [Ca^{2+}]_{d} in Eq. EC3 . An improved nonlinear analysis is beyond the scope of this paper.