The Journal of Neurophysiology Vol. 79 No. 3 March 1998, pp. 1549-1566
Copyright ©1998 by the American Physiological Society

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 Ca²⁺ conductance (g_Ca) and a Ca²⁺-dependent K⁺ conductance (g_AHP) 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 [Ca²⁺] transients are matched with recent dendritic [Ca²⁺] imaging measurements. Spike-frequency 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²⁺] 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²⁺] decay time constant _Ca, is assessed quantitatively. In addition, we demonstrate that the intracellular [Ca²⁺] 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 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 spike-adapting current I_AHP, as delayed inhibition through intracellular Ca²⁺ 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²⁺]-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 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 Ca²⁺ 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 Ca²⁺-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) = f_ss + (f₀  f_ss) exp(t/_adap), where f₀ is the initial firing rate, f_ss is the steady-state 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₀  f_ss)/f₀. 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 spike-frequency 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 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, Ca²⁺-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²⁺ that the intrinsic gating of I_AHP is rapid; its slow activation is thus attributable to the kinetics of the cytoplasmic calcium concentration [Ca²⁺] (Lancaster and Zucker 1994). Third, spike-evoked [Ca²⁺] 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²⁺] indicator dye like Fura-2 is also a [Ca²⁺] buffer, Helmchen et al. (1996) used increasingly low concentrations of Fura-2 and, by extrapolation to zero dye concentration, obtained measurements of putatively intrinsic [Ca²⁺] transient signals. Their estimated spike-evoked [Ca²⁺] 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, [Ca²⁺] 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) Ca²⁺ imaging measurements from pyramidal cells show spike-evoked [Ca²⁺] 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²⁺] influx is produced primarily by Ca²⁺ 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 [Ca²⁺]-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 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), I_Ca, and a calcium-dependent 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 current-balance equations

(1)

(2)

where C_m = 1 µF/cm² and I_L = g_L (V  V_L) is the leak current. Following Pinsky and Rinzel (1994), we express the current flows between the soma and dendrite [proportional to (V_s  V_d)] in microamperes per square centimeter, with the coupling conductance g_c = 2 mS/cm², 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/cm²) to the soma or by a random synaptic input I_syn of the -amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid type to the dendrite. I_syn =g_syn s(V  E_syn); the gating variable s obeys the equation ds/dt =(t)  s/_s, where (t) is a Poisson point-process with a rate ,E_syn = 0 mV, _s = 0.5 ms, and g_syn = 0.08 mS/cm² (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

(3)

The sodium current I_Na = g_Nam³(V)h(V  V_Na), where the fast activation variable is replaced by its steady state, m = _m/(_m + _m), _m = 0.1(V + 33)/{exp[0.1(V + 33)]  1}, _m =4 exp[(V + 58)/12], _h = 0.07 exp[(V + 50)/10], and_h = 1/{exp[0.1(V + 20)] + 1}. The delayed rectifier I_K =g_Kn⁴(V  V_K), where _n = 0.01(V + 34)/{exp[0.1(V +34)]  1}, and _n = 0.125 exp[(V + 44)/25]. The temperature factor _h = _n = 4.

The high-threshold calcium current I_Ca = g_Cam(V  V_Ca), where m is replaced by its steady-state m(V) = 1/{1 + exp[(V + 20)]/9} (Kay and Wong 1987). The voltage-independent, calcium-activated potassium current I_AHP = g_AHP[[Ca²⁺]/([Ca²⁺] + K_D)](V  V_K), with K_D = 30 µM. The intracellular calcium concentration [Ca²⁺] is assumed to be governed by a leaky-integrator (Helmchen et al. 1996; Tank et al. 1995; Traub 1982)

(4)

where  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  = 0.002 [in µM (msµA)¹cm²] so that the [Ca²⁺] 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 constant_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 I_AHP 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): g_L = 0.1, g_Na = 45,g_K = 18, dendritic g_Ca = 1, and g_AHP = 5 (in mS/cm²); 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.

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 adap,1 = 29.4 ms and 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 {t₁, t₂, t₃, . . . , 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)

(6)

The statistical interdependence between consecutive ISIs is measured by the coefficient of correlation (CC) (Perkel et al. 1967; Tuckwell 1988)

(7)

