INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX REFERENCES

Oscillations have long been observed in neuronal structures (Adrian and Matthews 1934) but their role, mechanisms, and interplay with single neuron biophysical characteristics have only recently been submitted to detailed scrutiny. Experiments have tested the response of neurons to oscillating current injection. Subthreshold resonance, in which the response of the induced oscillating voltage peaks at a preferred input frequency, has been found in inferior olive neurons (De Zeeuw et al. 1998; Lampl and Yarom 1993, 1997; Llinas and Yarom 1986), trigeminal root ganglion neurons (Puil et al. 1986), thalamic neurons (Hutcheon et al. 1994; Jahnsen and Karnup 1994; Puil et al. 1994), cortical neurons (Dickson et al. 2000; Gutfreund et al. 1995; Hutcheon et al. 1996b; Llinas et al. 1991), and both hippocampal CA1 pyramidal cells (Leung and Yu 1998; Pike et al. 2000) and interneurons (Pike et al. 2000). Many of these structures are known to support oscillations in vivo, suggesting an interplay between single-cell frequency preference and oscillations at the network level. Most of the recorded neurons show a single peak at a finite frequency in their voltage response. However, some interneurons of the hippocampus show a more complex response with a trough at low frequency followed by a peak at higher frequencies (Pike et al. 2000). Although a great deal of effort has been directed at understanding the input properties of resonant neurons, surprisingly little attention has been addressed to the effect of subthreshold resonance on the temporal properties of the firing rate. This is despite the common assumption that the presence of resonant neurons might provide a stabilizing influence on oscillations at the level of the network.

It is known from Hodgkin and Huxley (1952) and many studies since (see e.g., Gutfreund et al. 1995; Hutcheon et al. 1996a, 1994; Koch 1984; Mauro et al. 1970; Rinzel and Ermentrout 1989; White et al. 1995) that the resonance properties of neurons can be related to their ionic channel characteristics through a mathematical linearization of the corresponding conductance-based description. Several scenarios involving voltage-gated ionic currents have been shown to generate resonant behavior (for a review, see Hutcheon and Yarom 2000). Reduced two-variable descriptions have proven useful as a mathematical tool to study these and other neuronal properties (Gutfreund et al. 1995; Hutcheon et al. 1996a; Rinzel and Ermentrout 1989; White et al. 1995).

In the first part of this paper, a systematic classification of two-variable models is provided. The analysis highlights the possible types of subthreshold behavior associated with different neuronal characteristics. The results can be summarized in a graphical description. The change of membrane properties as the neuron is depolarized toward threshold is represented by trajectories crossing boundaries separating different types of behavior (e.g., passive from resonant). This description is illustrated with two conductance-based model neurons. More complex types of resonance cannot be described by a two-variable model. For this reason, a three-variable model, which exhibits a richer repertoire of behaviors, is also analyzed. A parameter region is identified with a suppression as well as a resonance in the impedance curve, a feature recently observed experimentally in hippocampal fast-spiking interneurons (Pike et al. 2000).

In the second part of the paper, the circumstances are examined in which a resonant neuron can communicate its subthreshold frequency preference through the dynamics of its firing rate. This property of resonant neurons manifests itself as a preferential amplification of input signals that are at the resonant frequency and requires an analysis of how the firing rate is modulated by an oscillatory current. To this end, the two-variable approach is extended to include spike emission, providing a generalized integrate-and-fire or GIF model. The model captures a wide range of subthreshold dynamics with a simplified firing and reset mechanism. The firing-rate dynamics of this model, as well as two specific conductance-based models which exhibit a subthreshold resonance, are studied in detail. The crucial role that noise plays in shaping the response is highlighted.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX REFERENCES

Glossary

	deviation of the membrane potential from the holding potential (mV).
C or CM	membrane capacity (nF).
g	effective leak conductance (µS).
w or w1	auxiliary variable characterizing the membrane dynamics (mV).
1	time scale of the dynamics of the w variable (ms).
g1	conductance measuring the membrane potential variation resulting from a change of w (µS).
	dimensionless parameter proportional to leak conductance g.
	dimensionless parameter proportional to conductance g1.
