INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORY SIMULATION STUDIES WHAT IS A SIGNIFICANT... DISCUSSION REFERENCES

The recording of single ion channel currents by the so-called patch clamp technique has deepened the understanding of the fundamental physiological mechanisms of a cell in the last two decades (Neher and Sakmann, 1976; Hamill et al., 1981; Sakmann and Neher, 1995). In steady state, ion channels permanently perform transitions among a number of unobserved states which are divided into two groups: the open and closed states. (Some channels posses sub-conductance levels, i.e., open states with different conductivities. All results in this paper generalize for this case.) In an open and in a closed state, the measured currents fluctuate around a certain conductance level which is the same for all open and all closed states, respectively. Thus, the observed ion current provides an aggregated image of the underlying process.

Aggregated Markov models are a suitable model class to describe the dynamics of the measured currents (Colquhoun and Hawkes, 1977, 1982; Colquhoun and Sigworth, 1995; Fredkin et al., 1983).

In thermodynamic equilibrium, the dynamics of the gating are subject to the principle of detailed balance. Therefore, a violation of detailed balance would indicate the presence of an external energy source. It has been observed for Cl channels from Torpedo electroplax that a transmembrane electrochemical gradient can cause such a non-equilibrium behavior of the Cl ion-channels (Richard and Miller, 1990).

Often, aggregated Markov models used to model realistic ion channels have one or more loops (Horn and Lange, 1983; Ball and Sansom, 1989; Bates et al., 1990; Vandenberg and Bezanilla, 1991). In these models the principle of detailed balance is equivalent to the condition that the products of the transition rates in clockwise and counterclockwise direction in each loop are the same (Kelly, 1979; Kijima, 1997).

In the following, we investigate the power of likelihood ratio tests to detect violations of detailed balance in two different gating schemes each with one loop: the simplest four-state model with one loop (see Fig. 1), which we refer to as the "loop model", and a six-state model in which the open and closed state alternate and in which the states are arranged in a circle (see Fig. 2). This model was used by Song and Magleby (1994) to investigate their method to test for microscopic reversibility based on visual inspection of estimated two dimensional dwell time histograms. Below we will refer to this model as the SM model.

View larger version (11K): [in this window] [in a new window]   FIGURE 1   Loop gating scheme: Identifiable gating scheme with one loop and four states. O, open state; C, closed state.

View larger version (8K): [in this window] [in a new window]   FIGURE 2   Gating scheme 4 from Song and Magleby (1994).

Both models serve as simplified examples for investigating the difficulties in testing for detailed balance in more realistic, but also more complicated aggregated Markov models. This is not a severe restriction because any more complicated aggregated Markov model with loops contains these simple models as submodels and, therefore, the larger model will exhibit an analogous behavior.

The loop model and the SM model have been chosen for this paper because they differ in the identifiability of their parameters. Whereas the transition rates in the SM model are always identifiable, the transition rates in the loop model loop are not identifiable if the open or the closed dwell times are equal, that is, the time constants characterizing the exponentially distributed dwell times in the open states or in the closed states are the same. Moreover, this non-identifiability severely influences the estimation of the transition rates in the case of nearly equal dwell times; the standard deviations of the estimated transition rates become extraordinarily large even for the number of data points typically recorded in experiments (Wagner et al., 1999). In the present paper we will show that this non-identifiability causes a drop in the power of tests for detailed balance if the open or the closed dwell times are almost equal.

The following section describes the theory of the non-identifiability for equal open times. Then we demonstrate in the Simulation Studies Section the effect of the non-identifiability for equal open times on the power of a likelihood ratio test by a simulation study. The number of data points is chosen to be in the order of magnitude of the number recorded in experiments.

	THEORY

TOP ABSTRACT INTRODUCTION THEORY SIMULATION STUDIES WHAT IS A SIGNIFICANT... DISCUSSION REFERENCES

Parameter estimation