which is between 1 and +1. The CC can be interpreted as follows. Suppose that we have an ISI return map, where t_n+1 is plotted against t_n. To assess how ISIs (t_n+1) depend on their preceding values (t_n), we calculate the conditional average of t_n+1 for each given t_n, t_n+1|t_n. If this function of t_n 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 Ca²⁺ accumulation, at a rate of 200 nM/spike [comparable with the [Ca²⁺] imaging measurements from proximal apical dendrites of cortical layer V pyramidal cells (Helmchen et al. 1996)]. The I_AHP increases with [Ca²⁺], 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²⁺] dynamics, when the spike-evoked [Ca²⁺] influx rate is balanced with the [Ca²⁺] decay rate. After the current pulse, there is a long-lasting AHP that mirrors the Ca²⁺ (hence the I_AHP) 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 I  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 f₀ and the final f_ss increases with the current intensity. As a result, the steady-state input-output relation is linearized by the I_AHP. We also computed the mean dendritic [Ca²⁺] as I was varied. Its steady-state plateau level depends linearly on f_ss (Fig. 1B, 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²⁺ 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  = 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/_adap), where A = f_ss, and B = f₀  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₀  f_ss)/f₀ = B/(A + B) = 57%. The [Ca²⁺] also follows an exponential time course with the same time constant _adap, the steady-state plateau is [Ca²⁺]_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. 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+]. , data obtained when [Ca2+] was varied as a parameter. , 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 Ca²⁺ evolves much more slowly than the membrane potential and the other channel gating variables of the model and that this slow [Ca²⁺](t) determines the adaptation time course. Thus we can first analyze the fast subsystem and the slow [Ca²⁺] 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 [Ca²⁺] as if it was a static parameter rather than a dynamical variable. This allows us to determine the functional [Ca²⁺] dependence of the firing frequency f([Ca²⁺]) as well as other voltage-dependent quantities like I_Ca, where x denotes an average of x over a typical ISI. The functional [Ca²⁺] dependence of f([Ca²⁺]) and of I_Ca was found to be approximately linear (Figs. 2, B and C)

(8)

where f₀ = 271 Hz is the initial firing frequency, and G_f = df/d[Ca²⁺] = 84 (in Hz/µM) is a negative-feedback "gain parameter" for the firing frequency. I_Ca0 = 28.8 µA/cm² is the initial I_Ca, G_cc = dI_Ca/d[Ca²⁺] = 10 (in µA/cm² µM) is a negative-feedback gain parameter for the [Ca²⁺]-dynamics (Ahmed et al., unpublished results), and G_cc (in 1/ms) is the rate at which the [Ca²⁺] influx is reduced by [Ca²⁺] itself.

In step 2, we turn to the [Ca²⁺] 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²⁺]. The resulting equation now only depends on [Ca²⁺]

(9)

with

(10)

Solving Eq. 9, we obtain

(11)

with [Ca²⁺]_ss = I_{Ca0 adap}. Note that the adaptation time constant _adap (Eq. 10) is always smaller than _Ca due to the presence of the negative feedback term G_cc. For instance, with  = 0.002 and G_cc = 10, G_cc = 0.02 is larger than 1/_Ca = 0.0125, hence contributes more to _adap. We have _adap = 30.8 ms, whereas _Ca = 80 ms; and [Ca²⁺]_ss = 1.77 µM.

Finally, in step 3, we insert the time course for [Ca²⁺] into the expression f = f₀  G_f[Ca²⁺], which yields the adaptation process for the firing frequency as function of time

(12)