w2	second auxiliary variable (mV).
g2	analogous to g1 for the second variable w2 (µS).
Iapp	total applied external current (nA).
Isyn	synaptic current (nA).
geo and gio	average excitatory and inhibitory total synaptic conductances (µS).
e and i	magnitude of excitatory and inhibitory synaptic noise (µS).
e and i	correlation timescales of the excitatory and inhibitory synaptic noise (ms).
IN	magnitude of the fluctuations of synaptic current (nA).
I0	constant (DC) current (nA).
I1	magnitude of oscillatory current (nA).
f	frequency of injected current (Hz).
V	strength of the synaptic noise as measured by the resulting amplitude of membrane potential fluctuations (mV).
Z(f)	cell impedance for an injected current of frequency f (M).
fR	resonant frequency corresponding to a maximum of the amplitude of Z(f) (Hz).
f0	natural frequency of the membrane potential damped oscillations (Hz).
Q	strength of the resonance peak (dimensionless).
r0	average spike rate (Hz).
r1(f)	magnitude of oscillatory component in spike rate induced by injected oscillatory current (Hz).
|A(f)|	signal gain (Hz/nA).
(f)	phase of oscillatory component in spike rate with respect to oscillatory current (deg).
	threshold for spike emission for the GIF model (defined in METHODS; mV)
r	membrane potential reset after spike emission for the GIF model (defined in METHODS; mV)

Linearization of conductance-based models

The starting point for the analysis in this paper is the conductance-based Hodgkin-Huxley formalism. The state of a neuron is described by a potential difference V across a membrane with a capacitance C_M, a set of trans-membrane currents I_mem (comprising the leak and various active ionic currents), a synaptic current I_syn (to be described in the following text) and an applied current I_app

(1)

The active ionic currents comprise both activation and inactivation variables x_k where k = 1, ..., N counts over all the variables that obey equations of the form

(2)

where both the relaxation times _k(V) and the steady-state values x_k,(V) are functions of the membrane voltage.

Below threshold for action potential generation, Eqs. 1 and 2 can be linearized around a holding voltage V* (see e.g., Koch 1999, chapter 10, and refs therein). For the sake of simplicity, the notation X* will be used to denote the quantity X evaluated at V = V*. Linearization of the Eq. set 1 and 2 allows for a direct categorization of the range of behavior that a neuron exhibits in its response to small input currents, for example, the response to an oscillating or square-pulse current considered here. The linearized equations will also provide the basis for a generalization of the IF model, to be described at the end of this section. The linear equations can be written in the following form

(3)

with  = V  V* being the deviation of the voltage from its steady-state value and g_M = (I_mem/V)* is the slope of the instantaneous I-V curve. The time-dependent variables

(4)

are proportional to the deviation of the activation or inactivation variables x_k from their steady-state values x^*_k = x_k, (V*) and are expressed in units of millivolts. The time constants _k correspond to those of the activation and inactivation variables evaluated at V* and the parameters

(5)

written in units of conductance, measure the strength of the effect that the variable x_k has on the voltage. Note that in the linear approximation, the dynamical variables w_k are no longer multiplied by a voltage-dependent term as they were in the original conductance-based description.

Activation or inactivation variables can be classified according to the sign of their corresponding parameter g_k. Examples of variables with g_k < 0 are the activation variables of Na⁺ and Ca²⁺ currents and inactivation variables of K⁺ currents. Examples of variables with g_k > 0 include inactivation variables of Na⁺, Ca²⁺ currents, activation variables of K⁺ currents, and the activation variable of the H current. For g_k > 0, the corresponding variable opposes voltage change (negative feedback), whereas g_k < 0 indicates that the variable amplifies voltage change (positive feedback). Previous modeling studies (reviewed in Hutcheon and Yarom 2000) have shown that a variable with g_k > 0 can create a subthreshold resonance (a resonant variable), whereas a variable with g_k < 0 can amplify an existing resonance (an amplifying variable).

Conductance-based models generally comprise many active ionic currents and are therefore described by a large number of activation or inactivation variables. Despite the simplification of linearity, such systems of equations can still be hard to handle analytically. However, it is often possible to reduce the number of variables to two or three, while still accurately modeling the behavior near the holding voltage. This can be achieved by considering that very fast variables (such as the activation variable of fast sodium channels) are instantaneous, by merging together variables with similar time constants, and by noting that very slow variables average over the voltage to provide a steady current. The resulting equations have the same form as Eq. 3 but with effective values C and g for the capacitance and leak respectively. The effective leak g can be zero or even negative, while the resting potential remains stable. Examples of the linearization method are given in the APPENDIX for two conductance-based models together with the further approximations leading to reductions in the number of variables to two or three.

