Estimation of ion channel kinetics from fluctuations of macroscopic currents

For single channel recordings, the maximum likelihood estimation (MLE) of kinetic rates and conductance is well established. A direct extrapolation of this method to macroscopic currents is computationally prohibitive: it scales as a power of the number of channels. An approximated MLE that ignored the local time correlation of the data has been shown to provide estimates of the kinetic parameters. In this article, an improved approximated MLE that takes into account the local time correlation is proposed. This method estimates the channel kinetics using both the time course and the random fluctuations of the macroscopic current generated by a homogeneous population of ion channels under white noise. It allows arbitrary kinetic models and stimulation protocols. The application of the proposed algorithm to simulated data from a simple three-state model on nonstationary conditions showed reliable estimates of all the kinetic constants, the conductance and the number of channels, and reliable values for the standard error of those estimates. Compared to the previous approximated MLE, it reduces by a factor of 10 the amount of data needed to secure a given accuracy and it can even determine the kinetic rates in macroscopic stationary conditions.     

Maximum likelihood estimation of ion channel kinetics from macroscopic currents

We describe a maximum likelihood method for direct estimation of rate constants from macroscopic ion channel data for kinetic models of arbitrary size and topology. The number of channels in the preparation, and the mean and standard deviation of the unitary current can be estimated, and a priori constraints can be imposed on rate constants. The method allows for arbitrary stimulation protocols, including stimuli with finite rise time, trains of ligand or voltage steps, and global fitting across different experimental conditions. The initial state occupancies can be optimized from the fit kinetics. Utilizing arbitrary stimulation protocols and using the mean and the variance of the current reduce or eliminate problems of model identifiability (Kienker, 1989). The algorithm is faster than a recent method that uses the full autocovariance matrix (Celentano and Hawkes, 2004), in part due to the analytical calculation of the likelihood gradients. We tested the method with simulated data and with real macroscopic currents from acetylcholine receptors, elicited in response to brief pulses of carbachol. Given appropriate stimulation protocols, our method chose a reasonable model size and topology.     

Use of the covariance matrix in directly fitting kinetic parameters: application to GABAA receptors

A new method of analysis is described that begins to explore the relationship between the phases of ion channel desensitization and the underlying states of the channel. The method, referred to as covariance fitting (CVF), couples Q-matrix calculations with a maximum likelihood algorithm to fit macroscopic desensitization data directly to kinetic models. Unlike conventional sum-of-squares minimization, CVF fits both the magnitude of the recorded current and the strength of the correlations between different time points. When applied to simulated data generated using various kinetic models with up to 11 free parameters, CVF leads to reasonable parameter estimates. Coupled with the likelihood ratio test, it accurately discriminates between models with different numbers of states, discriminates between most models with the same number but a different arrangement of states, and extracts meaningful information on the relationship between the desensitized states and the phases of macroscopic desensitization. When applied to GABA(A) receptor traces (outside out patches, alpha 1 beta 2 gamma 2S, 1 mM GABA, >2.5 s), a model with two open states and three desensitized states is favored. When applied to simulated data generated using a consensus model, CVF leads to reasonable parameter estimates and accurately discriminates between this and other models.     

Global parameter optimization for cardiac potassium channel gating models.

Quantitative ion channel model evaluation requires the estimation of voltage dependent rate constants. We have tested whether a unique set of rate constants can be reliably extracted from nonstationary macroscopic voltage clamp potassium current data. For many models, the rate constants derived independently at different membrane potentials are not unique. Therefore, our approach has been to use the exponential voltage dependence predicted from reaction rate theory (Stevens, C. F. 1978. Biophys. J. 22:295-306; Eyring, H., S. H. Lin, and S. M. Lin. 1980. Basic Chemical Kinetics. Wiley and Sons, New York) to couple the rate constants derived at different membrane potentials. This constrained the solution set of rate constants to only those that also obeyed this additional set of equations, which was sufficient to obtain a unique solution. We have tested this approach with data obtained from macroscopic delayed rectifier potassium channel currents in voltage-clamped guinea pig ventricular myocyte membranes. This potassium channel has relatively simple kinetics without an inactivation process and provided a convenient system to determine a globally optimized set of voltage-dependent rate constants for a Markov kinetic model. The ability of the fitting algorithm to extract rate constants from the macroscopic current data was tested using "data" synthesized from known rate constants. The simulated data sets were analyzed with the global fitting procedure and the fitted rate constants were compared with the rate constants used to generate the data. Monte Carlo methods were used to examine the accuracy of the estimated kinetic parameters. This global fitting approach provided a useful and convenient method for reliably extracting Markov rate constants from macroscopic voltage clamp data over a broad range of membrane potentials. The limitations of the method and the dependence on initial guesses are described.     

