Calcium Coding and Adaptive Temporal Computation in Cortical Pyramidal Neurons

Xiao-Jing Wang

Center for Complex Systems and Department of Physics, Brandeis University, Waltham, Massachusetts 02254

    ABSTRACT
Abstract
Introduction
Methods
Results
Discussion
References

Wang, Xiao-Jing. 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 spike-frequency adaptation in cortical pyramidal cells. Our model pyramidal neuron has two-compartments (a "soma" and a "dendrite") with a voltage-gated Ca2+ conductance (gCa) and a Ca2+-dependent K+ conductance (gAHP) located at the dendrite or at both compartments. Its frequency-current relations are comparable with data from cortical pyramidal cells, and the properties of spike-evoked intracellular [Ca2+] transients are matched with recent dendritic [Ca2+] imaging measurements. Spike-frequency adaptation in response to a current pulse is characterized by an adaptation time constant tau adap and percentage adaptation of spike frequency Fadap [% (peak - steady state)/peak]. We show how tau adap and Fadap can be derived in terms of the biophysical parameters of the neural membrane and [Ca2+] dynamics. Two simple, experimentally testable, relations between tau adap and Fadap are predicted. The dependence of tau adap and Fadap on current pulse intensity, electrotonic coupling between the two compartments, gAHP as well the [Ca2+] decay time constant tau Ca, is assessed quantitatively. In addition, we demonstrate that the intracellular [Ca2+] signal can encode the instantaneous neuronal firing rate and that the conductance-based model can be reduced to a simple calcium-model 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 gAHP is blocked, the dendritic gCa 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 spike-adapting current IAHP, as delayed inhibition through intracellular Ca2+ 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 [Ca2+]-dependent IAHP 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 frequency-dependent manner, while the Poisson input is totally uncorrelated in time. Possible functional implications of these results are discussed.

    INTRODUCTION
Abstract
Introduction
Methods
Results
Discussion
References

Cortical neurons display a large repertoire of voltage- and calcium-gated 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 input-output 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 Ca2+ 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 temporal-history-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).