The aggregated Markov models used to describe ion channel gating are parametrized by transition rates and not by transition probabilities. The estimation of transition rates in aggregated Markov models is burdened by the possibility that the transition rates might not be identifiable (Kienker, 1989). For instance, the maximum number of identifiable parameters in aggregated Markov models with two output levels for open and closed states is restricted to two times the number of open states times the number of closed states (Fredkin et al., 1983; Fredkin and Rice, 1986). Moreover, in aggregated Markov models, there is not only a maximum number of parameters that can be estimated from the recorded data but, under certain circumstances, the transition rates might not be identifiable, although the upper limit of parameters is not exceeded. For example, the transition rates in the loop model with equal open dwell times are not identifiable because for a given loop model with equal open dwell times, certain arbitrary small variations can be performed on the transition rates without changing the probability distribution of the observed data (Wagner et al., 1999).

The restriction of equal open times would rarely be put on analyses of experimental ion channel data. However, the non-identifiability for equal open times in the loop model has an effect on the estimation of the transition rates. The transition rates in aggregated Markov models are typically estimated by the maximum likelihood method (Horn and Lange, 1983; Ball and Sansom, 1989; Bates et al., 1990; Fredkin and Rice, 1992; Colquhoun et al., 1996; Qin et al., 1996, 1997; Ohno et al., 1996). For a large sample size, the covariance matrix of the estimated transition rates is given by the inverse Hessian matrix of the likelihood function evaluated at the maximum likelihood parameters (Bickel et al., 1998). Assuming that the maximum likelihood parameters correspond by chance to equal open times, then the Hessian matrix would be singular due to the non-identifiability for equal open times in the loop model, its condition number would be infinity (Golub and VanLoan, 1996). Since the condition number of a matrix depends continuously on its entries, the Hessian matrix is ill-conditioned if the maximum likelihood parameters correspond to nearly equal open times. Therefore, if the true model has almost equal open times and the maximum likelihood estimators converge to the true parameters, the estimated standard deviations of the parameters might become unexpectedly large even for a large amount of data (Wagner et al., 1999).

The transition rates in the SM model are always identifiable even if the open times are equal. The reason for this difference between the loop model and the SM model is given in the next section, in which a criterion for the non-identifiability will be derived.

Criterion for non-identifiability

Indistinguishable aggregated Markov models for ion channel gating are related by a so-called similarity transformation, i.e. the generator matrices Q and Q', respectively, of both models obey Q' = S¹ QS where the matrix S is of the form

(1)

and each row of S sums to 1 (Kienker, 1989). Let n_O and n_C denote the number of open and the number of closed states, respectively. The matrix S in general can be parameterized with n_O(n_O  1) + n_C(n_C  1) parameters because of the row normalization, e.g., for n_O = n_C = 2, the transformation matrix S has 4 degrees of freedom and a possible choice parameterization of S would be

(2)

We combine these parameters to a vector denoted by , for the example of Eq. 2 is (₁^O, ₂^O, ₁^C, ₂^C). In the following, we will investigate similarity transformations that are near to the identity transformation. Therefore, we adopt the convention that S() is the identity transformation for  = 0. Under this assumption, the condition that S needs to be invertible does not impose any further constraints on the maximal number of parameters of the transformation matrix S if the magnitude of is sufficiently small.

A given aggregated Markov model is described by its generator matrix Q. The gating scheme of this model is determined by the vanishing entries in the generator matrix. A similarity transformation leads to a new aggregated Markov model which is both statistically indistinguishable from the original one and compatible with the gating scheme of the original model if the vanishing entries of the original generator matrix Q are preserved in the transformed generator matrix Q' and the off-diagonal non-zero entries of Q' are not negative. For similarity transformations sufficiently near to the identity transformation, the requirement of non-negative entries is always fulfilled due to the continuity of the transformation.

