This Article Abstract Full Text (PDF) All Versions of this Article: biophysj.104.053256v1 88/4/2494    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Milescu, L. S. Articles by Sachs, F. Search for Related Content PubMed PubMed Citation Articles by Milescu, L. S. Articles by Sachs, F.

Maximum Likelihood Estimation of Ion Channel Kinetics from Macroscopic Currents

Lorin S. Milescu ^*, Gustav Akk  and Frederick Sachs ^*

^* Department of Physiology and Biophysics, State University of New York, Buffalo, New York; and Department of Anesthesiology, Washington University in St. Louis, St. Louis, Missouri

Correspondence: Address reprint requests to Frederick Sachs, E-mail: sachs{at}buffalo.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION THEORY AND ALGORITHM MATERIALS AND METHODS RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
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.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORY AND ALGORITHM MATERIALS AND METHODS RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
Ion channels control the flow of ions through cell membranes. These molecular gates exist in a continuum of conformations, but in practice these are reducible to a small number of kinetic states that live long enough to be resolved experimentally. Some of the states do not conduct ions (closed), whereas others are partially or fully conducting (open). "Instantaneous" transitions occur between these conformations, with the probability of transition described by a rate constant. The rate constant contains two intrinsic terms that describe the coupling of the rates to driving forces, such as ligand concentration, transmembrane voltage, tension, etc. (Hille, 2001).

The objective of kinetic analysis is to create a model—a hypothesis—with the minimum number of kinetic states and connections that can describe the data. Successful models are wonderful summaries of the data, and in that simplification help us to understand the structure and function of individual ion channels, and their mean behavior in cells. Hodgkin and Huxley provided the earliest prototype of this kind of analysis (Hodgkin and Huxley, 1952). The rate constants quantify the free energy difference between conformational states, providing insight into molecular structure and gating mechanisms. At the cellular level, the differential expression of ion channels, each with specific kinetic and conductance properties and modulation inputs, determine the static and dynamic membrane properties, and ultimately determine the network behavior of excitable cells.

Most properties of ion channels are best studied with single channel data (Neher and Sakmann, 1992), but technical difficulties often restrict recordings to samples containing multiple channels that cannot be resolved into unitary events. The statistical description of single channel data is well developed thanks to the work of many clever scientists (e.g., Ball and Rice, 1992; Ball and Sansom, 1989; Colquhoun and Hawkes, 1982, 1995a; Colquhoun and Sigworth, 1995; Horn and Lange, 1983; Kienker, 1989), and sophisticated analysis techniques based on Markov models have been published (e.g., Colquhoun et al., 1996; Qin et al., 1996, 2000). The statistical properties of macroscopic currents generated by ensembles of ion channels are well known, but the detailed kinetics are blurred by the overlap of many independent channels (Colquhoun and Hawkes, 1977). The time dependence (relaxation) of the current I_t of a model with N_S states, elicited by a step change in experimental conditions, is a sum of N_S – 1 exponentials (unless there are uncommon degenerate solutions), regardless of how states are connected (i.e., regardless of model topology; Colquhoun and Hawkes, 1995b):

(1)

The equilibrium current I_eq, the amplitudes A_i, and the time constants _i are not independent quantities because they are all functions of the same rate constants. A traditional method of macroscopic data analysis is to fit exponentials to nonstationary currents, with the main objective of estimating the time constants _i. However, the time constants thus obtained are estimated as independent parameters, without imposing the constraints determined by the particular topology of the model, or other constraints, such as loop balance. Hence, the topology of the model remains ambiguous. Rate constants cannot be easily calculated from time constants, and, without molecular rate constants, clues to molecular structure are obscure. Despite these shortcomings, exponential fitting remains the most widely used analytical method; partly for its intuitive feel, partly because of the illusion of model independence, and partly because few procedures for direct estimation of rate constants are available (e.g., d'Alcantara et al., 2002; Celentano and Hawkes, 2004).

A fundamental problem of kinetic modeling is the identifiability of the model, i.e., are there multiple topologies that are equally likely? In a now-classic example (Kienker, 1989), the model C-C-O (Closed-Closed-Open) cannot be distinguished from the model C-O-C, under equilibrium conditions. However, if one or more of the rates is stimulus dependent, they can be distinguished with nonstationary stimuli. Variable stimuli provide higher degrees of model discrimination.

We present here a maximum likelihood method for direct estimation of molecular rate constants from macroscopic ion channel data, using hidden Markov models. The algorithm can estimate the number of channels in the preparation, and the mean and standard deviation of a binary unitary current. The principle is to formulate the nonstationary macroscopic current I_t as a random variable with a Gaussian probability distribution, and to maximize the likelihood of the data with respect to the estimated parameters. For a given model, the mean µ_t and variance V_t of the macroscopic current I_t are functions of the rate constants, channel count, unitary conductance, and experimental protocol. These functions can be complicated, and involve calculation of the exponential of the rate matrix Q, as shown in Theory. The method is "direct" because the a priori assumptions about the topological and other constraints are explicitly specified in the model. We make two key assumptions: i), the ion channels are identical and independent; and ii), the experimental variables (such as ligand concentration, voltage, etc.) are not functions of channel activity. However, environmental variables need not be constant. We approximated continuously changing stimuli with a series of step functions. Although voltage clamp experiments generally satisfy these requirements, concentration and pressure jumps with significant rise time can be well approximated.