Spike-frequency adaptation that depends on a Ca2+-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 steady-state plateau level within hundreds of milliseconds. This phenomenon has been studied intensively in in vitro slice experiments (as is the case for all afore-cited references). Recently, Ahmed et al. (1993; B. Ahmed, C. Anderson, R. J. Douglas; K.A.C. Martin, unpublished results) observed and quantified spike-frequency 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) = fss + (f0 - fss) exp(-t/tau adap), where f0 is the initial firing rate, fss is the steady-state firing rate, and tau adap is an adaptation time constant. Thus this time course is characterized by two quantities: tau adap and the percentage adaptation of firing frequency Fadap = (f0 - fss)/f0. Ahmed et al. (1993; unpublished results) found that tau adap = 10-50 ms and Fadap = 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 spike-frequency adaptation temporal dynamics, which, in particular, yields analytical expressions for tau adap and Fadap in terms of the cellular biophysical parameters. We also explore possible implications of this phenomenon in the real-time input-output computation of cortical neurons. Compared with an early quantitative work modeling spike-frequency 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 spike-frequency adaptation phenomenon. First, it is well known that the spike-frequency adaptation is produced mainly by a voltage-independent, Ca2+-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 (IAHP) (Hotson and Prince 1980; Lancaster and Adams 1986; Schwindt et al. 1988). Second, it has been demonstrated by photolytic manipulation of Ca2+ that the intrinsic gating of IAHP is rapid; its slow activation is thus attributable to the kinetics of the cytoplasmic calcium concentration [Ca2+] (Lancaster and Zucker 1994). Third, spike-evoked [Ca2+] transients now can be measured by fluorescence imaging techniques (see Yuste and Tank 1996 for a review). Recently, to overcome the problem that a [Ca2+] indicator dye like Fura-2 is also a [Ca2+] buffer, Helmchen et al. (1996) used increasingly low concentrations of Fura-2 and, by extrapolation to zero dye concentration, obtained measurements of putatively intrinsic [Ca2+] transient signals. Their estimated spike-evoked [Ca2+] 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 (spike-evoked 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, spike-frequency 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 moment-to-moment neural computation in response to time-varying inputs, has not been emphasized. Through spike-frequency adaptation, [Ca2+] 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), spike-frequency 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
Abstract
Introduction
Methods
Results
Discussion
References

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 two-compartment model for three reasons. 1) Ca2+ imaging measurements from pyramidal cells show spike-evoked [Ca2+] 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 [Ca2+] influx is produced primarily by Ca2+ entry through voltage-gated channels (Miyakawa et al. 1992). In the present model, we focus mainly on spike-frequency adaptation that is caused by a dendritic [Ca2+]-dependent IAHP. 2) We wanted to see whether our theoretical analysis can be carried out even with two (or more) compartments. When IAHP is present both at soma and dendrite, we show that there are two "calcium modes" and the spike-frequency 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 high-threshold calcium current (L type), ICa, and a calcium-dependent potassium current IAHP. The somatic compartment contains spike generating currents (INa and IK) and possibly also ICa and IAHP. The somatic and dendritic membrane potentials Vs and Vd obey the following current-balance equations
<IT>C</IT><SUB>m</SUB><FR><NU>d<IT>V</IT><SUB>s</SUB></NU><DE>d<IT>t</IT></DE></FR>= −I<SUB>L</SUB><IT>− I</IT><SUB>Na</SUB><IT>− I</IT><SUB>K</SUB><IT>− I</IT><SUB>Ca</SUB><IT>− I</IT><SUB>AHP</SUB><IT>− <FR><NU>g<SUB><UP>c</UP></SUB></NU><DE>p<UP></UP></DE></FR>(V</IT><SUB>s</SUB><IT>− V</IT><SUB>d</SUB>) + <IT>I</IT> (1)
<IT>C</IT><SUB>m</SUB><FR><NU>d<IT>V</IT><SUB>d</SUB></NU><DE>d<IT>t</IT></DE></FR>= −<IT>I</IT><SUB>L</SUB><IT>− I</IT><SUB>Ca</SUB><IT>− I</IT><SUB>AHP</SUB><IT>− <FR><NU>g<SUB><UP>c</UP></SUB><UP></UP></NU><DE>(1 − p<UP>)</UP></DE></FR>(V</IT><SUB><IT>d</IT></SUB><IT>− V</IT><SUB>s</SUB>) − <IT>I</IT><SUB>syn</SUB> (2)
where Cm = 1 µF/cm2 and IL = gL (V - VL) is the leak current. Following Pinsky and Rinzel (1994), we express the current flows between the soma and dendrite [proportional to (Vs - Vd)] in microamperes per square centimeter, with the coupling conductance gc = 2 mS/cm2, and the parameter p = somatic area/total area = 0.5. Other, voltage-gated currents are described below. The cell is either excited by an injected current I (in µA/cm2) to the soma or by a random synaptic input Isyn of the alpha -amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid type to the dendrite. Isyn =gsyn s(V - Esyn); the gating variable s obeys the equation ds/dt =eta (t- s/tau s, where eta (t) is a Poisson point-process with a rate lambda ,Esyn = 0 mV, tau s = 0.5 ms, and gsyn = 0.08 mS/cm2 (if present).

The voltage-dependent currents are described by the Hodgkin-Huxley formalism (Hodgkin and Huxley 1952). Thus a gating variable x satisfies a first-order kinetics
<FR><NU>d<IT>x</IT></NU><DE>d<IT>t</IT></DE></FR>= φ<SUB><IT>x</IT></SUB>[α<SUB>x</SUB>(<IT>V</IT>)(1 − <IT>x</IT>) − β<SUB>x</SUB>(<IT>V</IT>)<IT>x</IT>] = φ<SUB>x</SUB>[<IT>x</IT><SUB>∞</SUB>(<IT>V</IT>) − <IT>x</IT>]/τ<SUB><IT>x</IT></SUB>(<IT>V</IT>) (3)
The sodium current INa = gNam3infinity (V)h(V - VNa), where the fast activation variable is replaced by its steady state, minfinity alpha m/(alpha m + beta m), alpha m = -0.1(V + 33)/{exp[-0.1(V + 33)] - 1}, beta m =4 exp[-(V + 58)/12], alpha h = 0.07 exp[-(V + 50)/10], andbeta h = 1/{exp[-0.1(V + 20)] + 1}. The delayed rectifier IK =gKn4(V - VK), where alpha n = -0.01(V + 34)/{exp[-0.1(V +34)] - 1}, and beta n = 0.125 exp[-(V + 44)/25]. The temperature factor phi h = phi n = 4.

The high-threshold calcium current ICa = gCaminfinity (V - VCa), where m is replaced by its steady-state minfinity (V) = 1/{1 + exp[-(V + 20)]/9} (Kay and Wong 1987). The voltage-independent, calcium-activated potassium current IAHP = gAHP[[Ca2+]/([Ca2+] + KD)](V - VK), with KD = 30 µM. The intracellular calcium concentration [Ca2+] is assumed to be governed by a leaky-integrator (Helmchen et al. 1996; Tank et al. 1995; Traub 1982)
<FR><NU>d[Ca<SUP>2+</SUP>]</NU><DE>d<IT>t</IT></DE></FR>= −α<IT>I</IT><SUB>Ca</SUB>− [Ca<SUP>2+</SUP>]/τ<SUB>Ca</SUB> (4)
where alpha  is proportional to S/V with S being the membrane area and V the volume immediately beneath the membrane (Yamada et al. 1989). We used alpha  = 0.002 [in µM (msµA)-1cm2] so that the [Ca2+] influx per spike is ~200 nM (Helmchen et al. 1996). The various extrusion and buffering mechanisms are described collectively by a first-order decay process with a time constanttau Ca = 80 ms (Helmchen et al. 1996; Markram et al. 1995; Svoboda et al. 1997).

The parameter values for the spike-generating currents and the IAHP were chosen so that the model displayed the initial as well as the adapted f-I curves similar to the data of McCormick et al. (1985) and Mason and Larkman (1990): gL = 0.1, gNa = 45,gK = 18, dendritic gCa = 1, and gAHP = 5 (in mS/cm2); VL = -65, VNa = +55, VK = -80, VCa = +120 (in mV). The somatic gCa and gAHP are zero except for Fig. 7. The resting state (with I = 0 and gsyn = 0) is at Vs = -64.8 and Vd = -64 mV.


View larger version (14K):
[in this window]
[in a new window]
 
FIG. 7. Two Ca2+ modes. With gCa and gAHP on both the somatic and dendritic compartments, the time course of spike-frequency adaptation is a sum of 2 exponentials, with tau adap,1 = 29.4 ms and tau adap,2 = 191 ms. Faster component (largely due to the dendritic gAHP) dominates (top). Dendritic [Ca2+](t) is not monotonic and displays a maximum at around tmax = 106 ms. Red curves: empirical fits.

The model was simulated on a Silicon Graphics Workstation, using a fourth-order Runge-Kutta 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 {t1, t2, t3, . . . , tN}. The interspike interval (ISI) histogram is tabulated, with a mean µ and a standard deviation sigma  given by
μ = <FR><NU>1</NU><DE><IT>N</IT></DE></FR><LIM><OP>∑</OP><LL><SUB><IT>n</IT>=1</SUB></LL><UL><SUP><IT>N</IT></SUP></UL></LIM><IT>t<SUB>n</SUB></IT>;
σ = <RAD><RCD><FR><NU>1</NU><DE><IT>N</IT></DE></FR><LIM><OP>∑</OP><LL><SUB><IT>n</IT>=1</SUB></LL><UL><SUP><IT>N</IT></SUP></UL></LIM>(<IT>t<SUB>n</SUB></IT> − μ)<SUP>2</SUP></RCD></RAD>

The ISI variability is quantified by the coefficient of variation (CV) (Softky and Koch 1993; Tuckwell 1988)
CV = σ/μ (6)
The statistical interdependence between consecutive ISIs is measured by the coefficient of correlation (CC) (Perkel et al. 1967; Tuckwell 1988)
CC = <FR><NU>⟨(<IT>t<SUB>n<UP>+1</UP></SUB></IT>− μ)(<IT>t<SUB>n</SUB></IT>− μ)⟩</NU><DE>σ<SUP>2</SUP></DE></FR>= <FR><NU><FR><NU>1</NU><DE><IT>N</IT>− 1</DE></FR><LIM><OP>∑</OP><LL><SUB><IT>n</IT>=1</SUB></LL><UL><SUP><IT>N</IT>−1</SUP></UL></LIM>(<IT>t<SUB>n<UP>+1</UP></SUB></IT>− μ)(<IT>t<SUB>n</SUB></IT>− μ)</NU><DE>σ<SUP>2</SUP></DE></FR> (7)
which is between -1 and +1. The CC can be interpreted as follows. Suppose that we have an ISI return map, where tn+1 is plotted against tn. To assess how ISIs (tn+1) depend on their preceding values (tn), we calculate the conditional average of tn+1 for each given tn, < tn+1|tn> . If this function of tn is linear, then its slope is the same as the CC (see APPENDIX A for a derivation of this statement).

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
Abstract
Introduction
Methods
Results
Discussion
References

Time course of spike-frequency adaptation

In response to a depolarizing current pulse, the model neuron initially fires at a high frequency, then adapts to a lower steady-state frequency (Fig. 1A). Spike-frequency adaptation is accompanied by a gradual increase of the fast spike AHP (from -53 to -57 mV, see Fig. 1A, inset). This firing pattern is in parallel with the time course of Ca2+ accumulation, at a rate of sime 200 nM/spike [comparable with the [Ca2+] imaging measurements from proximal apical dendrites of cortical layer V pyramidal cells (Helmchen et al. 1996)]. The IAHP increases with [Ca2+], hence the cell is gradually hyperpolarized and the firing frequency is decreased in time. In the steady state, an equilibrium is reached in the [Ca2+] dynamics, when the spike-evoked [Ca2+] influx rate is balanced with the [Ca2+] decay rate. After the current pulse, there is a long-lasting AHP that mirrors the Ca2+ (hence the IAHP) decay (Fig. 1A).


View larger version (22K):
[in this window]
[in a new window]
 
FIG. 1. Spike-frequency adaptation characteristics. A: an example of spike-frequency adaptation in response to a current pulse. Adaptation is accompanied by a gradual increase of the fast spike afterhyperpolarization (AHP; top, inset). Each action potential generates a [Ca2+] influx of ~200 nM (bottom, inset), and the adaptation time course follows that of [Ca2+] (hence IAHP) accumulation. Slow AHP after the spike firing mirrors the [Ca2+] decay process. B: 1st, 3rd, and 5th instantaneous firing rates and the steady-state firing rate vs. the applied current intensity (left). Initial f-I curves are nonlinear, but the steady-state f-I relation is essentially linear. Plateau [Ca2+] level is a linear function of the steady-state firing rate, with a slope of ~13 nM/Hz (right).

Frequency-current f-I curves are shown in Fig. 1B (left), for the initial first, third, and fifth interspike intervals, as well as the steady state. At the onset of repetitive firing (rheobase sime  0.5), the firing frequency starts at zero, through a homoclinic bifurcation of the saddle-node type (see also Crook et al. 1997). It is noticeable that the initial f-I curves are quite nonlinear, but the steady-state f-I relation is very close to linear, similar to regular spiking pyramidal neurons (cf. Figs. 8-9 in Mason and Larkman 1990). Intuitively, the adapting AHP-current provides a delayed negative feedback to the cell. It is larger at higher firing frequencies, thus the difference between the initial f0 and the final fss increases with the current intensity. As a result, the steady-state input-output relation is linearized by the IAHP. We also computed the mean dendritic [Ca2+] as I was varied. Its steady-state plateau level depends linearly on fss (Fig. 1B, right), with a slope sime  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 Ca2+ level encodes the neuronal firing activity (Helmchen et al. 1996; Johnston 1996).


View larger version (29K):
[in this window]
[in a new window]
 
FIG. 8. Adaptation to a Poisson synaptic input. A: an example of the membrane potential and [Ca2+] time course (top and middle top). Cell initially fires rapidly that appears as a burst of spikes, followed by a firing pattern that is uneven in time. Instantaneous firing rate (averaged over 100 trials) has a single exponential time course (middle bottom). ISI variability increases with decreasing firing rate (bottom). Red curves: theoretical predictions. B: linear dependence of f and < ICa> on [Ca2+]. C: ISI return map (tn+1 vs. tn). Blue curve: conditional average of tn+1 for each fixed tn, which is approximately linear with a negative slope kappa  = -0.3. (Dendritic gAHP = 8 mS/cm2.)

For the spike train in Fig. 1, the instantaneous firing rate f(t) is defined as the reciprocal of ISIs (Fig. 2A, top). Its time course can be well fitted by a mono-exponential curve, f(t) A + B exp(-t/tau adap), where A = fss, and B = f0 - fss. The empirical best fit is given by f(t) = 116 + 156 exp(-t/33). Thus tau adap = 33 ms, and the percentage adaptation Fadap = (f0 - fss)/f0 = B/(A B) = 57%. The [Ca2+] also follows an exponential time course with the same time constant tau adap, the steady-state plateau is [Ca2+]ss = 1.74 µM. Note that tau adap is much shorter than the decay time constant tau 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. 1-4. We shall see further that this description leads to a calcium-coding model of neuronal output, even when the input varies temporally.


View larger version (17K):
[in this window]
[in a new window]
 
FIG. 2. Theoretical derivation of the spike-frequency adaptation time course. A: instantaneous firing rate [f = 1/interspike interval (ISI)] and [Ca2+] as function of time. Red curves: linear theory predictions. Green curve: empirical best fit. Blue curve: computed using a square-root fit for the f-vs.-[Ca2+] relation (see B). B: neuronal firing rate as function of [Ca2+]. bullet , data obtained when [Ca2+] was varied as a parameter. square , data from A with f plotted against [Ca2+] averaged over each individual ISI. Two data sets yield a same curve, demonstrating that the functional dependence of f on [Ca2+] can be obtained with [Ca2+] varied as a parameter. Red curve: linear regression fit; blue curve: square-root fit. C: average ICa as function of [Ca2+], fitted by a straight line (red). In B and C, the dashed vertical lines indicate the plateau [Ca2+] level in A.

The fast-slow 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 Ca2+ evolves much more slowly than the membrane potential and the other channel gating variables of the model and that this slow [Ca2+](t) determines the adaptation time course. Thus we can first analyze the fast subsystem and the slow [Ca2+] dynamics separately then put them back together. This is done in three steps (see APPENDIX B for details). In step 1, the fast electrical subsystem is analyzed while considering the slow [Ca2+] as if it was a static parameter rather than a dynamical variable. This allows us to determine the functional [Ca2+] dependence of the firing frequency f([Ca2+]) as well as other voltage-dependent quantities like < ICa> , where < x> denotes an average of x over a typical ISI. The functional [Ca2+] dependence of f([Ca2+]) and of < ICa> was found to be approximately linear (Figs. 2, B and C)
<IT>f ≃ f</IT><SUB>0</SUB><IT>− G</IT><SUB><IT>f</IT></SUB>[Ca<SUP>2+</SUP>];
⟨<IT>I</IT><SUB>CA</SUB>⟩ ≃ ⟨<IT>I</IT><SUB>Ca</SUB>⟩<SUB>0</SUB>+ <IT>G</IT><SUB>cc</SUB>[Ca<SUP>2+</SUP>] (8)
where f0 = 271 Hz is the initial firing frequency, and Gf = -df/d[Ca2+] = 84 (in Hz/µM) is a negative-feedback "gain parameter" for the firing frequency. < ICa> 0 = -28.8 µA/cm2 is the initial < ICa> , Gcc = d< ICa> /d[Ca2+] = 10 (in µA/cm2 µM) is a negative-feedback gain parameter for the [Ca2+]-dynamics (Ahmed et al., unpublished results), and alpha Gcc (in 1/ms) is the rate at which the [Ca2+] influx is reduced by [Ca2+] itself.

In step 2, we turn to the [Ca2+] dynamics Eq. 4. Because it does not vary a lot within a sufficiently short ISI, ICa is substituted by < ICa> , which is a function of [Ca2+]. The resulting equation now only depends on [Ca2+]
<FR><NU>d[Ca<SUP>2+</SUP>]</NU><DE>d<IT>t</IT></DE></FR>= (−α ⟨<IT>I</IT><SUB>Ca</SUB>⟩<SUB>0</SUB>) − [Ca<SUP>2+</SUP>]/τ<SUB>adap</SUB> (9)
with
<FR><NU>1</NU><DE>τ<SUB>adap</SUB></DE></FR>= α<IT>G<SUB>cc</SUB></IT>+ <FR><NU>1</NU><DE>τ<SUB>Ca</SUB></DE></FR> (10)
Solving Eq. 9, we obtain
[Ca<SUP>2+</SUP>](<IT>t</IT>) = [Ca<SUP>2+</SUP>]<SUB>ss</SUB>[1 − exp(−<IT>t</IT>/τ<SUB>adap</SUB>)] (11)
with [Ca2+]ss = -alpha < ICa> 0 tau adap. Note that the adaptation time constant tau adap (Eq. 10) is always smaller than tau Ca due to the presence of the negative feedback term alpha Gcc. For instance, with alpha  = 0.002 and Gcc = 10, alpha Gcc = 0.02 is larger than 1/tau Ca = 0.0125, hence contributes more to tau adap. We have tau adap = 30.8 ms, whereas tau Ca = 80 ms; and [Ca2+]ss = 1.77 µM.

Finally, in step 3, we insert the time course for [Ca2+] into the expression f = f0 - Gf[Ca2+], which yields the adaptation process for the firing frequency as function of time
<IT>f</IT>(<IT>t</IT>) = <IT>f</IT><SUB>ss</SUB>+ (<IT>f</IT><SUB>0</SUB><IT>− f</IT><SUB>ss</SUB>)exp(−<IT>t</IT>/τ<SUB>adap</SUB>) = 122 + 149 exp(−<IT>t</IT>/30.8) (12)
with fss = f0 - Gf[Ca2+]ss. The theoretically predicted time course for spike-frequency adaptation is shown in Fig. 2A (red curve), which compares well with the empirical fit (green curve). However, there is a small discrepancy between the numerical and predicted fss. This is mainly due to the fact that f([Ca2+]) is not exactly linear. Indeed, f goes to zero at a critical value of [Ca2+sime  2.2 µM (when the IAHP becomes too strong) via a homoclinic bifurcation of the saddle-node type. The mathematical theory of such a bifurcation predicts that, near the bifurcation, f([Ca2+]) behaves as a square-root function of [Ca2+] (see Guckenheimer et al. 1997; Rinzel and Ermentrout 1987). We found that a square-root function fits well even the global f([Ca2+]): f([Ca2+]) = 265 <RAD><RCD>1 − 0.455[Ca<SUP>2+</SUP>]</RCD></RAD>, except near the bifurcation threshold (Fig. 2B, blue curve). Using this nonlinear expression for f([Ca2+]) and the same [Ca2+](t) as before, f(t) can now be very accurately predicted (Fig. 2A, blue curve). In the following, we shall limit ourselves to the linear approximation.

Note that the percentage adaptation
<IT>F</IT><SUB>adap</SUB>= (<IT>f</IT><SUB>0</SUB><IT>− f</IT><SUB>ss</SUB>)/<IT>f</IT><SUB>0</SUB>= (−α ⟨<IT>I</IT><SUB>Ca</SUB>⟩<SUB>0</SUB>/<IT>f</IT><SUB>0</SUB>)<IT>G</IT><SUB><IT>f</IT></SUB>τ<SUB>adap</SUB> (13)
which depends on the [Ca2+] dynamics and the IAHP only through the factor Gftau adap.

Calcium model of neuronal activity

We have described above how to reduce the biophysical membrane model of the neuron to a calcium-model, for a given applied current I
<IT>f</IT>(<IT>t</IT>) = [<IT>f</IT><SUB>0</SUB>(<IT>I</IT>) − <IT>G</IT><SUB><IT>f</IT></SUB>(<IT>I</IT>)[Ca<SUP>2+</SUP>]]<SUB>+</SUB>
d[Ca<SUP>2+</SUP>]/d<IT>t</IT>= −α(⟨<IT>I</IT><SUB>Ca</SUB>⟩<SUB>0</SUB>(<IT>I</IT>) + <IT>G</IT><SUB>cc</SUB>(<IT>I</IT>)[Ca<SUP>2+</SUP>]) − [Ca<SUP>2+</SUP>]/τ<SUB>Ca</SUB> (14)
where the firing rate is always positive by using the half-rectifying function [x]+ = x if x >=  0, and 0 otherwise. The dependence on the current intensity I of the four quantities f0, Gf, < ICa> 0, and Gcc are shown in Fig. 3A. We observe that tau adap is larger and Fadap is smaller with larger current intensities (see also Ahmed et al. unpublished results). This is because at higher I, spike width becomes slightly narrower, hence Ca2+ influx per spike is reduced. All four curves in Fig. 3A can be fitted reasonably well by logarithmic functions
<IT>f</IT><SUB>0</SUB>(<IT>I</IT>) = 9[ln (<IT>I</IT>/0.3)]<SUB>+</SUB><IT>G</IT><SUB><IT>f</IT></SUB>(<IT>I</IT>) = 59 − 15[ln (<IT>I</IT>/0.3)]<SUB>+</SUB>
⟨<IT>I</IT><SUB>Ca</SUB>⟩<SUB>0</SUB>(<IT>I</IT>) = −84[ln (<IT>I</IT>/0.3)]<SUB>+</SUB><IT>G</IT><SUB>cc</SUB>(<IT>I</IT>) = 551 − 143[ln (<IT>I</IT>/0.3)]<SUB>+</SUB> (15)
where I0 = 0.3 is the estimated rheobase (The actual rheobase is somewhat larger, I0 sime  0.5). Note that such a curve-fitting is purely empirical and is not based on theoretical grounds (see also Agin 1964; Koch et al. 1995).


View larger version (24K):
[in this window]
[in a new window]
 
FIG. 3. Dependence on applied current intensity I. A: f0, Gf, < ICa> 0, and Gcc as functions of I. , empirical fitting functions of Eq. 15. B: adaptation time constant tau adap increases, and the percentage adaptation Fadap decreases, with I. C: examples with 4 current intensities (indicated on the right), the smooth curves are theoretical predictions.

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 [Ca2+] 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 calcium-model of neuronal discharges to be valid only when the temporal change of I(t) is slower than both the individual ISIs and the [Ca2+] dynamics. Because of the adapting IAHP, the effective [Ca2+] time constant is given by tau adap, typically much shorter than tau Ca. Hence, one can expect that the [Ca2+] 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 [Ca2+] as function of time are shown in Fig. 4, middle (in blue), superimposed with the predictions from the calcium-model (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.


View larger version (22K):
[in this window]
[in a new window]
 
FIG. 4. Calcium coding of neuronal electrical activity. In response to a temporally varying input I(t) (bottom), the cell's firing (blue dots, middle top) and [Ca2+] time course (blue curve, middle bottom) are well predicted by the reduced calcium model Eqs. 14 and 15 (red curves).

Dependence on the electrotonic coupling gc

We next consider how the spike-frequency adaptation properties depend on the various biophysical parameters of the model. First, the spike-frequency adaptation is influenced strongly by the electrotonic coupling gc, which controls the two-way current flow between the somatic and dendritic compartments. An example is illustrated in Fig. 5A, with two different values of gc and a same current pulse to the soma. With a larger gc, there is greater current loss to the dendrite, hence the initial firing frequency is lower (Fig. 5A, left). On the other hand, the dendritic membrane potential repolarizes more rapidly after the somatic spike, the dendritic spike width is narrower, and the [Ca2+] influx per spike is reduced (Fig. 5A, right). This leads to a slower [Ca2+] accumulation, larger tau adap, and smaller percentage adaptation Fadap (Fig. 5B). As in the case of varying the applied current I (Fig. 3B), the two quantities tau adap and Fadap change in an antagonistic manner (faster adaptation time course is correlated with a greater degree of adaptation).


View larger version (17K):
[in this window]
[in a new window]
 
FIG. 5. Dependence on the electrotonic coupling gc between the 2 compartments. A: with an increased gc, the initial firing frequency is lower but the steady-state frequency remains the same (left). Reduced percentage adaptation is due to a narrower dendritic spike (because Vd repolarizes more rapidly at larger gc), hence a smaller [Ca2+] influx per spike (right). B: tau adap increases, and Fadap decreases, with gc. C: burst firing patterns with modified electrotonic properties (gc = 1.4, P = 0.3) and gCa = 0.5,gAHP = 18.

On the other hand, burst firing of spikes can be observed when the electrotonic coupling gc 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 soma-to-dendrite 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. 5C). This firing pattern can be viewed as an extremely strong form of spike-frequency 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 tau adap and Fadap

In addition to the neuronal electrotonic structure, spike-frequency adaptation depends also on the channel conductances gCa and gAHP, as well as the [Ca2+] kinetic parameters alpha  and tau Ca. These dependences were explored within the framework of our calcium-model. First, the initial firing rate f0 should not depend on any of these parameters. Second, because < ICa>  = < ICa> 0 + Gcc[Ca2+], < ICa> 0 and Gcc should be proportional to gCa. Third, both gain parameters Gf and Gcc must be proportional to gAHP. In summary, we can write
<IT>G<SUB>f</SUB></IT>= <IT><A><AC>G</AC><AC>¯</AC></A><SUB>f</SUB>g</IT><SUB>AHP</SUB>,
⟨<IT>I</IT><SUB>Ca</SUB>⟩<SUB>0</SUB>= ⟨<IT><A><AC>I</AC><AC>¯</AC></A></IT><SUB>Ca</SUB>⟩<SUB>0</SUB><IT>g</IT><SUB>Ca</SUB>,
<IT>G</IT><SUB>cc</SUB><IT>= <A><AC>G</AC><AC>¯</AC></A></IT><SUB>cc</SUB><IT>g</IT><SUB>Ca</SUB><IT>g</IT><SUB>AHP</SUB> (16)
Inserting these scaling relations into Eqs. 10 and 13 we have
1/τ<SUB>adap</SUB> = (α<IT>g</IT><SUB>Ca</SUB><IT>g</IT><SUB>AHP</SUB>)<IT><A><AC>G</AC><AC>¯</AC></A></IT><SUB>cc</SUB> + 1/τ<SUB>Ca</SUB>
<IT>F</IT><SUB>adap</SUB>= (α<IT>g</IT><SUB>Ca</SUB><IT>g</IT><SUB>AHP</SUB>)(−<IT><A><AC>G</AC><AC>¯</AC></A></IT><SUB><IT>f</IT></SUB>⟨<IT><A><AC>I</AC><AC>¯</AC></A></IT><SUB>Ca</SUB>⟩<SUB>0</SUB>/<IT>f</IT><SUB>0</SUB>)τ<SUB>adap</SUB> (17)
From Eq. 17 we can conclude that tau adap and Fadap depend only on tau Ca (a neuron's [Ca2+] extrusion and buffering properties) and on the combination alpha gCagAHP (the product of the spike-evoked [Ca2+]-influx size alpha gCa and the adaptation conductance gAHP). Given fixed alpha gCagAHP, as tau Ca is increased (Fig. 6A), the adaptation is slower (tau adap is larger), but the plateau [Ca2+] level is higher (cf. Eqs. 4 and 11), hence Fadap is larger. A plot of Fadap versus tau adap is linear with a positive slope (Fig. 6B). This is evident in Eq. 17, where f0, &Gmacr;f, and < ICa> 0 (computed with [Ca2+] as a static parameter) are all independent of tau Ca, hence Fadap is simply proportional to tau adap. The slope of the linear curve, however, depends on the input current intensity I.


View larger version (18K):
[in this window]
[in a new window]
 
FIG. 6. Dependence of tau adap and Fadap on tau Ca (A and B) and the combination alpha gCagAHP (C and D) (cf. Eq. 17). Three curves in each panel correspond to different applied current intensities. B: Fadap vs. tau adap from data in A is linear with a positive slope (which depends on I). D: Fadap vs. tau adap from data in C: is linear with a negative slope sime  -1/tau Ca. This is also true for data obtained with I (Fig. 3B) or gc (Fig. 5B) being varied as parameter.

By contrast, if tau Ca is held constant and alpha gCagAHP is increased (Fig. 6C), adaptation is faster (tau adap is smaller) and the percentage adaptation Fadap is larger. The plot of Fadap versus tau adap shows a linear relation with a negative slope (Fig. 6D). In Fig. 6D, we also plotted the Fadap versus tau adap simulation data, obtained when the input current intensity I is varied (Fig. 3B) and when the electrotonic coupling gc is varied (Fig. 5B). Quite surprisingly, in all cases, the Fadap-tau adap curve is linear with approximately the same slope sime  1/65, which is close to, but differs from, the reciprocal of the [Ca2+] decay time constant 1/tau Ca = 1/80.

To gain some insights about this general relation between Fadap and tau adap, let us suppose that the ratio between the initial value and the gain parameter is roughly the same for f and < ICa> , i.e., f0/Gf sime  (-< ICa> 0)/Gcc. With this assumption, we then can write -< ICa> 0 Gf/f0 sime  Gcc. Inserting this relation into Eq. 13, we have
<IT>F</IT><SUB>adap</SUB><IT>≃ αG</IT><SUB>cc</SUB>τ<SUB>adap</SUB> (18)
or because alpha Gcc = 1/tau adap - 1/tau Ca (Eq. 10), we have
<IT>F</IT><SUB>adap</SUB>≃ 1 − τ<SUB>adap</SUB>/τ<SUB>Ca</SUB> (19)
which is the desired relation between Fadap and tau adap. This theoretical prediction is directly testable by experimental measurements from cortical pyramidal neurons.

Two calcium modes

We have developed our calcium model of neuronal activity Eqs. 14 and 15 with the ion currents ICa and IAHP 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 [Ca2+] signaling and [Ca2+]-dependent spike-frequency adaptation are distributed across multiple compartments. For our two-compartment model, suppose that the distributions of gCa and gAHP are uniform at the soma and dendrite, gCa = 1 and gAHP = 5 (in mS/cm2) for both compartments. We then have two equations like Eq. 4 for somatic and dendritic [Ca2+], respectively
<FR><NU>d[Ca<SUP>2+</SUP>]<SUB>s</SUB></NU><DE>d<IT>t</IT></DE></FR> = −α<SUB>s</SUB><IT>I</IT><SUB>Ca,s</SUB> − [Ca<SUP>2+</SUP>]<SUB>s</SUB>/τ<SUB>Ca,s</SUB>
<FR><NU>d[Ca<SUP>2+</SUP>]<SUB>d</SUB></NU><DE>d<IT>t</IT></DE></FR>= −α<SUB>d</SUB><IT>I</IT><SUB>Ca,d</SUB>− [Ca<SUP>2+</SUP>]<SUB>d</SUB>/τ<SUB>Ca,d</SUB> (20)
Because the parameter alpha  is proportional to the surface:volume ratio, it should be much smaller for the soma than for the dendrite. Similarly, tau Ca is expected to be longer at the soma than at the dendrite. For instance, the spike-evoked [Ca2+] transient peak amplitude and decay rate (1/tau Ca) can be three times as large in some dendritic sites as in the soma (Schiller et al. 1995). With the simple two-compartment model, let us assume that for the dendritic compartment, alpha d = 0.002 and tau Ca,d = 80 ms, whereas for the somatic compartment we multiple the right-hand side of Eq. 4 by a factor of 1/3, so that alpha s = 0.000667 and tau Ca,s = 240 ms. As shown in Fig. 7, with two calcium modes, [Ca2+]s(t), [Ca2+]d(t), and f(t) can be empirically well fitted by a sum of two exponentials, with the first rapid phase followed by a second much slower phase (Fig. 7, red curves)
[Ca<SUP>2+</SUP>]<SUB>s</SUB>(<IT>t</IT>) = 0.74 − 0.3 exp(−<IT>t</IT>/τ<SUB>adap,1</SUB>) − 0.44 exp(−<IT>t</IT>/τ<SUB>adap,2</SUB>)
[Ca<SUP>2+</SUP>]<SUB>d</SUB>(<IT>t</IT>) = 1.13 − 1.63 exp(−<IT>t</IT>/τ<SUB>adap,1</SUB>) + 0.5 exp(−<IT>t</IT>/τ<SUB>adap,2</SUB>)
<IT>f</IT>(<IT>t</IT>) = 73 + 225 exp(−<IT>t</IT>/τ<SUB>adap,1</SUB>) + 16 exp(−<IT>t</IT>/τ<SUB>adap,2</SUB>) (21)
where tau adap,1 = 29.4 ms and tau adap,2 = 191 ms. None that the faster [Ca2+] mode (which is primarily of the dendritic origin) dominates the spike-frequency adaptation, while the slower mode is only 16/225 = 7% of the faster mode. Moreover, the time course of [Ca2+]d can be nonmonotonic (Fig. 7). In the first rapid phase of adaptation, [Ca2+]s is small, both [Ca2+]d and the firing rate converge to their plateau levels as if the somatic IAHP did not exist. But as [Ca2+]s (hence the somatic IAHP) slowly accumulates, the firing rate further decreases in the second phase of adaptation, which leads to a decay of [Ca2+]d to its actual final plateau, which is smaller than would be expected without a somatic IAHP. Hence, the nonmonotonic behavior of [Ca2+]d(t).

The theoretical analysis for the adaptation time course can be generalized in this case, but only approximately. As is detailed in APPENDIX C, our liner approach yields good estimates for the two adaptation time constants tau adap,1 and tau adap,2. However, the steady state plateaus and the actual time courses cannot be predicted accurately unless nonlinear interactions between the two [Ca2+] modes are taken into account.

Adaptation to stochastic Poisson input

So far, spike-frequency 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 spike-frequency adaptation to natural stimuli, we stimulated the response of our model neuron to random synaptic inputs generated by a Poisson process with rate lambda . As illustrated in Fig. 8A, 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 [Ca2+]. 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 trial-averaged firing rate f and calcium current ICa. As shown in Fig. 8B, both f and < ICa> are approximately linear with [Ca2+] when the latter is considered as a parameter
<IT>f = f</IT><SUB>0</SUB><IT>− G</IT><SUB><IT>f</IT></SUB>[Ca<SUP>2+</SUP>];
⟨<IT>I</IT><SUB>Ca</SUB>⟩ = ⟨<IT>I</IT><SUB>Ca</SUB>⟩<SUB>0</SUB>+ <IT>G</IT><SUB>cc</SUB>[Ca<SUP>2+</SUP>] (22)
where f0 = 213 Hz, Gf = 252 Hz/µM, < ICa> 0 = -22.6 µA/cm2, and Gcc = 27.5 µA/cm2 µM. Solving Eq. 4 with ICa substituted by < ICa> , we obtain
[Ca<SUP>2+</SUP>](<IT>t</IT>) = [Ca<SUP>2+</SUP>]<SUB>ss</SUB>(1 − exp(−<IT>t</IT>/τ<SUB>adap</SUB>))
<IT>f</IT>(<IT>t</IT>) = <IT>f</IT><SUB>ss</SUB>+ (<IT>f</IT><SUB>0</SUB><IT>− f</IT><SUB>ss</SUB>) exp(−<IT>t</IT>/τ<SUB>adap</SUB>) (23)
where tau adap = 14.8 ms, [Ca2+]ss = 0.67 µM, fss = 44.2 and f0 = 213 (in Hz). (Note, again, tau adap is much shorter than tau Ca = 80 ms.) These curves fit accurately with the simulation data (red curves in Fig. 8A).

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. 8A): 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 spike-adapting current IAHP. Statistically, if by chance the cell happens to fire rapidly in a short time window, significant [Ca2+] influx will be produced, the enhanced IAHP will hyperpolarize the cell, and the firing rate subsequently will be reduced. Conversely, during a time period of relative low firing, IAHP decreases as a result of the [Ca2+] decay, and the probability of firing increases. In short, the adaptation generates negative statistical temporal correlations between ISIs. This is shown in Fig. 8C by the ISI return map (computed in the steady state) where the (n + 1)th ISI is plotted against the nth ISI. Given ISIn, the average ISIn+1 is calculated (blue curve) and yields a linear relation with a slope kappa  = -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 (kappa  = -1) possible.

The steady-state input-output relation, i.e., the mean firing rate f versus the Poisson input rate lambda , is linear (Fig. 9A, red curve). We also computed the ISI CV and CC as function of lambda ; 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. 9B, red curve). The CC shows a much stronger frequency dependence with a large negative peak at f sime  20 Hz (Fig. 9C, red curve). This strong frequency dependence can be explained in terms of two time constants: tau Ca = 80 ms and tau adap sime  15 ms (for the sake of argument, we assume that tau adap remains about the same as lambda  is varied). The negative ISI correlation is produced by the spike-adapting current IAHP, and the effect is large only in a neuronal firing rate range such that the mean ISI (1/f) is between tau adap and tau Ca. This is because when the firing rate f is very low (f < 1/tau Ca ~ 10 Hz), [Ca2+] decays back to baseline between spikes so that IAHP cannot accumulate and kappa  sime  0. On the other hand, at very high firing rates (f > 1/tau adap ~ 70 Hz), no significant adaptation dynamics takes place during two consecutive ISIs, and their correlation should again be small.


View larger version (13K):
[in this window]
[in a new window]
 
FIG. 9. Dependence on the Poisson input rate lambda . A: average steady-state firing rate f as function of lambda . ISI variability (CV) (B) and ISI correlation (CC) (C) are plotted as function of f. Red: control (with dendritic gCa = 1 and gAHP = 8 mS/cm2), blue: gAHP = 0, green: gCa = 0. Note that when gAHP = 0 but gCa is present, the CV is larger than 1 at low firing frequencies (B). On the other hand, the CC is virtually 0, hence the ISIs are essentially uncorrelated in time when gAHP = 0 (C). In control condition, CC depends strongly on firing frequency and displays a negative peak at sime  20 Hz.

This result demonstrates that, when the input is stochastic, the [Ca2+]-dependent current IAHP strongly sculptures the firing patterns by producing negative correlations between ISIs in a frequency-dependent manner, even in the stationary state where the trial-averaged firing rate is constant.

Stochastic burst firing

We confirmed that the ISI correlation (CC) is essentially zero when the spike-frequency adaptation is absent by blocking either gAHP or gCa (Fig. 9C, blue and green curves). Interestingly, with gAHP = 0, the ISI variability (CV) is >1 at low frequencies, whereas it is always <1 if gCa = 0 (Fig. 9B). Note that the behavior with gAHP = 0 should be the same as the initial unadapted firing state when gAHP is not blocked (but IAHP = 0). Therefore, the f-lambda curves with or without gAHP blockade in Fig. 9A can be considered as the initial and the plateau input-output relations, respectively, under control conditions.

The firing is highly irregular (with CV > 1) when the gAHP 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. 10A). There is no significant ISI correlation in the absence of IAHP (Fig. 10B). The switching is not periodic but stochastic, and the ISI distribution is not bimodal (Fig. 10C). 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. 10C). 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. 10D, blue curve). By contrast, in an adapting state with gAHP not equal  0, powers at lower frequencies are strongly suppressed (Fig. 10D, red curve). Such high-pass filtering behavior is a common and important signature of temporally adapting dynamics.