Under these conditions, the constraint, to be compatible with the gating scheme of the original model, imposes m restrictions on the similarity matrix S where m is lower than or equal to the number of vanishing entries in the original generator matrix Q. For example, in the loop gating scheme with equal open times the submatrix Q_OO is proportional to the 2 × 2-identity matrix. Consequently, the requirement that the off-diagonal elements of Q_OO vanish is fulfilled by any similarity transformation. We combine the m entries in the transformed generator matrix Q' = S()¹QS() which are required to vanish, to a vector-valued function () of dimension m, in particular F() = . In the loop model with equal open times, a direct transition from open state O₂ to closed state C₁ and vice versa is not allowed, i.e., the entries q₂₃ and q₃₂ of the generator matrix Q must also vanish in the transformed generator matrix Q', and the function () is given by () = (q'₂₃(), q'₃₂()) in the example.

If the similarity transformation S() preserves the gating scheme, the function () satisfies () =  for all sufficiently small .

If m is strictly smaller than n_O(n_O  1) + n_C(n_C  1), the theorem of implicit functions can be applied to (). We split the parameter vector into a (n_O(n_O  1) + n_C(n_C  1)  m)-dimensional part and a m-dimensional part : () = (,). Under mild regularity conditions there exists a vector-valued function () of dimension m with (,()) = 0 for of sufficiently small magnitude. Under these conditions the continuous family of transformation matrices S() = S( = (,())) yields a continuous family of aggregated Markov processes which are both statistically indistinguishable from the original one and compatible with the gating scheme of the original model.

Thus, the transition rates in a given aggregated Markov model are not identifiable if the minimal number of restrictions m needed to preserve its gating scheme is smaller than the degrees of freedom n_O(n_O  1) + n_C(n_C  1) of the similarity transformation.

Theorem B of Fredkin and Rice (1986) derives an upper bound for the maximum number of identifiable parameters. In contrast to this result, we provide a checkable sufficient condition for non-identifiability which is applicable in cases where this upper bound is not exceeded, e.g., in the loop model with equal open times.

For equal open times in the SM model, the submatrix Q_OO is proportional to the 3 × 3 identity matrix, but the submatrices Q_OC, Q_CO, and Q_CC still contain 12 non-vanishing entries, which is equal to the number of degrees of freedom of a similarity transformation for a gating scheme with 3 open and 3 closed states. Thus, the assumptions of the above criterion are not fulfilled, and all transition rates are identifiable in the SM model for equal open times.

As a further illustration of the non-identifiability criterion, we investigate a slight variation of the loop gating scheme given in Fig. 3. According to Theorem B of Fredkin and Rice (1986), this gating scheme exceeds the maximum number of identifiable parameters. We can also derive the same result with the non-identifiability criterion. The degrees of freedom of the similarity transformation for this gating scheme are 2 × 1 + 2 × 1 = 4. Each of the 2 × 2 submatrices Q_OC, Q_CO contains 2 vanishing entries, however, these are not four independent restrictions to preserve the gating scheme, but only two. A similarity transformation cannot change the rank of the submatrices Q_OC, Q_CO; both have rank 1. Therefore, from the restriction that one of the vanishing entries in Q_OC and Q_CO, respectively, is preserved under similarity transformations, it follows that the other vanishing entry must also vanish under similarity transformations, otherwise the transformed submatrices Q_OC, Q_CO would have rank 2. Thus, the degrees of freedom of the similarity transformation for this gating scheme are 4 and the minimal number of restrictions to preserve the gating scheme is only 2, i.e., this gating scheme is not identifiable.

View larger version (11K): [in this window] [in a new window]   FIGURE 3   Variation of the loop gating scheme which is not non-identifiable.

This example indicates that, in general, the minimal number of restrictions m may not be trivial to determine because the ranks of the submatrices Q_CO and Q_OC are preserved by similarity transformations. We suppose that the investigation of infinitesimal similarity transformations is a successful strategy for determining m. In particular, the first order approximation of a similarity transformation is easy to compute even for more complicated models, and the requirement that the first order approximation of a similarity transformation is compatible with the given gating scheme gives a good hint for the minimal number of restrictions m.

Testing for detailed balance