The algorithm is quite general, and has a number of features to assist users in model building: i), It works with models of arbitrary size and topology because the structural equations are internally formulated in terms of matrices and vectors. Models can be quantitatively compared, because distinct topologies with different constraints typically have different likelihoods. The modulation of individual state transitions by driving forces can be naturally studied. ii), Constraints can be imposed on rate constants. These constraints may be a priori knowledge of rate values, imposition of identical ligand-binding sites, microscopic reversibility, etc. The constraints reduce the number of free parameters, permitting exploration of larger allosteric models. iii), The current is modeled as a random variable, i.e., both the mean and the variance are used, which tends to reduce the problem of identifiability, and permits a simultaneous estimate of the mean single channel conductance. iv), The algorithm can estimate the number of active channels in the patch, if the single channel amplitude is known or coestimated. v), It works with nonstationary currents elicited in response to arbitrary experimental protocols, allowing for stimuli with finite rise times. A single model with driven rates can be globally fit across multiple data sets, with arbitrary stimulation protocols, and data from the same or from different patches can be combined in a single analysis. vi), The algorithm can analyze a sequence of stimulus steps (e.g., voltage) without assuming that channels reach equilibrium between steps. This reduces the amount of data required to solve the kinetics, and improves identifiability. vii), The starting probabilities (the state occupancies during a conditioning stimulus) are calculated from the fit kinetics and the experimental conditions. Hence, channels need not be assumed to occupy an arbitrary starting state. viii), Finally, the method efficiently calculates the likelihood function. The algorithm is faster than curve fitting using traditional differential equation solvers (d'Alcantara et al., 2002). Because the program utilizes only the diagonal elements of the autocovariance matrix, it is much faster then a recent method that uses the full autocovariance matrix (Celentano and Hawkes, 2004), and hence can handle larger data sets (although with slightly compromised precision). The gradients of the likelihood function are calculated analytically, enabling use of a fast variable metric optimizer (Fletcher and Powell, 1963; Fletcher, 1981; Press et al., 1992), with quadratic convergence.

This article focuses on the algorithm, its basic properties, and the results of extensive testing on simulated data. We concentrated on testing two ligand-dependent models (C-C-O and C-O-C), as the minimum nontrivial kinetic schemes that capture the behavior of more realistic models. As pointed out above, these simple models have identifiability issues: they are not distinguishable with equilibrium data (Kienker, 1989). However, they can be resolved by applying a suitable stimulation protocol.

We evaluated the two most important properties of the algorithm: the accuracy and the precision, and studied the effect of experimental factors, such as the number of channels in the patch, the amount of data, the stimulation protocol, etc. Because the kinetics of ion channels are rarely as simple as our three-state models, and the noise models of experimental data may depart significantly from the Gaussian assumption, we tested the algorithm with actual macroscopic currents from cells expressing acetylcholine receptors (AChRs). The AChR has a relatively complex kinetics, but the equilibrium rates are known in detail from single channel analysis (e.g., Akk and Auerbach, 1999). We showed that our method can choose the "correct" model size and topology based on differences in the goodness of fit. Because some of the time constants of AChR kinetics are quite short relative to the rise time of concentration jumps practically obtained, we explored the effect of finite rise time. In future work, we will apply the method to voltage-gated ion channels. The algorithm was implemented in software and integrated in the QuB electrophysiology suite, available for free at www.qub.buffalo.edu.

	THEORY AND ALGORITHM

TOP ABSTRACT INTRODUCTION THEORY AND ALGORITHM MATERIALS AND METHODS RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
Kinetic model
We review here some basic results of ion channel kinetics, for an arbitrary Markov model with N_S conformational states. Rate constants are conventionally expressed as an N_S x N_S rate matrix Q (also termed the "generator" matrix; Colquhoun and Hawkes, 1995b). The Q matrix has the properties that each off-diagonal element q_ij is equal to the rate constant of the transition between states i j, and that each diagonal element q_ii is equal to the negative sum of the off-diagonal elements of row i, so that the sum of each row is zero. When there is no direct transition between states i j, q_ij = 0. In Eyring-type barrier models, the rate constant q_ij includes a preexponential factor multiplied by the concentration C_ij of some ligand, and an exponential factor multiplied by the magnitude of the driving force V_ij (such as voltage, pressure, etc.; Hille, 2001):

(2)

The expression in Eq. 2 is equivalent to the common formulation where

At any given time, a state probability vector P, of dimension N_S, represents the occupancy of the N_S conformational states. The Kolmogorov differential equation describes the dynamics of a Markov process (Elliott et al., 1995):

(3)

The superscript ^T denotes vector transposition. If P₀ is the state probability vector at time t₀ < t, then P_t is given by the Chapman-Kolmogorov solution:

(4)

where is an N_S x N_S matrix, known as the transition probability matrix. Each element a_ij is equal to the conditional probability that the channel is in state j at time t, given that it was in state i at time 0. If the eigenvalues of the Q matrix are all distinct (and most ion channel models have distinct eigenvalues, with one being zero and representing the equilibrium solution, and the others negative), the matrix exponential A_t can be efficiently obtained from Q by spectral decomposition (Moler and Van Loan, 1978) as The A_i values are the spectral matrices obtained from the eigenvectors of Q, and the _i values are the eigenvalues of Q.

Conductance model
To describe the conductance properties of an ensemble of N_C channels, we make the following simplifying assumptions: i), the baseline current (i.e., the current when all channels are closed) is a -correlated Gaussian random variable, with mean µ⁰ and variance V⁰; ii), the current passing through a single channel in state i is a -correlated random Gaussian variable, with mean µ⁰ + µ_i and variance V⁰ + V_i (Sigworth, 1985). By definition, µ_i = 0 and V_i = 0 for all the closed states. For two channels, one in state i and one in state j, the current is a Gaussian with mean µ⁰ + µ_i + µ_j and variance V⁰ + V_i + V_j. In general, the current given by an ensemble of N_C channels is the sum of Gaussian random variables, which is another Gaussian random variable. Hence, the conditional probability distribution of I_t, for a given probability vector P_t, has the following expression:

(5)

(6,7)

is the number of channels occupying state i at time t, and can be formulated as a function of state probabilities: where is the occupancy probability of state i at time t. Let µ and V be vectors of dimension N_S, with elements µ_i and V_i as defined. Equations 6 and 7 can now be written in a more compact, vectorial form, as follows:

(8,9)

Equations 8 and 9 give, respectively, the conditional mean and variance of the current I_t recorded from an ensemble of N_C channels, given a state probability vector P_t. To obtain the unconditional distribution of I_t, we must remember that P_t, as calculated with Eq. 4, is the mean of a probability distribution, and that state occupancies will have fluctuations around this mean. Hence, the recorded current has two sources of variance: i), baseline and intrinsic state noise, and ii), state fluctuations. The state fluctuations theoretically have a multinomial distribution (the extension of the binomial to N_S states), but we will approximate it with a multivariate Gaussian to avoid a complicated convolution of probability distributions. This approximation will introduce a bias for currents arising from small numbers of channels. The covariance matrix of state fluctuations for the ensemble of N_C channels, at time t, is defined as:

(10)

where is the N_S x N_S state covariance matrix, with elements and We obtain the unconditional distribution of I_t as follows:

(11)

(12,13)

The vector product represents the linear transformation of the state covariance matrix from the N_S x N_S state space into the unidimensional measurement space. To speed the computation, traces recorded in response to the same stimulation protocol can be averaged (assuming stationary kinetics), but to properly exploit the information contained in the sweep-to-sweep fluctuations, the expected variance V_t needs to be divided by the number of averaged traces. The above equations are derived for an arbitrary number of conductance classes. Hence, at least in principle, the method is suitable for models with substates (partially conducting states).

Starting state probabilities
The initial probability vector P₀ can be specified in two ways. One is to let the user assign values from a priori assumptions. This is commonly done with ligand-gated channels, when zero ligand concentration drives channels into an unliganded state, which is assigned a unity starting probability. With voltage gated channels, a conditioning potential is assumed to populate a specific state. However, this approach is not suitable for models with aggregated initial states (e.g., multiple unliganded states). The solution is to calculate the initial state probabilities as the equilibrium probabilities determined by the initial (conditioning) experimental conditions, which are assumed to last long enough for channels to reach equilibrium. The equilibrium assumption requires that P₀ satisfies the condition where is the rate matrix for which the equilibrium conditions hold. One way to calculate P₀ is the following (Colquhoun and Hawkes, 1995b):