TWO-VARIABLE SUBTHRESHOLD DYNAMICS. For the case of two variables, the neuron is described by the two equations

(6)

with four parameters, C, g, g₁, and ₁. However, expressing time in units of ₁, and dividing the voltage Eq. 6 by C makes it apparent that the model only depends on two dimensionless parameters  = g₁/C and  = g₁₁/C. The quantities  and  parameterize the behavior of the neuron near V* and can be considered as representing a point on a plane.  represents an effective leak, whereas  represents an effective coupling between the two variables.  measures the influence of the w variable on the membrane potential.

THREE-VARIABLE SUBTHRESHOLD DYNAMICS. The analysis is also extended to include a third variable. The subthreshold dynamics is then described by

(7)

where in this paper a restriction is made to g > 0. Four independent and dimensionless parameters g₁/g, g₂/g, ₁g/C, and ₂g/C are now needed to fully describe the model.

Models of spiking neurons

One of the major goals of the present paper is to investigate the effect of subthreshold resonance on the dynamics of spike emission. To this end, a simple spiking neuron model that exhibits subthreshold resonance, the generalized IF neuron, is introduced. To demonstrate that the general results derived for this simplified model carry over to more realistic neurons, two representative conductance-based models that are known from the literature to produce subthreshold resonance are also examined (Models I and II). In an attempt to cover the range of possible behaviors, the models are chosen to have different resonance mechanisms (hyperpolarization or depolarization activated currents) and also different resonant frequencies (near 10 and 50 Hz, respectively).

GENERALIZED IF NEURON. The IF model neuron provides a powerful tool for the understanding of neurons with passive membrane properties and is the standard component of large numerical simulations of recurrent networks. However, its passive subthreshold behavior cannot capture the phenomenon of resonance. An extension of the IF model, which captures the subthreshold behavior of the two-variable model with a simple spike mechanism, is therefore the first spiking neuron model to be introduced here. The generalized IF (GIF) neuron is obtained by supplementing Eq. 3 with a threshold for spike generation at  = , followed by a reset of the membrane voltage at  = _r (the auxiliary variables w_k considered here have a slower dynamics than the spike, and therefore it is not appropriate to reset them also). In the case g_k = 0 for all k, the voltage equation reduces to the IF model. With two variables, this model is similar to a model recently proposed by Izhikevich (2001). The two-variable GIF model subject to an applied current is described by Eq. 6 where the parameters C, g, g₁, and ₁ are kept fixed for the whole of the subthreshold regime. I_syn is the modeled synaptic current and I_app(t) represents the applied current, to be described in the following text. In this paper, the threshold is chosen to be at 50 mV, the rest (in absence of any input currents) at 70 mV, and the reset at 56 mV. Because  measures the deviation from rest  = V  V_rest, this corresponds to  = 20 mV and _r = 14 mV.

Conductance-based neurons

MODEL I. A NEURON WITH I_NA, I_K, AND I_H CURRENTS. The first model comprises a hyperpolarization-activated mixed cation current I_H and the Hodgkin-Huxley spike-generating currents. The form of the I_H current is taken from Spain et al. (1987) and comprises both fast f and slow s activation variables. The time scales of the two components are _f = 38 ms and _s = 319 ms and, as in Spain et al. (1987), taken to be voltage independent. The fast component has the greater contribution and determines the resonant frequency f_R, which is near 10 Hz at physiological temperatures. A detailed model description can be found in the APPENDIX.

MODEL II. A NEURON WITH I_NA, I_K, I_NAP, AND I_KS CURRENTS. In contrast to the I_H model defined in the preceding section, the second model neuron features two depolarization-activated currents: the slow potassium current I_Ks and the persistent sodium current I_NaP. In the language of Hutcheon and Yarom (2000), the I_Ks current generates the resonance and the I_NaP current amplifies its effect. A noninactivating form (Gutfreund et al. 1995) is used for the I_Ks with an activation time scale of _q = 6 ms (Wang 1993), giving a subthreshold resonance that is strongest at resonant frequencies ~35-55 Hz. Neurons with a resonance frequency or subthreshold oscillations at ~40 Hz are widespread (Pike et al. 2000; Puil et al. 1986), and the underlying mechanism is thought to sometimes involve the I_Ks and I_NaP currents (Llinas et al. 1991). Again, full details of this model are given in the APPENDIX.