View larger version (48K):
[in this window]
[in a new window]
 
FIG. 10. A: an example with gAHP = 0 and lambda  = 0.3 kHz, which shows high ISI variability (CV > 1). Membrane potential randomly switches in time between a resting (inactive) phase and a firing (active) phase. B: consecutive ISIs are not correlated in time. C: ISI distribution shows a long tail that spans several mean ISI (blue curve), in contrast with the case of Fig. 8 with adaptation (red curve). ISI distributions are plotted against ISI/< ISI> . D: power spectrum of the spike train has a large amplitude at low frequencies (blue curve), whereas the adaptation current IAHP filters out response power at low-frequencies (red curve).

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. 5C with Poisson inputs (Fig. 11). In this case, the membrane potential displayed a temporal mixture of single spikes and bursts of spikes (Fig. 11A). Unlike the bistable dynamics, bursting state showed a large negative correlation between ISIs due to the gAHP (Fig. 11B). Furthermore, the ISI distribution can be bimodal: a broad peak at typical ISI values between single spikes, whereas a second sharp peak atISI sime  3-5 ms corresponds to the short ISIs within a burst (Fig. 11C). The bimodality is most evident at low mean firing rates when the two maxima are well separated (Fig. 11C, blue and red curves). The bimodality also is reflected by the presence of a distinct peak in the power spectrum of spike trains (Fig. 11D). At higher mean firing rates, the broad peak moves closer to, and eventually overlaps with, the sharp peak of the ISI distribution (Fig. 11C, black curve), and the hump in the power spectrum disappears (Fig. 11D).