(14)

where R is a N_S x (N_S + 1) matrix, formed by adding a column of ones to the right-hand side of and u is a column vector of dimension N_S, with each element u_i = 1.

Likelihood function
We describe here the likelihood function that is to be maximized by the algorithm. The likelihood L is defined as the conditional probability of the macroscopic current trace I, recorded for time t = [0...T], given the model parameters :

(15)

We assume that the state probability vector P_t is a function of the initial state probabilities P₀, but does not depend on the state probabilities at any other time. Hence, the distribution of the current at any time t is a function of the initial state probabilities only, which is equivalent to saying that the data have no autocorrelation beyond that contained in the A matrix. Thus, the conditional probability of the entire current trace is equal to the product of the individual conditional probabilities:

(16)

where with the mean µ_t and the variance V_t as defined in Eqs. 12 and 13. The objective is to find the parameter set ^ML that maximizes the likelihood L:

(17)

The likelihood may be a very small, or a very large number, because it involves Gaussians, so to reduce numerical errors we actually maximize the logarithm of the likelihood function:

(18)

The parameters to be estimated are and (Eq. 2), and the channel count N_C (in the following, for the sake of clarity, we omit the derivations for conductance parameters). There is no mathematical constraint on the exponential factor but the preexponential factor and the channel count N_C must be positive numbers, leading to the following useful reparameterization:

(19,20)

The exponential formulation of the rates provides an implicit constraint that the rate constants are all positive, and hence physically realizable.

Parameter constraints
In addition to positivity of the preexponential factors, models may require other constraints. For example, if the two binding sites of the AChR are assumed to be identical and noninteracting, then the rates of the two consecutive binding steps should be constrained in a 2:1 ratio. Without such a constraint, the rates would be different; on the other hand, it is possible to test for independent sites using the likelihood. A priori constraints simplify otherwise complicated models. Each constraint, in addition to incorporating previously known results, reduces the number of free parameters by one. Typical constraints are: i), fix the pre- or exponential factor; ii), scale one pre- or exponential factor proportional to another; iii), enforce microscopic reversibility of cycles (the product of the rate constants going in one direction of a cyclic reaction is equal to the product of the rate constants going the other way). These constraints can be formulated as follows, where ct denotes the constrained value:

(21)

(22)

(23)

(24)

(25)

The constraint expressions above can be written in matrix form as a set of linear equations (Qin et al., 1996, 2000):

(26)

where r is the vector of dependent (constrained) variables, with entries for and u. M_c is the coefficient matrix and V_c is the constant vector. The dependent variables r can in turn be expressed as a function of the independent variables x (Qin et al., 2000):

(27)

where x is the vector of free parameters. The matrix A_c and the vector B_c can be determined from M_c and V_c using, e.g., the singular value decomposition (Qin et al., 2000). Other, nonlinear constraints could be included as well, by adding penalty terms to the likelihood function.

Computation of the likelihood function and its derivatives
We describe next the details of the computation of the likelihood function and its derivatives. The log likelihood LL of a macroscopic current, recorded for time t = [0...T] under constant experimental conditions, is calculated as follows:

(28)

(29)

The log-likelihood function in Eq. 29 reduces to the equally weighted (i.e., nonweighted) sum of squares if the constant terms and the variance-containing terms are eliminated.

The quantities of interest are the pre- and exponential factors and and the channel count N_C, but the optimizer maximizes LL with respect to the free parameters x. and N_C are subsequently obtained from x using Eqs. 19, 20, and 27. The derivative of LL (or SS) with respect to a free parameter x_k can be calculated with the chain rule as follows:

(30)

with and where i and j are indices over the off-diagonal entries in the Q matrix corresponding to nonzero rates in the kinetic model. The partial derivatives of and u with respect to x_k are obtained with Eq. 27. Let z denote either q_ij or N_C. Then the partial derivative of LL with respect to z is calculated as follows:

(31)

(32)

Similarly, The derivatives of µ_t and V_t with respect to z are:

(33)

(34)

The following results can be immediately derived: The derivative of the matrix with respect to q_ij (and similarly for ) is another matrix with elements:

(35)

The derivatives of the state probability vector P_t with respect to q_ij and are:

(36)

The following results can be immediately derived: The derivative of the initial state probability vector P₀ can be obtained using the following equality:

(37)

We skip the derivation details and give the final expression:

(38)

The derivative of the A_t matrix (Ball and Sansom, 1989; Najfeld and Havel, 1995) is:

(39)

(40)