Modeling the noisy synaptic input

The massive synaptic bombardment received by neurons in vivo represents a strong source of noise. Destexhe et al. (2001) provided evidence that an appropriate model of such a synaptic input is given by a fluctuating conductance with short correlation times related to the shapes of typical excitatory and inhibitory postsynaptic potentials. Hence, for the analysis of the firing rates of the two conductance-based models I and II, the noise is modeled as in Destexhe et al. (2001) by the equations

(8)

where _e, _i are delta-correlated Gaussian white-noise terms. Representative values of the correlation times _e = 3 ms and _i = 10 ms are used here. The reversal potentials are taken to be E_e = 0 mV and E_i = 75 mV. The average conductances g_eo, g_io and the noise amplitudes _e and _i can be varied to explore a range of input conditions.

MODELING NOISE FOR THE GIF NEURON. To have a simple model with a membrane-potential-independent subthreshold resonance, the synaptic inputs are modeled by a current comprising a direct drive I₀ and a white-noise source

(9)

where (t) is the delta-correlated Gaussian white-noise term with unit variance and I_N is the measure of the noise strength in nanoAmpere. The factor _N is introduced to preserve units, and throughout this paper, it is arbitrarily fixed at _N = 1 ms without affecting the generality of the results. The noise strength I_N can be related to a more intuitive measure: the SD of the membrane voltage (in the absence of the spiking mechanism). In the two-variable GIF model, the SD of the voltage _V takes the form

(10)

Response of the neuron to an oscillatory drive

SUBTHRESHOLD RESPONSE. To characterize the subthreshold response, an oscillating current of frequency f is used. The applied current and resulting voltage response to this current are given by

(11)

where both the phase difference (f) and the magnitude of the impedance

(12)

are functions of the driving frequency f. The existence of a peak in the |Z(f)| versus frequency curve provides the definition of subthreshold resonance. The impedance Z(f) can also be measured experimentally using a ZAP current (Puil et al. 1986).

FIRING-RATE RESPONSE. In the context of examining the interactions between membrane frequency preference (resonance) and network oscillations, it is of interest to examine how the instantaneous firing rate of a neuron responds to a sine-wave modulation in the background of a noisy synaptic current
A regime is considered where the noisy synaptic drive is sufficiently strong to cause the neuron to fire stochastically, at an average rate r₀. The weak sinusoidal component then causes a weak modulation of the firing rate that will be apparent over many trials, see Fig. 1. This quantity can also be thought of as the firing rate, averaged over a population of neurons each individually receiving a noisy drive but responding collectively to the same weak oscillatory component present in the firing rates of presynaptic neurons. The form of this population, or trial-averaged instantaneous rate is

(13)

The analogy with the subthreshold voltage form in Eq. 11 is clear. Similarly, the response r₁(f) is proportional to the strength of the modulatory current, leading to the introduction of the following quantity that measures the ability of a neuron to amplify a particular frequency

(14)

called the signal gain, see for example (Gerstner 2000). In the same way that a peak in the impedance |Z(f)| quantifies subthreshold resonance, the existence and position of the peak in the quantity |A(f)| will be the corresponding firing-rate measure of resonance.

View larger version (49K): [in this window] [in a new window]   Fig. 1. The response of the instantaneous firing rate to a weak sinusoidal input. The oscillatory current is applied on top of a noisy input that itself elicits firing at an average rate r0 (see top). By averaging over many realizations (shown in the raster plot), the modulation of the instantaneous firing rate can be computed, and the characteristics (amplitude and phase) of the induced sinusoidal component of the firing rate obtained (bottom). A bin width of 8 ms was used for illustrative purposes in this figure.

ANALYTICAL METHODS. The firing-rate response of the GIF neuron can be computed analytically in the limit of large ₁. Methods are sketched in the APPENDIX (see also Brunel et al. 2003).

NUMERICAL METHODS. The numerical analysis of the firing-rate response of the GIF and conductance-based models was performed using a stochastic second-order Runge-Kutta algorithm (Honeycutt 1992) with a time step of 10 and 20 µs, respectively. The amplitude of the modulatory current I₁ used in numerical measurements of the signal gain was varied until it was sufficiently small such that higher order nonlinear effects were negligible. The length of simulation time needed to get accurate measurements for each frequency point varied between 1,000 and 50,000 s, depending on the firing rate and the particular level of noise chosen in the input current. To estimate the firing rate modulation at a given frequency, the instantaneous firing rate is computed in bins of 1 ms. The resulting histogram, sketched in Fig. 1, is then fitted by a sinusoid with a frequency equal to that of the oscillatory input current.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX REFERENCES