View larger version (39K):
[in this window]
[in a new window]
 
FIG. 11. Burst firing pattern by Poisson input. A: in the bursting regime (with the same parameters as in Fig. 5C), the cell's response to random synaptic inputs typically show a temporal mixture of single spikes and bursts of spikes. From top to bottom: lambda  = 1.0, 1.5, and 2.5 kHz. B: ISI return map for top panel in A, showing a strong negative correlation between the consecutive ISIs (blue curve: conditional average of tn+1 for given tn). C: ISI distributions for the 3 examples in A. There is a bimodality at low rates (blue and red). An increased firing frequency shifts the second peak to the left, so that the 2 peaks overlap with each other (black). D: bimodality of the ISI distribution indicates a statistical repetition of the firing events in time, which is reflected in the power spectra of the spike train by the presence of a significant peak (blue and red). Peak disappears when the ISI distribution is no longer bimodal (black).

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 spike-frequency 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 IAHP 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 tau  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. 12C). The masking effect was produced by the residual [Ca2+], 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 tau  was not too large. Therefore, a significant amount of residual IAHP 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 [Ca2+] to decay further and leads to a gradually increasing time course of the neuronal response (Fig. 12B). The requirement for the relative timing of the two inputs can be quantified by plotting the ratio of the two peak responses f2/f1 versus tau . The recovery follows a decay time course like ~exp(- tau /tau Ca) (Fig. 12D). This is expected because it should be the same as the decay time course of [Ca2+] during the time tau , which, in the absence of spike-triggered influx, satisfies d[Ca2+]/dt = -[Ca2+]/tau Ca.