The derivative of Q (and R) with respect to q_ij (and ) is another matrix Q_ij', with the properties

Extension to arbitrary stimulation protocols
All the above calculations assume that the experimental conditions (i.e., the stimuli) do not change across the data set. We show next how to extend the method to arbitrary stimulation protocols. We assume that the applied stimulus is approximated by a sequence of S_C steps of constant value and arbitrary duration. Each step s starts and ends at times and respectively. We rewrite the expressions of the log likelihood and its derivative, to show the dependence of µ_t and V_t on the step s, as follows:

(41)

(42)

In this case, the log likelihood is a function of multiple Q matrices (in addition to Q^eq), one for each set of experimental conditions that change the ion channel response. Although for each value of the stimulus there is a different Q matrix, the set of free parameters remains the same, as it does not depend on experimental conditions. Accordingly, the derivative of the log-likelihood function LL with respect to a free parameter x_k has the following form:

(43)

where c is an index over all different experimental conditions, including the initial equilibrium. Let us consider the simpler case when the stimulus changes the kinetics, but not the unitary channel current, e.g., a change of ligand concentration. The expressions of µ_t,s and V_t,s become:

(44,45)

For an efficient computation of the state probability vector P_t,s and its derivatives, we derive recursive relations. The state vector at time t, contained in the step s, is calculated as follows:

(46)

where By definition, the vector at the start of step s is equal to the vector at the end of step s – 1:

(47)

where Equation 47 can be written recursively and initiated as follows:

(48-50)

The derivative of the state vector can also be calculated recursively and initiated as follows:

(51)

(52,53)

The derivative of A_t,s with respect to is equal to zero when the experimental conditions for step s have an index different than c, thus simplifying the calculations. A similar treatment can be derived for voltage stimuli that change both the kinetics and the unitary current.

Variance of the estimates
Our choice of optimization routine is the Davidon-Fletcher-Powell algorithm (Fletcher and Powell, 1963; Fletcher, 1981), which is a quasi-Newton method with approximate line search, with efficient quadratic convergence. We used the "dfpmin" implementation in Press et al. (1992), with some modifications (Qin et al., 2000). The optimizer provides variance estimates for the optimized, free parameters x. The variance of the pre- and exponential factors and and of the channel count N_C can be derived as follows (Qin et al., 2000):

(54)

(55)

(56)

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION THEORY AND ALGORITHM MATERIALS AND METHODS RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
Electrophysiology
The adult mouse muscle-type nicotinic AChRs (-, ß-, -, and -subunits) were transiently expressed in tsA 201 cells using a calcium phosphate-based transfection technique, as described (Akk, 2002). The cells were plated on coverslips in 35-mm dishes. At the time of experiments, one coverslip at a time was transferred into the recording chamber. The electrophysiological recordings were carried out using a piezo-driven liquid filament system to measure the macroscopic responses to 1 mM carbachol on outside-out patches (Heckmann and Pawlu, 2002). The system had a solution exchange time of <50 µs when measured with changes in the liquid junction potential with an open patch electrode, or <100 µs when measured with an electrode containing a patch. The bath solution contained (in mM): 140 NaCl, 5 KCl, 2 CaCl₂, 1 MgCl₂, 10 glucose, and 10 Hepes (pH 7.3). The pipette solution contained (in mM): 140 KCl, 5 MgCl₂, 5 EGTA, 10 glucose, and 10 Hepes (pH 7.3). The agonist was dissolved in the bath saline solution and the membrane potential was held at –40 mV. The experiments were carried out at room temperature. Currents were amplified with an Axopatch 200B amplifier (Axon Instruments, Union City, CA), low-pass filtered at 10 kHz, and digitized at 50 kHz. Traces of 30-ms duration were recorded at intervals of 1 s using the ISO2 software (MFK, Taunusstein, Germany). Averaged currents from 10 to 25 traces were used for kinetic analysis.

Simulation and data analysis
Stochastic and deterministic simulations were generated with the QuB software (www.qub.buffalo.edu), using the ligand-gated kinetic models from Fig. 1, referred to as C-O-C and C-C-O, respectively, and the three stimulation protocols from Fig. 2, referred to as P₁, P₂, and P₃, respectively. We call the parameter values used for simulation true parameters, as opposed to estimated parameters. For each simulation, 10,000 traces were generated. Simulations with other models and protocols are described in the Results. A sampling interval of 10 µs was used for simulations, except in the analysis of experimental data, for which we used 20 µs. The QuB program can generate a stochastic macroscopic simulation as the summation of N_C single-channel stochastic simulations, in conformity with Eqs. 6 and 7. Alternatively, it can make a deterministic macroscopic simulation by calculating the current according to Eq. 8. The starting state probabilities are calculated as the equilibrium probabilities under the conditioning stimuli. For stochastic data, the starting state for each of the N_C channels is chosen randomly, such that the probability of starting in state i is proportional to The estimation of rate constants and channel count, and the nonweighted exponential fitting were done with QuB.