with f_ss = f₀  G_f[Ca²⁺]_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 f_ss. This is mainly due to the fact that f([Ca²⁺]) is not exactly linear. Indeed, f goes to zero at a critical value of [Ca²⁺]  2.2 µM (when the I_AHP 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([Ca²⁺]) behaves as a square-root function of [Ca²⁺] (see Guckenheimer et al. 1997; Rinzel and Ermentrout 1987). We found that a square-root function fits well even the global f([Ca²⁺]): f([Ca²⁺]) = 265 , except near the bifurcation threshold (Fig. 2B, blue curve). Using this nonlinear expression for f([Ca²⁺]) and the same [Ca²⁺](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

(13)

which depends on the [Ca²⁺] dynamics and the I_AHP only through the factor G_fadap.

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

(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 f₀, G_f, I_Ca0, and G_cc are shown in Fig. 3A. We observe that _adap is larger and F_adap 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 Ca²⁺ influx per spike is reduced. All four curves in Fig. 3A can be fitted reasonably well by logarithmic functions

(15)

where I₀ = 0.3 is the estimated rheobase (The actual rheobase is somewhat larger, I₀  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, ICa0, and Gcc as functions of I. , empirical fitting functions of Eq. 15. B: adaptation time constant 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 [Ca²⁺] 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 [Ca²⁺] dynamics. Because of the adapting I_AHP, the effective [Ca²⁺] time constant is given by _adap, typically much shorter than _Ca. Hence, one can expect that the [Ca²⁺] 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²⁺] 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 g_c

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 g_c, 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 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. 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 [Ca²⁺] influx per spike is reduced (Fig. 5A, right). This leads to a slower [Ca²⁺] accumulation, larger _adap, and smaller percentage adaptation F_adap (Fig. 5B). As in the case of varying the applied current I (Fig. 3B), the two quantities _adap and F_adap 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: 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 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 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 _adap and F_adap

In addition to the neuronal electrotonic structure, spike-frequency adaptation depends also on the channel conductances g_Ca and g_AHP, as well as the [Ca²⁺] kinetic parameters  and _Ca. These dependences were explored within the framework of our calcium-model. First, the initial firing rate f₀ should not depend on any of these parameters. Second, because I_Ca = I_Ca0 + G_cc[Ca²⁺], I_Ca0 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

(16)

Inserting these scaling relations into Eqs. 10 and 13 we have

(17)

From Eq. 17 we can conclude that _adap and F_adap depend only on _Ca (a neuron's [Ca²⁺] extrusion and buffering properties) and on the combination g_Cag_AHP (the product of the spike-evoked [Ca²⁺]-influx size g_Ca and the adaptation conductance g_AHP). Given fixed g_Cag_AHP, as _Ca is increased (Fig. 6A), the adaptation is slower (_adap is larger), but the plateau [Ca²⁺] level is higher (cf. Eqs. 4 and 11), hence F_adap is larger. A plot of F_adap versus _adap is linear with a positive slope (Fig. 6B). This is evident in Eq. 17, where f₀, _f, and I_Ca0 (computed with [Ca²⁺] as a static parameter) are all independent of _Ca, hence F_adap is simply proportional to _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 adap and Fadap on Ca (A and B) and the combination gCagAHP (C and D) (cf. Eq. 17). Three curves in each panel correspond to different applied current intensities. B: Fadap vs. adap from data in A is linear with a positive slope (which depends on I). D: Fadap vs. adap from data in C: is linear with a negative slope  1/Ca. This is also true for data obtained with I (Fig. 3B) or gc (Fig. 5B) being varied as parameter.

By contrast, if _Ca is held constant and g_Cag_AHP is increased (Fig. 6C), 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. 6D). In Fig. 6D, we also plotted the F_adap versus _adap simulation data, obtained when the input current intensity I is varied (Fig. 3B) and when the electrotonic coupling g_c is varied (Fig. 5B). 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²⁺] 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₀/G_f  (I_Ca0)/G_cc. With this assumption, we then can write I_Ca0 G_f/f₀  G_cc. Inserting this relation into Eq. 13, we have

(18)

or because G_cc = 1/_adap  1/_Ca (Eq. 10), we have

(19)

which is the desired relation between F_adap and _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 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²⁺] signaling and [Ca²⁺]-dependent spike-frequency adaptation are distributed across multiple compartments. For our two-compartment 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²) for both compartments. We then have two equations like Eq. 4 for somatic and dendritic [Ca²⁺], respectively

(20)

Because the parameter  is proportional to the surface:volume ratio, it should be much smaller for the soma than for the dendrite. Similarly, _Ca is expected to be longer at the soma than at the dendrite. For instance, the spike-evoked [Ca²⁺] transient peak amplitude and decay rate (1/_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, _d = 0.002 and _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 _s = 0.000667 and _Ca,s = 240 ms. As shown in Fig. 7, with two calcium modes, [Ca²⁺]_s(t), [Ca²⁺]_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)

(21)

where _adap,1 = 29.4 ms and _adap,2 = 191 ms. None that the faster [Ca²⁺] 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 [Ca²⁺]_d can be nonmonotonic (Fig. 7). In the first rapid phase of adaptation, [Ca²⁺]_s is small, both [Ca²⁺]_d and the firing rate converge to their plateau levels as if the somatic I_AHP did not exist. But as [Ca²⁺]_s (hence the somatic I_AHP) slowly accumulates, the firing rate further decreases in the second phase of adaptation, which leads to a decay of [Ca²⁺]_d to its actual final plateau, which is smaller than would be expected without a somatic I_AHP. Hence, the nonmonotonic behavior of [Ca²⁺]_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 _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²⁺] 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 . 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 [Ca²⁺]. 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 I_Ca. As shown in Fig. 8B, both f and I_Ca are approximately linear with [Ca²⁺] when the latter is considered as a parameter

(22)

where f₀ = 213 Hz, G_f = 252 Hz/µM, I_Ca0 = 22.6 µA/cm², and G_cc = 27.5 µA/cm² µM. Solving Eq. 4 with I_Ca substituted by I_Ca, we obtain

(23)

where _adap = 14.8 ms, [Ca²⁺]_ss = 0.67 µM, f_ss = 44.2 and f₀ = 213 (in Hz). (Note, again, _adap is much shorter than _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 I_AHP. Statistically, if by chance the cell happens to fire rapidly in a short time window, significant [Ca²⁺] 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²⁺] 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 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 steady-state input-output relation, i.e., the mean firing rate f versus the Poisson input rate , is linear (Fig. 9A, red curve). We also computed the ISI CV and CC as function of