View larger version (22K):
[in this window]
[in a new window]
 
FIG. 12. Forward masking effect. A: 2 input pulses are presented to the cell model with a time separation of tau  ms. With increasing amplitude of the first (masking) input, the response to the second (test) input is dramatically reduced. B: an example with the spike rastergram (top), the average instantaneous firing frequency (middle top). [Ca2+] does not have time to decay back to its baseline during the time interval between the masking and test pulses (middle bottom), this residual [Ca2+] is at the origin of the forward masking phenomenon. C: peak response to the test input is negatively proportional to the amplitude of the masking input. D: masking effect decreases with tau  according to a recovery time course about exp (-tau /tau Ca), tau Ca = 200 ms. (Dendritic gAHP = 6 mS/cm2.)

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. 13A). As is shown in Fig. 13B, the input-output 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 [Ca2+] 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 f2/f1 versus that of inputs lambda 2/lambda 1, which yields a very nonlinear curve (Fig. 13C). 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. 13C). Therefore, in the presence of two or more inputs, the [Ca2+]-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.


View larger version (22K):
[in this window]
[in a new window]
 
FIG. 13. Cellular selective attention. A: 2 input pulse trains (blue and red) are presented to the cell model with slightly different amplitudes (lambda 1 and lambda 2, respectively). Response to the weaker input f2 is suppressed differentially as compared with f1 to the stronger input. B: f2 as function of lambda 2, for 3 different lambda 1 values. Presence of the masking input lambda 1 shifts the f-I curve to the right but does not change the input-output slope. C: response ratio f1/f2 is a highly nonlinear function of the input ratio lambda 1/lambda 2, demonstrating a selective inhibition of the response to the weaker input. Dendritic gAHP = 6 mS/cm2, tau Ca = 200 ms.

    DISCUSSION