View larger version (21K): [in this window] [in a new window]   FIGURE 1  The models used for simulation. (A) The transitions C1O2 (in model COC) and C1C2 (in model CCO) are ligand-binding steps. (B) An example of single-channel data simulated with the COC model. The single-channel current in either the Closed or the Open state is a random Gaussian variable.

View larger version (19K): [in this window] [in a new window]   FIGURE 2  The macroscopic response to different stimulation protocols. Macroscopic traces were simulated deterministically (smooth red lines) and stochastically (noisy traces) with the CCO and COC models, with NC = 1000. The starting equilibrium probabilities were (1, 0, 0), for protocol P1, and (0.571, 0.143, 0.286), for protocols P2 and P3. The units of the ligand concentration are arbitrary.

	RESULTS

TOP ABSTRACT INTRODUCTION THEORY AND ALGORITHM MATERIALS AND METHODS RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

With respect to the first question, Kienker (1989) showed that two models of different topology cannot be distinguished under equilibrium conditions if their Q matrices are related by a similarity transform. For example, the C-C-O and C-O-C models from Fig. 1 can describe the same equilibrium data with identical likelihood. If some three-state models cannot be distinguished, then we should be wary of more complicated models. However, the use of driving stimuli greatly enhances the ability to differentiate among models. For instance, when started in the leftmost state, the C-C-O model has a second-order latency, whereas the C-O-C has a first-order latency.

Can these two models be distinguished from macroscopic nonstationary data? The answer is yes. We fit the deterministic simulation of the C-O-C model in response to the P₁ stimulation protocol (Fig. 2 A), with the C-C-O model, and vice versa. The results are not shown, but we could not find any set of parameters that could generate the same current as the other model. The two models are therefore mathematically distinguishable with the simplest possible protocol, but how well can they be differentiated from noisy (possibly nonideal) data is a different issue.

With respect to the second question, a model can be classified (Cobelli and DiStefano, 1980) as a priori: i), "uniquely identifiable" if all its parameters have a unique solution; ii), "nonuniquely identifiable" if at least one of its parameters has a multiple but finite number of solutions, and all other parameters have unique solutions; iii) "nonidentifiable" if at least one of its parameters is undetermined, i.e., it has an infinite number of solutions.

We examined the identifiability of the C-C-O and C-O-C models. Macroscopic currents were deterministically simulated in response to the three driving stimuli from Fig. 2, with the rate constants shown in Fig. 1, and the number of channels N_C = 1000. We searched for other values of k₁₂, k₂₁, k₂₃, k₃₂, and N_C that gave the same average current as the true parameters (i.e., other solutions) by fixing one or more parameters and having the algorithm find the others, restarting several times with different initial values. Because the simulated data were noiseless, we minimized the sum of squares. Only the solutions with zero SS (within numerical precision) were accepted.

The results obtained with model C-C-O and protocol P₁ are shown in Fig. 3 A.1. The figure shows that, for a given N_C, for each value of k₂₁ there were two solutions for the other three rates. Hence, with protocol P₁, only two parameters can be uniquely estimated from model C-C-O. When repeated with different values of N_C, the experiment yielded similar results. Interestingly, the admissible domain of k₂₁ had an upper bound (341 s^–1). With protocol P₂, N_C was varied, and the four rate constants were again estimated (Fig. 3 A.2). For each value of N_C, there were again two solutions for the four rates. Hence, protocol P₂ is more informative, yielding three free parameters. Not surprisingly, the admissible domain of N_C had a lower bound (740), determined by the peak current amplitude and the unitary current.

View larger version (17K): [in this window] [in a new window]   FIGURE 3  The model identifiability depends on model topology and stimulation protocol. Macroscopic currents of identical mean (but not variance) can be obtained with different sets of the five parameters (four rate constants and NC). Macroscopic traces were simulated deterministically, as in Fig. 2. The solid and the dashed lines in panels A.1, A.2, and B.1 represent two possible solutions of the estimated parameters, for each value of the varied parameter (k21 and NC). In panel A.1, NC was fixed at its true value of 1000 (varying NC adds an extra dimension to the solution space). In panel B.2, solutions of the four rate constants (solid circles) exist only for two values of NC (1000 and 3333.333). With protocol P1, a maximum of two and three parameters could be uniquely estimated for the CCO (A.1) and the COC (B.1) models, respectively. Protocol P2 increased by one the maximum number of uniquely identifiable parameters (A.2 and B.2), whereas with protocol P3 all five parameters were uniquely identifiable (results not shown).