The principle of detailed balance imposes one constraint on the transition rates per loop, i.e., one transition rate per loop can be calculated from the others by the law of detailed balance. In the case of the loop model, one restriction is not sufficient to remove the non-identifiability in the loop model with equal open times, because the family of indistinguishable aggregated Markov models with the same equal open times has two degrees of freedom and the constraint imposed by the law of detailed balance lowers the degrees of freedom only by one. In particular, it follows that a loop model with equal open times, which satisfies the principle of detailed balance, is statistically not distinguishable from a model that violates this principle. Detailed balance, therefore, is not determinable in the loop model with equal open times. This does not apply for a loop model with different open times but which follows the principle of detailed balance. The transition rates in such a model are always identifiable. However, as in the case of parameter estimation, tests for detailed balance will be affected by the non-identifiability for equal open times if the true underlying model has almost equal open times. This will be investigated further in the following.

The hypotheses are as follows. Under the null hypothesis the principle of detailed balance is fulfilled and under the alternative it is violated:

(3)

(4)

We use the likelihood ratio as a test statistic because it has the following favorable asymptotic property under the null hypothesis (Cox and Hinkley, 1974):

(5)

The asymptotic normality of the maximum likelihood estimators is a prerequisite of this result. This was recently shown by Bickel et al. (1998). L_n(.) denotes the log likelihood function as a function of the parameters without the constraint of detailed balance, the parameters are all transition rates, the parameter vector is the maximum likelihood estimator, and L_n^DB(.) denotes the likelihood function with the constraint of detailed balance. In the case of the loop model and the SM model, the parameters are all transition rates but one which is calculated from the principle of detailed balance, the parameter vector is the maximum likelihood estimator under the restriction of detailed balance. The subscript n denotes the number of data points. There are efficient algorithms to calculate the maximum likelihood estimator from the measured data (Baum et al., 1970; Horn and Lange, 1983; Fredkin and Rice, 1992; Albertsen and Hansen, 1994; Michalek and Timmer, 1999).

As indicated in Eq. 5, the twofold log likelihood ratio is asymptotically ₁²-distributed under the null hypothesis. The ² distribution has one degree of freedom because the law of detailed balance imposes one constraint on the transition rates per loop.

In the following, the power of the likelihood ratio test to detect violations of the law of detailed balance will be investigated. The natural logarithm ln K of the ratio of products of the transition rates in clockwise and counterclockwise direction serves as a measure for the strength of the violation of detailed balance. A biophysical justification for this measure is given below. If the gating of an ion channel follows the principle of detailed balance, ln K vanishes, and it is positive after possibly inverting the ratio of transition rates for a gating, which violates the principle of detailed balance.

The detection of violations of detailed balance requires a precise estimation of ln K. For a loop model with almost equal open times, the transition rates cannot be estimated reliably even for a large number of data; consequently, this is also true for the estimation of ln K. So a loop model with nearly equal open times that follows the law of detailed balance is hardly distinguishable from a model which violates detailed balance. Thus, the power of the test for detailed balance will depend not only on the strength of violation of the null hypothesis expressed by ln K but also on the ratio of open times denoted by ₂/₁. The dependence of the power on ₂/₁ is the subject of the first simulation study in the next section.

	SIMULATION STUDIES

TOP ABSTRACT INTRODUCTION THEORY SIMULATION STUDIES WHAT IS A SIGNIFICANT... DISCUSSION REFERENCES

The power of a likelihood ratio test approaches 1 for the number of data points going to infinity (Cox and Hinkley, 1974). The purpose of the first simulation study described in this section is to show that the power to detect violations of detailed balance in the loop model is still quite low for a number of data points in the order of magnitude available from experiments if the ratio of the open times ₂/₁ is only small enough.

The data were simulated by the loop model with the following generator matrix:

(6)

All transition rates are given in Hz. We denote the entry of a generator matrix Q in the ith row and jth column by q_ij. The open dwell times are: ₁ = 10 ms and ₂. The shut dwell times are given by the inverses of the eigenvalues of the submatrix Q_cc. The ratio of the open times ₂/₁ varied from 2 to 14 and the natural logarithm ln K of the ratio of products of the transition rates in clockwise and counterclockwise direction is varied independently from 0.22 to 3.0. For given ₂/₁ and ln K the parameters ₂ and q₄₁ of the generator matrix are calculated by the following formulas:

(7)

(8)

For each ₂/₁ and ln K, we simulate 500 recordings of length 105 s with a sampling rate of 5 kHz (2¹⁹ data points) and estimate the transition rates by the maximum likelihood method twice, namely with and without the restriction of detailed balance, and calculate the twofold difference in the log likelihood functions evaluated at the maximum likelihood parameters according to Eq. 5. The maximization of the likelihood function is performed numerically by the EM algorithm (Michalek and Timmer, 1999) and a non-linear maximization routine based on a quasi-newton method (NAG, 1997). For the calculation of the first derivatives of likelihood function, we use Fisher's identity (Fisher, 1925; Jamshidian and Jennrich, 1997) and the "sinch" algorithm described by Najfeld and Havel (1995) to evaluate the derivatives of the matrix exponential.

Fig. 4 summarizes the results of these simulations. The probability of rejecting the null hypothesis of detailed balance against the ratio of open time constants for a test to the 5% level is shown for different values of ln K. Below a ratio of 2 of the open time ratio ₂/₁, a reliable numerical estimation of the transition rates is not possible (see Wagner et al., 1999, for the order of magnitude of the estimation errors in the transition rates). The power of the likelihood ratio test significantly drops for smaller values of the open time ratio, as expected, whereas for increasing values of ln K, that is for stronger violations of detailed balance, the power of the likelihood ratio test increases.

View larger version (22K): [in this window] [in a new window]   FIGURE 4   Loop model: The probability to reject the null hypothesis of detailed balance against the ratio of open time constants for a test to the 5% level is shown for different values of ln K of the products of the transition rates in clockwise and counterclockwise direction. The simulated data sets have an approximate length of 105 s with a sampling rate of 5 kHz (524288 data points). The topmost line corresponds to ln K  3.00. The error bars indicate the standard deviation of the estimated rejection probability.

In the second simulation study, we demonstrate that detailed balance is determinable in a SM model with equal open times. We use the following generator matrix according to Scheme 4 from Song and Magleby (1994):

(9)

We investigate the distribution of the test statistic for the following cases: The law of detailed balance holds (q₆₃ = 23) and it is violated (q₆₃ = 2) with a ln K  2.44. We simulate 1000 recordings of length 105 s with a sampling rate of 10 kHz (2²⁰ data points) and estimate the transition rates and the twofold difference of the log likelihood functions in same way as in the first simulation study. Fig. 5 shows the distribution of test statistic under the null hypothesis. As expected it follows a ₁² distribution. In Fig. 6 the distribution of twofold log likelihood ratio under the alternative with a ln K  2.44 is shown. In the following section, we will discuss that ln K  2.44 is already a strong violation of detailed balance. Therefore, the power to detect violations is almost 1 for recordings of length 105 s. So in contrast to the loop model, detailed balance is determinable for equal open times in the SM model.

View larger version (14K): [in this window] [in a new window]   FIGURE 5   SM model. The probability distribution of the twofold log likelihood ratio under the null hypothesis compared to the asymptotic 12 distribution. The simulated data sets have an approximate length of 105 s with a sampling rate of 10 kHz (1,048,576 data points).

View larger version (14K): [in this window] [in a new window]   FIGURE 6   SM model. The probability distribution of the twofold log likelihood ratio under the alternative with ln K = 2.44. The simulated data sets have an approximate length of 105 s with a sampling rate of 10 kHz (1,048,576 data points).

	WHAT IS A SIGNIFICANT VIOLATION OF DETAILED BALANCE?

TOP ABSTRACT INTRODUCTION THEORY SIMULATION STUDIES WHAT IS A SIGNIFICANT... DISCUSSION REFERENCES