Abstract
Introduction
Methods
Results
Discussion
References

The main findings of this work are twofold and are summarized in the following text.

Quantitative theory of spike-frequency adaptation and calcium coding

Cortical neurons and networks display many forms of adaptation to sensory inputs. Here we focused on the spike-frequency adaptation and developed a method to predict its time course (in response to a current pulse) characterized by two quantities: the adaptation time constant tau adap and the percentage adaptation Fadap. This procedure led to a reduction of the two-compartment conductance-based 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 tau adap for the [Ca2+] dynamics is much shorter than tau Ca due to the adaptation feedback. This result confirms and strengthens the suggestion (Helmchen et al. 1996) that the intracellular Ca2+ 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 firing-rate model of neuronal activity that is derived from a conductance-based model and that incorporates the spike-frequency adaptation, a main firing characteristic of regular spiking cortical pyramidal neurons. Such a calcium-based model can be generalized to a network of interconnected pyramidal neurons.

Our analysis permitted us to derive tau adap and Fadap in terms of the biophysical parameters of the membrane and intracellular Ca2+ dynamics. With plausible approximations, two relations between tau adap and Fadap were obtained. The first relation is given by Eq. 18: Fadap sime  alpha Gcctau adap, where alpha Gcc is a "negative feedback gain parameter" and is proportional to gCa and gAHP but independent of all [Ca2+] extrusion and buffering processes. The second relation is given by Eq. 19: Fadap sime  1 - tau adap/tau Ca, which predicts a negative proportionality between the two quantities as long as tau 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 IAHP = gAHP ([Ca2+]/([Ca2+] +KD)q × (V - EK), 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 tau adap and Fadap still held roughly with the slope comparable with, but larger than, 1/tau Ca. Another issue concerned other ionic currents such as the [Ca2+]-independent IM (Madison and Nicoll 1984) or/and the nonspecific cation IH (Lorenzon and Foehring 1992), which in some pyramidal neurons contribute to spike-frequency adaptation although to a much lesser extent than IAHP. In principle, their effects could be analyzed in the same way as we did for IAHP; the predicted relations between tau adap and Fadap 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 tau adap and Fadap while a current pulse is applied with different intensities or when the gAHP is blocked pharmacologically at various drug concentrations. These manipulations preserve tau Ca, and our theory predicts that Fadap and tau adap vary in such a way that they satisfy Fadap sime  1 - tau adap/tau Ca. Such measurements should provide an independent (if only crude) estimate for the [Ca2+] decay time tau Ca. Furthermore, because gAHP 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 gAHP could strongly regulate the spike-frequency adaptation dynamics while preserving the simple relation between tau adap and Fadap.

In addition, these relations could be used to assess how the biophysical mechanism for spike-frequency adaptation differs from cell to cell in neocortical circuits. Ahmed et al. (1993, 1997) have quantitatively described spike-frequency 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 (tau adap =11.5 ± 1.3 ms, and Fadap = 67 ± 3%) than deep layer cells(tau adap = 51.4 ± 6.4 ms and Fadap = 51 ± 5%). According to our theoretical results, this finding cannot be explained primarily by a difference in the [Ca2+] extrusion/buffering processes between superficial and deep layer neurons. If that was the case, the relation Fadap sime  alpha Gcctau adap would predict that a smaller tau adap should be correlated with a smaller Fadap, contrary to the observations.

Cellular selective attention and decorrelation

We investigated in detail the effect of spike-frequency adaptation to stochastic synaptic inputs of the Poisson type. In particular, we showed that the [Ca2+]-gated IAHP 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 IAHP 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 IAHP 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 IAHP (and IM) (McCormick and Williamsson 1989; Nicoll 1988), it is important to know how strong the spike-frequency 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 spike-frequency 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 IAHP), 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 IAHP is operative all the time and is responsible for generating negative temporal correlations between output ISIs in a frequency-dependent 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 spike-frequency adaptation. Computational implications of this observation will be explored in a separate work.

    ACKNOWLEDGEMENTS

  I thank G. Turrigiano for comments on the manuscript.

  This work was supported by National Institute of Mental Health Grant MH-53717-01), the Alfred P. Sloan Foundation, and the W. M. Keck Foundation.

    APPENDIX A: ISI RETURN MAP AND COEFFICIENT OF CORRELATION