Gating characteristics of a steeply voltage-dependent gap junction channel in rat Schwann cells

The gating properties of macroscopic and microscopic gap junctional currents were compared by applying the dual whole cell patch clamp technique to pairs of neonatal rat Schwann cells. In response to transjunctional voltage pulses (Vj), macroscopic gap junctional currents decayed exponentially with time constants ranging from < 1 to < 10 s before reaching steady-state levels. The relationship between normalized steady-state junctional conductance (Gss) and (Vj) was well described by a Boltzmann relationship with e-fold decay per 10.4 mV, representing an equivalent gating charge of 2.4. At Vj > 60 mV, Gss was virtually zero, a property that is unique among the gap junctions characterized to date. Determination of opening and closing rate constants for this process indicated that the voltage dependence of macroscopic conductance was governed predominantly by the closing rate constant. In 78% of the experiments, a single population of unitary junctional currents was detected corresponding to an unitary channel conductance of approximately 40 pS. The presence of only a limited number of junctional channels with identical unitary conductances made it possible to analyze their kinetics at the single channel level. Gating at the single channel level was further studied using a stochastic model to determine the open probability (Po) of individual channels in a multiple channel preparation. Po decreased with increasing Vj following a Boltzmann relationship similar to that describing the macroscopic Gss voltage dependence. These results indicate that, for Vj of a single polarity, the gating of the 40 pS gap junction channels expressed by Schwann cells can be described by a first order kinetic model of channel transitions between open and closed states.     

Inactivation viewed through single sodium channels

Recordings of the sodium current in tissue-cultured GH3 cells show that the rate of inactivation in whole cell and averaged single channel records is voltage dependent: tau h varied e-fold/approximately 26 mV. The source of this voltage dependence was investigated by examining the voltage dependence of individual rate constants, estimated by maximum likelihood analysis of single channel records, in a five-state kinetic model. The rate constant for inactivating from the open state, rather than closing, increased with depolarization, as did the probability that an open channel inactivates. The rate constant for closing from the open state had the opposite voltage dependence. Both rate constants contributed to the mean open time, which was not very voltage dependent. Both open time and burst duration were less than tau h for voltages up to -20 mV. The slowest time constant of activation, tau m, was measured from whole cell records, by fitting a single exponential either to tail currents or to activating currents in trypsin-treated cells, in which the inactivation was abolished. tau m was a bell-shaped function of voltage and had a voltage dependence similar to tau h at voltages more positive than -35 mV, but was smaller than tau h. At potentials more negative than about -10 mV, individual channels may open and close several times before inactivating. Therefore, averaged single channel records, which correspond with macroscopic current elicited by a depolarization, are best described by a convolution of the first latency density with the autocorrelation function rather than with 1 - (channel open time distribution). The voltage dependence of inactivation from the open state, in addition to that of the activation process, is a significant factor in determining the voltage dependence of macroscopic inactivation. Although the rates of activation and inactivation overlapped greatly, independent and coupled inactivation could not be statistically distinguished for two models examined. Although rates of activation affect the observed rate of inactivation at intermediate voltages, extrapolation of our estimates of rate constants suggests that at very depolarized voltages the activation process is so fast that it is an insignificant factor in the time course of inactivation. Prediction of gating currents shows that an inherently voltage-dependent inactivation process need not produce a conspicuous component in the gating current.     

Fractal methods to analyze ion channel kinetics

We describe the traditional nonfractal and the new fractal methods used to analyze the currents through ion channels in the cell membrane. We discuss the hidden assumptions used in these methods and how those assumptions lead to different interpretations of the same experimental data. The nonfractal methods assumed that channel proteins have a small number of discrete states separated by fixed energy barriers. The goal was to determine the parameters of the kinetic diagram, which are the number of states, the pathways between them, and the kinetic rate constants of those pathways. The discovery that these data have fractal characteristics suggested that fractal approaches might provide more appropriate tools to analyze and interpret these data. The fractal methods determine the characteristics of the data over a broad range of time scales and how those characteristics depend on the time scale at which they are measured. This is done by using a multiscale method to accurately determine the probability density function over many time scales and by determining how the effective kinetic rate constant, the probability of switching states, depends on the effective time scale at which it is measured. These fractal methods have led to new information about the physical properties of channel proteins in terms of the number of conformational substates, the distribution of energy barriers between those states, and how those energy barriers change with time. The new methods developed from the fractal paradigm shifted the analysis of channel data from determining the parameters of a kinetic diagram to determining the physical properties of channel proteins in terms of the distribution of energy barriers and/or their time dependence.     