Subthreshold properties of the membrane potential

The different classes of behavior in the subthreshold regime are examined first. Two- and three-variable models, with parameters directly related to measurable membrane properties, are used to classify the different types of response to standard test currents. It is shown that the types of behavior that the neuron can exhibit at different holding voltages can be conveniently presented in graphical form for both the two- and three-variable descriptions. The results are illustrated by two conductance-based models of spiking resonant neurons.

Subthreshold behavior of the two-variable model

The subthreshold behavior of the two-variable model is first classified with respect to stability and the response to standard test currents. Because the two-variable model is related to an underlying conductance-based description near a holding voltage V*, it is parameterized by an effective leak  = g₁/C and an effective coupling between the two variables  = g₁₁/C. The parameters  and  can be used to represent the behavior of the neuron by a string of points on a plane, as the neuron is depolarized or hyperpolarized by an injected current. The borders separating different types of behavior are obtained through the analysis of Eq. 6, presented in detail in the APPENDIX.

STABILITY. The classification of the subthreshold regime starts with the determination of the parameter region where the neuron remains stable at the holding potential (without e.g., subthreshold oscillations or spike emission). Analysis of the stability of the membrane potential of the two-variable model determines an unstable region shown in brown in Fig. 2. The region is bounded on one side by the vertical dashed line that signals the onset of spontaneous oscillations. On the other side, it is bounded by a diagonal line that corresponds to the total input conductance becoming zero, which can lead to spike emission. The rest of the analysis will focus on the stable region to the right of these two lines.

View larger version (34K): [in this window] [in a new window]   Fig. 2. The subthreshold behavior of two-variable models. Brown marks the unstable region. Insets: the qualitative response to the relevant test current in regions in the space of parameters  = g1/C (the effective leak) and  = g11/C (the effective coupling between the two variables). A: oscillating current injection, amplitude of the impedance. B: phase of the impedance. C: response to a square-pulse current. D: trajectories of the two conductance-based model neurons in parameter space as their holding potential is increased. Model I (full line): g1 = gf (see text and APPENDIX) parameterizes the effect of the fast component of the current IH on the voltage. The trajectory covers the range from 100 to 56.5 mV at which point the resting state is unstable and spikes are emitted. Model II (dotted line): g1 = gq (see text and APPENDIX) parameterizes the effect of the slow-potassium current IKs on the voltage. The trajectory covers the range 100 to 57 mV at which point spontaneous oscillations occur. In both cases, the black points are at 5-mV intervals with the last point plotted before destabilization occurring at 60 mV, and the arrows represent the direction of depolarizing membrane potential.

RESPONSE TO OSCILLATING CURRENT. The first experimental measure of subthreshold properties considered here is the voltage response to an oscillating input current. The magnitude and phase of this response measure the impedance of the neuronal membrane. A subthreshold resonance, signaled by the existence of peak in |Z(f)| at some nonzero frequency f_R, occurs in the whole of the green region of Fig. 2A (Hutcheon et al. 1996a). The line that bounds the region of the phase diagram in which resonance occurs starts at the point  = 1,  = 1 (the intersection between the two instability lines), and for large , it tends toward the axis  = 0. Thus in most of the stable region with  > 0 resonant behavior occurs. Positivity of  implies that the associated activation or inactivation variable is a resonant variable (Hutcheon and Yarom 2000).

PHASE RESPONSE TO OSCILLATING CURRENT. Another quantity of interest is the existence of a zero phase-lag in the membrane potential response at nonzero frequency, seen in cortical neurons (Gutfreund et al. 1995). Analysis of the phase difference (f) between the oscillating current and the voltage response, defined in Eq. 11, shows that a zero phase-lag exists for  > 1. This quantity also coincides with the existence of a maximum in the phase and a phase advance from the driving current. A second line  < 0 signifies the existence of a minimum in the phase, implying that at some frequencies the phase lag is greater than 90°. Taken together, these criteria divide the phase diagram into four regions, plotted in Fig. 2B. It should be noted that none of the lines separating the different qualitative responses of the phase correspond exactly to the presence of a resonance in the amplitude of the impedance.