To summarize, with protocols P₁ and P₂, two and three free parameters can be uniquely estimated for model C-C-O, respectively, and three and four for model C-O-C. With protocol P₃, both models become fully identifiable, with unique solutions for all five parameters (results not shown). In conclusion, model C-C-O is nonidentifiable with protocols P₁ and P₂, and uniquely identifiable with P₃. Model C-O-C is nonidentifiable with P₁, nonuniquely identifiable with P₂, and uniquely identifiable with P₃. Although our search approach was not designed to find all admissible solutions rigorously, it made it clear that even these simple models create identifiability issues, even with nonstationary stimuli. The identifiability can be improved by more complex stimulation protocols. C-C-O is more ambiguous than C-O-C for the same stimulation.

Next, we tested whether including the variance-dependent terms in the likelihood function improves identifiability. Macroscopic currents were stochastically simulated with the C-O-C model, in response to stimulation protocols P₁ and P₂, respectively. The number of channels N_C was fixed at different values and the four rate constants estimated using the nonweighted sum of squares or the log likelihood. The results shown in Fig. 4, A.1 and B.1, essentially restate the conclusion drawn from Fig. 3, i.e., that the parameters of the C-O-C model are not uniquely identifiable with either the P₁ or P₂ protocols. With P₁, the SS corresponding to the best fit is constant for all values of N_C (Fig. 4 A.1), or has two identical minima with P₂ (Fig. 4 B.1). In contrast, the LL function has a global maximum, centered at the true parameter values, with either protocol P₁ (Fig. 4 A.2) or P₂ (Fig. 4 B.2). In conclusion, analyzing the response to more complex stimulation protocols and adding the variance-containing terms to the mathematical model are different ways to increase the number of equations, and hence to improve identifiability.

View larger version (18K): [in this window] [in a new window]   FIGURE 4  The model identifiability improves if both the mean and the variance of the macroscopic current are utilized for parameter estimation. Macroscopic currents were stochastically simulated with the COC model, with NC = 1000, in response to protocols P1 and P2 (as in Fig. 2). NC was varied and the four rate constants were estimated from 100 simulated traces, using either the sum of squares or the log-likelihood method. The SS method uses only the mean, whereas the LL method uses both the mean and the variance. For easier comparison between LL and SS, the graphics show the negative SS.

View larger version (20K): [in this window] [in a new window]   FIGURE 5  Conductance parameters can be estimated from macroscopic currents with the likelihood function. Macroscopic currents were stochastically simulated with the COC model, with NC = 1000, in response to protocol P1 (as in Fig. 2). The mean µO and the standard deviation O of the single-channel Open current were varied and the four rate constants and NC were estimated with the LL method.

We started with a qualitative look at the effects of biological variability. The signal/noise ratio (SNR) of a macroscopic current is proportional to and to as follows from Eqs. 12 and 13:

(57)

where N_T is the number of averaged traces. The deviation between the fit and the current calculated using the true parameters is inversely proportional to the SNR (Fig. 6). Probabilistically, at low SNR the fit may differ substantially from the deterministically simulated current, but it matches well both the deterministic and the stochastic simulations at high SNR.

View larger version (31K): [in this window] [in a new window]   FIGURE 6  The fit of the stochastic simulation may differ substantially from the deterministic simulation, at low signal/noise ratio. Macroscopic currents were simulated stochastically (noisy traces) and deterministically (red lines) with the COC model, with NC = 10, 100, or 1000, in response to protocol P1. The stochastic simulations were fit (green lines) as individual traces or as averages of 10 or 100 traces. The fits are virtually the same for the LL (shown), SS, or exponential fitting.

View larger version (22K): [in this window] [in a new window]   FIGURE 7  The accuracy and precision of the estimates are functions of the signal/noise ratio. Macroscopic traces were simulated stochastically with the COC model, with NC = 10, 100, or 1000, in response to protocol P1. Rate constants (A), or rate constants and channel count (B) were obtained with the LL method, and time constants were obtained with exponential fitting (C). The average and the SE of the estimated parameters were calculated from 1000 fits of single traces or averages of multiple traces. The horizontal dotted lines in the average graphs mark the true parameter values. For easier comparison, the 10% SE is marked with a line in the standard error graphs.