A sodium channel gating model based on single channel, macroscopic ionic, and gating currents in the squid giant axon.

Sodium channel gating behavior was modeled with Markovian models fitted to currents from the cut-open squid giant axon in the absence of divalent cations. Optimum models were selected with maximum likelihood criteria using single-channel data, then models were refined and extended by simultaneous fitting of macroscopic ionic currents, ON and OFF gating currents, and single-channel first latency densities over a wide voltage range. Best models have five closed states before channel opening, with inactivation from at least one closed state as well as the open state. Forward activation rate constants increase with depolarization, and deactivation rate constants increase with hyperpolarization. Rates of inactivation from the open or closed states are generally slower than activation or deactivation rates and show little or no voltage dependence. Channels tend to reopen several times before inactivating. Macroscopic rates of activation and inactivation result from a combination of closed, open and inactivated state transitions. At negative potentials the time to first opening dominates the macroscopic current due to slow activation rates compared with deactivation rates: channels tend to reopen rarely, and often inactivate from closed states before they reopen. At more positive potentials, the time to first opening and burst duration together produce the macroscopic current.     

Estimating kinetic parameters for single channels with simulation. A general method that resolves the missed event problem and accounts for noise.

Analysis of currents recorded from single channels is complicated by the limited time resolution (filtering) of the data which can prevent the detection of brief intervals. Although a number of approaches have been used to correct for the undetected intervals (missed events) when identifying kinetic models and estimating parameters, none of them provide a general method which takes into account the true effects of noise and limited time resolution. This paper presents such a method. The approach is to use simulated single-channel currents to incorporate the true effects of filtering and noise on missed events and interval durations. The simulated currents are then analyzed in a manner identical to that used to analyze the experimental currents. An iterative search process using likelihood comparison of two-dimensional dwell-time distributions obtained from the simulated and experimental single-channel currents then allows the most likely rate constants to be determined. The large errors and false solutions that can result from the more typically applied assumptions of no noise and an absolute dead time (idealized filtering) are excluded by the iterative simulation method, and the correlation information contained in the two-dimensional distributions should increase the ability to distinguish among different gating mechanisms. The iterative simulation method is generally applicable to channels which typically open to a single conductance level. For these channels the method places no restrictions on the proposed gating mechanism or the form of the predicted dwell-time distributions.     

The calcium-independent transient outward potassium current in isolated ferret right ventricular myocytes. I. Basic characterization and kinetic analysis

Enzymatically isolated myocytes from ferret right ventricles (12-16 wk, male) were studied using the whole cell patch clamp technique. The macroscopic properties of a transient outward K+ current I(to) were quantified. I(to) is selective for K+, with a PNa/PK of 0.082. Activation of I(to) is a voltage-dependent process, with both activation and inactivation being independent of Na+ or Ca2+ influx. Steady-state inactivation is well described by a single Boltzmann relationship (V1/2 = -13.5 mV; k = 5.6 mV). Substantial inactivation can occur during a subthreshold depolarization without any measurable macroscopic current. Both development of and recovery from inactivation are well described by single exponential processes. Ensemble averages of single I(to) channel currents recorded in cell-attached patches reproduce macroscopic I(to) and indicate that inactivation is complete at depolarized potentials. The overall inactivation/recovery time constant curve has a bell-shaped potential dependence that peaks between -10 and -20 mV, with time constants (22 degrees C) ranging from 23 ms (-90 mV) to 304 ms (-10 mV). Steady-state activation displays a sigmoidal dependence on membrane potential, with a net aggregate half- activation potential of +22.5 mV. Activation kinetics (0 to +70 mV, 22 degrees C) are rapid, with I(to) peaking in approximately 5-15 ms at +50 mV. Experiments conducted at reduced temperatures (12 degrees C) demonstrate that activation occurs with a time delay. A nonlinear least- squares analysis indicates that three closed kinetic states are necessary and sufficient to model activation. Derived time constants of activation (22 degrees C) ranged from 10 ms (+10 mV) to 2 ms (+70 mV). Within the framework of Hodgkin-Huxley formalism, Ito gating can be described using an a3i formulation.     

New method of determining kinetic constants for immobilized enzymes in a packet-bed reactor system.

A new graphical method of determining the kinetic constants for immobilized-enzyme systems is proposed and illustrated with some examples. The advantages of using this method are that these kinetic constants can be determined accurately and conveniently from the conversion data of a packed-column enzyme-reactor system.     

Multiple kinetic states underlying macroscopic M-currents in bullfrog sympathetic neurons.