RESPONSE TO A SQUARE-PULSE CURRENT. The application of a small step change in applied current provides a different assessment of subthreshold membrane properties. The response of the neuron to such a current can be obtained explicitly (see APPENDIX for details). At the level of the two-variable description the neuron can exhibit three different types of response to the square-pulse current shown, the voltage-time profiles of which are shown in the insets of Fig. 2C:

Subthreshold behavior of the two conductance-based models

The diagram introduced in Fig. 2 allows a visualization of the trajectory of the neuron through the space of parameters (V*) and (V*) as the holding voltage V* is varied. This is shown for the two model neurons defined in METHODS and detailed in the APPENDIX. As the trajectory crosses different boundaries, so the neuronal response to input current will change qualitatively. The trajectories of the two model neurons, as parameterized by the changing effective leak g and coupling variable g₁ = g_f or g₁ = g_q, are plotted in Fig. 2D.

MODEL I. A NEURON WITH I_NA, I_K, AND I_H CURRENTS. The resonance curve of the neuron is plotted in Fig. 3A for a holding potential of 65 mV. The spike-generating currents are much faster than other time scales in the system and can be taken as instantaneous, reducing the full model to a three-variable description: the membrane potential and the two activation variables of the H current (see APPENDIX for details). The impedance curve of the reduced three-variable description is also shown in Fig. 3A. Comparison with the full model shows that this reduction is extremely accurate. A further approximation, that the slow variable of the I_H current averages to a steady value, provides the two-variable description. The behavior of the neuron is therefore classified by its leak conductance g and the effect of the I_H fast variable f, through the two dimensionless parameters  = g_f/C and  = g_ff/C. This two-variable description provides an excellent approximation of the original model for driving frequencies greater than 2 Hz as shown in Fig. 3A. At frequencies greater than 2 Hz, the dynamics of the slow variable of the H current is too slow to follow the voltage changes and therefore to affect the resonance curve.

View larger version (21K): [in this window] [in a new window]   Fig. 3. The frequency-dependent input impedance for the two model neurons held at 65 mV, showing the level of approximation between the full ( and ), three-variable () and two-variable (- - -) descriptions. A: model I. The three-variable model is obtained by taking the spike-generating currents to be instantaneous. The full and three-variable models agree closely for all frequencies plotted. The two-variable approximation is obtained by noting that the slow variable of the IH current averages to a steady value for frequencies greater than ~2 Hz. B: model II. The full noninactivating IKs model () and its two-variable approximation, obtained by taking the spike-generating currents to be fast. If an inactivation variable is included in the definition of the IKs current, an impedance profile () with a trough at 3 Hz as well as a resonance at 30 Hz is seen (see Three-variable model and the APPENDIX for details). The corresponding three-variable model obtained by taking the spike-generating currents to be fast, but retaining the activation and inactivation variables of the IKs current, provides a good approximation of the full model.

View larger version (16K): [in this window] [in a new window]   Fig. 4. The resonance frequency fR () and frequency of damped-oscillations f0 (- - -) of the 2 model neurons as a function of holding voltage. The corresponding Q values are also given that, following convention, measure the relative strength of the resonant peak Q = |Z(fR)|/|Z(0)|. A: model I. A resonance exists for most of the subthreshold regime, whereas damped oscillations occur only in a narrow voltage range near the firing threshold (at 56.4 mV). The resonance is strongest at 80 mV with fR = 10 Hz. B: model II. In this case, both resonance and damped oscillations exist above 72.5mV. The resonance strength increases as the membrane potential approaches the onset of spontaneous oscillations (indicated by  at 57.2 mV), above which the neuron fires periodically. The regions in which damped oscillations and resonance exist shown in these graphs can be compared with the trajectories in Fig. 2D.

