In central neurons, the threshold for spike initiation can depend on the stimulus and varies between cells and between recording sites in a given cell, but it is unclear what mechanisms underlie this variability. Properties of ionic channels are likely to play a role in threshold modulation. We examined in models the influence of Na channel activation, inactivation, slow voltage-gated channels and synaptic conductances on spike threshold. We propose a threshold equation which quantifies the contribution of all these mechanisms. It provides an instantaneous time-varying value of the threshold, which applies to neurons with fluctuating inputs. We deduce a differential equation for the threshold, similar to the equations of gating variables in the Hodgkin-Huxley formalism, which describes how the spike threshold varies with the membrane potential, depending on channel properties. We find that spike threshold depends logarithmically on Na channel density, and that Na channel inactivation and K channels can dynamically modulate it in an adaptive way: the threshold increases with membrane potential and after every action potential. Our equation was validated with simulations of a previously published multicompartemental model of spike initiation. Finally, we observed that threshold variability in models depends crucially on the shape of the Na activation function near spike initiation (about −55 mV), while its parameters are adjusted near half-activation voltage (about −30 mV), which might explain why many models exhibit little threshold variability, contrary to experimental observations. We conclude that ionic channels can account for large variations in spike threshold.
Spike initiation in neurons follows the all-or-none principle: a stereotypical action potential is produced and propagated when the neuron is sufficiently excited, while no spike is initiated below that threshold. The value of that threshold sets the firing rate and determines the way neurons compute, for example their coincidence detection properties [1], [2]. It is generally described as a voltage threshold: spikes are initiated when the neuron is depolarized above a critical value, when voltage-dependent sodium channels start to open. That biophysical mechanism is well understood since the studies of Hodgkin and Huxley in the squid giant axon [3] and subsequent modeling studies [4]–[7]. Recent findings have renewed the interest in the spike threshold. First, there is an intense ongoing debate about the origin of threshold variability observed in vivo [8]–[14]. In particular, it is unclear whether threshold variability is mainly due to experimental artifacts or molecular mechanisms, which might question the relevance of the Hodgkin-Huxley model for central neurons. Moreover, numerous experiments have shown that spike initiation does not only depend on the membrane potential but also on complex features of the inputs. For example, it depends on the preceding rate of depolarization [15]–[21] and on the preceding interspike intervals [12], [22]. Those properties are functionally important because they enhance the selectivity of neurons in several sensory modalities, in particular in audition [23], vision [24], and touch [21]. Developmental and learning studies have also shown that the threshold adapts to slow changes in input characteristics. This phenomenon is known as long-term plasticity of intrinsic excitability and may be involved in the regulation of cell firing, short term memory and learning [25]–[31]. The excitability threshold also varies with the distance to the soma in a given neuron and with cell type [2], [15], [32]–[35], which may explain functional differences. The modulation of cell excitability might be explained by the activation of voltage-gated potassium channel Kv1 [36]–[41], inactivation of voltage-gated sodium channels [15], [16], [19], [21], fluctuations in sodium channel gating [42], inhibitory synaptic conductance [43]–[45] and the site of spike initiation [14]. To understand the origin of spike threshold variability, we examined the role of several candidate mechanisms in biophysical neuron models: activation and inactivation of the sodium channel, slow voltage-gated channels (e. g. Kv1), synaptic conductances and the site of spike initiation. Our analysis is based on a simplification of the membrane equation near spike initiation and results in a simple formula for the spike threshold that quantifies the contribution of all those mechanisms. The threshold formula provides an instantaneous time-varying value which was found to agree well with numerical simulations of Hodgkin-Huxley type models driven by fluctuating inputs mimicking synaptic activity in vivo, and with simulations of a realistic multicompartmental model of spike initiation [54].
To derive the threshold equation, we made a quasi-static approximation, assuming that all mechanisms that modulate the threshold are slow processes (compared to the timescale of spike initiation). That threshold equation provides an instantaneous value of the spike threshold, as a function of modulating variables. Here we show how the dynamics of sodium inactivation, voltage-gated conductances and synaptic conductances translate into spike threshold dynamics.
Since Na channels are responsible for the generation of action potentials, the threshold is firstly determined by their activation characteristics. Activation curves for Na channels are well approximated by Boltzmann functions with similar slope factors (ka = 4–8 mV in neuronal channels). The threshold is linearly related to the half-activation value Va and logarithmically related to the maximum Na conductance ga. The threshold also depends logarithmically on the Na inactivation variable h, so that it increases with the membrane potential and with every emitted spike. The modulating effect of inactivation is most pronounced when the half-activation value Vi is lowest (i. e. , Na channels are partially inactivated at rest). Finally, the threshold depends logarithmically on the total conductance, which includes the leak conductance, voltage-gated conductances and synaptic conductances. In particular, Kv1 channels, which are expressed with high density at the spike initiation site [37], [39], [41], increase the threshold in an adaptive manner (the threshold increases with the membrane potential). This change in threshold occurs simultaneously with the effective membrane time constant, whereas threshold changes due to Na inactivation have no effect on the time constant, which might suggest a way to experimentally distinguish between the two effects. Indeed, the effective membrane time constant (as measured in vivo for example in Léger et al. , 2003 [73]) is (C is the membrane capacitance) while the effect of total conductance on spike threshold varies as, therefore as −. It is currently unclear whether threshold modulation is mainly due to Na inactivation or delayed-rectifier K currents. Our simulations with the multicompartmental model of spike initiation in pyramidal cells [54] suggest that the spike threshold is essentially determined by Na inactivation, but this may not be universally true. Recent experimental findings in hippocampal mossy fibers [74] suggest that delayed K+ currents are closed at spike initiation, which minimizes charge movements across the membrane and is thus more metabolically efficient. It emphasizes the fact that Na inactivation is a more metabolically efficient way to modulate spike threshold than K+ activation, since the former reduces charge transfer while the latter increases it. We have not considered the effect of channel noise, i. e. , fluctuations in Na channel gating [42], [75]–[78], which result in random threshold variations. Although dynamical equations of fluctuations in Na channel gating are well set [79], [80], they cannot be included in our theoretical framework because we neglected the time constant of Na activation (which leads to the exponential model). There are two additional sources of variability which are artefactual: the fact that the threshold is not measured at the site of spike initiation, and threshold measurement methods. The latter source is difficult to avoid in vivo because only spike onsets can be measured. The former one also seems technically very difficult to avoid in vivo, since spikes are initiated in the axon hillock, which is only a few microns large. Although the soma and AIS are virtually isopotential below threshold, experimentally measured values of threshold differ between the two sites [34] because, as we previously remarked, in vivo measurements correspond to spike onset rather than threshold and therefore take place after spike initiation, when the two sites are not isopotential anymore. This experimental difficulty may introduce artefactual variability in threshold measurements [14]. To derive the threshold equation, we made several simplifying assumptions. First, we assumed that Na activation is instantaneous. It is indeed significantly faster than all other time constants but not instantaneous. The approximation is legitimate as long as the effective membrane time constant in the membrane equation is small (, including all conductances), which is generally true before threshold. When Na channels open, the Na conductance dominates the total conductance and drastically reduces the effective time constant. Thus, we expect this approximation to be reasonable to predict spike initiation properties but not spike shape characteristics. Our second major assumption is a quasistatic approximation, i. e. , we assume that near spike initiation, all modulating variables and the input current can be considered as constant. In other words, we assume that the time constants (except that of Na activation) are larger than a few ms. This is clearly only a mathematically convenient approximation, but our predictions empirically agreed with numerical simulations. To investigate the role of Na inactivation, we also assumed that activation and inactivation are independent, which is a standard simplifying hypothesis (Hille, 2001). Although it is debatable [49], [56], it should remain valid in the case where activation and inactivation time constants are well separated. We also assumed that Na activation and inactivation curves were Boltzmann functions. Experimental data is indeed well fitted by Boltzmann functions, but the reported parameter values (Va, ka) correspond to fits on the entire voltage range, whereas we are interested in hyperpolarized voltage regions where the activation values are small. When only the relevant part of the experimental data is considered, different parameter values might be obtained. For example, when analyzing previously published biophysical models, we found that better results were obtained when Na activation curves, which were not exactly Boltzmann functions, were fitted in the spike initiation region (−60 to −40 mV) rather than on the entire voltage range (Fig. 9). We examined this issue in the biophysical model used in this paper (see Materials and Methods). The Na activation curve of this model seemed to be well fit to a Boltzmann function (Fig. 9A), however the fit was poor in the spike initiation zone (−60 to −40 mV, Fig. 9C) where activation is close to zero, which makes fit errors relatively larger. Although the slope factor is about 6 mV when the activation curve is fit over the entire voltage range, similar to experimental measurements [51], it is only half this value when fit in the spike initiation region (Fig. 9D), which explains why this model, as many other biophysical models, exhibits little threshold variability (since threshold modulation is proportional to). We calculated the slope factor as a function of the voltage region and we found that it varies between 1 and 6 mV (Fig. 9D). This finding motivates a reexamination of Na channel voltage-clamp data, focusing on the spike initiation region rather than on more depolarized regions, which are more relevant for spike shape than spike initiation. Fig. 10 addresses two potential difficulties. In experiments, activation curves are obtained by measuring the peak conductance after the clamp voltage is changed from an initial value V0 to a target value V, and normalizing over the entire range of target voltages. Thus, it assumes that inactivation is still at resting level h (V0) when the peak current is measured. This would not be the case if the inactivation time constant τh were close to the activation time constant τm. Fig. 10A shows the effect of this overlap on the measurement of ka with simulated voltage-clamp data, where is a Boltzmann function with ka = 6 mV. It appears that ka is overestimated if τh is very close to τm, up to 50% when the two time constants are equal (to 0. 3 ms in these simulations). However the error quickly decreases as τh increases (e. g. 10% error for τh = 1 ms). Another potential difficulty is the lack of data points in the relevant voltage range and the measurement noise, because currents are small. In Fig. 10B, we digitized an experimentally measured activation curve (black dots), where clamp voltages were spaced by 5 mV. A Boltzmann fit over the entire voltage range gave ka = 7. 2 mV while a fit over the hyperpolarized region V<−40 mV gave ka = 4 mV. However, the latter is not a reliable estimate because it corresponds to only 4 non-zero data points, which also seem to be corrupted by noise. Therefore it might be necessary to perform new measurements, specifically focusing on the spike initiation zone, perhaps with multiple measurements to reduce the measurement noise. Alternatively, ka could be measured with a phenomenological approach, using white noise injection in current clamp [22]. Another possible approach would be to directly fit the excitability model to the current-clamp response of a cell in which only Na channels would be expressed (perhaps with fluctuating inputs). Finally, our analysis relies on a single compartment model. In the compartmental model, we found that between spikes, the membrane potential was 1. 8±0. 6 mV more depolarized at the soma than at the AIS. This is small compared to the slopes of all activation and inactivation curves in this model (5–9 mV). This agrees with our analysis of the electrotonic length in the subthreshold range, which is much larger than the distance between the soma and the AIS, although very fast synaptic inputs or proximal axonal inhibition could produce larger voltage gradients. Thus, our analysis should remain valid if the compartment represents both the soma and initiation site (and also proximal dendrites). However, that approximation is not valid anymore after spike initiation (see below). Spikes look sharper in the soma than in the AIS, presumably because they are initiated in the AIS and back-propagated to the soma [10], [13], [14]. That property is also seen in numerical simulations of multicompartmental models [34], [70]. Yet, linear cable theory predicts the opposite property: the voltage at the soma is a low-pass filtered version of the voltage at the AIS, therefore spikes should look less sharp in the soma. Thus, increased sharpness must be due to active backpropagation of the action potential, which cannot be seen in a two compartment model (such as described in Text S1). From a theoretical point of view, the sharpening effect of backpropagation can be intuitively understood from the cable equation: It appears that the membrane equation is augmented by a diffusion term, which is positive and large in the rising phase of the action potential between the initiation site and the soma. Thus, for the same membrane potential V, the time derivative gets larger as this diffusion term increases, which sharpens action potentials. Sharpness can be measured in numerical simulations by plotting dV/dt vs. V in response to a suprathreshold DC current, and fitting it to an exponential model (). In the model of Hu et al. [54], we found that the slope factor, characterizing spike sharpness, was = 1. 6 mV in the AIS and only 0. 8 mV in the soma. This is in approximate agreement with empirical fits of exponential integrate-and-fire models to cortical neurons stimulated with fluctuating inputs, which reveal a surprisingly small slope factor, slightly above 1 mV [22]. Thus, in the multicompartmental model, active backpropagation did increase spike sharpness in the soma, but also in the AIS, since the slope factor was about twice smaller than predicted from fitting the Na activation curve to a Boltzmann function (3. 6 mV). This increased sharpness did not affect the magnitude of threshold modulation. In single-compartment models, sharpness of spikes and threshold modulation are determined by the same quantity, related to the sharpness of the Na activation curve (ka). It appears that this link does not hold anymore when active backpropagation is considered (in multicompartmental models). Thus, in the threshold equation, the modulating factor is indeed ka (from the Na activation curve) rather than (from spike sharpness, measured in the phase plot (dV/dt, V) ). This explains that Na inactivation can produce large threshold variability (10 mV in our simulations) even though spikes are very sharp.
We consider a single-compartment neuron model with voltage-gated sodium channels and other ion channels (voltage-gated or synaptic), driven by a current I. The membrane potential V is governed by the membrane equation: where C is the membrane capacitance, gL (resp. EL) is the leak conductance (resp. reversal potential), gi (resp. Ei) is the conductance (resp. reversal potential) of channel i, gNa (resp. ENa) the maximum conductance (resp. reversal potential) of sodium channels, PNa is the proportion of open Na channels and I is the input current. In this article, we used the following convention for conductances: lower case (g) for the total conductance over the surface of a compartment (typically in units of nS) and upper case (G) for conductances per unit area (in units of S/cm2). We assume that sodium channel activation and inactivation are independent, as in the Hodgkin-Huxley model [3], i. e. , , where Pa is the probability that activation gates are open and Pi is the probability that a channel is inactivated. Following the Hodgkin-Huxley formalism, we define. The steady-state activation curve can be empirically described as a Boltzmann function [51]: where is the half-activation voltage () and the activation slope factor (). We make the approximation that Na activation is instantaneous and we replace Pa by its equilibrium value, so that. With instantaneous activation, the sodium current is: Action potentials are initiated well below Va (about −30 mV, Angelino and Brenner, 2007 [51]), so that except during the spike. Similarly, ENa is very high (about 55 mV), so that is not very variable below threshold. We make the approximation and we obtain: where. This approximation is meaningful for spike initiation but not for spike shape. With a reset (ignoring inactivation and other ionic channels), we obtain the exponential integrate-and-fire model [48], which predicts the response of cortical neurons to somatic injection with good accuracy, in terms of spike timings [22], [81], [82]. In this model, VT is the voltage threshold for constant input currents I and ka (originally denoted ΔT) is the slope factor, which measures the sharpness of spikes: in the limit mV, the model becomes a standard integrate-and-fire model with threshold VT (although this is different in multicompartmental models, see Discussion). The resulting approximated model is thus: It is convenient to sum all conductances (except for the Na channel), which gives a simpler expression: where is the total conductance and E* is the effective reversal potential: Finally, the inactivation variable h can be inserted in the exponential function: whereis the voltage threshold if all other variables are constant, i. e. , it is such that F′ () = 0, where F is the current-voltage function. The effect of Na inactivation on the threshold can be seen in the exponential model above, neglecting other conductances (thus). Assuming that inactivation is slow compared to spike initiation (quasi-static approximation), the voltage threshold is now, and it changes with the inactivation variable h. We assume, as in the Hodgkin-Huxley model, that inactivation has first-order kinetics: The steady-state value of the threshold is thus. We differentiate the threshold equation with respect to time: We now express h as a function of using the invert relationships: and: If the threshold remains close to its steady-state value (), the equation simplifies to: with. The same method applies for voltage-gated conductances (e. g. Kv1). We compared our theoretical predictions with numerical simulations of a previously published point-conductance model with fluctuating synaptic inputs [83]. The membrane equation is: where and n are respectively the maximal conductance and the activation variable of thedelayed-rectifier potassium current, and and p are respectively the maximal conductance and the activation variable of the non-inactivating K current. All channel variables have standard Hodgkin-Huxley type dynamics. In Fig. 3A–C, only Na channel activation was considered, with instantaneous dynamics, i. e. , , h = 1, n = 0, p = 0, I = 0: In Fig. 3D, the threshold equation was used to calculate VT for the Na channel properties reported in Angelino and Brenner (2007) [51], since only the values of Va and ka were available. To evaluate our threshold equation with time-varying inputs (Figs. 2C, 5 and 6), we simulated the full conductance-based model with fluctuating synaptic conductances (same parameters as in Destexhe et al. , 2001 [83]). In Fig. 6, we shifted the voltage dependence of Na inactivation toward hyperpolarized potentials by −12. 5 mV so as to obtain more threshold variability. To measure the time-varying threshold, we used a similar method as one previously used in vitro by Reyes and Fetz [84], [61]. We simulated the model for 200 ms and measured the instantaneous threshold at regular time intervals T as follows. The model was simulated repeatedly with the same synaptic inputs (frozen noise). In each trial, the neuron was depolarized at time nT (only once per 100 ms run) to a voltage value between −51 mV and −38 mV. With T = 0. 6 ms and 65 voltage values, we ran 22,000 trials. The threshold at a given time is defined as the minimal voltage value above which a spike is elicited. The measured threshold was compared to the prediction obtained with the threshold equation (see Results), where VT and ka were obtained from a Boltzmann fit to the activation function over the range −51 mV to −38 mV, giving VT = −68 mV and ka = 3. 7 mV (Va = −30. 4 mV). The values of VT and ka depended on the fitting window (see Discussion and Fig. 9). In Fig. 7, the voltage dependence of Na inactivation was shifted by −20 mV to induce more threshold variability (giving Vi = −62 mV instead of −42 mV with the original parameter values) and the maximum Na conductance was multiplied by 3 (to keep threshold values in the same range). The standard deviations of synaptic conductances were also increased. In Fig. 8, we simulated a multicompartmental model of spike initiation recently published by Hu et al. (2009) [54], with fluctuating injected current modeled as an Ornstein-Uhlenbeck process (mean 0. 7 nA, standard deviation 0. 2 nA, time constant 10 ms). The model was otherwise unchanged. The spike threshold, both at the soma and AIS, is defined at the voltage value when dV/dt first exceeds 10 V. s−1 preceding a spike. In some panels (Fig. 8C, D, H), we extracted spikes from voltage traces by removing parts between spike onsets and 7 ms later. We estimated the activation properties of the Nav1. 6 channel, which is responsible for spike initiation in this model, by fitting a Boltzmann function to the activation curve () in the spike initiation zone (−60 mV to −40 mV), which gave Va = −33 mV and ka = 3. 6 mV. We then calculated the total maximal conductance of the Nav1. 6 channel over the AIS, by integrating the channel density over the surface of the AIS (using the morphology and channel density implemented in the published model code). We found gNa = 236 nS. Calculating the total leak conductance in this way was more difficult because leak channels were uniformly distributed on the whole morphology, including the dendrites, so that spatial attenuation should be taken into account. While this is theoretically possible using linear cable theory, we chose a simpler approach by directly measuring the membrane resistance at the soma with a DC current injection, and we found gL = 38 nS. With these values, the threshold equation predicted that the base threshold is VT = −55. 9 mV. The model had a slow K+ current (Im) with the same channel density as the leak channels. Therefore the maximum total conductance was estimated as gKm = gL = 38 nS. It also had a fast K+ current which was distributed inhomogeneously on the whole neuron morphology, including dendrites. We estimated its total maximum conductance as, where the effective dendritic area was estimated from the ratio of total leak conductance over leak channel density, i. e. , . We found gKv = 906 nS. We then calculated the theoretical threshold using these parameters and the instantaneous values of the relevant channel variables (h, nKm, nKv4). In Fig. 10A, we simulated a voltage clamp experiment in a simplified model with only Na channels, assuming the leak current was subtracted, where both activation and inactivation curves (and) were Boltzmann functions, with parameters Va = −30 mV, ka = 6 mV, Vi = −65 mV and ki = 6 mV. The activation and inactivation time constant were fixed (τm = 0. 3 ms and τh between 0. 3 and 3 ms). The conductance was measured at the current peak after the clamp voltage was switched from a fixed initial voltage V0 = −70 mV to a test voltage V, which was varied between −100 and 50 mV (the current was divided by V-ENa to obtain the conductance, and we assumed that ENa was known - in an experiment it would be obtained from a linear fit to the highest voltage region). The conductance was normalized by the maximal conductance over the tested voltage range and the resulting curve was fit to a Boltzmann function in the hyperpolarized region V<−40 mV. All simulations were written with the Brian simulator [85] on a standard desktop PC, except the simulation of the multicompartmental model of Hu et al. [54], for which we used Neuron.
Neurons communicate primarily with stereotypical electrical impulses, action potentials, which are fired when a threshold level of excitation is reached. This threshold varies between cells and over time as a function of previous stimulations, which has major functional implications on the integrative properties of neurons. Ionic channels are thought to play a central role in this modulation but the precise relationship between their properties and the threshold is unclear. We examined this relationship in biophysical models and derived a formula which quantifies the contribution of various mechanisms. The originality of our approach is that it provides an instantaneous time-varying value for the threshold, which applies to the highly fluctuating regimes characterizing neurons in vivo. In particular, two known ionic mechanisms were found to make the threshold adapt to the membrane potential, thus providing the cell with a form of gain control.
neuroscience/theoretical neuroscience
computational biology/computational neuroscience
A Threshold Equation for Action Potential Initiation