View larger version (28K): [in this window] [in a new window]   FIGURE 8  The probability distribution of parameter estimates is not Gaussian. Macroscopic traces with NC = 10, and single-channel traces were simulated stochastically with the COC model, in response to protocol P1. Macroscopic traces were fit individually or as averages of 10, and single-channel traces were analyzed individually or in groups of 10. (A) The distribution of rate constants and NC estimated from macroscopic traces with the LL method. Similar results (not shown) were obtained with NC fixed. (B) The distributions of rate constants estimated from single-channel traces, using the QuB software (www.qub.buffalo.edu). (C) The distribution of time constants estimated with exponential fitting. In all graphs, the dotted vertical lines mark the true parameter values.

View larger version (22K): [in this window] [in a new window]   FIGURE 9  The probability distribution of macroscopic data can be well approximated with a Gaussian, except for very low (or very high) pOpen and small NC. Ten-thousand macroscopic traces were simulated stochastically with the COC model, with NC = 10, 100, or 1000, in response to protocol P1. (A) Probability distribution of the macroscopic current, calculated from the 10,000 traces, for pOpen = 0.1 and 0.5. Overlapped is the fit with a Gaussian (red lines). (B) Standard deviation of the macroscopic current, measured from the 10,000 traces. Overlapped is the calculated standard deviation (red lines).

View larger version (25K): [in this window] [in a new window]   FIGURE 10  The bias of the estimates is caused by autocorrelation. Ten-thousand macroscopic traces were simulated stochastically with the COC model, with NC = 10, in response to protocol P1. The traces were further processed, by randomly reshuffling the data across traces, for each data point. The processed traces had the same current probability distribution, but lacked the autocorrelation of the raw data. The rate constants and the channel count were estimated with the LL method, from individual traces or from averages of multiple traces. (A) Examples of a simulated raw trace (top figure) and of a processed trace (bottom figure). Overlapped is the fit (green line), and the deterministic simulation (red line). (B) The average and the SE of the estimates obtained from the processed traces. (C) The distribution of the estimates. The dotted vertical lines mark the true parameter values.

View larger version (27K): [in this window] [in a new window]   FIGURE 11  Three different estimation methods indicate substantial data overfitting. Macroscopic traces were simulated stochastically with the COC model, with NC = 10, in response to protocol P1, and were fit individually or as averages of 10. Rate constants were estimated with the LL and the SS methods, and time constants were estimated with exponential fitting (ExpFit). The sums of squares between the data and the fit line (termed SSLL, SSSS, and SSExpFit, respectively), and between the data and the deterministic simulation (termed SSModel) were calculated. (A.1–A.3) The distribution of the SS for the three fitting methods. Overlapped is the distribution of SSModel (red lines). (B.1–B.3) The correlation between the SS of the three fitting methods and the SSModel. The graphics indicate that, for each fit, the residuals obtained with the best fit parameters are considerably smaller than the residuals obtained with the true parameters. The estimated parameters are therefore biased.

View larger version (31K): [in this window] [in a new window]   FIGURE 12  The parameter estimates are correlated. Macroscopic traces were simulated stochastically with the COC model, with NC = 10 or 1000, in response to protocols P1 or P3, and were fit individually or as averages of 10. Rate constants were estimated with the LL method. The red solid circles mark the true parameter values. Each graph is a two-dimensional histogram of the estimates for a given pair of two parameters, under a given set of conditions (NC, protocol, and number of averaged traces). There was strong correlation between, e.g., the k32 and k23 rate constants (row C or D), but little, if any, between, e.g., k23 and k12 (row C). Similar correlations were found between rates and NC (not shown). The correlations were not eliminated by a higher signal/noise ratio (compare rows B and C with row A), nor by a more complex stimulation protocol (compare row D with row C). Notice the change in scale between the first two and the last two rows of figures.

View larger version (18K): [in this window] [in a new window]   FIGURE 13  The accuracy and precision of the estimates are functions of model identifiability. Macroscopic traces were simulated stochastically with the COC model, with NC = 1000, in response to protocols P1, P2, and P3 (A, B, and C, respectively). Rate constants and NC (left-hand side figures of each panel), or rate constants, NC and the mean of the unitary current µO were estimated with the LL method. The average and the SE of the estimated parameters were calculated as in Fig. 7.

View larger version (22K): [in this window] [in a new window]   FIGURE 14  The kinetics of the acetylcholine receptor can be determined from the macroscopic response to brief pulses of 1 mM carbachol. The experimental procedure is described in Materials and Methods. (A) The average of 10 macroscopic traces. The green and the red lines are the fits with the CCCO and CCCOD models, respectively (C). The rate constants and NC were estimated from the average trace, using the SS method. For the single-channel amplitude we used the value found in single-channel experiments. The blue line is the fit with the CCCOD model, in which all rate