MODEL II. A NEURON WITH I_NA, I_K, I_KS, AND I_NAP CURRENTS. In a similar way to the previous case, the full conductance-based model can be reduced to a two-variable description by noting that the spike-generating currents are fast. This approximation is very accurate, as can be seen in Fig. 3B (the noninactivating I_Ks and two-variable profiles). The behavior of this neuron is classified by its leak g, the I_Ks coupling variable g_q, and its time constant _q through the dimensionless parameters  = g_q/C and  = g_qq/C. These quantities are calculated for a subthreshold voltage range of 100 to 57 mV, and the corresponding trajectory is plotted in Fig. 2D. In distinction to model I, the resonance here is driven by depolarization-activated currents. This can be seen in the vertical-moving trajectory and the increasing Q value as the line of onset of spontaneous oscillations is approached. This neuron features both a resonant current in the activation of I_Ks as well as an amplifying current I_NaP. The amplification effect is clearly seen in a comparison of the Q values of the resonance of model II with model I (which lacks an amplifying mechanism). An examination of Fig. 4B shows that the resonant frequency steadily increases as the line of onset of spontaneous oscillations is approached, a feature reminiscent of the "broad-frequency cells" in Llinas et al. (1991) that were also shown to feature a TTX-sensitive persistent sodium current as well as a delayed rectifier. The amplification is due to the fact that the I_NaP current decreases the effective leak (because the current amplifies voltage changes). Hence, the parameter  decreases as an effect of the activation of this current, and correspondingly the model moves toward the left in the diagram of Fig. 2D, in the direction of the line where spontaneous oscillations occur. On this line, the Q value diverges.

Subthreshold behavior of the three-variable model

Although two-variable models capture the properties of a broad class of neurons, Fig. 2A shows that they only provide two classes of impedance curves, either monotonously decreasing with increasing frequency (nonresonant case shown in white) or with a single peak at a preferred frequency (resonance case shown in green). Therefore more complex impedance curves cannot be described with only two variables. For this reason, the same classification described in the preceding text was performed for a model with three variables defined in Eq. 7. The model is now specified by four parameters: the conductance ratios g₁/g and g₂/g and the time constants ₁ and ₂ of the variables w₁ and w₂. Without loss of generality, the time constant of w₂ was chosen to be the faster variable, thus ₂ < ₁. When the faster variable is taken to be instantaneous, ₂ = 0, the model becomes equivalent to the two-variable model with  = (g + g₂)₁/C and  = g₁₁/C.

The classification of the three-variable model was done in terms of the presence or absence of damped oscillatory behavior in response to transient inputs and presence or absence of resonant behavior as defined by the peaks of the impedance profile. Figure 5 shows the different regions of interest in the (g₁/g) versus (g₂/g) plane for several values of the times constants ₁ and ₂ relative to the time constant  = C/g.

View larger version (27K): [in this window] [in a new window]   Fig. 5. Phase diagram of the three-variable model in the plane g1/g, g2/g, for different values of the time constants 1 and 2 relative to  = C/g (marked above each panel). Brown, unstable region; green, regions where resonance occur; dark green, resonance with a trough at a lower frequency; red lines, the boundaries of the regions in which damped oscillations occur in response to a current step. See text for more details.

STABILITY. As in the two-variable case, the equilibrium voltage may destabilize in two different ways: either the total conductance g₁ + g₂ + g becomes negative (Fig. 5, ), or spontaneous oscillations appear (Fig. 5, - - -). These two lines intersect at the point where the frequency of the spontaneous oscillations becomes zero. When the time constants are such that _{2 1}, the line where spontaneous oscillations arise becomes vertical. Thus as stated in the preceding text, in the limit ₂ = 0, the behavior of the two-variable model is recovered (compare Fig. 5B with Fig. 2).

RESONANCE AND DAMPED OSCILLATIONS. Again, similar but distinct regions are observed in which damped oscillations and resonant behavior occur. In the upper-left quadrant of each panel (where the slower variable is "resonant," g₁ > 0, while the faster variable is "amplifying," g₂ < 0), both damped oscillations and resonant behavior are present in most of the stable region. In the upper-right quadrant where both variables are "resonant," damped oscillations and resonant behavior are found almost throughout. In the lower-left quadrant, both variables are amplifying: the neuron exhibits neither damped oscillations nor resonance and destabilizes exclusively by the total conductance becoming zero.