M-current is a time- and voltage-dependent potassium current which is suppressible by muscarinic receptor activation. We have used curve fitting and noise analysis to determine if macroscopic M-currents deviate from a previously predicted simple two-state kinetic scheme. The M-current was best described by three kinetically distinct components: 'fast' (tau 0), 'intermediate' (tau 1) and 'slow' (tau 2) time constants. The 'fast' (tau 0) and 'intermediate' (tau 1) components were identified from the spectra of M-current noise at potentials positive to the cells' resting membrane potential. The 'intermediate' (tau 1) and 'slow' (tau 2) components were seen by curve fitting M-current deactivation currents. The 'intermediate' (tau 1) time constant was voltage dependent (decreasing e-fold in 23 mV), but voltage dependence of the 'fast' (tau 0) and 'slow' (tau 2) components was not obvious. All kinetic components were sensitive to muscarine, with the 'intermediate' (tau 1) and 'slow' (tau 2) being equally so. These data suggest that all components may derive from the same channel population, and that the M-channel may have at least four kinetic states.     

Optimal estimation of ion-channel kinetics from macroscopic currents

Markov modeling provides an effective approach for modeling ion channel kinetics. There are several search algorithms for global fitting of macroscopic or single-channel currents across different experimental conditions. Here we present a particle swarm optimization(PSO)-based approach which, when used in combination with golden section search (GSS), can fit macroscopic voltage responses with a high degree of accuracy (errors within 1%) and reasonable amount of calculation time (less than 10 hours for 20 free parameters) on a desktop computer. We also describe a method for initial value estimation of the model parameters, which appears to favor identification of global optimum and can further reduce the computational cost. The PSO-GSS algorithm is applicable for kinetic models of arbitrary topology and size and compatible with common stimulation protocols, which provides a convenient approach for establishing kinetic models at the macroscopic level.     

Gating of the squid sodium channel at positive potentials: II. Single channels reveal two open states.

Single-channel recordings from squid axon Na+ channels were made under conditions of reverse sodium gradient. In the range of potentials studied, +40-(+)120 mV, channels opened promptly after depolarization, closed and reopened several times during the pulse. In patches containing only one channel, the distributions of open dwell times showed two components showing the existence of a second open state. The ensemble average of single-channel records showed incomplete inactivation that became more pronounced at more positive potentials, showing that the maintained phase of the current is the result of only one type of sodium channel with two open states. Analysis of bursts indicated that the dwell times of the events at the onset of the depolarization are longer than those later in the pulse. The dwell open times of the first events could be fitted with a single exponential. This indicated that the channels open preferentially through the first open state, the access to the second open state happening subsequently. Maximum likelihood analysis was used to evaluate several possible kinetic schemes incorporating a second open state. The best model to fit the data from single channels, and consistent with the data from macroscopic and gating currents, has a second open state evolving from the inactivated state. A kinetic model is proposed that incorporates information obtained from dialyzed axons.     

Quantitative simulation of endplate currents at neuromuscular junctions based on the reaction of acetylcholine with acetylcholine receptor and acetylcholinesterase.

Two kinetic models are introduced which predict amplitudes and time-courses of endplate currents and miniature endplate currents at neuromuscular junctions, at both normal and acetylcholinesterase-inhibited endplates. Appropriate differential rate equations reflecting interactions of acetylcholine with acetylcholine receptor and with esterase, diffusion of acetylcholine both within and from the synaptic cleft, and cooperativity between receptor site occupancy and ion channel opening are solved. Acetylcholine release into the cleft is assumed to be instantaneous. The simpler homogeneous reaction space model accurately predicts decay phase time constants are inaccurate. The two-reaction space model predicts amplitudes and time constants within a factor of two of those observed experimentally. The simulations indicate that the amplitudes and time-courses are primarily determined by the chemical reaction rates that characterize acetylcholine interactions with receptor and esterase and that these interactions occur under nonequilibrium conditions. Approximately 50% of the total ion channels in the initial reaction space are predicted to be opened at the peak endplate current. The cooperative opening of ion channels by acetylcholine requires that acetylcholine be introduced into the cleft in discrete, concentrated elements. Virtually all the open channels are confined to the initial reaction space, although acetylcholine-bound receptor sites can be much more widely distributed.     

Fractal analysis of a voltage-dependent potassium channel from cultured mouse hippocampal neurons.

