- Email: [email protected]
Parametric survival analysis for gating kinetics of single potassium channels
Brain Research 770 Ž1997. 96–104
Research report
Parametric survival analysis for gating kinetics of single potassium channels Xuesi M. Shao
)
Department of Physiology, UCLA School of Medicine, Los Angeles, CA90095, USA Accepted 10 June 1997
Abstract The gating mechanism and the role of the S4 region in the activation of Shaker potassium channel was studied by statistical analysis on the wild type and mutant channels which have mutations in the S4 region. Single channel currents were recorded with the patch-clamp technique. The first latency time was analyzed with a parametric survival time regression model in which the generalized gamma distribution for the error term was assumed. Discrimination among Weibull, gamma, log-normal and exponential distributions was done with the likelihood ratio test. The three-parameter generalized gamma distribution was shown to be appropriate for the data set. The multiple regression function provides the statistical tests and the quantitative descriptions for the relationships between the mutations and the voltage dependence of the gating process. These results on statistical relations support the hypothesis that the S4 region plays an important role in sensing transmembrane voltage in the gating process, but the gating mechanism is not solely accounted for by the electrostatic interaction between the charged amino acids and the transmembrane voltage field. This work demonstrates that the survival analysis procedures can be useful tools for analyzing neurophysiological signals. q 1997 Elsevier Science B.V. Keywords: Shaker Kq channel; Single channel gating kinetics; Survival analysis; Shaker S4 region; Multiple regression analysis; Generalized gamma model
1. Introduction The Shaker Kq channel, a voltage gated potassium channel, is a membrane protein that controls the permeability of cells to Kq in response to changes in the transmembrane electric field. It is widely accepted that the channel consists of four identical subunits w20x and the amino acid sequence of the subunit has been identified w32x. The fourth transmembrane region ŽS4. of some voltage gated cation channels is hypothesized to be the voltage sensor which senses the voltage difference across the cell membrane and by way of its conformational change finally switches the channel open w10,11,14,16,17,19,21,23,31 x. The S4 sequences contain a repeated triplet motif of one positively charged amino acid Žarginine or lysine. plus two hydrophobic amino acids. The S4 sequence is considered to have a transmembrane disposition or to be buried in the low dielectric environment of the protein so that the charged amino acids can interact electrostatically with the )
Corresponding author. Systems Neurobiology Lab., Department of Physiological Science, University of California at Los Angeles, Los Angeles, CA90095-1527, USA. Fax: q1 Ž310. 206-9184; E-mail: [email protected] 0006-8993r97r$17.00 q 1997 Elsevier Science B.V. All rights reserved. PII S 0 0 0 6 - 8 9 9 3 Ž 9 7 . 0 0 7 6 0 - 9
transmembrane voltage to control the activity of the channel. In order to understand the role of the S4 region in channel activation and the gating mechanism of the Shaker Kq channel, mutations were made for a single amino acid residue in the S4 sequence. Two different mutant channels, R377K, in which the arginine at position 377 has been replaced by another positively charged amino acid, lysine, and R368Q, in which the arginine at position 368 has been replaced by the neutral amino acid, glutamine, have been made w29x. Single channel analysis w33,34x of the wild type Shaker channels has led to a picture of channel activation in which the channel undergoes a series of voltage-dependent Žcharge-moving. conformational changes between closed states on the way to opening. If the role of the S4 sequence in the activation process only involves the electrostatic interaction, we would expect that R368Q which reduces net positive charge in S4 will alter the voltage-dependent properties of the gating kinetics while R377K will not w18x. During the experiments we delivered repeated voltage pulses across the membrane, and single channel current was recorded w12x which is a random process with two values: 1 represents channel open, and 0 represents chan-
X.M. Shao r Brain Research 770 (1997) 96–104
nel closed. These data are a kind of survival data for which the time to the occurrence of some events is our major concern: e.g., the first latency, the time from starting of the voltage pulse to the first opening event; the channel open time, the time the channel stays open until it closes; the closed time; bursting time; etc. To analyze these conductance-state survival time data recorded with the above experimental scheme, we may face some difficulties: Ž1. frequently, data are censored due to observation time limits – the length of the voltage pulses – on data collection. Each observation ends at the end of each pulse but the channel events we are interested in may not be completed by then. For example, if we are interested in first latency, those records without opening during the entire pulses are considered ‘right censoring’. It does not mean there is no information available on the censored observations, but the information is partial. Dropping the censored data may cause bias and is a waste if the sample size is not large. Ž2. Sometimes we are interested in investigating a number of experimental variables simultaneously because most of the time more than one experimental variable influences the temporal properties of the channel event. There are possibly interactions among the experimental variables; therefore, descriptions for the relations between response variable and each explanatory variable individually would be inadequate. Ž3. Usually the distribution of this survival time is not normal and it may be very complicated under the Markov process assumption. A field of statistics, survival analysis, provides useful techniques for dealing with these types of problems. Many important developments in this field are quite recent. In this study the analysis was focused on the first latency of channel opening. The time of first latency reflects the series of conformational changes from closed states on the way to opening. I analyzed the latency data using the survival time regression procedure to determine whether the opening process is voltage dependent and how the mutations in the voltage sensor alter the properties of voltage-dependent gating of the channel. For the conventional curve fitting methods, establishing a random Markovian model and estimating the parameters for first latency is the major concern. But in this study I address a different aspect of the process; I establish a regression model of which the deterministic component describes the relationship between several experimental treatments and the response variable. The regression function provides quantitative descriptions of the relations, and the statistical significance of the relations can be tested as well. The coefficients of the interaction terms of the voltage and the mutant types yield interesting information about the hypothesized mechanism of the S4 region in sensing transmembrane voltages. This paper shows the potential usefulness of survival analysis methodology in the field of neurophysiology. The capabilities of survival analysis in dealing with censoring and non-normality problems make it a powerful technique for analyzing the dwell times of
97
the voltage gated ion channels or other event time physiological signals.
2. Methods 2.1. Electrophysiological recording and data collection Wild type and mutant channels were expressed in Xenopus oocytes w29x. Patch-clamp recordings were made on oocytes bathed in a ‘depolarizing solution’ Žcontaining 98 mM KCl, 3 mM MgCl 2 , 1 mM EGTA and 10 mM HEPES pH 7.2., in which the cell membrane potential should be close to 0 mV. The actual transmembrane potentials were close to the clamp voltages applied in the cell-attached mode. Patch electrodes were pulled to a tip diameter of about 1 m m and coated with Sylgard. Electrodes were filled with frog Ringer’s solution consisting of 115 mM NaCl, 2.5 mM KCl, 1.8 mM CaCl 2 and 10 mM HEPES ŽpH 7.2.. The patch potential was nulled just before seal formation. Single channel currents were recorded from cell-attached patches w12x using an Axopatch-1D with integrating headstage ŽAxon Instruments.. Most of the capacitative current was compensated using the capacitance compensation circuit of the patch-clamp amplifier. Voltage-clamp protocols were controlled, and data were collected using pCLAMP v.5.5.1 software and a TL-1 Labmaster interface ŽAxon Instruments. with an 80386-based computer. In each trial, step depolarizations of 80 ms duration were delivered from a holding potential of y80 mV. The data were filtered at a cutoff frequency of 5 kHz Žy3 dB. with an eight-pole low pass Bessel filter ŽFrequency Devices. and then digitized at a sampling interval of 50 m s. The data were then refiltered at 1500 Hz using the low pass digital Gaussian filter in the pCLAMP program. Linear leak currents and uncompensated capacity transients were digitally subtracted using templates made by averaging several blank records. Single channel data were idealized using a 50% amplitude criterion to detect opening and closing transitions w5x. The idealized single channel data were sorted and a ‘dead time’ correction was applied w4x. The data in this study were obtained from six patches from different oocytes injected with wild type Shaker B RNA, five patches from the mutant R377K and seven patches from the mutant R368Q. Each patch contained only one active channel. The data were collected under step voltages y20, 0, 20, 40, 60, 80 mV and the sample size of the data set was 8829 of which 2083 were right censored. 2.2. Statistical analysis The response variable of each trial is the time duration from the beginning of a voltage clamp pulse to the first
X.M. Shao r Brain Research 770 (1997) 96–104
98
opening event. If no open event was recorded during the 80 ms voltage pulse period, we considered it to be right censoring. The two explanatory variables are the value of the clamp voltage pulse and the type of channel. Regression analysis expresses two essential components of a statistical relation: Ž1. a tendency of the response variable to vary with the explanatory variables in a systematic fashion; Ž2. a scattering of points around the curve of statistical relationship which has a certain probability distribution. With the assumption that the effect of the explanatory variables on the survival time distribution is multiplicative Žwe will show this in Section 3., a general regression model, referred to as an ‘accelerated failure time model’ w28x, is given by: ln Ž T . s xX b q s ln Ž T0 . . In this equation, ln signifies natural logarithm, T is the survival time, lnŽT0 . is the error term in which T0 has a specified distribution and s is a scale parameter. x is a column vector of explanatory variables or covariates Žand xX its transpose., and b is a column vector of unknown regression parameters. In our model, we have ln Ž T . s b 0q b 1Vq b 2 RKq b 3 RQq b 12 RK)V q b 13 RQ)V q s ln Ž T0 .
Ž 1.
where V is the clamped membrane voltage, RK and RQ are indicator variables RK s RQ s
½ ½
1, 0,
if the channel is R377K otherwise
1, 0,
if the channel is R368Q otherwise
and RK)V and RQ)V are the interaction terms accounting for the experimental variables whose effects may interact, or in other words, for the possible effect of mutations on the voltage dependence of channel first latency. The regression function separates the contributions of different experimental variables and their possible interaction to the logarithm survival time and provides a coefficient to quantify the contribution of each of them. For the probability distribution of the random component, or error term, there are four widely used survival time models: exponential, Weibull, gamma and log-normal models. It is difficult to discriminate between them and rather little has been written on how to choose which model is appropriate w1,6–9x. In addition, the maximum likelihood estimation for the parameters of the gamma distribution is very difficult. By extending the generalized gamma model w24,30x, all these models can be embraced by a single parametric family. The generalized gamma density function is f Ž t. s
g Ž trt .
g ky1
tG Ž k .
exp Ž y Ž trt .
t ) 0, t ) 0, g ) 0, k ) 0.
g
. Ž 2.
As special cases of it, the exponential distribution has g s k s 1, the Weibull has k s 1, the gamma has g s 1 and the log-normal has k ™ `. We can discriminate the four competing models by ordinary parametric inference w9,22x. LIFEREG is a standard procedure in the statistical package SAS for parametric survival time analysis. If T0 is assigned to be gamma in LIFEREG procedure w28x, it is assumed to have the generalized gamma distribution with probability density function Žp.d.f.. g Ž t0 . s
d Ž t 0drd 2 .
1r d 2
t 0 G Ž 1rd 2 .
exp Ž yt 0drd 2 . ,
t 0 ) 0, d ) 0
Ž 3. which is a special case of Eq. Ž2. with k s 1rd 2 , g s d and t s Ž d 2 .1r d w13x. Here reparameterization has been done for computational convenience Žsee Appendix A.. The method of maximum likelihood is used in the regression process. Since we have censoring in the survival time data, a term that takes care of the censoring is included in the log-likelihood function: L Ž u . s Ýi ln f Ž y i < x i ,u . q Ý j ln S Ž Ymax < x j ,u .
Ž 4.
where the first sum is concerned with trials in which the first latencies t i Ž y i s lnŽ t i .. are uncensored. It involves f, the density function given by Eq. Ž14. in Appendix A, which depends on the vector of parameters us Ž k, s , b. and, corresponding to the ith observation y i , the vector of explanatory variables x i ; in this case the parameter m in Eq. Ž14. is taken as xXb. The second sum is concerned with trials in which no opening is observed before the end of the pulse, and the duration is Tmax Ž Ymax s lnŽTmax .., so-called censored times. It involves the survival function, defined as `
S Ž y < x j ,u . s P Ž Yj ) y < x j ,u . s
Hy f Ž y < x ,u . d y. j
Ž 5.
Maximizing the likelihood LŽu. with respect to u yields estimates of regression coefficients for the explanatory variables and estimates of the parameters for the probability distribution of the error term d and s from which k and g of the generalized gamma distribution can be derived. Discrimination among the four different models can be done by testing if k s 1, g s 1 or k ™ ` with the maximum likelihood ratio test. Let the parameter vector us Žf,q .. If we want to test the hypothesis that the parameter q s q0 , take R s y2 L Ž fˆ Ž q0 . ,q0 . y L Ž uˆ.
Ž 6.
where uˆ is the maximum likelihood estimator Žm.l.e.. of u and fˆŽ q0 . is the m.l.e. of f with q held at q0 . R follows an approximate chi-squared distribution with 1 degree of freedom.
X.M. Shao r Brain Research 770 (1997) 96–104
99
Fig. 1. Single channel currents were recorded in the cell-attached configuration from a wild type Shaker Kq channel. The top trace is the depolarization pulse which was delivered to the patch of membrane through the microelectrode. The pulses were from a holding potential of y80 mV to q60, q20 or y20 as indicated in traces 2, 3 and 4. First latency T is defined as the time from the beginning of the pulse to the first opening of the channel. ti is an observation of T.
Fig. 2. Histograms of first latencies of sample records from a R368Q channel. A: the abscissa is linear. The voltage pulse was 20 mV. B: the abscissa is logarithmic. The voltages are as indicated.
100
X.M. Shao r Brain Research 770 (1997) 96–104
Upon establishing the regression model, we should examine how good the model is, namely whether the assumptions about the systematic part and the assumed probability distribution for the random part of the model are appropriate. Residual analysis to assess the goodness of fit is essential. The residual was defined as e i s Ž y i y xXi b . rs
Ž 7.
where y i is lnŽ t i .. The residual plot, e i vs. xXib, is usually used to assess the adequacy of the regression function to see whether the functional part of the model is misspecified or whether there is a systematic lack of fit. If the systematic structure of the data is already taken into account during the regression process the plot should show that the residuals spread evenly around zero. A probability plot is a useful graphical tool for comparing the empirical distribution of the data and a theoretical distribution or, in other words, assessing how good the fit between the data and the theoretical distribution is. It is done by plotting the probability integral transform of the ordered observations F Ž x Žj. . against the points 1rŽ N q 1., ... jrŽ N q 1. ... NrŽ N q 1.. If the assumed theoretical distribution is correct, the plot should look roughly linear w25x. For regression analysis, residuals defined by Eq. Ž7. are often considered as behaving approximately like a random sample of size N from the underlying distribution of the error term. Thus the probability plot of residuals F Ž eŽ j. . Õs. jr Ž N q 1 .
Ž 8.
provides a check for the assumed distribution about the error term. F Ž eŽj. . is the cumulative distribution function evaluated at the point of eŽj. , where eŽj. is the ordered version of e i so that eŽ1. F eŽ2. F ... F eŽ N . . The censored data points were not included when I constructed the probability plot. N is the sample size, the number of the uncensored observations. The cumulative distribution func-
Table 1 SAS printout of LIFEREG procedure Variable
df
Estimate
S.E.M.
x2
P)x
Intercept V RK RQ RK)V RQ)V Scale s Shape d
1 1 1 1 1 1 1 1
1.462 y0.0126 0.422 0.268 y0.0056 y0.0052 1.313 y1.894
0.0376 0.0011 0.0480 0.0466 0.00144 0.00128 0.0169
1514.07 141.87 77.31 33.02 15.12 16.71
0.0001 0.0001 0.0001 0.0001 0.0001 0.0001
Log-likelihood for gammasy14905.07. df, degree of freedom. The data used in this regression analysis were obtained from six patches from different oocytes injected with wild type Shaker B RNA, five patches from the mutant R377K and seven patches from the mutant R368Q. Each patch contained only one active channel. The data were collected under step voltages y20, 0, 20, 40, 60, 80 mV and the sample size of the data set was 8829 of which 2083 were right censored.
tion is defined as F Ž e . s H0e f Ž w .dw, where f is given by ˆ The Eq. Ž14. with k equal to the estimated value k. coefficient of correlation for F w eŽj. x and jrŽ N q 1. was also calculated.
3. Results An example of random openings of a channel under voltage clamp pulses and different first latencies is shown in Fig. 1. Fig. 2A is a histogram of first latency for a R368Q channel clamp to 20 mV. It shows that the distribution of the first latency is highly skewed to the left with most latencies occurring within a short time. If we take logarithm for the x-axis ŽFig. 2B., the histograms look more symmetric and we can see that the position parameter
Fig. 3. Likelihood ratio statistics Ž2)log-likelihood. for gamma parameters. The dashed lines indicate the value of the chi-squared distribution with 1 df at the 0.05 significance level. A: likelihood ratio statistics from fitting the data with the shape parameter k fixed at different values around the m.l.e. B: likelihood ratio statistics from fitting the data with the scale parameter g fixed at different values around the m.l.e.
X.M. Shao r Brain Research 770 (1997) 96–104
shifts to the left as the clamp voltage increases. This suggests that the accelerated failure time model which assumes that the effect of explanatory variables on the event-time distribution is multiplicative would be appropriate. By fitting the model ŽEq. Ž1.., we obtained the following regression function: ln Ž T . s 1.46y0.013Vq0.42RKq0.27RQ y 0.0056RK)V y 0.0052RQ)V
Ž 9.
Table 1 is the SAS printout giving the estimates of regression coefficients for explanatory variables and estimates of the parameters for the probability distribution of the error term. It also gives standard error and P-value for each
101
estimate. The m.l.e. of the gamma scale parameter s is 1.313 and the shape parameter d is y1.894. From these parameters, the gamma distribution parameters g s 1.443 and k s 0.279 were found Žsee reparameterization in Appendix A.. Table 1 also shows that all the regression coefficients are statistically significant Ž P-values - 0.0001. They are significant even with the consideration of simultaneous inferences with the Bonferroni approach w27x.. The regression coefficients for V, RQ and RK indicate that the voltage and the two mutations have effects on the first latency time. The minus sign of the coefficient for V shows that the slope of the voltage dependence of the latency time is negative. The values of the coefficients for
Fig. 4. A: probability plot for the log residuals Žsee the definitions in Eq. Ž7. and Eq. Ž8. in Section 2.. The ordinate is the estimated distribution function of the ordered residuals. B: residual plot for the fitted model Žsee Eq. Ž7. and Section 2..
102
X.M. Shao r Brain Research 770 (1997) 96–104
the terms RQ and RK in the regression function ŽEq. Ž9.. are large suggesting that two mutations have strong effects on latency time and the difference between these effects is small. And, surprisingly, the coefficient for R377K is larger than that for R368Q. The significant coefficients for the two interaction terms RK)V and RQ)V indicate that the voltage and the mutations affecting the first latency time do interact. In other words, the two mutations affect the voltage dependence of the opening process. Possibly, the charge moving process has been altered by the two mutations in the S4 region. The minus signs of the coefficients for the two interaction terms RK)V and RQ)V show that the mutations decrease the slope of the voltage dependence. Since the slope is negative, the mutations actually increase the absolute value of the slope of voltage dependence. From the values of the coefficients of the two interaction terms, we can see that the effects of R377K and R368Q on voltage dependence are similar. It implies that the activation mechanism cannot be explained solely by the electrostatic interaction. To evaluate the regression model, we consider the two essential components, the distribution of the random error term and the systematic dependence of the response variable on the explanatory variables, separately, as follows. We assessed the distribution of the data using the likelihood ratio tests with H 0 : g s 1 vs. H 1: g / 1. We did regressions by fixing the scale parameter at different values corresponding to the values of g shown in Fig. 3B. We then got the log-likelihood ratio statistics from each regression. Fig. 3B established by this method shows that H 0 should be rejected at the 5% significance level. Fig. 3A is constructed with the same procedure by fixing the shape parameter at different values. It shows that H 0 : k s 1 and k ™ ` should both be rejected. These results indicate that the data are not Weibull, gamma, exponential or log-normal distributed. From this result a question is naturally raised: are the data really following the three-parameter generalized gamma distribution? In other words, is the empirical distribution of the data sufficiently well described by a generalized gamma distribution? A probability plot for the residuals was made ŽFig. 4A.. The straightness of the estimated generalized gamma distribution of the ordered residuals vs. the probability is satisfactory and the coefficient of correlation is 0.996. There is no large departure from a straight line for most of the plot. To assess the adequacy of the regression function, residuals were computed by Eq. Ž7. in Section 2. The plot of the residuals vs. the fitted values ŽFig. 4B. shows the residuals spread evenly around zero. The upper bound of the residuals seems to decrease as the fitted value increases, which may be due to the fact that our censoring is always at 80 ms. We certainly do not have any observation bigger than 80. The y i are always less than lnŽ80. s 4.38, so for the larger fitted value xXb, Ž4.38 y xXb.rs would be smaller. Thus, the distribution of the residuals here does
not suggest a deviation from the regression function and the regression model is adequate and acceptable.
4. Discussion Parametric survival analysis has been applied to the first latency for single channel opening of wild type and mutant Shaker Kq channels. The application of the survival time regression procedure to single channel behavior is novel. By introducing a multiple regression model, we reveal some properties of gating kinetics which are difficult to detect and quantitatively measure by classical single channel analysis methods w29x. The effects of the experimental variables Žvoltage and mutant types. on the statistical properties of the first latency were summarized by a single regression function as Eq. Ž9.. Each different effect is quantified with the corresponding coefficient. The alternative survival time models, Weibull, exponential, ordinary gamma and log-normal, were evaluated by parametric inference under the generalized gamma model and they were all rejected by likelihood ratio tests. The goodness of fit and the adequacy of the regression function were assessed by the probability plot and residual plot. The regression results suggest: Ž1. First latency is voltage dependent. The minus sign of the regression coefficient for voltage indicates that the channel takes a shorter time to make the transition from channel closed to open under higher voltages. It implies that in the process of going from closed states to opening, in other words, the series of conformational changes that occur before the channel opens involve charge movement. Ž2. The two mutations have significant effects on the first latency. The positive regression coefficients for types R377K and R368Q tell us how much longer the mutant channels take from closed to open as compared to the wild type channel. The values of the coefficients show that R377K has a stronger effect on the first latency than R368Q rather than the other way around as expected w18x. Since the mutations in the S4 region alter the voltage sensitive process of channel activation in some way, it supports the hypothesis that the S4 region is a part of the voltage sensor w16,17,21,23,29,31x. Ž3. The regression coefficients for the two interaction terms show that the mutations R377K and R368Q have a significant effect on the voltage dependence of first latency. The latency times of the two mutants as a function of voltages have different slopes from the wild type channel, but the slope changes for the two mutants are similar. It implies that the mutations change the charge moving process but the voltage sensitive gating does not solely depend on the charge contents of the S4 region. This result supports the idea that the mechanism of activation is not due solely to an electrostatic interaction between the S4 positively charged amino acids and the membrane potential. These positive residues are not functionally equivalent
X.M. Shao r Brain Research 770 (1997) 96–104
w21x. Other properties of the residues apart from the charge, e.g., the length of the side chains or the interaction between these residues and the residues of other parts of the channel, may contribute to the gating mechanism of the channel as well. Ž4. The probability distribution of first latency may be well modeled by the three-parameter generalized gamma distribution. The Shaker Kq channel is thought to be a tetramer comprised of four identical subunits w20x and all four subunits have to be in open state for the channel to be open. The channel is closed if any one subunit is in closed state. If the simple kinetic model a
Close ° Open
Ž 10 .
b
for each subunit is correct and if the transition rates a and b are constants Ždo not vary with time., we would expect that the time it stays in any one state is exponentially distributed with parameter a or b respectively w4x. This means that the first latency from closed state to open state for each subunit Žthis is just for a single subunit from closed state to open state and it is not observable until the whole channel is open. should follow an exponential distribution with mean 1ra . If we assume the latency times for the subunits have a common distribution and the subunits gate independently, and if we further assume that the data obtained at high voltages where the contribution of the backward rate b values to first latency could be negligible w33x, the distribution of the time to channel opening would be the distribution of the maximum of the four-order statistics of exponential w26x. Then the cumulative distribution function would be 4
F Ž t . s Ž 1 y exp Ž ya t . . .
Ž 11 .
From this, we have the p.d.f.: f Ž t . s 4a exp Ž ya t . y12 a exp Ž y2 a t . q 12 a exp Ž y3 a t . i y 4a exp Ž y4a t . .
Ž 12 .
If the backward rate b values are not negligible, the probability expression becomes more complicated and has to be expressed in matrix form w4x. In addition to this complexity, recent experiments have shown that functional or structural interactions among subunits during the transitions from closed to open exist w2,33,34x. Furthermore, models like Eq. Ž10. for each subunit may be oversimplified. Each subunit may undergo more than one transition from closed state to open state w3x. Unfortunately, increasing the number of exponentials to account for these different closed states makes the estimation of parameters unreliable since the curve fitting would involve a large number of parameters. More fundamentally, whether the transitions in the gating process are Markovian is still controversial in channel biophysics w15x. From the regression analysis point of view, the three-parameter gamma distribution is flexible
103
in dealing with this complexity. We see from the residual analysis that the generalized gamma model is an acceptable approximation for the distribution. Our conclusions based on parametric survival analysis are qualitatively consistent with those obtained from conventional methods w29x. If we plot the empirical cumulative distribution functions of first latency for different mutants at different voltages, and plot the medians of first latency for different mutants at different voltages, we can show that first latency is voltage dependent and the two mutations have effects on the first latency and on the voltage dependence of the first latency. It is not clear, however, which effects are stronger and how strong they are. If we do multi-exponential curve fittings for different mutants at different voltages, we get an estimate for each combination of the mutant and voltage, but the interactions between the functional relations between first latency and the mutations, and between first latency and voltage are unclear. The parametric survival analysis method suggested in this paper is different from the conventional methods in the following aspects: Ž1. it provides quantitative descriptions of the relations between the response variable and each experimental treatment and the interactions between experimental variables; Ž2. it provides significance tests for those relations; Ž3. it deals with the censoring problems and thus it results in less bias in parameter estimation; Ž4. when the number of experimental treatments is large or when these experimental variables vary simultaneously, these multiple regression models can be extended very easily and can be standardized with recent powerful statistical analysis packages. It provides a simpler way to deal with large amounts of data before we have a sound biophysical understanding of the channel gating process to establish more analytical models. A number of different signals from nervous systems, e.g., the dwell time of single channel records, the interspike intervals of single neuronal activity and spontaneous or miniature postsynaptic events, excitatory postsynaptic potentials ŽEPSP. or inhibitory postsynaptic potentials ŽIPSP., can be described by certain kinds of survival time models. The application of survival analysis methodology, as illustrated in this paper, can be extended in the field of neuroscience.
Acknowledgements The author wishes to thank the referees for very helpful comments which led to substantial improvements, Dr. Diane Papazian for her support of this work and allowing the author to use the resources of her lab, Dr. Donald Guthrie for his advice on this project, and Mr. Calvin Cliff for his comments on the manuscripts. This work was supported by NIH Grant GM43459.
X.M. Shao r Brain Research 770 (1997) 96–104
104
Appendix A For computational reasons, reparameterization is done as follows. We consider W s lnŽT0 ., which has p.d.f. f Ž w. s
kk
'k G Ž k . exp
'k w y k )exp Ž wr'k .
.
Ž 13 .
Let Y s lnŽT . and xXbs m , then W s Ž Y y m .rs and Y has p.d.f. f Ž y. s
kk
'k sG Ž k . exp
'k Ž y y m .
rs y k )exp Ž Ž y y m . r Ž s'k . .
Ž 14 .
and the corresponding distribution of T s expŽ Y . has p.d.f. of the form of Eq. Ž2., with g s 1rŽ s 6k . and t s expŽ m .rk 1r g or, equivalently, s s 1rŽg 6k . and m s lnŽt . q lnŽ k .rg . References w1x L.J. Bain, M. Engelhardt, Probability of correct selection of Weibull versus gamma based on likelihood ratio, Commun. Stat. Theory Methods A9 Ž1980. 375–381. w2x F. Bezanilla, E. Perozo, D.M. Papazian, E. Stefani, Molecular basis of gating charge immobilization in Shaker potassium channels, Science 254 Ž1991. 679–683. w3x F. Bezanilla, E. Perozo, E. Stefani, Gating of Shaker Kq channels. II. The component of gating currents and a model of channel activation, Biophys. J. 66 Ž1994. 1011–1021. w4x D. Colquhoun, A.G. Hawkes, The principles of the stochastic interpretation of ion-channel mechanisms, in: B. Sakmann, E. Neher ŽEds.., Single Channel Recording, Plenum Press, New York, 1983, pp. 135–174. w5x D. Colquhoun, F.J. Sigworth, Fitting and statistical analysis of single channel records, in: B. Sakmann, E. Neher ŽEds.., Single Channel Recording, Plenum Press, New York, 1983, pp. 191–263. w6x D.R. Cox, Tests of separate families of hypotheses, Proceedings of the 4th Berkeley Symposium, vol. 1, 1961, pp. 105–123. w7x D.R. Cox, Further results on tests of separate families of hypotheses, J. R. Stat. Soc. B24 Ž1962. 406–424. w8x R. Dumonceaux, C.E. Antle, Discrimination between the log-normal and the Weibull distributions, Technometrics 15 Ž1973. 923–926. w9x V.T. Farewell, R.L. Prentice, A study of distributional shape in life testing, Technometrics 19 Ž1977. 69–75. w10x R.E. Greenblatt, E. Blatt, M. Montal, The structure of the voltagesensitive sodium channel. Inferences derived from computer-aided analysis of the Electrophorus electricus channel primary structure, FEBS Lett. 193 Ž1985. 125–134. w11x H.R. Guy, P. Seetharamulu, Molecular model of the action potential sodium channel, Proc. Natl. Acad. Sci. USA 83 Ž1986. 508–512. w12x O.P. Hamill, A. Marty, E. Neher, B. Sakmann, F.J. Sigworth, Improved patch-clamp techniques for high resolution current recording from cells and cell-free membrane patches, Pflugers Arch. 391 ¨ Ž1981. 85–100. w13x J.F. Lawless, The generalized gamma distribution, in: Statistical Models and Methods for Lifetime Data, John Wiley and Sons, New York, 1982, pp. 237–252.
w14x H.R. Leuchtag, Long-range interactions, voltage sensitivity, and ion conduction in S4 segments of excitable channels, Biophys. J. 66 Ž1994. 217–224. w15x L.S. Liebovitch, J.M. Sullivan, Fractal analysis of a voltage-dependent potassium channel from cultured mouse hippocampal neurons, Biophys. J. 52 Ž1987. 979–988. w16x E.R. Liman, P. Hess, F. Weaver, G. Koren, Voltage-sensing residues in the S4 region of a mammalian Kq channel, Nature 353 Ž1991. 752–756. w17x D.E. Logothetis, S. Movahedi, C. Satler, K. Lindpaintner, B. NadalGinard, Incremental reductions of positive charge within the S4 region of a voltage-gated Kq channel result in corresponding decreases in gating charge, Neuron 8 Ž1992. 531–540. w18x D.E. Logothetis, B.F. Kammen, K. Lindpaintner, D. Bisbas, B. Nadal-Ginard, Gating charge differences between two voltage-gated Kq channels are due to the specific charge content of their respective S4 regions, Neuron 10 Ž1993. 1121–1129. w19x G.A. Lopez, Y.N. Jan, L.Y. Jan, Hydrophobic substitution mutations in the S4 sequence alter voltage-dependent gating in Shaker Kq channels, Neuron 7 Ž1991. 327–336. w20x R. MacKinnon, Determination of the subunit stoichiometry of a voltage-activated potassium channel, Nature 350 Ž1991. 232–235. w21x D.M. Papazian, L.C. Timpe, Y.N. Jan, L.Y. Jan, Alteration of voltage-dependence of Shaker potassium channel by mutations in the S4 sequence, Nature 349 Ž1991. 305–310. w22x V.B. Parr, J.T. Webster, A method for discriminating between failure density functions used in reliability predictions, Technometrics 7 Ž1965. 1–10. w23x E. Perozo, L. Santacruz-Toloza, E. Stefani, F. Bezanilla, D.M. Papazian, S4 mutations alter gating currents of Shaker Kq channels, Biophys. J. 66 Ž1994. 345–354. w24x R.L. Prentice, A log-gamma model and its maximum likelihood estimation, Biometrika 61 Ž1974. 539–544. w25x J.A. Rice, Probability plot, in: Mathematical Statistics and Data Analysis, Wadsworth and Brook, Pacific Grove, CA, 1988, pp. 292–299. w26x J.A. Rice, Extrema and order statistics, in: Mathematical Statistics and Data Analysis, Wadsworth and Brook, Pacific Grove, CA, 1988, pp. 94–97. w27x J.A. Rice, The problem of multiple comparisons, in: Mathematical Statistics and Data Analysis, Wadsworth and Brook, Pacific Grove, CA, 1988, pp. 404–408. w28x SAS Institute Inc., The LIFEREG procedure, in: SASrSTAT User’s Guide, vol. 2, version 6, 4th ed., SAS Institute Inc., Cary, NC, 1989, pp. 997–1017. w29x X.M. Shao, D.M. Papazian, S4 mutations alter the single-channel gating kinetics of Shaker Kq channels, Neuron 11 Ž1993. 343–352. w30x E.W. Stacy, A generalization of the gamma distribution, Ann. Math. Stat. 33 Ž1962. 1187–1192. w31x W. Stuhmer, F. Conti, H. Suzuki, X. Wang, M. Noda, N. Yahagi, H. Kubo, S. Numa, Structural parts involved in activation and inactivation of the sodium channel, Nature 339 Ž1989. 597–603. w32x B.L. Tempel, D.M. Papazian, T.L. Schwarz, Y.N. Jan, L.Y. Jan, Sequence of a probable potassium channel component encoded at Shaker locus of Drosophila, Science 237 Ž1987. 770–775. w33x W.N. Zagotta, R.W. Aldrich, Voltage-dependent gating of Shaker A-type potassium channels in Drosophila muscle, J. Gen. Physiol. 95 Ž1990. 29–60. w34x W.N. Zagotta, T. Hoshi, R.W. Aldrich, Gating of single Shaker potassium channels in Drosophila muscle and in Xenopus oocytes injected with Shaker mRNA, Proc. Natl. Acad. Sci. USA. 86 Ž1989. 7243–7247.
Research report
Parametric survival analysis for gating kinetics of single potassium channels Xuesi M. Shao
)
Department of Physiology, UCLA School of Medicine, Los Angeles, CA90095, USA Accepted 10 June 1997
Abstract The gating mechanism and the role of the S4 region in the activation of Shaker potassium channel was studied by statistical analysis on the wild type and mutant channels which have mutations in the S4 region. Single channel currents were recorded with the patch-clamp technique. The first latency time was analyzed with a parametric survival time regression model in which the generalized gamma distribution for the error term was assumed. Discrimination among Weibull, gamma, log-normal and exponential distributions was done with the likelihood ratio test. The three-parameter generalized gamma distribution was shown to be appropriate for the data set. The multiple regression function provides the statistical tests and the quantitative descriptions for the relationships between the mutations and the voltage dependence of the gating process. These results on statistical relations support the hypothesis that the S4 region plays an important role in sensing transmembrane voltage in the gating process, but the gating mechanism is not solely accounted for by the electrostatic interaction between the charged amino acids and the transmembrane voltage field. This work demonstrates that the survival analysis procedures can be useful tools for analyzing neurophysiological signals. q 1997 Elsevier Science B.V. Keywords: Shaker Kq channel; Single channel gating kinetics; Survival analysis; Shaker S4 region; Multiple regression analysis; Generalized gamma model
1. Introduction The Shaker Kq channel, a voltage gated potassium channel, is a membrane protein that controls the permeability of cells to Kq in response to changes in the transmembrane electric field. It is widely accepted that the channel consists of four identical subunits w20x and the amino acid sequence of the subunit has been identified w32x. The fourth transmembrane region ŽS4. of some voltage gated cation channels is hypothesized to be the voltage sensor which senses the voltage difference across the cell membrane and by way of its conformational change finally switches the channel open w10,11,14,16,17,19,21,23,31 x. The S4 sequences contain a repeated triplet motif of one positively charged amino acid Žarginine or lysine. plus two hydrophobic amino acids. The S4 sequence is considered to have a transmembrane disposition or to be buried in the low dielectric environment of the protein so that the charged amino acids can interact electrostatically with the )
Corresponding author. Systems Neurobiology Lab., Department of Physiological Science, University of California at Los Angeles, Los Angeles, CA90095-1527, USA. Fax: q1 Ž310. 206-9184; E-mail: [email protected] 0006-8993r97r$17.00 q 1997 Elsevier Science B.V. All rights reserved. PII S 0 0 0 6 - 8 9 9 3 Ž 9 7 . 0 0 7 6 0 - 9
transmembrane voltage to control the activity of the channel. In order to understand the role of the S4 region in channel activation and the gating mechanism of the Shaker Kq channel, mutations were made for a single amino acid residue in the S4 sequence. Two different mutant channels, R377K, in which the arginine at position 377 has been replaced by another positively charged amino acid, lysine, and R368Q, in which the arginine at position 368 has been replaced by the neutral amino acid, glutamine, have been made w29x. Single channel analysis w33,34x of the wild type Shaker channels has led to a picture of channel activation in which the channel undergoes a series of voltage-dependent Žcharge-moving. conformational changes between closed states on the way to opening. If the role of the S4 sequence in the activation process only involves the electrostatic interaction, we would expect that R368Q which reduces net positive charge in S4 will alter the voltage-dependent properties of the gating kinetics while R377K will not w18x. During the experiments we delivered repeated voltage pulses across the membrane, and single channel current was recorded w12x which is a random process with two values: 1 represents channel open, and 0 represents chan-
X.M. Shao r Brain Research 770 (1997) 96–104
nel closed. These data are a kind of survival data for which the time to the occurrence of some events is our major concern: e.g., the first latency, the time from starting of the voltage pulse to the first opening event; the channel open time, the time the channel stays open until it closes; the closed time; bursting time; etc. To analyze these conductance-state survival time data recorded with the above experimental scheme, we may face some difficulties: Ž1. frequently, data are censored due to observation time limits – the length of the voltage pulses – on data collection. Each observation ends at the end of each pulse but the channel events we are interested in may not be completed by then. For example, if we are interested in first latency, those records without opening during the entire pulses are considered ‘right censoring’. It does not mean there is no information available on the censored observations, but the information is partial. Dropping the censored data may cause bias and is a waste if the sample size is not large. Ž2. Sometimes we are interested in investigating a number of experimental variables simultaneously because most of the time more than one experimental variable influences the temporal properties of the channel event. There are possibly interactions among the experimental variables; therefore, descriptions for the relations between response variable and each explanatory variable individually would be inadequate. Ž3. Usually the distribution of this survival time is not normal and it may be very complicated under the Markov process assumption. A field of statistics, survival analysis, provides useful techniques for dealing with these types of problems. Many important developments in this field are quite recent. In this study the analysis was focused on the first latency of channel opening. The time of first latency reflects the series of conformational changes from closed states on the way to opening. I analyzed the latency data using the survival time regression procedure to determine whether the opening process is voltage dependent and how the mutations in the voltage sensor alter the properties of voltage-dependent gating of the channel. For the conventional curve fitting methods, establishing a random Markovian model and estimating the parameters for first latency is the major concern. But in this study I address a different aspect of the process; I establish a regression model of which the deterministic component describes the relationship between several experimental treatments and the response variable. The regression function provides quantitative descriptions of the relations, and the statistical significance of the relations can be tested as well. The coefficients of the interaction terms of the voltage and the mutant types yield interesting information about the hypothesized mechanism of the S4 region in sensing transmembrane voltages. This paper shows the potential usefulness of survival analysis methodology in the field of neurophysiology. The capabilities of survival analysis in dealing with censoring and non-normality problems make it a powerful technique for analyzing the dwell times of
97
the voltage gated ion channels or other event time physiological signals.
2. Methods 2.1. Electrophysiological recording and data collection Wild type and mutant channels were expressed in Xenopus oocytes w29x. Patch-clamp recordings were made on oocytes bathed in a ‘depolarizing solution’ Žcontaining 98 mM KCl, 3 mM MgCl 2 , 1 mM EGTA and 10 mM HEPES pH 7.2., in which the cell membrane potential should be close to 0 mV. The actual transmembrane potentials were close to the clamp voltages applied in the cell-attached mode. Patch electrodes were pulled to a tip diameter of about 1 m m and coated with Sylgard. Electrodes were filled with frog Ringer’s solution consisting of 115 mM NaCl, 2.5 mM KCl, 1.8 mM CaCl 2 and 10 mM HEPES ŽpH 7.2.. The patch potential was nulled just before seal formation. Single channel currents were recorded from cell-attached patches w12x using an Axopatch-1D with integrating headstage ŽAxon Instruments.. Most of the capacitative current was compensated using the capacitance compensation circuit of the patch-clamp amplifier. Voltage-clamp protocols were controlled, and data were collected using pCLAMP v.5.5.1 software and a TL-1 Labmaster interface ŽAxon Instruments. with an 80386-based computer. In each trial, step depolarizations of 80 ms duration were delivered from a holding potential of y80 mV. The data were filtered at a cutoff frequency of 5 kHz Žy3 dB. with an eight-pole low pass Bessel filter ŽFrequency Devices. and then digitized at a sampling interval of 50 m s. The data were then refiltered at 1500 Hz using the low pass digital Gaussian filter in the pCLAMP program. Linear leak currents and uncompensated capacity transients were digitally subtracted using templates made by averaging several blank records. Single channel data were idealized using a 50% amplitude criterion to detect opening and closing transitions w5x. The idealized single channel data were sorted and a ‘dead time’ correction was applied w4x. The data in this study were obtained from six patches from different oocytes injected with wild type Shaker B RNA, five patches from the mutant R377K and seven patches from the mutant R368Q. Each patch contained only one active channel. The data were collected under step voltages y20, 0, 20, 40, 60, 80 mV and the sample size of the data set was 8829 of which 2083 were right censored. 2.2. Statistical analysis The response variable of each trial is the time duration from the beginning of a voltage clamp pulse to the first
X.M. Shao r Brain Research 770 (1997) 96–104
98
opening event. If no open event was recorded during the 80 ms voltage pulse period, we considered it to be right censoring. The two explanatory variables are the value of the clamp voltage pulse and the type of channel. Regression analysis expresses two essential components of a statistical relation: Ž1. a tendency of the response variable to vary with the explanatory variables in a systematic fashion; Ž2. a scattering of points around the curve of statistical relationship which has a certain probability distribution. With the assumption that the effect of the explanatory variables on the survival time distribution is multiplicative Žwe will show this in Section 3., a general regression model, referred to as an ‘accelerated failure time model’ w28x, is given by: ln Ž T . s xX b q s ln Ž T0 . . In this equation, ln signifies natural logarithm, T is the survival time, lnŽT0 . is the error term in which T0 has a specified distribution and s is a scale parameter. x is a column vector of explanatory variables or covariates Žand xX its transpose., and b is a column vector of unknown regression parameters. In our model, we have ln Ž T . s b 0q b 1Vq b 2 RKq b 3 RQq b 12 RK)V q b 13 RQ)V q s ln Ž T0 .
Ž 1.
where V is the clamped membrane voltage, RK and RQ are indicator variables RK s RQ s
½ ½
1, 0,
if the channel is R377K otherwise
1, 0,
if the channel is R368Q otherwise
and RK)V and RQ)V are the interaction terms accounting for the experimental variables whose effects may interact, or in other words, for the possible effect of mutations on the voltage dependence of channel first latency. The regression function separates the contributions of different experimental variables and their possible interaction to the logarithm survival time and provides a coefficient to quantify the contribution of each of them. For the probability distribution of the random component, or error term, there are four widely used survival time models: exponential, Weibull, gamma and log-normal models. It is difficult to discriminate between them and rather little has been written on how to choose which model is appropriate w1,6–9x. In addition, the maximum likelihood estimation for the parameters of the gamma distribution is very difficult. By extending the generalized gamma model w24,30x, all these models can be embraced by a single parametric family. The generalized gamma density function is f Ž t. s
g Ž trt .
g ky1
tG Ž k .
exp Ž y Ž trt .
t ) 0, t ) 0, g ) 0, k ) 0.
g
. Ž 2.
As special cases of it, the exponential distribution has g s k s 1, the Weibull has k s 1, the gamma has g s 1 and the log-normal has k ™ `. We can discriminate the four competing models by ordinary parametric inference w9,22x. LIFEREG is a standard procedure in the statistical package SAS for parametric survival time analysis. If T0 is assigned to be gamma in LIFEREG procedure w28x, it is assumed to have the generalized gamma distribution with probability density function Žp.d.f.. g Ž t0 . s
d Ž t 0drd 2 .
1r d 2
t 0 G Ž 1rd 2 .
exp Ž yt 0drd 2 . ,
t 0 ) 0, d ) 0
Ž 3. which is a special case of Eq. Ž2. with k s 1rd 2 , g s d and t s Ž d 2 .1r d w13x. Here reparameterization has been done for computational convenience Žsee Appendix A.. The method of maximum likelihood is used in the regression process. Since we have censoring in the survival time data, a term that takes care of the censoring is included in the log-likelihood function: L Ž u . s Ýi ln f Ž y i < x i ,u . q Ý j ln S Ž Ymax < x j ,u .
Ž 4.
where the first sum is concerned with trials in which the first latencies t i Ž y i s lnŽ t i .. are uncensored. It involves f, the density function given by Eq. Ž14. in Appendix A, which depends on the vector of parameters us Ž k, s , b. and, corresponding to the ith observation y i , the vector of explanatory variables x i ; in this case the parameter m in Eq. Ž14. is taken as xXb. The second sum is concerned with trials in which no opening is observed before the end of the pulse, and the duration is Tmax Ž Ymax s lnŽTmax .., so-called censored times. It involves the survival function, defined as `
S Ž y < x j ,u . s P Ž Yj ) y < x j ,u . s
Hy f Ž y < x ,u . d y. j
Ž 5.
Maximizing the likelihood LŽu. with respect to u yields estimates of regression coefficients for the explanatory variables and estimates of the parameters for the probability distribution of the error term d and s from which k and g of the generalized gamma distribution can be derived. Discrimination among the four different models can be done by testing if k s 1, g s 1 or k ™ ` with the maximum likelihood ratio test. Let the parameter vector us Žf,q .. If we want to test the hypothesis that the parameter q s q0 , take R s y2 L Ž fˆ Ž q0 . ,q0 . y L Ž uˆ.
Ž 6.
where uˆ is the maximum likelihood estimator Žm.l.e.. of u and fˆŽ q0 . is the m.l.e. of f with q held at q0 . R follows an approximate chi-squared distribution with 1 degree of freedom.
X.M. Shao r Brain Research 770 (1997) 96–104
99
Fig. 1. Single channel currents were recorded in the cell-attached configuration from a wild type Shaker Kq channel. The top trace is the depolarization pulse which was delivered to the patch of membrane through the microelectrode. The pulses were from a holding potential of y80 mV to q60, q20 or y20 as indicated in traces 2, 3 and 4. First latency T is defined as the time from the beginning of the pulse to the first opening of the channel. ti is an observation of T.
Fig. 2. Histograms of first latencies of sample records from a R368Q channel. A: the abscissa is linear. The voltage pulse was 20 mV. B: the abscissa is logarithmic. The voltages are as indicated.
100
X.M. Shao r Brain Research 770 (1997) 96–104
Upon establishing the regression model, we should examine how good the model is, namely whether the assumptions about the systematic part and the assumed probability distribution for the random part of the model are appropriate. Residual analysis to assess the goodness of fit is essential. The residual was defined as e i s Ž y i y xXi b . rs
Ž 7.
where y i is lnŽ t i .. The residual plot, e i vs. xXib, is usually used to assess the adequacy of the regression function to see whether the functional part of the model is misspecified or whether there is a systematic lack of fit. If the systematic structure of the data is already taken into account during the regression process the plot should show that the residuals spread evenly around zero. A probability plot is a useful graphical tool for comparing the empirical distribution of the data and a theoretical distribution or, in other words, assessing how good the fit between the data and the theoretical distribution is. It is done by plotting the probability integral transform of the ordered observations F Ž x Žj. . against the points 1rŽ N q 1., ... jrŽ N q 1. ... NrŽ N q 1.. If the assumed theoretical distribution is correct, the plot should look roughly linear w25x. For regression analysis, residuals defined by Eq. Ž7. are often considered as behaving approximately like a random sample of size N from the underlying distribution of the error term. Thus the probability plot of residuals F Ž eŽ j. . Õs. jr Ž N q 1 .
Ž 8.
provides a check for the assumed distribution about the error term. F Ž eŽj. . is the cumulative distribution function evaluated at the point of eŽj. , where eŽj. is the ordered version of e i so that eŽ1. F eŽ2. F ... F eŽ N . . The censored data points were not included when I constructed the probability plot. N is the sample size, the number of the uncensored observations. The cumulative distribution func-
Table 1 SAS printout of LIFEREG procedure Variable
df
Estimate
S.E.M.
x2
P)x
Intercept V RK RQ RK)V RQ)V Scale s Shape d
1 1 1 1 1 1 1 1
1.462 y0.0126 0.422 0.268 y0.0056 y0.0052 1.313 y1.894
0.0376 0.0011 0.0480 0.0466 0.00144 0.00128 0.0169
1514.07 141.87 77.31 33.02 15.12 16.71
0.0001 0.0001 0.0001 0.0001 0.0001 0.0001
Log-likelihood for gammasy14905.07. df, degree of freedom. The data used in this regression analysis were obtained from six patches from different oocytes injected with wild type Shaker B RNA, five patches from the mutant R377K and seven patches from the mutant R368Q. Each patch contained only one active channel. The data were collected under step voltages y20, 0, 20, 40, 60, 80 mV and the sample size of the data set was 8829 of which 2083 were right censored.
tion is defined as F Ž e . s H0e f Ž w .dw, where f is given by ˆ The Eq. Ž14. with k equal to the estimated value k. coefficient of correlation for F w eŽj. x and jrŽ N q 1. was also calculated.
3. Results An example of random openings of a channel under voltage clamp pulses and different first latencies is shown in Fig. 1. Fig. 2A is a histogram of first latency for a R368Q channel clamp to 20 mV. It shows that the distribution of the first latency is highly skewed to the left with most latencies occurring within a short time. If we take logarithm for the x-axis ŽFig. 2B., the histograms look more symmetric and we can see that the position parameter
Fig. 3. Likelihood ratio statistics Ž2)log-likelihood. for gamma parameters. The dashed lines indicate the value of the chi-squared distribution with 1 df at the 0.05 significance level. A: likelihood ratio statistics from fitting the data with the shape parameter k fixed at different values around the m.l.e. B: likelihood ratio statistics from fitting the data with the scale parameter g fixed at different values around the m.l.e.
X.M. Shao r Brain Research 770 (1997) 96–104
shifts to the left as the clamp voltage increases. This suggests that the accelerated failure time model which assumes that the effect of explanatory variables on the event-time distribution is multiplicative would be appropriate. By fitting the model ŽEq. Ž1.., we obtained the following regression function: ln Ž T . s 1.46y0.013Vq0.42RKq0.27RQ y 0.0056RK)V y 0.0052RQ)V
Ž 9.
Table 1 is the SAS printout giving the estimates of regression coefficients for explanatory variables and estimates of the parameters for the probability distribution of the error term. It also gives standard error and P-value for each
101
estimate. The m.l.e. of the gamma scale parameter s is 1.313 and the shape parameter d is y1.894. From these parameters, the gamma distribution parameters g s 1.443 and k s 0.279 were found Žsee reparameterization in Appendix A.. Table 1 also shows that all the regression coefficients are statistically significant Ž P-values - 0.0001. They are significant even with the consideration of simultaneous inferences with the Bonferroni approach w27x.. The regression coefficients for V, RQ and RK indicate that the voltage and the two mutations have effects on the first latency time. The minus sign of the coefficient for V shows that the slope of the voltage dependence of the latency time is negative. The values of the coefficients for
Fig. 4. A: probability plot for the log residuals Žsee the definitions in Eq. Ž7. and Eq. Ž8. in Section 2.. The ordinate is the estimated distribution function of the ordered residuals. B: residual plot for the fitted model Žsee Eq. Ž7. and Section 2..
102
X.M. Shao r Brain Research 770 (1997) 96–104
the terms RQ and RK in the regression function ŽEq. Ž9.. are large suggesting that two mutations have strong effects on latency time and the difference between these effects is small. And, surprisingly, the coefficient for R377K is larger than that for R368Q. The significant coefficients for the two interaction terms RK)V and RQ)V indicate that the voltage and the mutations affecting the first latency time do interact. In other words, the two mutations affect the voltage dependence of the opening process. Possibly, the charge moving process has been altered by the two mutations in the S4 region. The minus signs of the coefficients for the two interaction terms RK)V and RQ)V show that the mutations decrease the slope of the voltage dependence. Since the slope is negative, the mutations actually increase the absolute value of the slope of voltage dependence. From the values of the coefficients of the two interaction terms, we can see that the effects of R377K and R368Q on voltage dependence are similar. It implies that the activation mechanism cannot be explained solely by the electrostatic interaction. To evaluate the regression model, we consider the two essential components, the distribution of the random error term and the systematic dependence of the response variable on the explanatory variables, separately, as follows. We assessed the distribution of the data using the likelihood ratio tests with H 0 : g s 1 vs. H 1: g / 1. We did regressions by fixing the scale parameter at different values corresponding to the values of g shown in Fig. 3B. We then got the log-likelihood ratio statistics from each regression. Fig. 3B established by this method shows that H 0 should be rejected at the 5% significance level. Fig. 3A is constructed with the same procedure by fixing the shape parameter at different values. It shows that H 0 : k s 1 and k ™ ` should both be rejected. These results indicate that the data are not Weibull, gamma, exponential or log-normal distributed. From this result a question is naturally raised: are the data really following the three-parameter generalized gamma distribution? In other words, is the empirical distribution of the data sufficiently well described by a generalized gamma distribution? A probability plot for the residuals was made ŽFig. 4A.. The straightness of the estimated generalized gamma distribution of the ordered residuals vs. the probability is satisfactory and the coefficient of correlation is 0.996. There is no large departure from a straight line for most of the plot. To assess the adequacy of the regression function, residuals were computed by Eq. Ž7. in Section 2. The plot of the residuals vs. the fitted values ŽFig. 4B. shows the residuals spread evenly around zero. The upper bound of the residuals seems to decrease as the fitted value increases, which may be due to the fact that our censoring is always at 80 ms. We certainly do not have any observation bigger than 80. The y i are always less than lnŽ80. s 4.38, so for the larger fitted value xXb, Ž4.38 y xXb.rs would be smaller. Thus, the distribution of the residuals here does
not suggest a deviation from the regression function and the regression model is adequate and acceptable.
4. Discussion Parametric survival analysis has been applied to the first latency for single channel opening of wild type and mutant Shaker Kq channels. The application of the survival time regression procedure to single channel behavior is novel. By introducing a multiple regression model, we reveal some properties of gating kinetics which are difficult to detect and quantitatively measure by classical single channel analysis methods w29x. The effects of the experimental variables Žvoltage and mutant types. on the statistical properties of the first latency were summarized by a single regression function as Eq. Ž9.. Each different effect is quantified with the corresponding coefficient. The alternative survival time models, Weibull, exponential, ordinary gamma and log-normal, were evaluated by parametric inference under the generalized gamma model and they were all rejected by likelihood ratio tests. The goodness of fit and the adequacy of the regression function were assessed by the probability plot and residual plot. The regression results suggest: Ž1. First latency is voltage dependent. The minus sign of the regression coefficient for voltage indicates that the channel takes a shorter time to make the transition from channel closed to open under higher voltages. It implies that in the process of going from closed states to opening, in other words, the series of conformational changes that occur before the channel opens involve charge movement. Ž2. The two mutations have significant effects on the first latency. The positive regression coefficients for types R377K and R368Q tell us how much longer the mutant channels take from closed to open as compared to the wild type channel. The values of the coefficients show that R377K has a stronger effect on the first latency than R368Q rather than the other way around as expected w18x. Since the mutations in the S4 region alter the voltage sensitive process of channel activation in some way, it supports the hypothesis that the S4 region is a part of the voltage sensor w16,17,21,23,29,31x. Ž3. The regression coefficients for the two interaction terms show that the mutations R377K and R368Q have a significant effect on the voltage dependence of first latency. The latency times of the two mutants as a function of voltages have different slopes from the wild type channel, but the slope changes for the two mutants are similar. It implies that the mutations change the charge moving process but the voltage sensitive gating does not solely depend on the charge contents of the S4 region. This result supports the idea that the mechanism of activation is not due solely to an electrostatic interaction between the S4 positively charged amino acids and the membrane potential. These positive residues are not functionally equivalent
X.M. Shao r Brain Research 770 (1997) 96–104
w21x. Other properties of the residues apart from the charge, e.g., the length of the side chains or the interaction between these residues and the residues of other parts of the channel, may contribute to the gating mechanism of the channel as well. Ž4. The probability distribution of first latency may be well modeled by the three-parameter generalized gamma distribution. The Shaker Kq channel is thought to be a tetramer comprised of four identical subunits w20x and all four subunits have to be in open state for the channel to be open. The channel is closed if any one subunit is in closed state. If the simple kinetic model a
Close ° Open
Ž 10 .
b
for each subunit is correct and if the transition rates a and b are constants Ždo not vary with time., we would expect that the time it stays in any one state is exponentially distributed with parameter a or b respectively w4x. This means that the first latency from closed state to open state for each subunit Žthis is just for a single subunit from closed state to open state and it is not observable until the whole channel is open. should follow an exponential distribution with mean 1ra . If we assume the latency times for the subunits have a common distribution and the subunits gate independently, and if we further assume that the data obtained at high voltages where the contribution of the backward rate b values to first latency could be negligible w33x, the distribution of the time to channel opening would be the distribution of the maximum of the four-order statistics of exponential w26x. Then the cumulative distribution function would be 4
F Ž t . s Ž 1 y exp Ž ya t . . .
Ž 11 .
From this, we have the p.d.f.: f Ž t . s 4a exp Ž ya t . y12 a exp Ž y2 a t . q 12 a exp Ž y3 a t . i y 4a exp Ž y4a t . .
Ž 12 .
If the backward rate b values are not negligible, the probability expression becomes more complicated and has to be expressed in matrix form w4x. In addition to this complexity, recent experiments have shown that functional or structural interactions among subunits during the transitions from closed to open exist w2,33,34x. Furthermore, models like Eq. Ž10. for each subunit may be oversimplified. Each subunit may undergo more than one transition from closed state to open state w3x. Unfortunately, increasing the number of exponentials to account for these different closed states makes the estimation of parameters unreliable since the curve fitting would involve a large number of parameters. More fundamentally, whether the transitions in the gating process are Markovian is still controversial in channel biophysics w15x. From the regression analysis point of view, the three-parameter gamma distribution is flexible
103
in dealing with this complexity. We see from the residual analysis that the generalized gamma model is an acceptable approximation for the distribution. Our conclusions based on parametric survival analysis are qualitatively consistent with those obtained from conventional methods w29x. If we plot the empirical cumulative distribution functions of first latency for different mutants at different voltages, and plot the medians of first latency for different mutants at different voltages, we can show that first latency is voltage dependent and the two mutations have effects on the first latency and on the voltage dependence of the first latency. It is not clear, however, which effects are stronger and how strong they are. If we do multi-exponential curve fittings for different mutants at different voltages, we get an estimate for each combination of the mutant and voltage, but the interactions between the functional relations between first latency and the mutations, and between first latency and voltage are unclear. The parametric survival analysis method suggested in this paper is different from the conventional methods in the following aspects: Ž1. it provides quantitative descriptions of the relations between the response variable and each experimental treatment and the interactions between experimental variables; Ž2. it provides significance tests for those relations; Ž3. it deals with the censoring problems and thus it results in less bias in parameter estimation; Ž4. when the number of experimental treatments is large or when these experimental variables vary simultaneously, these multiple regression models can be extended very easily and can be standardized with recent powerful statistical analysis packages. It provides a simpler way to deal with large amounts of data before we have a sound biophysical understanding of the channel gating process to establish more analytical models. A number of different signals from nervous systems, e.g., the dwell time of single channel records, the interspike intervals of single neuronal activity and spontaneous or miniature postsynaptic events, excitatory postsynaptic potentials ŽEPSP. or inhibitory postsynaptic potentials ŽIPSP., can be described by certain kinds of survival time models. The application of survival analysis methodology, as illustrated in this paper, can be extended in the field of neuroscience.
Acknowledgements The author wishes to thank the referees for very helpful comments which led to substantial improvements, Dr. Diane Papazian for her support of this work and allowing the author to use the resources of her lab, Dr. Donald Guthrie for his advice on this project, and Mr. Calvin Cliff for his comments on the manuscripts. This work was supported by NIH Grant GM43459.
X.M. Shao r Brain Research 770 (1997) 96–104
104
Appendix A For computational reasons, reparameterization is done as follows. We consider W s lnŽT0 ., which has p.d.f. f Ž w. s
kk
'k G Ž k . exp
'k w y k )exp Ž wr'k .
.
Ž 13 .
Let Y s lnŽT . and xXbs m , then W s Ž Y y m .rs and Y has p.d.f. f Ž y. s
kk
'k sG Ž k . exp
'k Ž y y m .
rs y k )exp Ž Ž y y m . r Ž s'k . .
Ž 14 .
and the corresponding distribution of T s expŽ Y . has p.d.f. of the form of Eq. Ž2., with g s 1rŽ s 6k . and t s expŽ m .rk 1r g or, equivalently, s s 1rŽg 6k . and m s lnŽt . q lnŽ k .rg . References w1x L.J. Bain, M. Engelhardt, Probability of correct selection of Weibull versus gamma based on likelihood ratio, Commun. Stat. Theory Methods A9 Ž1980. 375–381. w2x F. Bezanilla, E. Perozo, D.M. Papazian, E. Stefani, Molecular basis of gating charge immobilization in Shaker potassium channels, Science 254 Ž1991. 679–683. w3x F. Bezanilla, E. Perozo, E. Stefani, Gating of Shaker Kq channels. II. The component of gating currents and a model of channel activation, Biophys. J. 66 Ž1994. 1011–1021. w4x D. Colquhoun, A.G. Hawkes, The principles of the stochastic interpretation of ion-channel mechanisms, in: B. Sakmann, E. Neher ŽEds.., Single Channel Recording, Plenum Press, New York, 1983, pp. 135–174. w5x D. Colquhoun, F.J. Sigworth, Fitting and statistical analysis of single channel records, in: B. Sakmann, E. Neher ŽEds.., Single Channel Recording, Plenum Press, New York, 1983, pp. 191–263. w6x D.R. Cox, Tests of separate families of hypotheses, Proceedings of the 4th Berkeley Symposium, vol. 1, 1961, pp. 105–123. w7x D.R. Cox, Further results on tests of separate families of hypotheses, J. R. Stat. Soc. B24 Ž1962. 406–424. w8x R. Dumonceaux, C.E. Antle, Discrimination between the log-normal and the Weibull distributions, Technometrics 15 Ž1973. 923–926. w9x V.T. Farewell, R.L. Prentice, A study of distributional shape in life testing, Technometrics 19 Ž1977. 69–75. w10x R.E. Greenblatt, E. Blatt, M. Montal, The structure of the voltagesensitive sodium channel. Inferences derived from computer-aided analysis of the Electrophorus electricus channel primary structure, FEBS Lett. 193 Ž1985. 125–134. w11x H.R. Guy, P. Seetharamulu, Molecular model of the action potential sodium channel, Proc. Natl. Acad. Sci. USA 83 Ž1986. 508–512. w12x O.P. Hamill, A. Marty, E. Neher, B. Sakmann, F.J. Sigworth, Improved patch-clamp techniques for high resolution current recording from cells and cell-free membrane patches, Pflugers Arch. 391 ¨ Ž1981. 85–100. w13x J.F. Lawless, The generalized gamma distribution, in: Statistical Models and Methods for Lifetime Data, John Wiley and Sons, New York, 1982, pp. 237–252.
w14x H.R. Leuchtag, Long-range interactions, voltage sensitivity, and ion conduction in S4 segments of excitable channels, Biophys. J. 66 Ž1994. 217–224. w15x L.S. Liebovitch, J.M. Sullivan, Fractal analysis of a voltage-dependent potassium channel from cultured mouse hippocampal neurons, Biophys. J. 52 Ž1987. 979–988. w16x E.R. Liman, P. Hess, F. Weaver, G. Koren, Voltage-sensing residues in the S4 region of a mammalian Kq channel, Nature 353 Ž1991. 752–756. w17x D.E. Logothetis, S. Movahedi, C. Satler, K. Lindpaintner, B. NadalGinard, Incremental reductions of positive charge within the S4 region of a voltage-gated Kq channel result in corresponding decreases in gating charge, Neuron 8 Ž1992. 531–540. w18x D.E. Logothetis, B.F. Kammen, K. Lindpaintner, D. Bisbas, B. Nadal-Ginard, Gating charge differences between two voltage-gated Kq channels are due to the specific charge content of their respective S4 regions, Neuron 10 Ž1993. 1121–1129. w19x G.A. Lopez, Y.N. Jan, L.Y. Jan, Hydrophobic substitution mutations in the S4 sequence alter voltage-dependent gating in Shaker Kq channels, Neuron 7 Ž1991. 327–336. w20x R. MacKinnon, Determination of the subunit stoichiometry of a voltage-activated potassium channel, Nature 350 Ž1991. 232–235. w21x D.M. Papazian, L.C. Timpe, Y.N. Jan, L.Y. Jan, Alteration of voltage-dependence of Shaker potassium channel by mutations in the S4 sequence, Nature 349 Ž1991. 305–310. w22x V.B. Parr, J.T. Webster, A method for discriminating between failure density functions used in reliability predictions, Technometrics 7 Ž1965. 1–10. w23x E. Perozo, L. Santacruz-Toloza, E. Stefani, F. Bezanilla, D.M. Papazian, S4 mutations alter gating currents of Shaker Kq channels, Biophys. J. 66 Ž1994. 345–354. w24x R.L. Prentice, A log-gamma model and its maximum likelihood estimation, Biometrika 61 Ž1974. 539–544. w25x J.A. Rice, Probability plot, in: Mathematical Statistics and Data Analysis, Wadsworth and Brook, Pacific Grove, CA, 1988, pp. 292–299. w26x J.A. Rice, Extrema and order statistics, in: Mathematical Statistics and Data Analysis, Wadsworth and Brook, Pacific Grove, CA, 1988, pp. 94–97. w27x J.A. Rice, The problem of multiple comparisons, in: Mathematical Statistics and Data Analysis, Wadsworth and Brook, Pacific Grove, CA, 1988, pp. 404–408. w28x SAS Institute Inc., The LIFEREG procedure, in: SASrSTAT User’s Guide, vol. 2, version 6, 4th ed., SAS Institute Inc., Cary, NC, 1989, pp. 997–1017. w29x X.M. Shao, D.M. Papazian, S4 mutations alter the single-channel gating kinetics of Shaker Kq channels, Neuron 11 Ž1993. 343–352. w30x E.W. Stacy, A generalization of the gamma distribution, Ann. Math. Stat. 33 Ž1962. 1187–1192. w31x W. Stuhmer, F. Conti, H. Suzuki, X. Wang, M. Noda, N. Yahagi, H. Kubo, S. Numa, Structural parts involved in activation and inactivation of the sodium channel, Nature 339 Ž1989. 597–603. w32x B.L. Tempel, D.M. Papazian, T.L. Schwarz, Y.N. Jan, L.Y. Jan, Sequence of a probable potassium channel component encoded at Shaker locus of Drosophila, Science 237 Ž1987. 770–775. w33x W.N. Zagotta, R.W. Aldrich, Voltage-dependent gating of Shaker A-type potassium channels in Drosophila muscle, J. Gen. Physiol. 95 Ž1990. 29–60. w34x W.N. Zagotta, T. Hoshi, R.W. Aldrich, Gating of single Shaker potassium channels in Drosophila muscle and in Xenopus oocytes injected with Shaker mRNA, Proc. Natl. Acad. Sci. USA. 86 Ž1989. 7243–7247.