APPEARANCE OF A TROUGH. In the bottom-right quadrant (where the slower variable is now amplifying g₁ < 0 and the faster variable is resonant g₂ > 0), a qualitatively new phenomenon is observed, depending on the values of the time constants. When the two variables ₁ and ₂ are slower than the time constant  = C/g, the amplitude of the impedance has a local minimum or trough at a finite frequency (indicating a suppression of the membrane response at that frequency) followed by a resonant peak at higher frequency. One example of such subthreshold dynamics is a neuron with an inactivating potassium current with a relatively large activation time constant ₂ and a much larger inactivation time constant ₁ with overlapping steady-state activation and inactivation functions (a window current). In fact, a model of the I_Ks current that also includes an inactivation variable (see APPENDIX for details) gives exactly this effect. The frequency-impedance curve for such a neuron is plotted in Fig. 3B and shows a close similarity to the experimentally measured impedance curve of the fast-spiking interneurons measured in Pike et al. (2000). A two- and three-variable reduction of this full conductance-based model are also plotted for comparison. Another possibility would be to have the two variables implemented in two active persistent currents: a potassium current with activation time constant ₂ and a sodium or a calcium current with slower activation time constant ₁.

Firing-rate resonance

In this section, the effect of the subthreshold resonance on the dynamics of the firing rate is investigated. The aim of the analysis is to determine when a small oscillatory component in the synaptic inputs of a given neuron will be amplified in its output and how this depends on the subthreshold properties of the considered neuron.

The current used to model the synaptic bombardment such a neuron would experience in vivo comprises a noisy hyperpolarizing or depolarizing drive as well as a weak sinusoidal component of frequency f. The signal gain A(f) defined in Eq. 14 and illustrated in Fig. 1 measures the strength of the temporal modulation of the instantaneous firing rate induced by the oscillating current. It is the firing-rate analog of the impedance Z(f), and it is the existence of a peak in the signal gain that categorizes the amplification of frequencies in the outgoing spike train of resonant neurons. As will be shown, the noise inherent in biological networks is an important factor in determining the frequency that is maximally amplified.

The range of behavior is first examined by the use of a GIF model neuron. These results are then illustrated by two conductance-based models with spike generating currents and also an I_H current or I_NaP and I_Ks currents.

Firing-rate resonance in the GIF model neuron

In the previous section, the subthreshold behavior of a general two-variable model was analyzed in detail. As described in METHODS, a simple spike mechanism (threshold and reset) can be added to the two-variable description to produce a generalization of the IF neuron. This provides the simplest mathematical description of a spiking neuron with resonant subthreshold dynamics and allows a direct link to be made between the subthreshold characteristics and the statistical properties of the outgoing spike train. In spite of its simplicity, the GIF model provides a good approximation to more complete descriptions of neurons as will be seen in the next section.

The model examined here is parameterized by C = 0.5 nF, g = 0.025 µS, g₁ = 0.025 µS, and ₁ = 100 ms, giving a subthreshold resonance frequency f_R near 5 Hz. The signal gain A(f) was examined as a function of frequency for different respective strengths of the constant I₀ and noisy I_N components of the injected current, defined in Eq. 9. It should be noted that the sinusoidal component I₁ is always taken to be weak in the present work.

The fact that the neuron is induced to fire at a frequency r₀ by the applied current implies that there are now two distinct and independent frequency scales: the subthreshold resonant frequency f_R, controlled by the subthreshold dynamics of the membrane potential, and the background firing frequency r₀, which is controlled by the characteristics of the externally applied noisy current. It is useful to distinguish situations depending on whether r₀ is lower or greater than f_R.

View larger version (39K): [in this window] [in a new window]   Fig. 6. Effect of noise strength on the amplification of a frequency f in the firing rate of a resonant neuron. Two cases are considered: low noise (red) and high noise (blue). The particular generalized integrate-and-fire (GIF) model neuron (C = 0.5 nF, g = 0.025 µS, g1 = 0.025 µS, and 1 = 100 ms) has a resonance near 5 Hz, and in all cases, the firing rate is kept at r0 = 20 Hz. A and B: simulation of a 3-s current injection protocol. 0-1 s: injected noisy current with no oscillatory component. 1-2 s: addition of a 5 Hz (fR) sinusoidal component. 2-3 s: sinusoidal component at 20 Hz (r0). A: low-noise case (I0 = 0.95 nA, I1 = 0.024 nA, and IN = 0.11 nA) in which the neuron fires regularly. The amplification is greatest at f = r0 = 20 Hz. B: high-noise case (I0 = 0.78 nA, I1 = 0.059 nA, and IN = 0.55 nA). The strongest amplification is at the resonance frequency f = fR = 5 Hz. C: signal gain amplitude |A(f)| vs. frequency. D: phase of the signal gain (f) vs. frequency. Both low (red curves) and high (blue cu