The kinetics of ion channels have been widely modeled as a Markov process. In these models it is assumed that the channel protein has a small number of discrete conformational states and the kinetic rate constants connecting these states are constant. In the alternative fractal model the spontaneous fluctuations of the channel protein at many different time scales are represented by a kinetic rate constant k = At1-D, where A is the kinetic setpoint and D the fractal dimension. Single-channel currents were recorded at 146 mM external K+ from an inwardly rectifying, 120 pS, K+ selective, voltage-sensitive channel in cultured mouse hippocampal neurons. The kinetics of these channels were found to be statistically self-similar at different time scales as predicted by the fractal model. The fractal dimensions were approximately 2 for the closed times and approximately 1 for the open times and did not depend on voltage. For both the open and closed times the logarithm of the kinetic setpoint was found to be proportional to the applied voltage, which indicates that the gating of this channel involves the net inward movement of approximately one negative charge when this channel opens. Thus, the open and closed times and the voltage dependence of the gating of this channel are well described by the fractal model.     

The Estimation of Rapid Rate Constants from Current-Amplitude Frequency Distributions of Single-Channel Recordings

A method is described for estimating rapid rate constants from the distributions of current amplitude observed in single-channel electrical recordings. It has the advantages over previous, similar approaches that it can accommodate both multistate kinetic models and adjustable filtering of the data using an 8-pole Bessel filter. The method is conceptually straightforward: the observed distributions of current amplitude are compared with theoretical distributions derived by combining several simplifying assumptions about the underlying stochastic process with a model of the filter and electrical noise. Parameters are estimated by approximate maximum likelihood. The method was used successfully to estimate rate constants for both a simple two-state kinetic model (the transitions between open and closed states during the rapid gating of an outward-rectifying K⁺-selective channel in the plasma membrane of Acetabularia) and a complex multistate kinetic model (the blockade of the maxi cation channel in the plasma membrane of rye roots by verapamil). For the two-state model, parameters were estimated well, provided that they were not too fast or too slow in relation to the sampling rate. In the three-state model the precision of estimates depended in a complex way on the values of all rate parameters in the model. Received: 4 October 1996/Revised: 2 September 1997     

Block and inactivation of sodium channels in nerve by amino acid derivatives. II. Dependence on temperature and drug concentration.

The interaction of n-propylguanidinium (nPG) with sodium channels has been further characterized. From experiments at varying temperatures, the Q10 for the sodium current decay time constant in the two [Na+] gradients is 2.6-2.9 independent of drug. Testing several nPG concentrations we find that peak sodium current declines sharply with [nPG] at all levels, but the decay time constant approaches an asymptote above 4 mM. No "hooks" in sodium tail currents are seen. If the sodium current is allowed to decay completely before repolarization no tail current is observed. We have developed a kinetic model in which nPG acts at a single site within the sodium channel. Reaction of nPG with its receptor requires two steps. Fitting the temperature data shows that the first step involves diffusion of the drug to the site and close association with it. The second step may include molecular reorganization of the complex. The rate constants for the reaction are all simple exponential functions of voltage. Using them, the model successfully predicts decay time constants and peak currents, and their dependence on potential, [Na+] gradient, temperature, and nPG concentration. The results are consistent with the idea that an arginine residue may be closely associated with inactivation.     

Glycine-dependent activation of NMDA receptors

N-methyl-d-aspartate (NMDA) receptors are the only neurotransmitter receptors whose activation requires two distinct agonists. Heterotetramers of two GluN1 and two GluN2 subunits, NMDA receptors are broadly distributed in the central nervous system, where they mediate excitatory currents in response to synaptic glutamate release. Pore opening depends on the concurrent presence of glycine, which modulates the amplitude and time course of the glutamate-elicited response. Gating schemes for fully glutamate- and glycine-bound NMDA receptors have been described in sufficient detail to bridge the gap between microscopic and macroscopic receptor behaviors; for several receptor isoforms, these schemes include glutamate-binding steps. We examined currents recorded from cell-attached patches containing one GluN1/GluN2A receptor in the presence of several glycine-site agonists and used kinetic modeling of these data to develop reaction schemes that include explicit glycine-binding steps. Based on the ability to match a series of experimentally observed macroscopic behaviors, we propose a model for activation of the glutamate-bound NMDA receptor by glycine that predicts apparent negative agonist cooperativity and glycine-dependent desensitization in the absence of changes in microscopic binding or desensitization rate constants. These results complete the basic steps of an NMDA receptor reaction scheme for the GluN1/GluN2A isoform and prompt a reevaluation of how glycine controls NMDA receptor activation. We anticipate that our model will provide a useful quantitative instrument to further probe mechanisms and structure–function relationships of NMDA receptors and to better understand the physiological and pathological implications of endogenous fluctuations in extracellular glycine concentrations.