Abstract
We examined the kinetic properties of voltage-gated Na+ channels and their contribution to the repetitive spiking activity of medullary raphé neurons, which exhibit slow pacemaking and strong spiking adaptation. The study is based on a combination of whole-cell patch-clamp, modeling and real-time computation. Na+ currents were recorded from neurons in brain slices obtained from male and female neonatal rats, using voltage-clamp protocols designed to reduce space-clamp artifacts and to emphasize functionally relevant kinetic features. A detailed kinetic model was formulated to explain the broad range of transient and stationary voltage-dependent properties exhibited by Na+ currents. The model was tested by injecting via dynamic clamp a model-based current as a substitute for the native TTX-sensitive Na+ currents, which were pharmacologically blocked. The model-based current reproduced well the native spike shape and spiking frequency. The dynamics of Na+ channels during repetitive spiking were indirectly examined through this model. By comparing the spiking activities generated with different kinetic models in dynamic-clamp experiments, we determined that state-dependent slow inactivation contributes significantly to spiking adaptation. Through real-time manipulation of the model-based current, we established that suprathreshold Na+ current mainly controls spike shape, whereas subthreshold Na+ current modulates spiking frequency and contributes to the pacemaking mechanism. Since the model-based current was injected in the soma, the results also suggest that somatic Na+ channels are sufficient to establish the essential spiking properties of raphé neurons in vitro.
Introduction
Neurons can spike with only two voltage-gated conductances: sodium (Nav) channels and potassium (Kv) channels, as demonstrated by Hodgkin and Huxley (1952). Although many other channels cooperate to implement a multitude of firing patterns (Bean, 2007), Nav channels play a central role. Apart from generating spikes, Nav channels are involved in pacemaking (Bennett et al., 2000; Taddese and Bean, 2002), fast spiking (Khaliq et al., 2003), spike adaptation (Fleidervish et al., 1996; Blair and Bean, 2003), spike timing (Vervaeke et al., 2006), and creating complex firing patterns (Williams and Stuart, 1999; Del Negro et al., 2002; Swensen and Bean, 2003). Not surprisingly, defects in neuronal Nav channels have been associated with severe neurological disorders, such as epilepsy and seizures, ataxia, or increased sensitivity to pain (Catterall et al., 2008).
A powerful tool for studying the physiological and pathological function of Nav channels in individual neurons is dynamic clamp (Sharp et al., 1993). The principle is to pharmacologically block Na+ currents, and then functionally replace them with an injected current, dynamically calculated on the basis of a Nav channel kinetic model (Milescu et al., 2008). To a first approximation, the neuron makes no distinction between the native Na+ currents and the model-based current, as long as the model is accurate. Thus, the sensitivity of the firing pattern to Nav channel properties and the contribution of Na+ currents to spiking can be easily studied, albeit indirectly, by varying the properties of the model and manipulating the model-based current in real-time. This hybrid experimental-computational approach, where only Nav channels must be modeled against a background of uncharacterized conductances, is potentially more accurate than mathematical simulations, where in principle every current must be well described, and more versatile than pharmacological approaches, which are limited by the available drugs.
In this study, we use a combination of whole-cell patch-clamp recording, modeling and real-time computation to examine the kinetic properties of Nav channels and their contribution to the repetitive spiking activity of medullary raphé neurons. These neurons exhibit regular and slow pacemaking, and strong spiking adaptation (Jacobs and Azmitia, 1992; Ptak et al., 2009), behaviors that potentially involve Nav channels. To test this hypothesis, we first search for a Nav channel model that can explain the variety of experimentally observed kinetic properties, including activity-dependent inactivation. Then, we validate the model by checking whether a model-based current can reproduce the native spiking pattern. Using this model as a substitute for the real channels, we examine the activity of Nav channels during repetitive spiking at different frequencies. Next, we test whether Nav channel slow inactivation contributes to spiking adaptation, by comparing different models. Finally, we examine the differential effects of subthreshold and suprathreshold Na+ currents on the spiking pattern, and determine the role of subthreshold Na+ current in pacemaking, via real-time manipulation of the injected current.
Materials and Methods
Animal procedures.
All animal procedures were approved by the NINDS Animal Care and Use Committee.
Slices.
In vitro medullary slice preparations containing midline raphé neurons were obtained from neonatal (postnatal days 0–4) Sprague Dawley rats (Koshiya and Smith, 1999). The medulla was dissected in artificial CSF (aCSF) containing the following (in mm): 124 NaCl, 25 NaHCO3, 3 KCl, 1.5 CaCl2, 1.0 MgSO4, 0.5 NaH2PO4 and 30 d-glucose, equilibrated with 95% O2 and 5% CO2 (pH 7.4 ± 0.05 at room temperature). Transverse slices (300–400 μm thick) containing nucleus raphé obscurus and hypoglossal nerve rootlets were transferred to the recording chamber and superfused with aCSF at room temperature.
Neuron identification.
For voltage-clamp (VC) experiments, raphé neurons were selected based on their location in the slice, adjacent to the midline. For current-clamp (CC) or dynamic-clamp (DC) experiments, the neurons were further selected based on their spiking pattern: regular and slow pacemaking (1–3 Hz) with broad action potentials (3–6 ms) (Ptak et al., 2009); intrinsic pacemaking behavior was verified by the presence of spiking activity in cell-attached mode. All experiments were done under infrared-differential interference contrast visualization.
Solutions.
For VC, pipettes were filled with a solution containing the following (in mm): 70 Cs-gluconate (prepared from CsOH and gluconic acid), 30 Na-gluconate, 10 TEA-Cl, 5 4-AP, 10 EGTA, 1 CaCl2, 10 HEPES, 4 Mg-ATP, 0.3 Na3-GTP, 10 Na2-phosphocreatine, pH 7.4 adjusted with CsOH (285 ± 5 mOsm/L). The elevated internal Na+ concentration (50 mm) reduced the size of Na+ currents, and improved voltage and space clamp. For CC or DC, pipettes were filled with (in mm): 125 K-gluconate, 4 NaCl, 10 EGTA, 1 CaCl2, 10 HEPES, 4 Mg-ATP, 0.3 Na3-GTP, 4 Na2-phosphocreatine, pH 7.4 adjusted with KOH (285 ± 5 mOsm/L).
Pharmacology.
For VC experiments, CdCl2 (200 μm) and 7-nitro-2,3-dioxo-1,4-dihydroquinoxaline-6-carbonitrile (CNQX; 20 μm) were added to the superfusing aCSF to block Ca2+ currents and inhibit synaptic transmission. CNQX (20 μm) was also used in some CC and DC experiments. To block Na+ currents, TTX (1 μm) was added to the superfusing aCSF. All reagents were purchased from Sigma-RBI.
Electrophysiology.
Electrodes (5–7 MΩ) were pulled from borosilicate glass and were coated with Sylgard to reduce capacitive transients. Pipette capacitance was compensated 100% in VC, and ≈70% in CC and DC. For DC and off-line analysis, the membrane capacitance (Cm) was approximated as the value used for compensation in VC (≈20 pF). Series resistance (Rs) was typically 9–15 MΩ. Cells with Rs > 15 MΩ were discarded. In VC experiments, Rs was compensated 80% (2 μs response time), and the compensation was readjusted before each voltage-clamp protocol. In CC and DC experiments, Rs was compensated 100% and periodically readjusted. A measured liquid junction potential of ≈10 mV for the K+-based and ≈8 mV for the Cs+-based solutions was corrected on-line. Somatic whole-cell recordings were obtained with a HEKA EPC10 Double patch-clamp amplifier. For VC, data were low-pass filtered at 40 kHz and digitally sampled at 100 kHz, using the amplifier's built-in digitizer, controlled by the Pulse software from HEKA. For CC and DC, the signal representing membrane voltage was digitally sampled at 50 kHz (open-bandwidth), using National Instruments data acquisition cards (PCI 6052E or PCI 6251) and NI-DAQmx 8.1 or newer driver, controlled by the QuB software (www.qub.buffalo.edu).
Voltage-clamp experiments.
Voltage-clamp protocols (Fig. 1) were constructed and applied with the Pulse program. In all protocols except that in Figure 1C, leak currents were subtracted using the P/n procedure (n = 2–10). Some protocols (Fig. 1A–C) were repeated under bath-applied TTX, to better isolate by subtraction the TTX-sensitive currents, particularly the smaller components. In all protocols, the intersweep interval was 6 s at −80 mV, which was necessary for complete recovery from inactivation of Na+ currents. Since most experiments were relatively long, recording stability and Na+ current stationarity during a protocol were carefully monitored, whenever possible, and unstable recordings were discarded.
The voltage dependence of activation and the time course of activation and fast inactivation were examined in 25 neurons (Fig. 1A). Stationarity was tested with a brief voltage pulse (1 ms at −15 mV), preceding each voltage step by 20 ms. Space clamp was improved with a brief subthreshold voltage prepulse (8 ms at −40 mV) (Milescu et al., 2010). Three representative TTX-subtracted datasets were selected for analysis. To construct the activation curve, the peak current at each voltage was baseline-subtracted and converted to conductance, assuming a linear I–V relationship and a reversal voltage of +33 mV, and then normalized to the conductance value at 0 mV (Fig. 2A, activation). The three datasets were individually fitted with a third-power Boltzmann function f = (1/(1 + exp((Vh − V) × s)))3, where f is the fraction of activation, Vh is the half-activation voltage and s is the slope factor. The average parameters are (mean ± SD): Vh = −34.2 ± 1.2 mV, s = 0.157 ± 0.021 mV−1. To obtain the time course of activation and fast inactivation, the three sets of recordings were averaged, and the result was normalized point-by-point to the largest negative current (Fig. 2B).
Deactivation was examined in 24 neurons (Fig. 1B). Sixteen of these were selected for analysis, consisting in measuring the time of tail current decay from 80% to 40% of the peak value, at each deactivation voltage (Fig. 2C). The 80–40% time of decay is approximately equal to the time constant but can be calculated more quickly during model fitting. The current evoked by the activation step (0.3 ms at −15 mV) was used as a measure of current stationarity. Space clamp was improved with a brief subthreshold voltage prepulse (8 ms at −40 mV) (Milescu et al., 2010).
The time course of slow inactivation was examined in 8 neurons. Satisfactory space clamp, signaled by the absence of uncontrolled current spikes, was obtained in only one of the examined neurons (Fig. 1C) but we were able to determine that the Na+ currents recorded from this neuron were representative. The TTX-subtracted current was converted to conductance and then normalized to the peak conductance at 0 mV (Fig. 2F).
Availability was examined in 22 neurons (Fig. 1D). Recordings were considered stable if the peak current evoked by the test pulse at −15 mV that follows the conditioning at −85 mV was within 10% of the peak current evoked by the conditioning step at −15 mV. Representative recordings (n = 9, 13 and 12 for the 50 ms, 300 ms and 2.5 s inactivation steps, respectively) were selected for analysis. The peak current at each voltage was baseline-subtracted and normalized to the current value at −85 mV. The data were individually fitted with a Boltzmann function (Fig. 2A, availability), with average parameters (mean ± SD): Vh = −45.0 ± 1.0 mV, s = −0.192 ± 0.018 mV−1 (50 ms); Vh = −48.6 ± 1.2 mV, s = −0.193 ± 0.011 mV−1 (300 ms); Vh = −50.2 ± 1.2 mV, s = −0.202 ± 0.018 mV−1 (2.5 s).
Subthreshold inactivation was examined in 14 neurons (Fig. 1E). For each tested voltage, the available fraction was calculated as the ratio between the baseline-subtracted test peak current at a given inactivation time and the baseline-subtracted test peak current at time 0 (Fig. 2D). To monitor current stationarity, the inactivation intervals were alternated in the protocol, and recordings with evidence of nonmonotonous decay or space-clamp artifacts were discarded. At each voltage, between 1 and 7 datasets were selected for analysis.
Recovery from inactivation was examined in 27 neurons (Fig. 1F). For each tested voltage, the available fraction was calculated for each recovery time, as the ratio between the baseline-subtracted test peak current and the baseline-subtracted conditioning peak current (Fig. 2E). At each voltage, between 6 and 18 datasets were selected for analysis. Consistency between development of inactivation and recovery from inactivation is confirmed by the convergence of the two types of curves to the same stationary values at −60 and −50 mV; those values also match the availability curves plotted in Figure 2A. The current evoked by the conditioning step (5 ms at −15 mV) was used as a measure of current stationarity.
Entry into slow inactivation was examined in 6 neurons (Fig. 1G). Current stationarity was determined as for subthreshold inactivation, and unstable recordings were discarded. One representative dataset was selected for analysis. The fraction that remains available after a given inactivation time was calculated as the ratio between the baseline-subtracted peak test current and the baseline-subtracted peak test current for zero inactivation time (Fig. 2G).
The development of inactivation within the interspike interval was examined in 26 neurons (Fig. 1H). The available fraction was calculated as the ratio between the baseline-subtracted peak test current and the baseline-subtracted peak conditioning current (Fig. 2H). The activity-dependent inactivation was examined in 6 neurons (Fig. 1I). The available fraction versus pulse index was calculated separately for each interpulse duration and voltage, as the ratio between the baseline-subtracted peak currents evoked by the ith pulse and the first pulse (Fig. 2I).
Current-clamp and dynamic-clamp experiments.
Approximately 100 neurons were examined in total, and all were patched on the soma. At least 5 neurons were tested with each CC or DC protocol. All CC and DC experiments were performed using the dynamic-clamp functionality of the QuB software (Milescu et al., 2008), which runs in Windows and uses National Instruments data acquisition cards. QuB features a scripting language that allows the design of custom protocols (VC, CC or DC) that combine real-time model computation with data processing, visualization, and recording. We used computers featuring a dual quad-core 3.0 GHz Intel Xeon processor, or a dual 2.6 GHz AMD Opteron processor. The dynamic clamp was run at 50 kHz, with a measured time step jitter demonstrating good real-time performance (coefficient of variation 0.2–2%, depending on the processor).
The ion channel kinetic models tested in this study were state models, and were integrated using a fast and stable method (Milescu et al., 2008). Briefly, the vector Pt, containing the state occupancies of a channel model, is advanced in time using the equation Pt + dt = Pt × eQ(V) × dt, where Q(V) is the rate matrix (Colquhoun and Hawkes, 1995) at voltage V, t is the time, and dt is the integration time step (20 μs). The quantity eQ(V) × dt was precalculated over a voltage range (−80… +60 mV, 0.2 mV increment) and stored in a look-up table, and recalculated each time a kinetic parameter was changed.
The total current ITot [pA] injected in a DC or CC experiment had the expression ITot = IInj + IHold − INa. The quantity IInj is a bias current, used to polarize the cell in current step protocols. IHold is a voltage-clamping current, with expression IHold = Cm × GHold × (Vm − VHold), where Vm [mV] is the measured membrane voltage, Cm [pF] is the estimated membrane capacitance, VHold is the holding voltage, and GHold is a scaling factor (1 nS/pF). We used this mechanism to effectively clamp the intersweep membrane potential, by setting GHold = 1 during the intersweep interval (5 s at −70 mV), and GHold = 0 in the rest. Unlike passive injection of bias current, this feedback method was less subject to cellular heterogeneity and noise.
INa is the current predicted by the Nav channel model, with expression INa = f × GNa × Cm × POpen × (Vm − VNa), where GNa is the density of Nav channels [nS/pF], POpen is the instantaneous open probability of the Nav channel model, and VNa is the Na+ current reversal potential (≈62 mV). The f value is a conditional expression, equal to fLow if Vm < VThr, and equal to fHigh if Vm ≥ VThr, where VThr is a voltage threshold parameter [mV], and fLow and fHigh are two adimensional factors. The rationale of using this conditional expression is to be able to separately manipulate the amplitude of the subthreshold and suprathreshold components of INa, defined as the currents flowing below and above a given voltage threshold VThr, respectively.
To reduce data throughput, only Vm and INa were continuously recorded in dynamic-clamp experiments. All parameters (including IInj, GHold, VHold, VThr, fLow, fHigh, GNa and kinetic parameters) were recorded at the beginning of an experiment and whenever they were changed. Any other voltage- and time-dependent continuously varying quantity was calculated off-line, by simulating the model in response to the recorded voltage waveform. We verified that the holding intervals (5 s at −70 mV) preceding each sweep in a DC protocol were sufficiently long to allow the solution calculated off-line to converge to the on-line (real-time) solution. In some experiments requiring an instantaneous change of kinetics, two models were run at the same time but only one generated current (see Fig. 7A).
Data analysis and statistics.
We used custom scripts in QuB to process the raw VC data (Fig. 1), to extract and prepare a set of data for kinetic modeling (Fig. 2), and to process CC and DC data and extract information (e.g., spiking frequency vs time). Statistics and curve fitting (e.g., activation or inactivation curves) were done with the Prism 4.0 software. The quantities shown in Figure 4 (see below) were obtained by exponential fits to either experimental or simulated data.
Kinetic modeling.
Rate constants were modeled as exponential functions of voltage: kij = kij0 × ekij1 × V, where kij [ms−1] is the rate constant of the transition between states i and j, kij0 [ms−1] is the rate at zero depolarization, kij1 [mV−1] is the voltage sensitivity factor, and V [mV] is the voltage. Parameter constraints, including allosteric factors and microscopic reversibility, were enforced as previously described, using the singular value decomposition technique (Milescu et al., 2005). Kinetic parameters were optimized with QuB as previously described (Milescu et al., 2005), with some modifications to allow for direct fitting of fractional data (e.g., the curves describing the activation or the recovery from inactivation), in addition to macroscopic current time course, and to allow for optimization of allosteric factors.
The model fitting algorithm minimized a cost function that was calculated as the sum of square errors between the data and the simulated predictions of the model. These simulations were calculated in response to the same voltage protocols as were used to record the experimental data, including holding intervals and test pulses. The components of the cost function, each corresponding to one of the quantities plotted in Figure 2, were normalized to 1 or −1, with the exception of those in Figure 2, C and F, and were empirically weighted to balance the goodness of fit across them. The cost function was minimized with a modified version of a fast variable metric optimizer with quadratic convergence (Fletcher and Powell, 1963; Press et al., 1992). The gradients of the cost function were calculated numerically. For efficiency, the fitting algorithm was coded for parallel computation, taking advantage of the currently available multiprocessor computers. Depending on model size and number of free parameters, the optimization reached convergence in <50–100 iterations, generally taking on the order of 1 h on a dual quad-core 3.0 GHz Intel Xeon processor, running Windows XP Pro 32 bit SP3. Each optimization was restarted several times with different initial parameters, to avoid local minima.
Results
Voltage-clamp experiments for kinetic analysis
Comprehensive information about Nav channel kinetics was collected from a suite of voltage-clamp experiments (Fig. 1). The logic of these protocols is to force the channels to visit as many conformational states as possible, via different trajectories. For example, during a voltage step from a negative to a positive potential and back, the channels change state from closed to open to inactivated, and then from inactivated to noninactivated. In principle, any such trajectory depends on all the rate constants of a kinetic model, but each pathway is predominantly influenced by a certain combination of rates. All these voltage-dependent transients, together with initial and final state occupancies, can be combined to extract rate constants and to formulate model structure (Milescu et al., 2005).
Voltage-clamp experiments for kinetic analysis of Na+ currents. Whole-cell recordings of Na+ currents from medullary raphé neurons in rat brainstem slices (postnatal days 0–4). A, Voltage-dependent activation, and time course of activation and fast inactivation. B, Time course of deactivation. C, Slow inactivation of persistent current. D, Voltage-dependent availability. E, Subthreshold inactivation. F, Recovery from inactivation. G, Entry into slow inactivation. H, Dynamic availability in the interspike interval. I, Activity-dependent slow inactivation. In A and B, a voltage prepulse was introduced (see insets) to inactivate out-of-control axonal Na+ currents. Solid red lines in A–C, and dotted red lines in E–G are exponential or biexponential fits, with time constants τ (Fig. 4A).
Thus, voltage sensor activation/deactivation and pore opening/closing were probed with step protocols as in Figure 1, A and B. Availability was resolved with inactivating steps of 50 ms, 300 ms and 2.5 s (Fig. 1D), matching the spiking frequency range of raphé neurons (1–25 Hz). The transient properties of inactivation were obtained from the time course of current decay (Fig. 1A,C), and from the time course of availability (Fig. 1E,F). The time course of entry into slow inactivation was examined as in Figure 1G. The development of inactivation within the interspike interval was probed with a protocol resembling the action potential and the interspike depolarization of a raphé neuron firing at ∼3 Hz (Fig. 1H). This synthetic protocol was used instead of action potential waveforms to reduce fitting time. A pulse protocol resembling trains of action potentials at different frequencies (interpulse intervals of 500, 200 and 80 ms, corresponding to frequencies of 2, 5 and 12.5 Hz) provided crucial information about activity-dependent inactivation (Fig. 1I).
It is generally very difficult to obtain good recordings of Na+ currents from intact neurons in brain slices, because of poor space clamp (Rall and Segev, 1985). Moreover, the large amplitude of Na+ currents and the greater access resistance typical for these in situ recordings may cause substantial voltage-clamp errors. To reduce these problems, we decreased Na+ current amplitude with elevated internal Na+ concentration (50 mm). Furthermore, we introduced in some protocols (Fig. 1A,B) a brief voltage prepulse (8 ms at −40 mV) that eliminates the out-of-control current spikes characteristic of slice recordings and overall improves space clamp. The prepulse intentionally triggers axonal spikes and thus transiently inactivates axonal Nav channels but leaves somatic channels largely unperturbed (Milescu et al., 2010).
The voltage prepulse has a time-limited effect (20–50 ms), and it could not be introduced in the slow inactivation protocol (Fig. 1C). In this case, only data with minimal evidence of artifacts were selected. Imperfect space clamp was less of a concern with all other protocols, where information is contained in the relative amplitude of a peak current. By simulations with a compartmental model, as in (Milescu et al., 2010), we determined that in this case the error was only a slight shift in half-inactivation voltage (1–2 mV) and a small bias in time constants (10–20%; results not shown). Overall, all protocols were optimized to reduce the magnitude of space-clamp artifacts, and only the cleanest recordings were selected for analysis.
Nav channel kinetic model
From the recordings illustrated in Figure 1, a set of data was prepared for model-based fitting, as shown in Figure 2 (black symbols and traces). The fitting procedure was based on a maximum likelihood method for direct estimation of rate constants from macroscopic ion channel currents, for kinetic models of arbitrary size and topology (Milescu et al., 2005). The objective was to find a kinetic model that could be an accurate functional substitute for the Nav channels in the cell. Although we sought a biophysically realistic model, a more detailed analysis was not practical, considering both the inherent restrictions of somatic recordings in slices (imperfect space clamp, limited resolution, mixture of channel types, etc.) and the potentially very complex Nav channel kinetics. Nevertheless, we were able to formulate a model that globally explains a variety of physiologically relevant data.
Model fitting results. Data (A–I, black symbols and traces) were prepared from the voltage-clamp experiments shown in Figure 1, and were simultaneously fitted. Model II (Fig. 3B) fits well (green curves) all the data except the availability curves during slow or fast repetitive spiking (I). Model III (Fig. 3C) fits well all the data (red curves). A, n = 3 for activation, n = 9, 13, 12 for inactivation with 50 ms, 300 ms and 2.5 s steps, respectively; B, average of 3 datasets; C, n = 16; D, n = 1 … 7; E, n = 6 . . 18; F, n = 1; G, n = 1, repeated three times; H, n = 26; I, n = 6. All error bars are SD.
We started from the allosteric state model (model I in Fig. 3A) proposed by Kuo and Bean (1994). Unlike the original model, however, all rates were allowed to be voltage-dependent. With this modification, the experimentally observed voltage sensitivity of inactivation is ascribed not only to the coupling of inactivation rates to the voltage-dependent activation pathway (encoded by the a and b factors in model I) but also to the intrinsic voltage sensitivity of inactivation rates (αh, βh, αho, and βho). Alternative parameterizations with voltage-independent inactivation rates were not flexible enough to explain our data. A simpler model, similar to the Hodgkin-Huxley formulation (Hodgkin and Huxley, 1952) but with four rather than three activation and one inactivation particles (m4h) and all rates exponential functions of voltage, was also inadequate.
Nav channel kinetic models. State models fitted to the data shown in Figure 2. A, Model I (Kuo and Bean, 1994) cannot explain slow inactivation. B, Model II accounts for slow inactivation with one extra state (I13), which provides a secondary pathway for inactivation (O6 → I13), with fast entry and slow exit, briefly competing with open state fast inactivation (O6 → I12). C, Model III has two gating modes, assembled from two model II schemes, differing only in their slow inactivation properties. Channels reside mostly in mode I during slow spiking, with a fraction switching to mode II during fast spiking. The fits with model II and model III are shown in Figure 2 (green and red curves, respectively). In all models, rates were exponential functions of voltage (k = k0 × ek1 × V). Numerical values are given in supplemental Table 1, available at www.jneurosci.org as supplemental material.
Model I could not at all account for the slow inactivation detected in the data (e.g., the slow components in Fig. 1E,F). Within the structural constraints of this model, no set of kinetic parameters could be found to support both fast and slow inactivation properties. However, by trial and error we discovered that slow inactivation could be modeled reasonably well (Fig. 2, green traces) with just one additional state. This kinetic scheme (Fig. 3B, model II) predicts that, during brief depolarizations resembling spikes, fast inactivation follows two pathways: from O6 to I12 (≈80%) and from O6 to I13 (≈20%). We note that the 80/20% ratio is valid when all channels are initially noninactivated but this ratio will change during repetitive pulse application. Upon repolarization, the fraction residing in states I12 and I11 recovers quickly, while the fraction residing in I13 recovers slowly, which explains the slow component apparent in the recovery time course (Figs. 1F, 2E).
Although in model II the slow inactivated state I13 is connected to the open state O6, we also tried other configurations and obtained qualitatively similar fits (data not shown). Thus, state I13 could be connected to either state C5 or I12. Also, an additional inactivated state connected to both C5 and I13 slightly improved the fits. The idea is that at least one state (I13) is necessary to introduce an additional inactivation pathway, characterized by fast entry and slow exit, briefly competing with the standard fast inactivation (C5 → O6 → I12) during channel activation. Although we cannot exclude other, possibly more complex configurations, the scheme with I13 connected to O6 is parsimonious, fits sufficiently well the data, and resembles the model proposed for the resurgent current (Raman and Bean, 2001). Moreover, the connectivity of state I13 is similar to that found by an automated model search procedure, applied to Nav channels in CA1 neurons (Menon et al., 2009), which exhibit inactivation properties similar to raphé neurons.
Model II fits well the voltage-clamp data (Fig. 2, green traces) but has a functionally relevant shortcoming: it has too much sensitivity to spiking frequency, in the sense that availability is predicted to decay too little or too much at low or high frequencies, respectively, compared with the data (Fig. 2I). Although we tried to expand the topology by intuitively adding states or connections, no simple model could be derived that explained better the pulse train data (Fig. 2I) without failing the other data. A possible solution was model III (Fig. 3C), which contains two interconnected kinetic modes, each identical to model II. The two modes differ only in the O6I13 and O19I26 rates.
The logic of model III is that during slow spiking channels reside mostly in mode 1, which by itself has slow inactivation qualitatively similar to model II. At high spiking frequencies, a fraction of the channels shifts occupancy to mode 2, which has no slow inactivation (state I26 was included mainly for symmetry, although it also improved the fit). This modal shift takes place when the open state O6 is transiently occupied during an action potential, which explains its dependence on spiking frequency (the fractional residence in state O6 is proportional to the spiking rate). Model III fits well all the data (Fig. 2, red traces), and has correct sensitivity to spiking frequency (Fig. 2I). Moreover, although large, this model is still reasonably parsimonious, as most parameters are identical between the two gating modes.
The kinetic parameters corresponding to the best fits are listed in supplemental Table 1, available at www.jneurosci.org as supplemental material, as found by our automated fitting algorithm. We note that the large amount of data and the inherent molecular and cellular heterogeneity and noise have made optimization difficult. Under these conditions, the search algorithm converged robustly and greatly minimized the error, up to a point where the parameter surface became shallow, with nonunique solutions. Once this level was reached, the search had to be stopped. At this point, 10–20% changes in the estimated parameters had small effects on the goodness of fit.
From fitting these data, we learned that all of the voltage-clamp protocols shown in Figure 1 were necessary to obtain a well constrained and functionally correct model. The data shown in Figure 1, A, B, and D, constrain basic kinetic properties (activation, deactivation and fast inactivation), those in 1C and 1E–H constrain fast and slow inactivation, while the data in 1I constrain activity-dependent inactivation. Although the presence of a mixture of Nav channel types in these data cannot be excluded, the experimentally observed properties are also fairly consistent with a kinetically homogenous population.
A summary of Nav channel kinetic properties
Model III matches the time constants detected in the data (Fig. 4A). The overlap between different time constants along the voltage scale is evidence of data and model self-consistency. The idea is that the time constants characterizing the transition of a channel between two state occupancy distributions depend on rate constants and voltage but not on initial state occupancy. For example, inactivation at a given voltage and recovery from inactivation at the same voltage should proceed exponentially with a mixture of identical time constants but different amplitudes. This identity is more apparent for those time constants that are well separated from the rest, such as the slow component of inactivation or recovery from inactivation.
Summary of Nav channel kinetic properties. A, Stationary and transient kinetic properties extracted from the experimental data shown in Figure 2 (solid symbols) or data simulated with model III (dashed lines or clear symbols). The gray boxes denote the voltage and time ranges of action potential (AP) and interspike intervals (ISI). The plotted quantities are: time-to-peak (tP); time constants of deactivation (τf,D), inactivation (τf,I, τm,I, τs,I), recovery from inactivation (τf,RI, τs,RI), subthreshold inactivation (τf,SI), and entry into slow inactivation (τf,ES, τs,ES; only one value, denoted by solid or clear symbols); stationary activation and availability curves (gray symbols represent data, dotted lines are Boltzmann fits). All error bars are SEM. B, Current-clamp recording illustrating the spiking cycle of a raphé neuron. Also shown are the voltage and time ranges of stationary and transient activation and inactivation, as predicted by the model.
In Figure 4A, the time constant of slow recovery (τs,RI) overlaps with subthreshold slow inactivation (τs,SI), slow inactivation of the persistent current (τs,I), and slow entry into slow inactivation (τs,ES). Remarkably, they are all the result of a single transition: O6I13 (we know this because model II has qualitatively similar properties). The same overlap exists between fast recovery (τf,RI), subthreshold fast inactivation (τf,SI), and suprathreshold inactivation (τm,I and τf,I). Interestingly, the time constant of fast inactivation changes more than one order of magnitude from −50 mV (τf,SI or τf,RI) to −30 mV (τf,I). This abrupt change likely reflects the switch at ≈ −40 mV from closed state inactivation (slower) to open state inactivation (faster).
Two components were detected in suprathreshold inactivation (τf,I and τm,I) but the model predicts only one (τf,I). Although the slower component could be a genuine kinetic property, it could also be an artifact, and we decided not to model it (but see Irvine et al., 1999). The time-to-peak (tP) values predicted by the model are slightly off from the experimental values. However, we have not tried to improve the fit in this regard, because the first few hundred μs after a voltage jump are the most strongly affected by resolution artifacts, and data are generally not reliable. At very negative voltages (<−80 mV), tail currents (Fig. 1B) are so fast that they are most likely subject to low-pass filtering, which means that in reality deactivation may be faster (smaller τf,D values). For clarity, we left out the slower component detected in the deactivation protocol (τm,D in Fig. 1B), which could also be an artifact. The activity of Nav channels during the spiking cycle of a raphé neuron can be qualitatively predicted from these kinetic properties, as illustrated in Figure 4B.
Functional model verification with dynamic clamp
Although model III (Fig. 3C) explains well the voltage-clamp data (Fig. 2, red plots), the critical test would be a functional validation: can a model-based, dynamic clamp-injected current replace the native Na+ current in a neuron, and induce similar firing properties? For this test, we compared the spiking patterns between control conditions (“control”), and after Nav channels were pharmacologically blocked with TTX and a virtual Na+ conductance was introduced via dynamic clamp, based on model III (“model”). The conductance value GNa was chosen so as to match the maximum slope of the action potential upstroke between control and model. In most of the examined cells, this value was ≈10 nS/pF.
Figure 5 illustrates typical results. The overall spiking pattern and the shape of the action potential are similar between control and model, as shown for levels of injected current IInj = 5 pA (Fig. 5A,B) and IInj = 50 pA (Fig. 5C,D). A striking feature is the time-dependent change in action potential shape, which is also reproduced well by model. The instantaneous spiking frequency (calculated as the inverse of the interspike interval) is strongly time- and injected current-dependent (Fig. 5E). For Iinj < 20 pA the frequency decays slowly but a faster component emerges when Iinj ≥ 20 pA. At 100 pA and more, the cells quickly enter depolarization block (data not shown). The time course of frequency decay and the current dependence are similar between control and model. The frequency versus current relationship is linear and steepest for the first interspike interval but becomes nonlinear and shallow during sustained spiking (Fig. 5F).
Functional model validation with dynamic clamp. The model-based current is a good substitute for the native Na+ currents, as demonstrated by a comparison between native spiking (control), and after TTX-block of Na+ currents and dynamic-clamp injection of current based on model III (model). All recordings are from the same neuron. A, Spiking pattern during 10 s constant current injection (Iinj = 5 pA), preceded by 5 s voltage clamp at −70 mV. B, Action potential shape. The black, blue and gray traces highlight, the first, last, and intervening spikes in the train, respectively. C, D, Same as in A and B, but with IInj = 50 pA. E, Instantaneous spiking frequency versus time, for Iinj = 5, 20, and 50 pA. The solid lines are exponential fits. F, Spiking frequency versus injected current, for the first spike and at 1, 5 and 10 s from the beginning of current injection. The frequency values are averages over 1… 5 spikes. The dotted lines are polynomial fits. Error bars are SD.
There are also some differences between control and model. Most visibly, the action potential starts more abruptly in control than in model (Fig. 5B,D). However, this is not a failure of the kinetic model but a technical limitation of dynamic clamp. Normally, in brain slice preparations action potentials are initiated in the axon and back-propagate to the soma, causing the abrupt onset (Yu et al., 2008). Whereas real Nav channels are compartmentalized within soma, dendrites, and axon, the model-based current is injected strictly in the soma, which becomes the place of action potential initiation. In fact, this configuration resembles the case of dissociated neurons that have lost their axon and exhibit similarly smoothly rising action potentials. In a previous study (Milescu et al., 2008), we have shown that the back-propagation and the abrupt action potential onset can be reproduced in a dynamic-clamp experiment by including a virtual axonal compartment. Nevertheless, it is remarkable that basic spiking properties can be reproduced so well by a soma-injected model-based current.
Another difference between control and model resides in the time course of spiking frequency, particularly during fast spiking (Fig. 5E,F). In addition to the causes discussed above, one potential explanation for this mismatch is the long time elapsed between control and model experiments (a few minutes), during which other currents as well as whole-cell recording properties may change. Additionally, the activity-dependent behavior of the model, which potentially modulates spiking frequency, is mostly based on the information contained in the voltage-clamp protocol shown in Figure 1I. It is possible that the pulse trains that we have used for fitting were too short, and the model is not perfectly constrained in terms of long-term activity. However, this aspect could be easily adjusted through the transitions connecting the two gating modes.
Nav channel activity during repetitive spiking
The functional accuracy of model III suggests that it would be possible to inspect the dynamics of Nav channels during the spiking cycle, indirectly, through the activity of the model in a dynamic-clamp experiment. The quantities of interest are the state occupancies and the output current (INa), all of which are accessible. However, since there is no guarantee that the structure of the model is correct, we preferred to examine channel activity in terms of more general quantities. These voltage- and time-dependent variables are: the fraction of noninactivated channels available for spiking (calculated as the sum occupancy of noninactivated states), the occupancy of the open state(s), the occupancy of the slow inactivated state(s), and the current INa. Also, it would be interesting to compare the dynamic channel availability with the equilibrium availability, which is voltage-dependent but time-independent, and is available off-line as the availability curve, as plotted in Figure 2A for 2.5 s inactivating steps.
During slow, spontaneous spiking (Iinj = 0), the peak values reached by the model-based current INa during spikes exhibit a slow decay (Fig. 6A,C, second panel from top). Surprisingly, the subthreshold INa flowing in the interspike interval has the opposite trend (third panel from top). The fraction of available channels (bottom panel, trace labeled “Available”) follows the same trend as the peak INa, and is matched by an incremental increase in slow inactivation (“Slow inact.”), which reaches >30%. During spikes, the maximum open state occupancy (“Open”) reaches only 15–20%. It is interesting to note that both the open probability and the generated current reach their maximum values around −10 mV, half-way between spike threshold and spike peak.
Dynamics of Nav channels during repetitive firing. Dynamic-clamp experiment with the same neuron and current step protocol as in Figure 5. A, Time course plots of membrane voltage, current injected based on model III (INa) with subthreshold range expanded below, fraction of dynamically available channels (Available), occupancy of slow inactivated states (Slow inact.), and open probability (Open), under spontaneous spiking (Iinj = 0 pA). B, Same as in A but with Iinj = 50 pA. Also plotted is the equilibrium availability (Available eq.), identical to the 2.5 s availability curve in Figure 2A. C, One spiking cycle expanded from A (a). D, The first three spiking cycles expanded from B (b).
During fast spiking (Iinj = 50 pA), the state occupancies and the current INa undergo greater changes (Fig. 6B,D). The subthreshold INa reaches significantly greater values, and then decays slowly. Apparently, the slow inactivated state I13 acts as a spiking frequency-dependent reservoir of current, which is turned on at the more negative potentials in the interspike interval. However, this subthreshold INa predicted by model III is small in absolute terms (only a few pA) and we were not able to determine whether it exists in reality. It is quite remarkable that, although maximum availability drops to ≈30% and maximum open state occupancy is <5%, the neuron is still capable of spiking.
The relationship between the dynamic availability and the equilibrium availability demonstrates the importance of kinetics to channel function. During spontaneous spiking, the predicted equilibrium availability drops to very small levels during the few milliseconds before the action potential onset (Fig. 6C, bottom panel, trace labeled “Available eq.”). Thus, if inactivation were faster at these voltages, the channels would inactivate before they had a chance to open and generate current. However, the dynamic availability is still ≈40% at the spike onset (“Available”), after which it decays quickly but allowing INa sufficient time to generate the action potential. This discrepancy between equilibrium and dynamic prespike availabilities is enhanced by the difference between closed state and open state inactivation rates, which differ by a factor of ≈10 (Fig. 4A, compare τf,RI at −50 mV vs τf,I at −30 mV).
After an action potential is completed, with very little INa flowing during the downstroke, the dynamic availability is regenerated quickly, reaching a voltage-dependent maximum in ≈5 ms (Fig. 6C, Available). This maximum is only ≈65%, limited mainly by the occupancy of the slow inactivated states, which remains at ≈25%. During fast spiking, the maximum availability is further limited by the insufficient recovery time (Fig. 6D). After the first spike, the maximum is only ≈35%, increasing a little as the frequency drops, and then slowly decreasing as slow inactivation accumulates.
To explain why the maximum open probability is so low (<0.15 during spontaneous spiking, and <0.05 during fast spiking) we should note that the maximum Popen that can be reached during voltage steps is ≈0.42, comparable to other estimates of 0.5–0.6 (Horn and Vandenberg, 1984; Kimitsuki et al., 1990); without inactivation, this value would be >0.9. Furthermore, the prespike availability is <0.4 (Fig. 6C), because the half-inactivation voltage is constrained to provide maximum availability at prespike potentials (≈−50 mV) and minimum availability above the spike threshold (>−30 mV), and the steepness of the availability curve is limited by physical constraints. Thus, if the channel were subjected to an instantaneous voltage depolarization starting from the spike threshold, the maximum Popen could not be >0.42 × 0.4 ≈ 0.17. In reality, considering the finite slope of action potentials, Popenmax ≈ 0.12.
Nav channel slow inactivation causes spiking adaptation
During sustained firing, most raphé neurons exhibit strong spike shape and frequency adaptation (Fig. 5, control), paralleled by a progressive reduction in Nav channel availability during pulse train voltage-clamp protocols (Fig. 1I). This kinetic behavior is reproduced by model III (Fig. 2I, red fits), which also causes similar spiking adaptation in dynamic-clamp experiments (Fig. 5, model). Therefore, the logical question is whether Nav channel slow inactivation is the (only) cause for adaptation. If true, then a dynamic clamp-injected current based on model I, which does not exhibit slow inactivation, should not result in adaptation.
This dynamic-clamp test, performed with the two models that differ in their slow inactivation properties (model III and model I), demonstrates that the spiking patterns are strikingly different (Fig. 7A, top two panels). Most visibly, the peak voltage decays slowly and steadily with model III but remains virtually constant with model I; the trough voltage is slightly more elevated with model III. These differences are emphasized by the model switch experiment (Fig. 7A, bottom panel), where model III was replaced after 5 s with model I. Interestingly, spike shape adaptation occurs with either model (Fig. 7B). However, as in natural spiking (Fig. 5B,D, control), the action potentials produced by model III (Fig. 7B, top panel) broaden via a progressive reduction in both upward and downward slopes. In contrast, only the downward slope changes (somewhat less) with model I (Fig. 7B, bottom panel). The upward and downward slopes are proportional to inward Na+ and outward K+ currents, respectively. Thus, in control and model III experiments, both inward and outward peak currents decay during fast spiking, whereas with model I only the outward current changes (Fig. 7D). We note that the peak inward current predicted by model III decays to slightly lower levels than the native current in control (Fig. 7D, model III vs control plots), which is the result of imperfectly modeling the long-term behavior of the channels in response to repetitive pulses, as discussed earlier.
Nav channel slow inactivation causes spiking adaptation. Dynamic-clamp experiment with the same neuron and current step protocol as in Figure 5. A, Time course of membrane voltage under Iinj = 20 pA, with current injected based on model III, model I, or model III replaced after 5 s with model I. B, Broadening of action potential shape with model III or model I. The upward and downward slopes are proportional to inward Na+ and outward K+ currents, respectively. C, Instantaneous spiking frequency versus time, for Iinj = 5, 20, and 50 pA. For the model III/model I switch experiment, the frequency plot is shown only for Iinj = 50 pA. D, Time plots of the peak outward and inward currents flowing during the action potential, for Iinj = 50 pA. Solid lines are exponential fits.
Interestingly, frequency adaptation occurred with either model III or model I (Fig. 7C). However, the frequency decays faster with model III (Fig. 7C, black vs red plots). At Iinj = 50 pA, the time constants of the fast and slow components are (mean ± SD): 50 ± 10 ms and 6.2 ± 0.7 s for model III, and 105 ± 25 ms and 22.1 ± 20.6 s for model I. In the model switch experiment, the spiking frequency (Fig. 7C, green plot) jumps upon transition from model III to model I, which maintains a larger fraction of available channels, but the rates of decay during each phase are virtually identical to the rates obtained separately with each corresponding model.
Together, all of these results suggest that Nav channel slow inactivation is responsible for some but not all aspects of spiking adaptation. Thus, the spike amplitude and the broadening of the action potential upstroke depend entirely on Nav channel slow inactivation, whereas spiking frequency adaptation is enhanced by slow inactivation but occurs to some degree even without it.
Differential effects of suprathreshold and subthreshold Na+ currents
Nav channels are active throughout the entire spiking cycle, as seen in Figure 6, and may affect action potential shape and spiking frequency through multiple mechanisms. Subthreshold Na+ current (INasub), flowing during the interspike interval, may add up to the net current to increase membrane depolarization rate, and thus spiking frequency. Suprathreshold Na+ current (INasupra), flowing during the spike, can modify the action potential waveform and thus may change the activation and inactivation of other voltage-dependent currents. In turn, these other currents exert their own effects. Slow inactivation brings further complication, as INasub and INasupra may exhibit opposite trends during repetitive firing, as shown in Figure 6.
To understand how INasub and INasupra modulate firing properties, and to distinguish between their effects, we designed a simple dynamic-clamp experiment, which allows separate manipulation of the amplitudes of the two current components. First, we choose a voltage threshold VThr to separate the two components, based on the approximate point where the action potential starts to rise abruptly. In the case of our experiments with raphé neurons, VThr ≈ −42 mV. Then, the model-based injected current is multiplied by a factor f that depends on the measured membrane voltage Vm. For example, if the subthreshold current alone is to be amplified, we make f > 1 when the measured voltage Vm ≥ VThr, and f = 1 when Vm < VThr. In this way, the magnitude of the model-based current injected with dynamic clamp can be instantaneously modulated, according to the value of the measured voltage. Effectively, this procedure can be used to divide the injected current in functional components, for example sub- versus suprathreshold current, which can be differentially amplified, reduced, or eliminated, to test their individual roles. In principle, the procedure can be further expanded to take into account the phase of the spiking cycle, for example action potential upstroke versus downstroke.
Suprathreshold Na+ current controls action potential shape
We first examined the effects of the suprathreshold Na+ current. For easier understanding, the results are presented as if the suprathreshold Nav channel conductance GNasupra were varied in the range 2–50 nS/pF (in fact, the f factor was varied between 0.2 and 5), while the subthreshold conductance GNasub was maintained constant at 10 nS/pF. The conductance was expressed as density (nS/pF), for easier comparison with other studies, where the membrane capacitance value may be different. The results presented in Figure 8A demonstrate that the shape of the action potential waveform depends strongly on the magnitude of the suprathreshold Na+ current. Thus, the parameters most affected are action potential height and upward slope (Fig. 8B). When the suprathreshold current is very large (GNasupra = 50 nS/pF), the spike height approaches the reversal voltage for Na+ (≈60 mV).
Suprathreshold Na+ current controls spike shape. Dynamic-clamp experiment where the current injected based on model III was divided in subthreshold and suprathreshold components, according to a voltage threshold (VThr). A, Spiking pattern of a representative neuron, where the suprathreshold conductance GNasupra was varied in the range 2… 50 nS/pF (when Vm ≥ VThr), with constant subthreshold conductance GNasub = 10 nS/pF (when Vm < VThr). B, Action potential shapes with different GNasupra values. C, Net current (calculated as −Cm × dVm/dt, shown as solid lines) and model-based current (INa, shown as dotted lines) flowing during the spike.
While the model-based injected current INa is directly accessible during a dynamic-clamp experiment, the net membrane current can be calculated with the formula Im = −Cm × dVm/dt, where Cm is the estimated membrane capacitance, Vm is the measured voltage, and t is the time. The minus sign reflects the convention that a negative (inward) current depolarizes the membrane. As shown in Figure 8C (solid noisy traces), the direction of the net current Im is inward (depolarizing) during the action potential upstroke, and outward (repolarizing) during the downstroke. During the inward phase, Im is proportional with INa (Fig. 8C, dotted smooth traces), but slightly greater than INa for small GNasupra values (Fig. 8C, blue traces), probably due to activation of Ca2+ currents. For large GNasupra values, the potentials reached during the action potential peak are substantially more positive, causing a stronger activation of K+ currents, as reflected by the larger net current Im flowing during the outward phase (Fig. 8C, red traces). In turn, this greater outward current causes faster repolarization and enhanced afterhyperpolarization (Fig. 8B).
The relationship between the peak INa and GNasupra is approximately linear (Fig. 8C, inset, dotted plot). This relationship is similar for the inward peak of the net current Im (Fig. 8C, inset, solid plot) but somewhat distorted, presumably by the flow of inward Ca2+ currents at the low end, and by the flow of outward K+ currents at the high end. This approximately linear function implicitly applies to the relationship between the maximum depolarization rate of the action potential and the peak suprathreshold Na+ current.
The subthreshold Na+ current controls spiking frequency
The effect of the subthreshold Na+ current was tested in a similar experiment, where the subthreshold conductance GNasub was varied in the range of 2–50 nS/pF, while the suprathreshold conductance GNasupra was maintained constant at 10 nS/pF. As shown in Figure 9A, the firing pattern depends strongly on the subthreshold INa. The action potential shape is virtually unaffected by INasub, but the spike threshold is changed by up to ≈7 mV in this example (Fig. 9B). The model-predicted subthreshold INa increases proportionally to GNasub, from a fraction of a pA to ≈6 pA (Fig. 9C). Most importantly, increasing the subthreshold INa can shorten the interspike interval to half, doubling the spiking frequency (Fig. 9D, solid plot).
Subthreshold Na+ current controls spiking frequency. Dynamic-clamp experiment where the current injected based on model III was divided in subthreshold and suprathreshold components, according to a voltage threshold (VThr). A, Spiking pattern of the same neuron as in Figure 8, where the subthreshold conductance GNasub was varied in the range 2… 50 nS/pF (when Vm < VThr), with constant suprathreshold conductance GNasupra = 10 nS/pF (when Vm ≥ VThr). B, Action potential shapes with different GNasub values. C, Subthreshold model-based current (INa) flowing during the spiking cycle. D, Average spiking frequency as a function of subthreshold conductance GNasub (solid symbols), or suprathreshold conductance GNasupra (open symbols; the data are taken from the experiment shown in Fig. 8). The values are mean and SD calculated for each conductance step. The lines are polynomial fits.
Interestingly, analysis of the previous experiment (Fig. 8) reveals that the suprathreshold INa also affects spiking frequency (Fig. 9D, dotted plot), although not as strongly. Since in that experiment the suprathreshold INa was not manipulated within the interspike interval, its effects on spiking frequency must be indirectly exerted by modulation of other currents, via changes in action potential shape. Together, these results demonstrate the complexity of the relationships between sub- and suprathreshold Na+ currents and the firing pattern.
Pacemaking and subthreshold INa
It might seem disproportionate that a 25-fold increase in GNasub (from 2 to 50 nS/pF) would cause only a two- or threefold increase in spiking frequency, as shown in Figure 9D (solid plots). To understand this apparent disparity, we compared the net current Im (calculated as −Cm × dVm/dt) with the model-based injected current INa, as obtained in the previous experiment described in Figure 9. Thus, Im and INa are plotted versus voltage in Figure 10A, in a phase-plane representation that depicts both the voltage and the time dependence of the two currents. The five trajectories of Im shown in the figure are the average action potential waveforms calculated for GNasub = 2… 50 nS/pF; the single INa trajectory is for GNasub = 10 nS/pF.
Pacemaking and subthreshold Na+ current. A, Phase-plane representation of the average net current Im (calculated as −Cm × dVm/dt, low-pass filtered at 5 kHz) and model-based current INa (green trace) during the spiking cycle (time direction denoted by arrow), for the dynamic-clamp experiment described in Figure 8. Only the INa corresponding to GNasub = 10 nS/pF is shown. B, The subthreshold range expanded from A. The net current is shown low-pass filtered at 0.2 kHz. C, D, Same analysis as in A and B, for the dynamic-clamp experiment described in Figure 5, for Iinj = 50 pA. The blue, red and gray trajectories of the net current correspond to the first, last and intervening spikes during the current injection, respectively. The model-based current INa is shown only for the first and last spikes. E, F, Dynamic-clamp experiment where the model-based current INa was zeroed below a voltage threshold (VThr), which was varied in the range −80… −49 mV. Spontaneous spiking stopped when VThr ≥ −49 mV.
During the action potential upstroke, Im is approximately equal to INa. Interestingly, it appears that a larger subthreshold GNa causes a slightly greater suprathreshold net current Im to flow during the upstroke (Fig. 10A, blue vs red curves). This counterintuitive effect occurs because the spike threshold is lower when GNasub is larger (Fig. 9B), and the more negative voltage increases prespike Nav channel availability, which in turn increases INa and implicitly Im. Surprisingly, a closer examination of the subthreshold range reveals that, when GNasub was increased from 2 to 50 nS/pF, Im changed by only ≈2 pA in the −65… −55 mV range, even though INasub was incremented by ≈6 pA (Fig. 10B). This suggests that the subthreshold INa does not simply add to the net current Im in the interspike interval.
To understand how the net subthreshold current Im responds to external stimuli, we examined the data obtained from the experiment described in Figure 5, where up to 50 pA of bias current was injected into the neuron. The net current Im corresponding to a current injection of IInj = 50 pA is plotted versus voltage in Figure 10C, in a phase-plane representation. The trajectories corresponding to the first and the last action potential waveforms are highlighted in blue and red, respectively. At the beginning of the current injection (before the first spike), the value of the subthreshold Im (≈−50 pA) matches the value of Iinj (with opposite sign) but then it decays slowly to less than −10 pA (Fig. 10D). According to these results, the membrane responds to excitatory input with a counteracting, slowly activating hyperpolarizing current (or, alternatively, a slowly inactivating depolarizing current), which explains why spiking frequency has lower than expected sensitivity to either subthreshold INa or bias current Iinj.
Under spontaneous spiking (no excitatory input, Iinj = 0 pA), the subthreshold INa is significantly smaller than the net current Im during most of the interspike interval (Fig. 10B). This leads to the prediction that, in a certain voltage range, INa is not the current controlling depolarization and may even be redundant. We tested this hypothesis in a dynamic-clamp experiment, where the model-based current INa was zeroed whenever the membrane potential was more negative than a certain threshold (INa = 0 if Vm < VThr).
In the example shown in Figure 10E, no visible change in the spiking pattern occurred when VThr was more negative than −50 mV (black section in the voltage trace). However, spiking became irregular and slower when VThr = −50 mV (blue section), and completely stopped for VThr = −49 mV (red section). The phase-plane trajectories of Im and INa are plotted in Figure 10F, for VThr = −51, −50 and −49 mV. The trajectories corresponding to VThr = −51 mV show that, even though INa (Fig. 10F, green plots) was zeroed below −51 mV, the net current Im (black plot) has a finite value and is sufficient to drive the membrane potential past VThr, at which point the subthreshold INa takes over and eventually a spike is generated.
In this example, the net current Im becomes null when the membrane potential Vm ≈ −51 mV, as shown by the trajectory corresponding to VThr = −50 mV (blue plot). Ideally, the membrane voltage would remain stationary under these conditions. In reality, the membrane potential resides around −51 mV for relatively long periods of time (Fig. 10E, blue section; Fig. 10F, blue plot) but occasional fluctuations are sufficiently large to drive the voltage above this level. When VThr = −49 mV, these fluctuations are not sufficient, and the membrane potential is trapped at ≈−51 mV. According to these results, subthreshold INa is only required to depolarize the membrane from −50 mV up to the spike threshold of ≈−40 mV (Fig. 10F). Under −50 mV, it appears that INa is not necessary for spontaneous spiking in raphé neurons.
Discussion
The role of voltage-gated ion channels in action potential generation was explained by Hodgkin and Huxley, in their seminal analysis of spiking in the squid giant axon (Hodgkin and Huxley, 1952). Since then, our understanding of ion channel gating and neuronal firing mechanisms has continued to grow. In this study, we have examined the kinetic properties and the role played by voltage-gated Na+ channels during repetitive firing, in a mammalian central neuron that exhibits slow pacemaking and strong spiking adaptation. In the following, we summarize and discuss our findings.
Voltage-clamp data and kinetic modeling
We characterized the kinetic properties of Na+ currents in medullary raphé neurons, from a suite of voltage-clamp experiments in vitro (Fig. 1). To make the results more physiologically relevant, we recorded Na+ currents from intact neurons in brainstem slices, taking precautions (Milescu et al., 2010) to reduce voltage-clamp and space-clamp artifacts (Rall and Segev, 1985). From the voltage-clamp data, we formulated a detailed kinetic model (Fig. 3C, model III) that explains the transient and stationary voltage-dependent properties of Na+ currents (Figs. 2, 4), over a broad voltage range (−120 to +20 mV) and a time scale spanning six orders of magnitude (50 μs to 5 s). This quantitative analysis was made possible by recent advances in kinetic modeling algorithms and software that allowed us to quickly construct and fit very large state models from comprehensive datasets (Milescu et al., 2005).
The proposed model III can be regarded as a hierarchical structure. (1) The core, represented by model I (Fig. 3A), describes basic activation and inactivation kinetics with minimally acceptable complexity (Bean, 1981; Patlak, 1991; Kuo and Bean, 1994; Cha et al., 1999; Bezanilla, 2000; Raman and Bean, 2001; Armstrong, 2006). (2) One additional state (Fig. 3B, model II) parsimoniously explains a slow inactivation process with fast entry and slow exit. A similar construct has been found by an automated topology-search algorithm (Menon et al., 2009) and the same topology has been proposed to model the resurgent Na+ current (Raman and Bean, 2001). (3) Finally, model III (Fig. 3C) assembles two model II units to explain activity-dependent inactivation (Fig. 2I) through modal gating. Although further study is necessary, modal gating is an intriguing hypothesis, supported by evidence in Nav (Nilius et al., 1987; Zhou et al., 1991; Alzheimer et al., 1993; Brown et al., 1994), Kv (Michaelevski et al., 2003) and Cav channels (Delcour et al., 1993; Fellin et al., 2004; Luvisetto et al., 2004), and NMDA receptors (Legendre and Westbrook, 1991; Popescu and Auerbach, 2003).
Testing the computational model in real neurons
Model III fits well the voltage-clamp data (Fig. 2) but the ultimate reality check was testing the model in real neurons, with dynamic clamp. For this, a model-based current was injected in neurons as a substitute for the native TTX-sensitive Na+ currents, which were pharmacologically blocked. The model-based current was capable of generating spikes of correct shape, and reproduced the relationship between spiking frequency and excitatory input current, including spontaneous firing and spike shape and frequency adaptation (Fig. 5). It is interesting to mention that a technical limitation of dynamic clamp is that the model-based current can only be injected in the soma, which means that axonal or dendritic conductances cannot be tested using this approach. Nevertheless, a model-based current injected in the soma is apparently sufficient to restore the essential spiking properties of raphé neurons.
Until recently, it was not possible to run large ion channel state models in dynamic-clamp experiments, because no real-time integration method was fast and precise enough. This was a particular difficulty for Nav channels, which have fast and steeply voltage-dependent kinetics, and operate on time scales that require a resolution of 10–20 μs. Thus, most dynamic-clamp studies involving Nav channels have been limited to using Hodgkin-Huxley models (Vervaeke et al., 2006), which are much less demanding computationally but cannot easily capture the full kinetic complexity of Na+ currents. Our present study was made possible by new real-time integration methods and dynamic-clamp software (Milescu et al., 2008) that permitted us to accurately run models with 26 states at 50–100 kHz.
Nav channel kinetics and neuronal firing mechanisms
By substituting the model-based current for the real Na+ current, we gained insight into the activity of Nav channels during repetitive spiking (Fig. 6). This has always been possible with computer simulations but now we have the added confidence that the model is not only based on comprehensive data, it is also functionally accurate in real neurons. Using this approach, the relationship between Nav channel kinetics and neuronal firing mechanisms could be more directly and precisely examined. Thus, we learned here that Nav channels are optimally operating throughout the range of spiking frequencies characteristic of raphé neurons (1–25 Hz), in a way that maximizes the amount of current available for spike generation (Fig. 6).
As a direct result of model structure (model I, Fig. 3A), closed state inactivation is substantially slower than open state inactivation (Fig. 4). Furthermore, subthreshold inactivation can be strongly voltage dependent, despite the relatively weak voltage sensitivity intrinsic to inactivation rates. Three important imperatives are thus satisfied: slow inactivation at the spike onset (time constant of 10–20 ms) to ensure maximum current availability and sharp spikes, very fast inactivation at peak voltages (submillisecond) for good spike efficiency (Carter and Bean, 2009), and fast recovery during afterhyperpolarization (milliseconds) to sustain relatively fast spiking (20–30 Hz), albeit transiently.
Slow inactivation and spiking adaptation
Nav channels in raphé and other neurons (Mickus et al., 1999) exhibit a type of slow inactivation characterized by fast entry at spike potentials, and slow recovery at subthreshold potentials (Figs. 1F, 2E). The functional result is activity-dependent, incremental accumulation into slow inactivation with each spike (Fig. 6), accompanied by reduction in effective availability. Recent studies link this behavior to proteins homologous to fibroblast growth factors (Goldfarb et al., 2007). The structure of model II implies that these proteins compete with fast inactivation for the open pore structure (O6I13 vs O6I12 transition). This interaction itself may be frequency-dependent, as phenomenologically described by model III.
Activity-dependent slow inactivation limits the fraction of available Nav channels to ≈70% during spontaneous spiking and ≈50% during fast spiking (Fig. 6). Hence, a third of the 70% effectively available channels can be functionally silenced, making slow inactivation a suitable candidate for Nav channel-dependent spiking adaptation. This hypothesis (Fleidervish et al., 1996) cannot be easily verified pharmacologically, without selective drugs. Here, we simply compared the spiking patterns generated by different Nav channel models, in a dynamic-clamp experiment. In this way, we could unambiguously demonstrate that Nav channel slow inactivation does indeed cause spiking adaptation in raphé neurons (Fig. 7), although in cooperation with other factors.
Na+ currents control spike shape and frequency and are involved in pacemaking
A truly remarkable power of dynamic clamp is the ability to run a simulation of a computational model in a real neuron. Not only can a model-based current be used as a substitute for a native current to test the relationship between ion channel kinetics and neuronal firing, but this model-based current can be manipulated in real-time to dissect the functional roles of individual current components during the different phases of the spiking cycle. Here, we took this approach to examine the differential contributions of subthreshold and suprathreshold Na+ currents to action potential shape, spiking frequency, and spontaneous spiking. These results would have been difficult to obtain by pharmacological means.
For these experiments, we defined the subthreshold and suprathreshold Na+ currents simply as the currents flowing below and above the spike threshold (≈−40 mV). These two components were separately amplified or reduced, via a voltage-dependent factor that was introduced in the model-based current equation. Thus, we found that the shape of the action potential depends strongly on the amount of suprathreshold current (Fig. 8). The reason is that a larger Na+ current causes faster depolarization, allowing spikes to reach higher peak voltages. In turn, these greater potentials cause increased activation of voltage-gated K+ channels, which results in faster repolarization. In contrast, we found that the subthreshold Na+ current has virtually no effect on spike shape but controls spiking frequency (Fig. 9).
Finally, we tested the contribution of subthreshold Na+ current to the spontaneous spiking of raphé neurons. The pacemaking mechanism of raphé neurons is currently not understood but it seemed logical that subthreshold Na+ current would be involved, considering its effect on spiking frequency. Taking a similar approach, the model-based current was zeroed below a voltage threshold, which was varied to determine the operational range of the subthreshold Na+ current. Our results (Fig. 10) show that currents other than Na+ depolarize the membrane from the afterhyperpolarization trough (≈−70 mV) up to −50 mV, possibly similar to the subthreshold currents found in other slow pacemaking neurons (Khaliq and Bean, 2010). Moreover, the net current flowing in this subthreshold range features some homeostatic mechanism that is involved in spiking adaptation, which remains to be established. Nevertheless, our experiments demonstrate unambiguously that the Na+ current is strictly required above ≈−50 mV, to depolarize the membrane toward the spike threshold (≈−40 mV) and generate spontaneous spiking in raphé neurons.
In conclusion, this study demonstrates that the relationships between the kinetic properties of voltage-gated sodium channels and the firing properties of pacemaker neurons can be quantitatively studied using a combination of experiment and recently developed computational tools. The analysis and methodology presented here can be applied to other types of voltage-gated channels and neurons, to advance our understanding of neuronal excitability.
Footnotes
This work was supported by the Intramural Research Program of the National Institutes of Health (NIH), National Institute of Neurological Disorders and Stroke. We thank our colleagues for discussing the experiments and the manuscript, particularly Drs. Mirela Milescu and Shai Silberberg (NIH, Bethesda, MD), Andrew Plested (Max Delbrück Center for Molecular Medicine, Berlin, Germany), and Joel Tabak (Florida State University, Tallahassee, FL). Dr. Ruli Zhang (NIH, Bethesda, MD) helped with experimental procedures. We are grateful to Dr. Bruce Bean (Harvard Medical School, Boston, MA) for suggestions and support.
- Correspondence should be addressed to Lorin S. Milescu, Department of Neurobiology, Harvard Medical School, 220 Longwood Avenue, Goldenson Building, Room 420, Boston, MA 02115. Lorin_Milescu{at}hms.harvard.edu