For a neuronal spike train converted into a sequence of ISIs {t1, t2, t3, . . ., tN} (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 < tn+1|tn> versus tn. Let p(tn+1, tn) be the joint probability for tn+1 and tn. Then by definition
⟨(<IT>t<SUB>n<UP>+1</UP></SUB></IT>− μ)(<IT>t<SUB>n</SUB></IT>− μ)⟩ = <LIM><OP>∫</OP></LIM><LIM><OP>∫</OP></LIM>d<IT>t<SUB>n<UP>+1</UP></SUB></IT>d<IT>t<SUB>n</SUB></IT>(<IT>t<SUB>n<UP>+1</UP></SUB></IT>− μ)(<IT>t<SUB>n</SUB></IT>− μ)<IT>p</IT>(<IT>t<SUB>n<UP>+1</UP></SUB></IT>, <IT>t<SUB>n</SUB></IT>) (A1)
We can write p(tn+1, tn) = p(tn)p(tn+1|tn), where p(tn+1|tn) is the conditional probability for tn+1 given tn. Suppose that the conditional average of tn+1 given a fixed tn is linear with tn (for examples, see Figs. 8C and 11B)
⟨<IT>t<SUB>n<UP>+1</UP></SUB></IT>‖<IT>t<SUB>n</SUB></IT>⟩ = <LIM><OP>∫</OP></LIM><IT>t<SUB>n<UP>+1</UP></SUB>p</IT>(<IT>t<SUB>n<UP>+1</UP></SUB></IT>‖<IT>t<SUB>n</SUB></IT>)d<IT>t<SUB>n<UP>+1</UP></SUB></IT>≃ κ<IT>t<SUB>n</SUB></IT>+ κ<SUB>0</SUB> (A2)
where kappa  is the slope and kappa 0 is a constant. Then
⟨(<IT>t</IT><SUB>n+1</SUB>− μ)(<IT>t</IT><SUB><IT>n</IT></SUB>− μ)⟩
= <LIM><OP>∫</OP></LIM>d<IT>t<SUB>n</SUB></IT>(<IT>t<SUB>n</SUB></IT>− μ)<IT>p</IT>(<IT>t<SUB>n</SUB></IT>)&cjs0358;<LIM><OP>∫</OP></LIM>d<IT>t<SUB>n<UP>+1</UP></SUB></IT>(<IT>t<SUB>n<UP>+1</UP></SUB></IT>− μ)<IT>p</IT>(<IT>t<SUB>n<UP>+1</UP></SUB></IT>‖<IT>t<SUB>n</SUB></IT>)<SUB>+1</SUB>&cjs0359;
= <LIM><OP>∫</OP></LIM>d<IT>t<SUB>n</SUB></IT>(⟨<IT>t<SUB>n<UP>+1</UP></SUB></IT>‖<IT>t<SUB>n</SUB></IT>⟩ − μ)(<IT>t<SUB>n</SUB></IT>− μ)<IT>p</IT>(<IT>t<SUB>n</SUB></IT>)
≃ κ <LIM><OP>∫</OP></LIM>d<IT>t<SUB>n</SUB></IT>(<IT>t<SUB>n</SUB></IT>− μ)<SUP>2</SUP><IT>p</IT>(<IT>t</IT><SUB><IT>n</IT></SUB>) = κσ<SUP>2</SUP> (A3)
Therefore, CC = kappa . In other words, if the conditional average of tn+1 is plotted as function of its preceding value tn, the slope kappa  is identical to the ISI coefficient of correlation.

    APPENDIX B: DERIVATION OF SPIKE-FREQUENCY 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. 1-3, where [Ca2+] is considered as a parameter rather than a variable. For a constant [Ca2+], the IAHP = gAHP[[Ca2+]/([Ca2+] + KD)](V - EK) is a "passive" outward current. Given I, the larger is [Ca2+], so is IAHP, and the less the neuron should fire. We simulated the membrane Eqs. 1-3 with different [Ca2+] values. For each [Ca2+] value, the firing rate f was computed as well as the average < ICa> over an ISI; they are plotted in Fig. 2, B and C, respectively. These curves are nonlinear, in fact f = 0 for [Ca2+] > 2.2 µM. However, if we restrict ourselves to the range of [Ca2+] between 0 and [Ca2+]ss = 1.74 µM of Fig. 2A, the two curves in Fig. 2, B and C, can be approximated well by straight lines. This yields
<IT>f ≃ f</IT><SUB>0</SUB>− <IT>G</IT><SUB><IT>f</IT></SUB>[Ca<SUP>2+</SUP>]
⟨<IT>I</IT><SUB>Ca</SUB>⟩ ≃ ⟨<IT>I</IT><SUB>Ca</SUB>⟩<SUB>0</SUB>+ <IT>G</IT><SUB>cc</SUB>[Ca<SUP>2+</SUP>] (B1)
with f0 = 271, Gf = 84, < ICa> 0 = -28.8, and Gcc = 10 (Fig. 2, B and C, red solid lines).

In step 2, we turn our attention to the slowly evolving dynamics of [Ca2+], Eq. 4. Because [Ca2+] varies slowly, we substitute ICa by its average < ICa> which is a function of [Ca2+] itself (Eq. B1). Therefore, we have
<FR><NU>d[Ca<SUP>2+</SUP>]</NU><DE>d<IT>t</IT></DE></FR> ≃ −α⟨<IT>I</IT><SUB>Ca</SUB>⟩ − [Ca<SUP>2+</SUP>]/τ<SUB>Ca</SUB>
= − α(⟨<IT>I</IT><SUB>Ca</SUB>⟩<SUB>0</SUB>+ <IT>G</IT><SUB>cc</SUB>[Ca<SUP>2+</SUP>]) − [Ca<SUP>2+</SUP>]/τ<SUB>Ca</SUB>
= −α⟨<IT>I</IT><SUB>Ca</SUB>⟩<SUB>0</SUB>− (α<IT>G</IT><SUB>cc</SUB>+ 1/τ<SUB>Ca</SUB>)[Ca<SUP>2+</SUP>]
= −α⟨<IT>I</IT><SUB>Ca</SUB>⟩<SUB>0</SUB>− [Ca<SUP>2+</SUP>]/τ<SUB>adap</SUB> (B2)
with
<FR><NU>1</NU><DE>τ<SUB>adap</SUB></DE></FR>= α<IT>G</IT><SUB>cc</SUB>+ <FR><NU>1</NU><DE>τ<SUB>Ca</SUB></DE></FR> (B3)
Solving the averaged equation, we have
[Ca<SUP>2+</SUP>](<IT>t</IT>) = [Ca<SUP>2+</SUP>]<SUB>ss</SUB>(1 − exp(−<IT>t</IT>/τ<SUB>adap</SUB>)),
with
[Ca<SUP>2+</SUP>]<SUB>ss</SUB>= −α⟨<IT>I</IT><SUB>Ca</SUB>⟩<SUB>0</SUB>τ<SUB>adap</SUB> (B4)
Using the numerical values of the parameters in this expression, we obtain [Ca2+]ss = 1.77 µM and tau adap = 30.8 ms. These predicted values are close to those from numerical simulations of the whole model system ([Ca2+]ss = 1.74 µM, tau adap = 33 ms, respectively). Finally, in step 3, we insert the time course for [Ca2+] into the expression f = f0 - Gf[Ca2+], which yields
<IT>f</IT>(<IT>t</IT>) = <IT>f</IT><SUB>0</SUB><IT> − G</IT><SUB><IT>f</IT></SUB>[Ca<SUP>2+</SUP>](<IT>t</IT>) = <IT>f</IT><SUB>0</SUB><IT> − G</IT><SUB><IT>f</IT></SUB>[Ca<SUP>2+</SUP>]<SUB>ss</SUB>(1 − exp(−<IT>t</IT>/τ<SUB>adap</SUB>))
= <IT>f</IT><SUB>0</SUB><IT>− G</IT><SUB><IT>f</IT></SUB>[Ca<SUP>2+</SUP>]<SUB>ss</SUB>+ <IT>G</IT><SUB><IT>f</IT></SUB>[Ca<SUP>2+</SUP>]<SUB>ss</SUB>exp(−<IT>t</IT>/τ<SUB>adap</SUB>)
= <IT>f</IT><SUB>ss</SUB>+ (<IT>f</IT><SUB>0</SUB><IT>− f</IT><SUB>ss</SUB>) exp(−<IT>t</IT>/τ<SUB>adap</SUB>) (B5)
where fss = f0 - Gf[Ca2+]ss. Using the numerical values for f0 and Gf, we obtain f0 = 271 Hz and fss = 122 Hz, compared with numerical values of f0 = 272 Hz and fss = 116 Hz. The slight overestimate of fss is due to a small nonlinearity of the f-[Ca2+] curve (Fig. 2B).

    APPENDIX C: ADAPTATION TIME CONSTANTS BY TWO COUPLED CALCIUM MODES

With ICa and IAHP in both somatic and dendritic compartments, we attempted to generalize our fast-slow variable analysis to predict the biphasic spike-frequency adaptation time course. In step 1, we hold [Ca2+]s and [Ca2+]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 ICa,s and ICa,d on each of the [Ca2+]s and [Ca2+]d separately. For instance, if we consider the effect of the somatic IAHP on spike-frequency adaptation, we let [Ca2+]d = 0 while [Ca2+]s is varied, we have (Fig. B1A)
<IT>f = f</IT><SUB>0</SUB><IT> − G</IT><SUB><IT>f</IT>,s</SUB>[Ca<SUP>2+</SUP>]<SUB>s</SUB>
⟨<IT>I</IT><SUB>Ca,s</SUB>⟩ = ⟨<IT>I</IT><SUB>Ca,s</SUB>⟩<SUB>0</SUB> + <IT>G</IT><SUB>cc,ss</SUB>[Ca<SUP>2+</SUP>]<SUB>s</SUB>
⟨<IT>I</IT><SUB>Ca,d</SUB>⟩ = ⟨<IT>I</IT><SUB>Ca,d</SUB>⟩<SUB>0</SUB>+ <IT>G</IT><SUB>cc,ds</SUB>[Ca<SUP>2+</SUP>]<SUB>s</SUB> (C1)
with f0 = 272, < ICa,s> 0 = -17.2, < ICa,d> 0 = -28.3, Gf,d = 81.4, Gcc,ss = 5, and Gcc,ds = 8.5. On the other hand, the effect of the dendritic IAHP is considered with [Ca2+]s = 0 while [Ca2+]d is varied as parameter, we have (Fig. B1B)
<IT>f = f</IT><SUB>0</SUB><IT> − G</IT><SUB><IT>f</IT>,d</SUB>[Ca<SUP>2+</SUP>]<SUB>d</SUB>
⟨<IT>I</IT><SUB>Ca,s</SUB>⟩ = ⟨<IT>I</IT><SUB>Ca,s</SUB>⟩<SUB>0</SUB> + <IT>G</IT><SUB>cc,sd</SUB>[Ca<SUP>2+</SUP>]<SUB>d</SUB>
⟨<IT>I</IT><SUB>Ca,d</SUB>⟩ = ⟨<IT>I</IT><SUB>Ca,d</SUB>⟩<SUB>0</SUB>+ G<SUB>cc,dd</SUB>[Ca<SUP>2+</SUP>]<SUB>d</SUB> (C2)
with Gf,d = 75, Gcc,sd = 4.6, and Gcc,dd = 9. 


View larger version (12K):
[in this window]
[in a new window]
 
FIG. B1. Average firing rate f, < ICa,s> and < /ICa,d> as functions of [Ca2+]s (A) and [Ca2+]d (B) for the fast-slow variable analysis. , linear approximations.

The combined effect of both IAHP,s and IAHP,d is expressed as
<IT>f = f</IT><SUB>0</SUB><IT> − G</IT><SUB><IT>f</IT>,s</SUB>[Ca<SUP>2+</SUP>]<SUB>s</SUB> − <IT>G</IT><SUB><IT>f</IT>,d</SUB>[Ca<SUP>2+</SUP>]<SUB>d</SUB>
⟨<IT>I</IT><SUB>Ca,s</SUB>⟩ = ⟨<IT>I</IT><SUB>Ca,s</SUB>⟩<SUB>0</SUB><IT> + G</IT><SUB>cc,ss</SUB>[Ca<SUP>2+</SUP>]<SUB>s</SUB><IT> + G</IT><SUB>cc,sd</SUB>[Ca<SUP>2+</SUP>]<SUB>d</SUB>
⟨<IT>I</IT><SUB>Ca,d</SUB>⟩ = ⟨<IT>I</IT><SUB>Ca,d</SUB>⟩<SUB>0</SUB>+ <IT>G</IT><SUB>cc,ds</SUB>[Ca<SUP>2+</SUP>]<SUB>s</SUB><IT>+ G</IT><SUB>cc,dd</SUB>[Ca<SUP>2+</SUP>]<SUB>d</SUB> (C3)
In step 2, we substitute ICa,s and ICa,d by their averages Eq. C3, which yields two-coupled linear equations
<FR><NU>d[Ca<SUP>2+</SUP>]<SUB>s</SUB></NU><DE>d<IT>t</IT></DE></FR> = <IT>A</IT><SUB>c,s</SUB><IT> − B</IT><SUB>ss</SUB>[Ca<SUP>2+</SUP>]<SUB>s</SUB> − <IT>B</IT><SUB>sd</SUB>[Ca<SUP>2+</SUP>]<SUB>d</SUB>
<FR><NU>d[Ca<SUP>2+</SUP>]<SUB>d</SUB></NU><DE>d<IT>t</IT></DE></FR>= <IT>A</IT><SUB>c,d</SUB><IT>− B</IT><SUB>ds</SUB>[Ca<SUP>2+</SUP>]<SUB>s</SUB><IT>− B</IT><SUB>dd</SUB>[Ca<SUP>2+</SUP>]<SUB>d</SUB> (C4)
where Ac,s = -alpha s< ICa,s> 0, Ac,d = -alpha d< ICa,d> 0; Bss = alpha sGcc,ss + 1/tau Ca,s, Bsd = alpha sGcc,sd, Bds = alpha dGcc,ds, Bdd = alpha dGcc,dd + 1/tau Ca,d.

The solution of Eq. C4 has the form
[Ca<SUP>2+</SUP>]<SUB>s</SUB>(<IT>t</IT>) = [Ca<SUP>2+</SUP>]<SUB>s,ss</SUB> + <IT>c</IT><SUB>s1</SUB>exp(−<IT>t</IT>/τ<SUB>adap,1</SUB>) + <IT>c</IT><SUB>s2</SUB>exp(−<IT>t</IT>/τ<SUB>adap,2</SUB>)
[Ca<SUP>2+</SUP>]<SUB>d</SUB>(<IT>t</IT>) = [Ca<SUP>2+</SUP>]<SUB>d,ss</SUB>+ <IT>c</IT><SUB>d1</SUB>exp(−<IT>t</IT>/τ<SUB>adap,1</SUB>) + <IT>c</IT><SUB>d2</SUB>exp(−<IT>t</IT>/τ<SUB>adap,2</SUB>) (C5)
where the time constants are given by the eigenvalues of the matrix
− &cjs0358;<AR><R><C><IT>B</IT><SUB>ss</SUB></C><C><IT>B</IT><SUB>sd</SUB><IT></IT></C></R><R><C><IT>B</IT><SUB>ds</SUB><IT></IT></C><C><IT>B</IT><SUB>dd</SUB></C></R></AR>&cjs0359; (C6)

The fact that cd1 and cd2 have different signs implies that the time course [Ca2+]d(t) is not monotonic and exhibits a maximum. Indeed, we can find the time tmax when the maximum occurs by letting d[Ca2+]d/dt = 0, which yields
− <FR><NU><IT>c</IT><SUB>d1</SUB></NU><DE>τ<SUB>adap,1</SUB></DE></FR>exp(−<IT>t</IT><SUB>max</SUB>/τ<SUB>adap,1</SUB>) − <FR><NU><IT>c</IT><SUB>d2</SUB></NU><DE>τ<SUB>adap,2</SUB></DE></FR>exp(−<IT>t</IT><SUB>max</SUB>/τ<SUB>adap,2</SUB>) = 0 (C7)
or, solving this algebraic equation,
<IT>t</IT><SUB>max</SUB>= <FR><NU>τ<SUB>adap,1</SUB>τ<SUB>adap,2</SUB></NU><DE>τ<SUB>adap,2</SUB>− τ<SUB>adap,1</SUB></DE></FR>ln &cjs0358;<FR><NU>−<IT>c</IT><SUB>d1</SUB>τ<SUB>adap,2</SUB></NU><DE><IT>c</IT><SUB>d2</SUB>τ<SUB>adap,1</SUB></DE></FR>&cjs0359; = 112 ms (C8)
which compared well with the numerical value of the peak time (= 106 ms).

Finally, in step 3, by inserting the solutions for [Ca2+]s(t) and [Ca2+]d(t) into the expression of the firing rate (Eq. C3), we obtain
<IT>f</IT>(<IT>t</IT>) = <IT>f</IT><SUB>ss</SUB><IT>+ b</IT><SUB>1</SUB>exp(−<IT>t</IT>/τ<SUB>adap,1</SUB>) + <IT>b</IT><SUB>2</SUB>exp(−<IT>t</IT>/τ<SUB>adap,2</SUB>) (C9)
with fss = 93, b1 = 153.6, b2 = 23.7 (in Hz).

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 steady-state 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 cross-product terms like [Ca2+]s × [Ca2+]d in Eq. C3. An improved nonlinear analysis is beyond the scope of this paper.

    FOOTNOTES

  Address reprint requests to X.-J. Wang.

  Received 8 July 1997; accepted in final form 3 December 1997 .

    REFERENCES
Abstract
Introduction
Methods
Results
Discussion
References

0022-3077/98 $5.00 Copyright ©1998 The American Physiological Society