If the gating of an ion channel violates the principle of detailed balance in steady state, the ion channel must be driven by an external energy source. The typical amount of energy needed to influence the gating of an ion channel is given by the amount of work for a change in the geometrical conformation of the channel protein. Since subunits of a channel protein are often either polarized or carry some elementary charges, the work for a conformational change of the channel protein, e.g., from state i to state j, is associated with an activation energy E_ij. These energies are roughly bounded by the amount of work needed to push an elementary charge e₀ against the membrane potential:

(10)

with a membrane potential of 70 mV. Due to the law of Arrhenius, the transition rates q_ij are proportional to
We exemplify the calculations for the loop model. It is analogous in the SM model. If the gating obeys the principle of detailed balance, the total activation energy for one turn in the loop model must be the same in the clockwise and counterclockwise directions:

(11)

Assuming that the proportionality factor in the law of Arrhenius is approximately the same for all transitions between the states of the channel protein, Eq. 11 is equivalent to the following condition on the transition rates:

(12)

Thus, in the case of a violation of detailed balance, kT ln K does not vanish any more, and it is a measure for the difference in activation energies in clockwise and counterclockwise directions. Under the assumption that E₀ (Eq. 10) approximately bounds the work needed to influence the gating of an ion channel, a value of ln K in the order of magnitude of 2 would indicate a significant violation of the principle of detailed balance (at room temperature (T = 300 K) and a membrane potential of 70 mV):

(13)

In the first simulation study in the previous section, we investigated the power of the likelihood ratio test in the given range of ln K for significant violations of the principle of detailed balance in the loop model. The number of data points in the simulation study is of the order of magnitude typically available in experiments. The power drops towards smaller ratios of the open times below any upper bound for errors of second kind that is acceptable in practice. Therefore, we expect that violations of the principle of detailed balance might only be detected reliably in experiments where the dwell times are at least an order of magnitude different or with a very large amount of data.

	DISCUSSION

TOP ABSTRACT INTRODUCTION THEORY SIMULATION STUDIES WHAT IS A SIGNIFICANT... DISCUSSION REFERENCES

The SM model and the loop model differ in the identifiability of their transition rates for equal open times. In particular, the loop model exemplifies that the power to detect deviations from detailed balance depends not only on the strength of the violation of detailed balance, but also on the true transition rate itself, namely the ratio of the open times. This property of the loop model is due to the non-identifiability of its transition rates for equal open times. In contrast to the loop model, the SM model does not suffer from this non-identifiability because this model has many states, but only a few allowed transition between in the states and all states are part of the loop. Therefore, enough vanishing entries in the generator matrix need to be preserved under similarity transformations. In gating schemes for practical purposes, however, typically only a few states form a loop, so that these models have gating schemes like the loop model as submodels (see Vandenberg and Bezanilla, 1991), for some examples of physiologically relevant models with loops. Thus, we expect that effects like a lack of power must be taken into account for the analysis of measured ion channel data with realistic gating schemes.

Song and Magleby (1994) present a method to detect violations of the principle of detailed balance by comparing the estimated two-dimensional distributions in forward and backward direction of adjacent open and closed dwell time distribution. This method is based mainly on the visual inspection of estimated histograms of two-dimensional dwell time distributions which requires very long observations and a good signal-to-noise ratio due to the missed event problem (Blatz and Magleby, 1986; Ball et al., 1993). Likelihood ratio testing, presented in this paper, provides a statistically very reliable alternative approach for the following reasons: the asymptotic distribution of the log likelihood ratio under the null hypothesis is known to be exactly a ₁² distribution. The true finite sample size distribution of the test statistic deviates from the ₁² distribution, but the magnitude of this deviation is already small for typical sample sizes in patch clamp experiments as the magnitude is determined by the deviation of the finite sample size distribution of the maximum likelihood estimators for the transitions rates from normality. These maximum likelihood estimators already follow a normal distribution for rather small sample sizes (see Wagner et al., 1999, for a simulation study). Moreover, likelihood ratio testing can easily be extended to cases where the application of hidden Markov models is more appropriate for the analysis of measured ion channel data (Chung et al., 1990; Michalek et al., 1999).