[51 ]
STATIONARYSINGLE-CHANNELANALYSIS
[51 ] S t a t i o n a r y S i n g l e - C h a n n e l
729
Analysis
By MEYER B. JACKSON Introduction The opening and closing of ion channels in artificial and cell membranes generates a signal of membrane current that varies in a stepwise manner with time. The stochastic manner in which this signal varies depends on the underlying channel-gating mechanisms. Thus, analysis of records of sin#e-channel current can provide much insight into the mechanisms by which channels open and close. Because channel gating is a stochastic process, the methods of analysis are statistical. This chapter describes a number of useful methods of kinetic analysis of stationary single-channel data. In a stationary system, the properties of the ensemble do not change with time; the mean probabilities of occupancy of specific configurations remain constant. Such a system must satisfy a condition of detailed balance, in which the frequency of transitions toward a particular configuration is equal to the frequency of transitions out. In nonstationary systems, detailed balance is not satisfied, and the mean probabilities of occupancy of a particular configuration change with time. The practical implication is that if the data are stationary, one simply can collect the data continuously for analysis. Nonstationary conditions are often achieved by perturbing the system, most commonly with changes in voltage. Methods for the analysis of nonstationary sin#e-channel data are presented elsewhere (see [52], this volume). Sin#e-channel kinetic analysis can be divided into three distinct steps. The first step involves the theoretical derivation of predictions of kinetic models. Examples of such theoretical predictions include the probability distribution function for open times and the joint probability density function for the lifetimes of two consecutive events. The second step of single-channel kinetic analysis is to extract something from real data that can be compared to theoretical predictions. These operations generally start with raw sin#e-channel data and generate open time and closed time histograms, or more complicated distillates of channel behavior. One thus obtains an experimental result in a form that can be directly compared with a theoretical prediction. This comparison, which most often involves some form of curve fit, is the third and final step of sin#e-channel kinetic analysis. The insights into channel gating processes are obtained through a comparison of distillates of sin#e-channel data with the theoretical predictions of models. METHODS IN ENZYMOIX~Y, VOL. 207
Colaytisht © 1992 by Academic Prem, Inc. All ri$1mof relm)duetion in any form reserved.
730
DATA STORAGE AND ANALYSIS
[51]
Single-channel analysis is very software intensive, and a large number of computer programs have been developed in many laboratories and by many companies. Many software packages are commercially available. The number of programs is so great, and the purposes so varied, that it is not possible to evaluate and compare them. The emphasis is rather on the specific tasks for which the programs are written. It is hoped that after reading this chapter, the reader will know what to expect from the software, how to use it intelligently, and what are the possibilities for developing programs for more specialized needs.
S u m m a r y of Kinetic T h e o r y The most widely used theories of single-channel kinetics treat channel gating as a finite Markov process. The complex trajectory of a channel through its protein conformational space is modeled as a series of jumps or transitions between states within which the conductance is relatively constant. Such models for channel gating then consist of a finite number of states that are either open, closed, or have an intermediate conductance. Each state can then interconvert to another with a finite rate. This is represented in a kinetic scheme by an arrow between the two states that can interconvert. When a rate is zero in both directions, then the two states can be considered unconnected. Within this framework, the relevant questions about a channel-gating mechanism are the following: (l) How many states are there with a given conductance? (2) How are they connected? (3) What are the rate constants for each step? This description defines the term "state," which is used extensively in the ensuing discussion. Relative to the time scales accessible in a singlechannel experiment (> 20/~sec), the transitions between states are instantaneous, whereas the states themselves are stable enough to observe. According to the Markov assumption, the probability of a particular transition depends only on which state the system is in. This model for channel gating corresponds to a specific physical picture which at the present time appears to be consistent with essentially all experiments on channels. One of the most useful results of the kinetic theory of channel gating provides a method of estimating the number of states with a particular conductance. Consider a membrane system with a single channel having n open states and m closed states. The open time and closed time distributions are represented by sums of n and m exponentials, respectively:
Po(t) = ~ x, eU~
(1 a)
[51]
STATIONARYSINGLE-CHANNELANALYSIS
7 31
m
Pc(t) = ~ Yj e vj
(lb)
where Po(t) and Pc(t) represent the probabilities of open times or closed times lasting a time t or longer.t.2 The uppercase symbols are used for the probability distribution, whereas the probability density function, obtained by differentiating the distribution, is denoted by lowercase letters. The parameters Pi, vj, x~, and yj depend on the rate constants and on the details of how the open and closed states are connected in the gating scheme. However, the form of Eqs. (1 a) and (1 b) does not depend on these details and is independent of the parameters and connectivity of the gating scheme. In effect, when one determines the number of exponentials that are needed to fit an observed distribution, one has "counted" the number of states with a particular conductance. More precisely, this method provides a lower bound to the number of such states, since one is always faced with the possibility that some of the terms of Eqs. (la) and (lb) are too small, too rapid, or too slow to detect. It is also possible that the system has two time constants that are too similar to be resolved. For even three closed or open states the values of the time constants and coefficients in Eqs. (la) and (lb) can be very complicated functions of the rate constants for transitions between states. Thus, for complex gating mechanisms it is often very difficult to go beyond state counting to a quantitative determination of specific rate constants. With complex gating behavior, it is often so arduous to determine the rate constants that, depending on the questions at hand, it may be prudent to focus on simpler quantities such as mean open time, mean closed time, or mean burst duration. Equations (la) and (lb) are examples of one-dimensional probability distributions; they are functions of a single experimental variable. One can also make use of higher dimensional probability distributions and density functions. For a closed time probability density function pc(to) and an open time probability density function po(to), where tc and to are independent (no correlations), the joint, two-dimensional probability density function for observing both to and tc is the product po(to)pc(tc). If the values of to and tc are correlated, then the joint probability density function for consecutive open and closed times is no longer simply the product of the two one-dimensional probability density functions. This is a qualitative difference, which can be used to estimate the number of distinct gating pathways and to make distinctions between these pathways. In general, correlation analI D. Colquhoun and A. G. Hawkes, Proc. R. Soc. London, Set. B 199, 231 (1981). 2 D. Colquhoun and A. G. Hawkes, in "Single-Channel Recording" (B. Sakmann and E. Neher, eds.), p. 135. Plenum, New York, 1983.
732
DATA STORAGE AND ANALYSIS
[51 ]
ysis is used more to answer questions about the basic topology or connectivity of a channel-gating mechanism, rather than to estimate the values of specific rate constants. One important example of correlation analysis is highly analogous with state counting described above, in that it provides an estimate of the minimum number of gating pathways of a channel. This arises from the theoretical result that the number of terms in the covariance function of two event lifetimes as a function of the number of intervening time intervals is equal to one less than the number of distinct pathways between the open and closed states. The number of pathways is defined as the minimum of the number of closed states that can open and open states that can close. 3,4 Thus, when one determines the number of exponentials that are needed to fit an observed covariance function, one has "counted" (or, more precisely, estimated a lower bound to) the number of gating pathways. Correlation and covariance functions can be derived with higher dimensions (i.e., with more time intervals). However, there is no additional information about the underlying channel-gating mechanism in these higher dimensional functions. 3 Thus, they have never been used in singlechannel analysis. As the number of states of a model increases, the number of kinetic parameters and the number of topological possibilities for connecting the states increase. Singie-channel analysis is limited in the number of parameters that can be determined, and in the power to discriminate models, a,4 Nevertheless, a great deal of information can be obtained purely from statistical analysis of single-channel records. To go beyond the limits circumscribed by single-channel kinetic theory, it is necessary to combine the stationary kinetic analysis with other approaches, such as nonstationary kinetic analysis and dose-response behavior. Idealization of Single-Channel Data The first step in single-channel kinetic analysis is the idealization of the data. Raw single-channel data represent the superposition of background noise and stepwise changes in channel current. When the data have been idealized, the background noise has been eliminated and only the stepwise changes caused by channel opening and closing remain. With an ideal channel record, no information is lost by replacing the current versus time signal with a two-column list consisting of a time and a current for each 3D. R. Fredkin, M. Montal, and J. A. Rice, in "Proceedingsof the BerkeleyConferencein Honor of Jerzy Neyman and Jack Kiefer" (L. M. Leeam and R. A. Olshen eds.), Vol. 1, p. 269. Wadsworth,Belmont, California, 1985. 4D. Colquhoun and A. G. Hawkes,Proc. R. Soc. London B 230, 15 (1987).
[51 ]
STATIONARY SINGLE-CHANNEL ANALYSIS
733
interval at a particular current level. Idealization is then the reduction of a single-channel data record to such a list. The ease with which single-channel data are idealized is critically dependent on the ratio of the single-channel current amplitude to the background noise. For noisy data, no ideal idealization method exists. Idealization can, and often is, done by hand, or with a computer in a semiautomatic, event-by-event mode. In this mode the computer will search the data for threshold crossings and then display the event, indicating guesses for the precise time of the transition and the levels. A person running the program is then given the opportunity to accept, reject, or modify the event and the initial guesses. Examination of channel current data event by event is essential with data from a new type of experiment where the basic properties are unknown. When one attempts an event-byevent analysis, even of fairly good data, one encounters ambiguous events. There is always a chance that a random background fluctuation will occur that is comparable in size to the single-channel current. For poor data where the channel current is only 3 or 4 times the peak-to-peak background noise, such fluctuations will occur frequently enough to make idealization very difficult. Ambiguities arise more often for brief intervals. 5 Analyzing records with thousands of events in an event-by-event mode is simply too time consuming to be useful if a systematic study is undertaken. A fully automated computer program for idealization is then required. Such programs are available from Axon Instruments (Foster City, CA) and Instrutech (Elmont, NY), as well as from other companies. Most programs allow the user to toggle between two modes, one of which is a semiautomated, event-by-event analysis by the operator. The other mode is fully automated. This toggling is an essential feature, because it is necessary to check the performance of these programs. The lists of numbers produced by a fully automated idealization program can be very deceptive. Glitches in data, such as those caused by power surges, seal breakdown, or bubbles in the perfusion system, will produce crossings of the threshold used to signify channel-gating transitions. Such events are easily recognized as spurious by the human eye, owing to their nonstepwise character. In an event-by-event analysis, an investigator can almost always avoid such mistakes, but to most computer programs, all threshold crossings look the same. Another source of false transitions has been described by Colquhoun. 6 If the records contain subconductance states (or for that matter infrequent appearances of another smaller channel) that are near the threshold for D. Colquhoun and F. J. Sigworth, in "Single-Channel Recording" (B. Sakmann and E. Neher, eds.), p. 191. Plenum, New York, 1983. 6 D. Colquhoun, Trends Pharmacol. Sci. 9, 157, (1988).
734
D A T A STORAGE A N D ANALYSIS
[51]
event detection, a single long dwell time at this level will produce a large number of very brief events in the output list. To avoid these and other errors, the raw data should be spot checked, and the output lists should be examined for anomalous behavior. The performance of idealization programs should be compared for different values of the event detection parameters. The amplitude and lifetime distributions should be checked visually with every data set for unusual peaks or components. An anomalous excess of points somewhere in the distribution can almost always be traced to a problem in the data, and editing these events will improve the quality of analysis. As Colquhoun succinctly put it, "single-channel analysis costs time, ''6 and the temptation to turn a software package loose, and publish the results without careful inspection and monitoring along the way, should be resisted. Because idealization is the stage of analysis where the noise level will have its primary impact, the bandwidth and sampling frequency should be seIected to optimize the quality of idealization. The Nyquist criterion of digitization at twice the corner frequency of the low-pass filter is inappropriate for single-channel analysis. 5 Most investigators have determined empirically that one should digitize at 8 to l0 times the corner frequency of the filter. At an early stage in a study, a data set should be analyzed with different bandwidths and digitization frequencies to see where problems develop. When a channel exhibits gating activity on a rapid time scale, the bandwidth must be widened to observe these processes, but if the noise increases enough to impede the idealization, then better data form the best solution. Curve-Fitting Lifetime Distributions Once data have been idealized it is generally straightforward to produce a list of all the open times or all the closed times. At this point one is ready to find a sum of exponentials that best describes the relative frequencies of different event lifetimes. One can either fit to a histogram or work with the list of lifetimes directly using likelihood maximization. The first goal of this procedure is to estimate the number of exponential terms needed to describe the open time and closed time data, and thus obtain a count of the minimum number of states. When the number of open or closed states is one or two, it may also be useful to determine the parameters to relate them to specific rate constants of the channel-gating process.
Fitting to Histograms From a list of state durations, one can construct a histogram, N(ti, ti+~), where ti and tj+~ are times that mark the boundaries of the bins. When the bins are narrow, then p(t)(tl+~ - ti) becomes a viable theoretical estimate of
[ 51 ]
STATIONARY SINGLE-CHANNEL ANALYSIS
735
N(ti, ti+l). The sum of squares difference between these two quantities can can then be used as an index of the quality of the fit. Fitting then proceeds by minimizing the sum of squares with respect to the parameters ofp(t). A function minimization program can then be used to carry out this optimization. In most programs, the sum of the coefficients in the probability density function Ix,. or yj of Eqs. (1 a) and (lb)] is constrained to be one, to reduce the number of free parameters. Minimizing ~a has also been suggested for curve fitting. 5 A likelihood maximization method can also be used for fitting to histograms (see below). A Z2 goodness-of-fit test can be used to evaluate the fit and the improvement in the fit resulting from additional exponential terms. The number of exponential components in a distribution is estimated by comparing the goodness-of-fit for different numbers of terms. With a small data set, one has fewer events per bin, and the statistical fluctuations in N become large. The histogram will then appear noisy, the curve fitting will not yield good parameter estimates, and low amplitude exponential terms will be difficult to detect. If an attempt is made to widen bins to remedy this problem, then p(t)(t,+l - ti) will be a poor approximation for the predicted value of N, and this expression should be replaced by the integral ofp(t) over each bin, or P(O - P(t~+1)- When bins are too wide, there will be too few data points to determine the fitting parameters accurately. Thus, fitting to histograms is problematic with small data sets. Likelihood Maximization For n event lifetimes, ti, and a probability density function p(t), the likelihood is
L = f i P(ti)
(2)
A curve fit proceeds by variation of the parameters ofp(t) to maximize L. 5,7 In practice, the method of likelihood maximization is carded out with the logarithm of the likelihood, because L itself exhibits such extreme values that its computation can be difficult. The likelihood ratio test s,9 can then be used to estimate the significance of the improvement of a fit by an additional exponential term. This method does not involve binning data, as does the sum of squares method with histograms. It is thus advantageous with small data sets. The time required to compute L increases with the size of the data set, and the scatter of N in histograms becomes smaller, so that with large data sets it is 7 R. Horn and K. Lange, Biophys. J. 43, 207 (1983). 8 C. R. Rao, "Linear Statistical Inference and Its Applications." Wiley, New York, 1973. 9 M. B. Jackson, Biophys. J. 49, 663 (1986).
736
DATA STORAGE AND ANALYSIS
[51]
better to use histograms. However, the method of likelihood maximization can still be used with histograms. The above likelihood function [Eq. (2)] then becomes ~° i
L = I-[ p(t,)',
(3)
where ni is the number of events in bin L This can be computed much more rapidly for a very large event lifetime list than it can using Eq. (2) and is thus a convenient function to use in the fitting of binned histogram data. With large data sets maximum likelihood and least-squares techniques become equivalent, and there are no compelling reasons for preferring one method over the other. Out-of-Range Events
There will always be some events which are too brief to observe. There will also be events with lifetimes right at the limit of detection, such that the frequency with which such events occur is not accurately determined. In addition, it is often important to exclude from analysis the longest times, which can have anomalous effects on curve fitting. When using the sum of squares method, the correction for out-of-range events is trivial. One simply does not include the N values for bins that are out of range. For likelihood maximization, this procedure is not correct; simply ignoring such events is equivalent to assigning zero values to the experimentally observed frequencies. To correct for out-of-range events, one must fit to a probability density function which includes the condition that the interval is in the specified range: Thus the likelihood for a single event of duration t becomes L = p(t)
p(t) dt = p(t)/[P(t,) - P(tl)]
(4)
where ts and tl are the short- and long-time cutoffs, respectively. When the adaptation of likelihood maximization to binned data [Eq. (3)] is corrected in this fashion, one obtains the following index for the quality of fit)°,l~: ~ P ( t , ) - P(t,+,)~ ln(L)---- ~ n/In i p - - ~ s ) - - - ~ l) j
(5)
Many investigators increase the bin width exponentially to compensate for fewer events at longer durations (tt + 1 = Oti, where t~ is approximately 1.25). This provides constant binwidth on a logarithmic time axis. For event lifetime lists of more than a few thousand, this implementation of l0 F. J. Sigworth and S. M. Sine, Biophys. J. 52, 1047 (1987). H O. B. McManus, A. L. Blatz, and K. L. Magleby, PfluegersArch. 410, 530 (1987).
[51 ]
STATIONARYSINGLE-CHANNELANALYSIS
737
likelihood maximization [Eq. (5)] is at present the method of choice in fitting to a sum of exponentials. Careful attention must be given to the precise value of the short-time cutoff in these expressions [Eqs. (4) and (5)]. When lifetimes are determined as multiples of a sampling frequency, then the cutoff must be n + 0,5 times the sampling frequency,9 where n is usually an integer between 3 and 6. This form of cutoff is necessary because the true lifetimes are effectively rounded up or down by carrying out the idealization procedure on digitized data. The true cutoffis then at the rounding point. When the idealization procedure employs an interpolation procedure for estimating the lifetimes more accurately, then there is no longer such a restriction on the cutoff time. However, it is then important to evaluate the interpolation procedure carefully to ascertain that it interpolates uniformly. Digitization and binning can introduce another error called the sampling promotion error when the fastest time constant approaches the bin width or digitization time interval.t t,~2 Additional Corrections for Brief Events The above discussion presented a correction for the primary error produced by out-of-range events in the curve fitting of probability density functions. However, intervals that are too brief to detect can lead to additional errors, and there has been much analysis of how missed intervals influence single-channel analysis. The most commonly encountered problem is that ifa brief open (closed) time is missed, then the two adjacent closed (open) times appear as one closed (open) event with a lifetime equal to the sum of the three event lifetimes. A simple and effective correction can be used for missed closures when there is a single exponential distribution of open times, t3,t4 The open times will remain exponentially distributed despite the missed brief closures, but with a lengthened apparent time constant. The first step in the correction is to estimate the number of missed brief closures from the zero-time intercept of the best-fitting closed time distribution. The corrected time constant, z ~ , is then related to the observed value, ~o~, by the expression z~or = ~o~ No~/(No~ + Nm~a)
(6)
where NoB is the number of observed open times and N=ma is the number of missed brief closures. Treatment of the analogous case for closed times interrupted by brief openings is essentially the same. ~2s. Sineand J. H. Steinbach,J. Physiol. (London) 373, 129 (1986). t3 E. Neher,J. Physiol. (London) 339, 663 (1983). t4 F. Sachs,J. Neil, and N. Barkakati,PfluegersArch. 395, 331 (1982).
738
DATA STORAGE AND ANALYSIS
[51]
For multiple exponential distributions the problem becomes more complicated. Equation (6) is still correct for the mean open time, but it is difficult to extend this correction to each time constant of the open time distribution. If one can demonstrate by correlation analysis (discussed below) the extent to which the brief closures are associated with the different open states, then Eq. (6) can be implemented by replacing No~ with the number of observed openings corresponding to a particular component, Zo~ with the time constant of that component, and N . ~ d with the number of missed closures estimated to be associated with that component. The general problem of missed events has been treated systematically by Roux and Sauve.~5 Blatz and Magleby ~6have also presented some useful results.
Plotting Exponential Distributions In a simple linear plot of a lifetime distribution, it is often difficult to recognize different exponential components (Fig. 1A). For this reason transforms have been developed to facilitate visual inspection. Semilogarithmic plotting provides some improvement by transforming an exponential to a line. However, when the time constants differ by more than a factor of 100 it is difficult to use a linear time axis. A particularly useful means of displaying such data is in a modified log plot (Fig. IC), in which bin size increases exponentially with time, and ni is not sealed to the size of the bin. I°,1~ Further transforming the y axis by taking the square root l° produces a variance-stabilized plot, that is, the deviation or scatter in the data becomes independent of the y value (Fig. 1D).
Fitting to Data from System with Many Channels The theory employed above is strictly valid only for a system with one channel? ,2 When a system has many channels this theory is no longer valid. When channel openings are well spaced from one another, and simultaneous openings are very infrequent, it has been assumed) 7 and can be shown, 18 that the one-channel kinetic theory provides a reasonable approximation to certain aspects of the gating process. The problems arising in the analysis of a many channel system can be illustrated with a few examples. When examining a closed time interval in a system known to have more than one channel, there is no way of knowing whether the channel that closed to begin the interval is the same 15 B. Roux and R. Sanve, Biophys. J. 48, 149 (1985). 16 A. L. Blatz and K. L. Magleby, Biophys, J. 49, 967 (1986). 17 D. Colquhoun and B. Sakmann, J. Physiol. (London) 369, 501 (1985). is M. B. Jackson, Biophys. J. 47, 129 (1985).
[51 ]
STATIONARY SINGLE-CHANNEL ANALYSIS ,
,
i
,
,
,m.
,
,
.
,
739
,
_1 30Q
lO 5
A
'
250
104
. . 200
~o 3
150
B
I0 2
~' ~oo io I 50 0
. I 0
'T. I,Y::7";S-(:-"-,:--'!r=--; .... 10
20
30
40
SO
60
70
80
90
i
100
-5 10
10- ¢
~0-3
Time (rnsec)
T~me
500
10-I
10
0
i
10
1
10
2
10
(sec)
300 25O 200
C
2SO ¢
10-2
D
m 150
200
100 ~ " 150 100 50 o
I
-5 10
-10 4
i
I
-3 1Q
-2 10
i
-1 10
Time (~ec)
i
10
i
0
10
I
1
10 2
I
10
-5 I0
-4 10
-3 10
-10 2
-101
10
0
10
I
I
10
2
I0
llme (se~)
FIG. 1. Four representations of a dwell time distribution with two exponential components (from Ref. 10). A total of 5120 random numbers were generated according to a distribution with time constants of 10 msec (70% of the events) and 100 msec (30%) and binned for display as histograms in the lower part of each figure. The theoretical probability density functions for each component (dashed curves) and their sum (continuous curve) are superimposed with each curve. In each part of the figure the upper panel shows the absolute values of the deviation of the height of each bin from the theoretical curve, with dashed curves showing the expectation value of the standard deviation for each bin. The upper panels were plotted with expansion factors of 2.1, 5.4, 3.1, and 4.9, respectively. (A) Linear histogram. Events were collected into bins of 1 msec width and plotted on a linear scale. The 100-msec component has a very small amplitude in this plot. (B) Log~log display with variable width (logarithmic) binding. The number of entries in each bin was divided by the bin width to obtain a probability density in events/second, which is plotted on the ordinate. (C) Direct display of logarithmic histogram. Events were collected into bins of width 6x = 0.2 and superimposed on the histogram as the sum of two functions, as in Eq. (11) of Ref. 10. (D) Square-root ordinate display of a logarithmic histogram, as in (C). Note that the scatter about the theoretical curve is constant throughout the display.
740
DATA STORAGE AND ANALYSIS
[51]
channel that opened to end the interval. This produces a systematic error in closed time distributions. When two channel openings overlap, there is no way of assigning the open time. Neglecting such events produces a systematic error in open time distributions. For a system with N identical and independent channels, with a singlechannel closed time probability density function of the form Xm vjyje-Vl [the derivative of Eq. (1 b)], the closed time probability density function is ~8 P~rc(t) = (:~ vyj e-V1)[:~(y./vj) e - ' d ~-1 + ( N - 1)(:~ Y1 e-"J)2[:~(Y'/Vi) e-,i],~-2
(7)
When the number of channels is known, Eq. (7) can be fitted to data by likelihood maximization, as has been demonstrated for a stretch-activated channel. 19With a single-channel open time probability density function of :~/zixi e--~'~[the derivative of Eq. (1 a)], the probability density function for isolated channel openings is is
t,'(t) -- r. x,l,,
[ (y/vj) (y/vj) j
(8)
where a prime is used to distinguish the probability for isolated openings from the true single-channel open-time probability density function. A fit with Eq. (8) requires input of parameters obtained from a fit to closed times, and knowledge of N. When N is difficult to determine, as is generally the case when N is large, then one can approximate N - 1 as N, which leads to p' ( t ) ~- ('2 x,l~, e-.~)[:~(z/coj) e--'o~]/~(z./o~j)
(9)
where the parameters o)j and zj are taken from a fit of :~ (ojzj e-",t to the closed times from the same data set) ° The results of such a fit are shown in Fig. 2. When the frequency of channel opening is low, there is little deviation of the best fitting one-channel open time distribution from the observed histogram of isolated open times. With data from a patch where the channel opening frequency is high, the observed histogram of isolated open times departs substantially from the best-fitting distribution, indicating that the distribution of isolated openings is a poor representation of the kinetics of closure of a single channel. 2° The similarity between the two functions (the solid curves of Fig. 2A,B) obtained by fitting Eq. (9) indi~9F. Guhary and F. Sachs, J. PhysioL (London) 363, 119 (1985). 20 M. B. Jackson, £ Physiol. (London) 397, 555 (1988).
[51 ]
lOOO
STATIONARYSINGLE-CHANNELANALYSIS A
741
B
N 100
I0 [
2
l 4
i 6 t
(msec)
I
/
i
8
I0
12
i
I
I
I
I
I
L
2
4
6
8
lO
12
d
t (msec)
FIG. 2. Comparison of cumulative distributions of lifetimes of isolated openings (doUed curve) with the best-fitting single-channel open time probability density function (smooth curve), obtained by curve fitting with Eq. (8) (from Ref. 20). In (A) the level of channel activity was low, and there were almost never simultaneous channel openings. In (B) the level of channel activity was high; most channel openings were interrupted by other channel openings, and isolated openings were relatively rare. The data were obtained from cultured mouse muscle using 2 # M carbachol. The difference in channel activity between the two patches reflects primarily differences in receptor density. Note that whereas the isolated open time distributions are very similar. The parameters for the best fits in (A) were #t = 6.0 m see-I, #2 = 0.24 msec-t, x~ = 0.26 for open times, with a frequency of channel opening of 0.029 msec- L In (B) the parameters were #~ = 5.5 msec-~, g2 = 0.021 msec-t, and xt = 0.30, with a frequency of channel opening of 0.21 msec- ~.
cates that with both sparse and crowded channel data, the same kinetics operates at the single-channel level. When fitting to these equations using likelihood maximization, it is important to incorporate the out-of-range event correction in an analogous manner to that described for Eq. (4). The integrations are lengthy but straightforward. Equation (9) can be used to analyze open times in very crowded data sets, provided that at least several hundred isolated openings are observed. When the frequency of channel opening is so high that the slowest time constant of the observed dosed time distribution is shorter than the slowest time constant of the single-channel open time distribution, then the parameters determined by fitting Eq. (9) will have large standard errors. 2° It will then be necessary to carry out such fits on many data sets to obtain accurate estimates of parameters.
742
DATA STORAGE AND ANALYSIS
[51]
Burst Analysis
Many channels exhibit "bursting behavior," in which openings occur in clusters. The openings within such a burst are separated by short closed time intervals. A great deal of jargon has arisen around this phenomenon. Terms such as "flickers, .... gaps," and "nachschlags" refer to various types of events in a "bursty" channel record. In the case of chemically activated channels, the bursts of openings are generally thought to be repeated openings of the same liganded receptor-channel complex. This mechanistic interpretation has provided a framework for a useful simplified approach to the kinetic analysis. Treating ligand binding and channel gating as separate steps of a kinetic scheme gives the following model: R + nA .-~- RA,,~'~R*A,~
(10)
where A represents the ligand, R the receptor, and R* indicates that the channel is open. A gap or flicker is a sojourn in state RA, between two openings. The mean burst duration, tb, mean gap duration, tg, and mean number of gaps per burst, n8, can be derived for this model and solved for the kinetic parameters to give2! a = (n 8 + 1)/(t b
-
-
nets)
k_~ = I/[2ts(n s + 1)1 fl = 2 k _ l n s
(1 1)
It then remains to estimate tb, ns, and ts from the data. The mean gap duration is best determined as the time constant of the fast component of the best-fitting closed time distribution. By inspecting a closed time distribution, one can define an operational cutoff time and assign all shorter closed intervals as gaps. One can then reconstruct bursts and determine burst duration and number of gaps. Although the assignment of all closures shorter than a cutoff as gaps may seem arbitrary, when there are two time constants of a closed time distribution differing by more than a factor of 10, then the number of incorrect assignments becomes insignificant. This type of burst analysis has been used in the study of both acetylcholine~7and ;,-aminobutyric acid (GABA)22 receptor channel data. Burst analysis is particularly useful with data showing frequently occurring subconductance states, or with which idealization is difficult owing to a low signal-to-noise ratio. By including subeonductance states as part of a burst, one can avoid the problem of how to idealize these events, which would require treating the subconductances as separate states. 2~ D. Colquhoun and A. G. Hawkes, Philos. Trans. R. Soc. London, Ser. B300, 1 (1982). 22 j. Bormann and D. E. Clapham, Proc. Natl. Acad. Sci. U.S.A. 82, 2168 (1985).
[51 ]
STATIONARY SINGLE-CHANNEL ANALYSIS
743
Burst analysis has the advantage of being simple and the disadvantage of being too simple. Kinetic models very different from Eq. (10) can generate bursting behavior. Equation (1 1) may not be valid for these models. If there are other pathways leading away from state RAn, such as a desensitization step, 23 then what is determined as k_ 1 will actually be the sum of k_ t and the rates of the other steps leading away from RA n. If there are multiple open states, but only one state that can gate, then burst analysis will still work. If there are three or more components in the closed time distribution, then it is very difficult to decide where to place the cutoff in closed times to classify gaps. If there is an open channel blocking process, then tg will no longer be a function solely of k_l and a. If experiments have established the validity of the model embodied in Eq. (10), then burst analysis can be carried out on the system under a wide variety of conditions to see how the rate constants a, t , and k_~ are affected by experimental manipulations. Correlation Analysis Correlations are evident in much single-channel data. Bursts of long openings may be separated by very brief closures, whereas brief openings are often separated by long closed times. To detect these correlations in event lifetimes, one employs the linear correlation coefficient or covariance function, which are, respectively (t~tb - t~tb)/aaab and t#b - t~tb. In general, it is easy to write a program that will accept as input an idealized data set and compute these functions for pairs of events related in any way. The estimation of the statistical significance of a linear correlation coefficient is a standard part of linear regression analysis. Alternatively, events can be classified as long and short with reference to a cutoff time, and a two-by-two contingency test can be used to evaluate the significance of the correlation. 24 The simple demonstration of statistically significant correlations allows one to make an important statement about the gating mechanism: there is more than one gating pathway. This can most easily be established by determining the correlation coefficient of adjacent open and closed times, or for successive openings separated by a single closed time. ~ For short closed times, there is a high probability that the two successive open times are of the same channel. Thus, focusing on the open intervals separated by relatively brief closures improves the detection of significant correlations. 23 j. Dudel, C. Franke, and H. Hatt, Biophys. J. 57, 533 0990). 24 M. B. Jackson, B. S. Wong, C. E. Morris, H. Lecar, and C. N. Christian, Biophys. J. 42, 109 0983).
744
D A T A S T O R A G E A N D ANALYSIS
A
[51]
1.o .8 .6 .4 .2 0
B
16 12 8
4
m 0 ~ range of closed 0.075- 0.125- 0 . 5 2 5 - 1.725time intervals 0.125 0.:525 1.725 21
@
//
• -
2141
>41
42
507
(msec) # poirs
128
89
69
150
FIO. 3. Single-channel currents through nicotinic receptor channels in cultured embryonic mouse muscle were activated by the agonist suberyldicholine and analyzed for correlations in successive open times (from Ref. 25). Pairs of openings were separated on the basis of the intervening closure lifetime (abscissa). The probability that a particular opening arose from the rapid or slow exponential component was calculated for a sum of two exponentials fitted to all open times. This probability, denoted as PF, is plotted in (A) versus the closed time between pairs of openings. Openings adjacent to long closed times have higher values of PF, meaning they are more likely to be brief. Openings adjacent to brief closed times have lower values of PF, meaning they are more likely to be long. This indicates an inverse correlation between open times and the adjacent closed time. (B) Correlation between successive open times was tested with a two-by-two contingency test. Z2 reflects the deviation from random mixing of long duration and short duration consecutive open times. A high Z2 means that pairs of consecutive openings are more often similar in open time; members of each exponential component arc segregated. A low)~ suggests that long and short open times are mixed randomly. Pairs of openings separated by dosed times between 0.125 and 1.725 msec are highly correlated in open time. This indicates that there arc at least two distinct gating pathways. Pairs of openings separated by long open times are uncorrelated, as expected for the gating of independent channels. Pairs of openings separated by very brief closures are uncorrelated, suggesting that most of these brief closures gate only to a long duration open state. This is consistent with the analysis shown in (A).
[5 l ]
STATIONARY SINGLE-CHANNEL ANALYSIS
745
A measure of the correlation for successive openings of acetylcholine receptor channels is displayed for different values of the intervening dosed time~ (Fig. 325). Significant correlations for intermediate closed times (Fig. 3B) show that there is more than one gating pathway. The same result can be inferred from Fig. 3A, which shows that long duration openings are more likely to be adjacent to short duration closures. Figure 3A thus shows that one of the gating pathways is between a short lifetime dosed state and a long lifetime open state. Similar results have been obtained for a Ca 2÷activated K + channel. ~ Figure 3B suggests that closed states with very short durations have access to only one gating pathway but that dosed states with intermediate durations can gate by two pathways to either long duration or short duration open states. With long duration dosed times there are no significant correlations between the durations of the two openings, as expected, since such openings are likely to be openings of different channels. When large numbers of events are available (typically more than 10,000), the condition on the lifetime of a neighboring event can be exploited much more effectively. The quantitative dependence of the lifetime distribution can be examined. This method has been termed adjacent state analysis, and it has been used to discriminate between very complex models in the gating of a C1- channel from rat skeletal muscle.27 As noted above, the covariance as a function of n is a sum of exponenrials, where n is the number of events separating two events for which the covariance is being computed. The number of exponential terms in the covariance function is one less than the minimum number of gating pathways. 3,4 Although specifying the number of dosed states, open states, and gating pathways is often deemed substantial progress toward understanding the gating mechanism, it generally does not identify a unique model. Correlation analysis can be pushed further to discriminate between models. The most general way to obtain additional information starts with a derivation of the two-dimensional probability density function for adjacent open and dosed times for each candidate model with the requisite number of open states, closed states, and gating pathways. 3 For some of the candidate models, there may be no difference in the form of the two-dimensional probability density function, and discrimination between these models cannot be achieved solely by the analysis of stationary single-channel data) ,4 The two-dimensional probability density functions can be curve fitted, using likelihood maximization, to the relevant list of pairs of dwell 25 M. B. Jackson, Adv. Neurol. 44, 171 (1986). 26 O. B. McManus, A. L. Blatz, and K. L. Magleby, Nature (London) 317, 625 (1985). 27 A. L. Blatz and K. L. Magleby, J. Physiol. (London) 410, 561 (1989).
746
DATA STORAGE AND ANALYSIS
[52]
times, prepared from the idealized data record. 3,2s The fits provided by the different models can then be compared and ranked according to a goodness-of-fit criterion. 28 These methods have been developed in detail and applied to the glutamate receptor channel of locust muscle. 29
Acknowledgments I thank Larry Truss¢ll for critical comments on the manuscript. 2s F. G. Ball and M. S. P. Sansom, Proc. R. Soc. London B 236, 385 (1989). 29 S. E. Bates, M. S. P. Sansom, F. G. Ball, R. L. Ramsey, and P. N. R. Usherwood, Biophys. J. 58, 219 (1990).
[52]
Analysis of Nonstationary
Single-Channel
Currents
By F. J. SIGWORTH and J. Z n o u Introduction A nonstationary process is a random process whose behavior changes with time. In this chapter we consider the analysis of channel activity as a particular kind of nonstationary process. A channel is prepared in a certain state or configuration of states, for example, by holding the membrane at a negative potential. Then a sudden perturbation is applied, for instance, a step depolarization, causing the channel to relax from this configuration, passing through various states on its way to a new extuilibdum configuration; the time course of this channel activity is what is analyzed. Some kinds of channels must be studied in pulsed experiments of this kind. Sodium channels and A-current potassium channels inactivate very strongly, so that one has little hope of observing much channel activity unless one gives pulsed depolarizations. Similarly, some neurotransmitter receptors desensitize rapidly, making pulsed application of agonist necessary. Analyzing the results of such pulsed experiments involves extra difficulty compared to continuous recordings: there are artifacts from the stimulus to deal with, and special consideration has to be given to dwell times immediately after the stimulus. However a great deal of extra information can be gained from such data. For example, if the channel starts in a particular closed state So, the time before the channel first opens (the first latency) can be analyzed to give a description of the kinetic steps between So and the open states of the channel. METHODS IN ENZYMOLOGY,VOL. 207
Cop~d~t © 1992by AcademicPress,Inc. Allr~htsof reproductionin any formr~c~'rved.
STATIONARYSINGLE-CHANNELANALYSIS
[51 ] S t a t i o n a r y S i n g l e - C h a n n e l
729
Analysis
By MEYER B. JACKSON Introduction The opening and closing of ion channels in artificial and cell membranes generates a signal of membrane current that varies in a stepwise manner with time. The stochastic manner in which this signal varies depends on the underlying channel-gating mechanisms. Thus, analysis of records of sin#e-channel current can provide much insight into the mechanisms by which channels open and close. Because channel gating is a stochastic process, the methods of analysis are statistical. This chapter describes a number of useful methods of kinetic analysis of stationary single-channel data. In a stationary system, the properties of the ensemble do not change with time; the mean probabilities of occupancy of specific configurations remain constant. Such a system must satisfy a condition of detailed balance, in which the frequency of transitions toward a particular configuration is equal to the frequency of transitions out. In nonstationary systems, detailed balance is not satisfied, and the mean probabilities of occupancy of a particular configuration change with time. The practical implication is that if the data are stationary, one simply can collect the data continuously for analysis. Nonstationary conditions are often achieved by perturbing the system, most commonly with changes in voltage. Methods for the analysis of nonstationary sin#e-channel data are presented elsewhere (see [52], this volume). Sin#e-channel kinetic analysis can be divided into three distinct steps. The first step involves the theoretical derivation of predictions of kinetic models. Examples of such theoretical predictions include the probability distribution function for open times and the joint probability density function for the lifetimes of two consecutive events. The second step of single-channel kinetic analysis is to extract something from real data that can be compared to theoretical predictions. These operations generally start with raw sin#e-channel data and generate open time and closed time histograms, or more complicated distillates of channel behavior. One thus obtains an experimental result in a form that can be directly compared with a theoretical prediction. This comparison, which most often involves some form of curve fit, is the third and final step of sin#e-channel kinetic analysis. The insights into channel gating processes are obtained through a comparison of distillates of sin#e-channel data with the theoretical predictions of models. METHODS IN ENZYMOIX~Y, VOL. 207
Colaytisht © 1992 by Academic Prem, Inc. All ri$1mof relm)duetion in any form reserved.
730
DATA STORAGE AND ANALYSIS
[51]
Single-channel analysis is very software intensive, and a large number of computer programs have been developed in many laboratories and by many companies. Many software packages are commercially available. The number of programs is so great, and the purposes so varied, that it is not possible to evaluate and compare them. The emphasis is rather on the specific tasks for which the programs are written. It is hoped that after reading this chapter, the reader will know what to expect from the software, how to use it intelligently, and what are the possibilities for developing programs for more specialized needs.
S u m m a r y of Kinetic T h e o r y The most widely used theories of single-channel kinetics treat channel gating as a finite Markov process. The complex trajectory of a channel through its protein conformational space is modeled as a series of jumps or transitions between states within which the conductance is relatively constant. Such models for channel gating then consist of a finite number of states that are either open, closed, or have an intermediate conductance. Each state can then interconvert to another with a finite rate. This is represented in a kinetic scheme by an arrow between the two states that can interconvert. When a rate is zero in both directions, then the two states can be considered unconnected. Within this framework, the relevant questions about a channel-gating mechanism are the following: (l) How many states are there with a given conductance? (2) How are they connected? (3) What are the rate constants for each step? This description defines the term "state," which is used extensively in the ensuing discussion. Relative to the time scales accessible in a singlechannel experiment (> 20/~sec), the transitions between states are instantaneous, whereas the states themselves are stable enough to observe. According to the Markov assumption, the probability of a particular transition depends only on which state the system is in. This model for channel gating corresponds to a specific physical picture which at the present time appears to be consistent with essentially all experiments on channels. One of the most useful results of the kinetic theory of channel gating provides a method of estimating the number of states with a particular conductance. Consider a membrane system with a single channel having n open states and m closed states. The open time and closed time distributions are represented by sums of n and m exponentials, respectively:
Po(t) = ~ x, eU~
(1 a)
[51]
STATIONARYSINGLE-CHANNELANALYSIS
7 31
m
Pc(t) = ~ Yj e vj
(lb)
where Po(t) and Pc(t) represent the probabilities of open times or closed times lasting a time t or longer.t.2 The uppercase symbols are used for the probability distribution, whereas the probability density function, obtained by differentiating the distribution, is denoted by lowercase letters. The parameters Pi, vj, x~, and yj depend on the rate constants and on the details of how the open and closed states are connected in the gating scheme. However, the form of Eqs. (1 a) and (1 b) does not depend on these details and is independent of the parameters and connectivity of the gating scheme. In effect, when one determines the number of exponentials that are needed to fit an observed distribution, one has "counted" the number of states with a particular conductance. More precisely, this method provides a lower bound to the number of such states, since one is always faced with the possibility that some of the terms of Eqs. (la) and (lb) are too small, too rapid, or too slow to detect. It is also possible that the system has two time constants that are too similar to be resolved. For even three closed or open states the values of the time constants and coefficients in Eqs. (la) and (lb) can be very complicated functions of the rate constants for transitions between states. Thus, for complex gating mechanisms it is often very difficult to go beyond state counting to a quantitative determination of specific rate constants. With complex gating behavior, it is often so arduous to determine the rate constants that, depending on the questions at hand, it may be prudent to focus on simpler quantities such as mean open time, mean closed time, or mean burst duration. Equations (la) and (lb) are examples of one-dimensional probability distributions; they are functions of a single experimental variable. One can also make use of higher dimensional probability distributions and density functions. For a closed time probability density function pc(to) and an open time probability density function po(to), where tc and to are independent (no correlations), the joint, two-dimensional probability density function for observing both to and tc is the product po(to)pc(tc). If the values of to and tc are correlated, then the joint probability density function for consecutive open and closed times is no longer simply the product of the two one-dimensional probability density functions. This is a qualitative difference, which can be used to estimate the number of distinct gating pathways and to make distinctions between these pathways. In general, correlation analI D. Colquhoun and A. G. Hawkes, Proc. R. Soc. London, Set. B 199, 231 (1981). 2 D. Colquhoun and A. G. Hawkes, in "Single-Channel Recording" (B. Sakmann and E. Neher, eds.), p. 135. Plenum, New York, 1983.
732
DATA STORAGE AND ANALYSIS
[51 ]
ysis is used more to answer questions about the basic topology or connectivity of a channel-gating mechanism, rather than to estimate the values of specific rate constants. One important example of correlation analysis is highly analogous with state counting described above, in that it provides an estimate of the minimum number of gating pathways of a channel. This arises from the theoretical result that the number of terms in the covariance function of two event lifetimes as a function of the number of intervening time intervals is equal to one less than the number of distinct pathways between the open and closed states. The number of pathways is defined as the minimum of the number of closed states that can open and open states that can close. 3,4 Thus, when one determines the number of exponentials that are needed to fit an observed covariance function, one has "counted" (or, more precisely, estimated a lower bound to) the number of gating pathways. Correlation and covariance functions can be derived with higher dimensions (i.e., with more time intervals). However, there is no additional information about the underlying channel-gating mechanism in these higher dimensional functions. 3 Thus, they have never been used in singlechannel analysis. As the number of states of a model increases, the number of kinetic parameters and the number of topological possibilities for connecting the states increase. Singie-channel analysis is limited in the number of parameters that can be determined, and in the power to discriminate models, a,4 Nevertheless, a great deal of information can be obtained purely from statistical analysis of single-channel records. To go beyond the limits circumscribed by single-channel kinetic theory, it is necessary to combine the stationary kinetic analysis with other approaches, such as nonstationary kinetic analysis and dose-response behavior. Idealization of Single-Channel Data The first step in single-channel kinetic analysis is the idealization of the data. Raw single-channel data represent the superposition of background noise and stepwise changes in channel current. When the data have been idealized, the background noise has been eliminated and only the stepwise changes caused by channel opening and closing remain. With an ideal channel record, no information is lost by replacing the current versus time signal with a two-column list consisting of a time and a current for each 3D. R. Fredkin, M. Montal, and J. A. Rice, in "Proceedingsof the BerkeleyConferencein Honor of Jerzy Neyman and Jack Kiefer" (L. M. Leeam and R. A. Olshen eds.), Vol. 1, p. 269. Wadsworth,Belmont, California, 1985. 4D. Colquhoun and A. G. Hawkes,Proc. R. Soc. London B 230, 15 (1987).
[51 ]
STATIONARY SINGLE-CHANNEL ANALYSIS
733
interval at a particular current level. Idealization is then the reduction of a single-channel data record to such a list. The ease with which single-channel data are idealized is critically dependent on the ratio of the single-channel current amplitude to the background noise. For noisy data, no ideal idealization method exists. Idealization can, and often is, done by hand, or with a computer in a semiautomatic, event-by-event mode. In this mode the computer will search the data for threshold crossings and then display the event, indicating guesses for the precise time of the transition and the levels. A person running the program is then given the opportunity to accept, reject, or modify the event and the initial guesses. Examination of channel current data event by event is essential with data from a new type of experiment where the basic properties are unknown. When one attempts an event-byevent analysis, even of fairly good data, one encounters ambiguous events. There is always a chance that a random background fluctuation will occur that is comparable in size to the single-channel current. For poor data where the channel current is only 3 or 4 times the peak-to-peak background noise, such fluctuations will occur frequently enough to make idealization very difficult. Ambiguities arise more often for brief intervals. 5 Analyzing records with thousands of events in an event-by-event mode is simply too time consuming to be useful if a systematic study is undertaken. A fully automated computer program for idealization is then required. Such programs are available from Axon Instruments (Foster City, CA) and Instrutech (Elmont, NY), as well as from other companies. Most programs allow the user to toggle between two modes, one of which is a semiautomated, event-by-event analysis by the operator. The other mode is fully automated. This toggling is an essential feature, because it is necessary to check the performance of these programs. The lists of numbers produced by a fully automated idealization program can be very deceptive. Glitches in data, such as those caused by power surges, seal breakdown, or bubbles in the perfusion system, will produce crossings of the threshold used to signify channel-gating transitions. Such events are easily recognized as spurious by the human eye, owing to their nonstepwise character. In an event-by-event analysis, an investigator can almost always avoid such mistakes, but to most computer programs, all threshold crossings look the same. Another source of false transitions has been described by Colquhoun. 6 If the records contain subconductance states (or for that matter infrequent appearances of another smaller channel) that are near the threshold for D. Colquhoun and F. J. Sigworth, in "Single-Channel Recording" (B. Sakmann and E. Neher, eds.), p. 191. Plenum, New York, 1983. 6 D. Colquhoun, Trends Pharmacol. Sci. 9, 157, (1988).
734
D A T A STORAGE A N D ANALYSIS
[51]
event detection, a single long dwell time at this level will produce a large number of very brief events in the output list. To avoid these and other errors, the raw data should be spot checked, and the output lists should be examined for anomalous behavior. The performance of idealization programs should be compared for different values of the event detection parameters. The amplitude and lifetime distributions should be checked visually with every data set for unusual peaks or components. An anomalous excess of points somewhere in the distribution can almost always be traced to a problem in the data, and editing these events will improve the quality of analysis. As Colquhoun succinctly put it, "single-channel analysis costs time, ''6 and the temptation to turn a software package loose, and publish the results without careful inspection and monitoring along the way, should be resisted. Because idealization is the stage of analysis where the noise level will have its primary impact, the bandwidth and sampling frequency should be seIected to optimize the quality of idealization. The Nyquist criterion of digitization at twice the corner frequency of the low-pass filter is inappropriate for single-channel analysis. 5 Most investigators have determined empirically that one should digitize at 8 to l0 times the corner frequency of the filter. At an early stage in a study, a data set should be analyzed with different bandwidths and digitization frequencies to see where problems develop. When a channel exhibits gating activity on a rapid time scale, the bandwidth must be widened to observe these processes, but if the noise increases enough to impede the idealization, then better data form the best solution. Curve-Fitting Lifetime Distributions Once data have been idealized it is generally straightforward to produce a list of all the open times or all the closed times. At this point one is ready to find a sum of exponentials that best describes the relative frequencies of different event lifetimes. One can either fit to a histogram or work with the list of lifetimes directly using likelihood maximization. The first goal of this procedure is to estimate the number of exponential terms needed to describe the open time and closed time data, and thus obtain a count of the minimum number of states. When the number of open or closed states is one or two, it may also be useful to determine the parameters to relate them to specific rate constants of the channel-gating process.
Fitting to Histograms From a list of state durations, one can construct a histogram, N(ti, ti+~), where ti and tj+~ are times that mark the boundaries of the bins. When the bins are narrow, then p(t)(tl+~ - ti) becomes a viable theoretical estimate of
[ 51 ]
STATIONARY SINGLE-CHANNEL ANALYSIS
735
N(ti, ti+l). The sum of squares difference between these two quantities can can then be used as an index of the quality of the fit. Fitting then proceeds by minimizing the sum of squares with respect to the parameters ofp(t). A function minimization program can then be used to carry out this optimization. In most programs, the sum of the coefficients in the probability density function Ix,. or yj of Eqs. (1 a) and (lb)] is constrained to be one, to reduce the number of free parameters. Minimizing ~a has also been suggested for curve fitting. 5 A likelihood maximization method can also be used for fitting to histograms (see below). A Z2 goodness-of-fit test can be used to evaluate the fit and the improvement in the fit resulting from additional exponential terms. The number of exponential components in a distribution is estimated by comparing the goodness-of-fit for different numbers of terms. With a small data set, one has fewer events per bin, and the statistical fluctuations in N become large. The histogram will then appear noisy, the curve fitting will not yield good parameter estimates, and low amplitude exponential terms will be difficult to detect. If an attempt is made to widen bins to remedy this problem, then p(t)(t,+l - ti) will be a poor approximation for the predicted value of N, and this expression should be replaced by the integral ofp(t) over each bin, or P(O - P(t~+1)- When bins are too wide, there will be too few data points to determine the fitting parameters accurately. Thus, fitting to histograms is problematic with small data sets. Likelihood Maximization For n event lifetimes, ti, and a probability density function p(t), the likelihood is
L = f i P(ti)
(2)
A curve fit proceeds by variation of the parameters ofp(t) to maximize L. 5,7 In practice, the method of likelihood maximization is carded out with the logarithm of the likelihood, because L itself exhibits such extreme values that its computation can be difficult. The likelihood ratio test s,9 can then be used to estimate the significance of the improvement of a fit by an additional exponential term. This method does not involve binning data, as does the sum of squares method with histograms. It is thus advantageous with small data sets. The time required to compute L increases with the size of the data set, and the scatter of N in histograms becomes smaller, so that with large data sets it is 7 R. Horn and K. Lange, Biophys. J. 43, 207 (1983). 8 C. R. Rao, "Linear Statistical Inference and Its Applications." Wiley, New York, 1973. 9 M. B. Jackson, Biophys. J. 49, 663 (1986).
736
DATA STORAGE AND ANALYSIS
[51]
better to use histograms. However, the method of likelihood maximization can still be used with histograms. The above likelihood function [Eq. (2)] then becomes ~° i
L = I-[ p(t,)',
(3)
where ni is the number of events in bin L This can be computed much more rapidly for a very large event lifetime list than it can using Eq. (2) and is thus a convenient function to use in the fitting of binned histogram data. With large data sets maximum likelihood and least-squares techniques become equivalent, and there are no compelling reasons for preferring one method over the other. Out-of-Range Events
There will always be some events which are too brief to observe. There will also be events with lifetimes right at the limit of detection, such that the frequency with which such events occur is not accurately determined. In addition, it is often important to exclude from analysis the longest times, which can have anomalous effects on curve fitting. When using the sum of squares method, the correction for out-of-range events is trivial. One simply does not include the N values for bins that are out of range. For likelihood maximization, this procedure is not correct; simply ignoring such events is equivalent to assigning zero values to the experimentally observed frequencies. To correct for out-of-range events, one must fit to a probability density function which includes the condition that the interval is in the specified range: Thus the likelihood for a single event of duration t becomes L = p(t)
p(t) dt = p(t)/[P(t,) - P(tl)]
(4)
where ts and tl are the short- and long-time cutoffs, respectively. When the adaptation of likelihood maximization to binned data [Eq. (3)] is corrected in this fashion, one obtains the following index for the quality of fit)°,l~: ~ P ( t , ) - P(t,+,)~ ln(L)---- ~ n/In i p - - ~ s ) - - - ~ l) j
(5)
Many investigators increase the bin width exponentially to compensate for fewer events at longer durations (tt + 1 = Oti, where t~ is approximately 1.25). This provides constant binwidth on a logarithmic time axis. For event lifetime lists of more than a few thousand, this implementation of l0 F. J. Sigworth and S. M. Sine, Biophys. J. 52, 1047 (1987). H O. B. McManus, A. L. Blatz, and K. L. Magleby, PfluegersArch. 410, 530 (1987).
[51 ]
STATIONARYSINGLE-CHANNELANALYSIS
737
likelihood maximization [Eq. (5)] is at present the method of choice in fitting to a sum of exponentials. Careful attention must be given to the precise value of the short-time cutoff in these expressions [Eqs. (4) and (5)]. When lifetimes are determined as multiples of a sampling frequency, then the cutoff must be n + 0,5 times the sampling frequency,9 where n is usually an integer between 3 and 6. This form of cutoff is necessary because the true lifetimes are effectively rounded up or down by carrying out the idealization procedure on digitized data. The true cutoffis then at the rounding point. When the idealization procedure employs an interpolation procedure for estimating the lifetimes more accurately, then there is no longer such a restriction on the cutoff time. However, it is then important to evaluate the interpolation procedure carefully to ascertain that it interpolates uniformly. Digitization and binning can introduce another error called the sampling promotion error when the fastest time constant approaches the bin width or digitization time interval.t t,~2 Additional Corrections for Brief Events The above discussion presented a correction for the primary error produced by out-of-range events in the curve fitting of probability density functions. However, intervals that are too brief to detect can lead to additional errors, and there has been much analysis of how missed intervals influence single-channel analysis. The most commonly encountered problem is that ifa brief open (closed) time is missed, then the two adjacent closed (open) times appear as one closed (open) event with a lifetime equal to the sum of the three event lifetimes. A simple and effective correction can be used for missed closures when there is a single exponential distribution of open times, t3,t4 The open times will remain exponentially distributed despite the missed brief closures, but with a lengthened apparent time constant. The first step in the correction is to estimate the number of missed brief closures from the zero-time intercept of the best-fitting closed time distribution. The corrected time constant, z ~ , is then related to the observed value, ~o~, by the expression z~or = ~o~ No~/(No~ + Nm~a)
(6)
where NoB is the number of observed open times and N=ma is the number of missed brief closures. Treatment of the analogous case for closed times interrupted by brief openings is essentially the same. ~2s. Sineand J. H. Steinbach,J. Physiol. (London) 373, 129 (1986). t3 E. Neher,J. Physiol. (London) 339, 663 (1983). t4 F. Sachs,J. Neil, and N. Barkakati,PfluegersArch. 395, 331 (1982).
738
DATA STORAGE AND ANALYSIS
[51]
For multiple exponential distributions the problem becomes more complicated. Equation (6) is still correct for the mean open time, but it is difficult to extend this correction to each time constant of the open time distribution. If one can demonstrate by correlation analysis (discussed below) the extent to which the brief closures are associated with the different open states, then Eq. (6) can be implemented by replacing No~ with the number of observed openings corresponding to a particular component, Zo~ with the time constant of that component, and N . ~ d with the number of missed closures estimated to be associated with that component. The general problem of missed events has been treated systematically by Roux and Sauve.~5 Blatz and Magleby ~6have also presented some useful results.
Plotting Exponential Distributions In a simple linear plot of a lifetime distribution, it is often difficult to recognize different exponential components (Fig. 1A). For this reason transforms have been developed to facilitate visual inspection. Semilogarithmic plotting provides some improvement by transforming an exponential to a line. However, when the time constants differ by more than a factor of 100 it is difficult to use a linear time axis. A particularly useful means of displaying such data is in a modified log plot (Fig. IC), in which bin size increases exponentially with time, and ni is not sealed to the size of the bin. I°,1~ Further transforming the y axis by taking the square root l° produces a variance-stabilized plot, that is, the deviation or scatter in the data becomes independent of the y value (Fig. 1D).
Fitting to Data from System with Many Channels The theory employed above is strictly valid only for a system with one channel? ,2 When a system has many channels this theory is no longer valid. When channel openings are well spaced from one another, and simultaneous openings are very infrequent, it has been assumed) 7 and can be shown, 18 that the one-channel kinetic theory provides a reasonable approximation to certain aspects of the gating process. The problems arising in the analysis of a many channel system can be illustrated with a few examples. When examining a closed time interval in a system known to have more than one channel, there is no way of knowing whether the channel that closed to begin the interval is the same 15 B. Roux and R. Sanve, Biophys. J. 48, 149 (1985). 16 A. L. Blatz and K. L. Magleby, Biophys, J. 49, 967 (1986). 17 D. Colquhoun and B. Sakmann, J. Physiol. (London) 369, 501 (1985). is M. B. Jackson, Biophys. J. 47, 129 (1985).
[51 ]
STATIONARY SINGLE-CHANNEL ANALYSIS ,
,
i
,
,
,m.
,
,
.
,
739
,
_1 30Q
lO 5
A
'
250
104
. . 200
~o 3
150
B
I0 2
~' ~oo io I 50 0
. I 0
'T. I,Y::7";S-(:-"-,:--'!r=--; .... 10
20
30
40
SO
60
70
80
90
i
100
-5 10
10- ¢
~0-3
Time (rnsec)
T~me
500
10-I
10
0
i
10
1
10
2
10
(sec)
300 25O 200
C
2SO ¢
10-2
D
m 150
200
100 ~ " 150 100 50 o
I
-5 10
-10 4
i
I
-3 1Q
-2 10
i
-1 10
Time (~ec)
i
10
i
0
10
I
1
10 2
I
10
-5 I0
-4 10
-3 10
-10 2
-101
10
0
10
I
I
10
2
I0
llme (se~)
FIG. 1. Four representations of a dwell time distribution with two exponential components (from Ref. 10). A total of 5120 random numbers were generated according to a distribution with time constants of 10 msec (70% of the events) and 100 msec (30%) and binned for display as histograms in the lower part of each figure. The theoretical probability density functions for each component (dashed curves) and their sum (continuous curve) are superimposed with each curve. In each part of the figure the upper panel shows the absolute values of the deviation of the height of each bin from the theoretical curve, with dashed curves showing the expectation value of the standard deviation for each bin. The upper panels were plotted with expansion factors of 2.1, 5.4, 3.1, and 4.9, respectively. (A) Linear histogram. Events were collected into bins of 1 msec width and plotted on a linear scale. The 100-msec component has a very small amplitude in this plot. (B) Log~log display with variable width (logarithmic) binding. The number of entries in each bin was divided by the bin width to obtain a probability density in events/second, which is plotted on the ordinate. (C) Direct display of logarithmic histogram. Events were collected into bins of width 6x = 0.2 and superimposed on the histogram as the sum of two functions, as in Eq. (11) of Ref. 10. (D) Square-root ordinate display of a logarithmic histogram, as in (C). Note that the scatter about the theoretical curve is constant throughout the display.
740
DATA STORAGE AND ANALYSIS
[51]
channel that opened to end the interval. This produces a systematic error in closed time distributions. When two channel openings overlap, there is no way of assigning the open time. Neglecting such events produces a systematic error in open time distributions. For a system with N identical and independent channels, with a singlechannel closed time probability density function of the form Xm vjyje-Vl [the derivative of Eq. (1 b)], the closed time probability density function is ~8 P~rc(t) = (:~ vyj e-V1)[:~(y./vj) e - ' d ~-1 + ( N - 1)(:~ Y1 e-"J)2[:~(Y'/Vi) e-,i],~-2
(7)
When the number of channels is known, Eq. (7) can be fitted to data by likelihood maximization, as has been demonstrated for a stretch-activated channel. 19With a single-channel open time probability density function of :~/zixi e--~'~[the derivative of Eq. (1 a)], the probability density function for isolated channel openings is is
t,'(t) -- r. x,l,,
[ (y/vj) (y/vj) j
(8)
where a prime is used to distinguish the probability for isolated openings from the true single-channel open-time probability density function. A fit with Eq. (8) requires input of parameters obtained from a fit to closed times, and knowledge of N. When N is difficult to determine, as is generally the case when N is large, then one can approximate N - 1 as N, which leads to p' ( t ) ~- ('2 x,l~, e-.~)[:~(z/coj) e--'o~]/~(z./o~j)
(9)
where the parameters o)j and zj are taken from a fit of :~ (ojzj e-",t to the closed times from the same data set) ° The results of such a fit are shown in Fig. 2. When the frequency of channel opening is low, there is little deviation of the best fitting one-channel open time distribution from the observed histogram of isolated open times. With data from a patch where the channel opening frequency is high, the observed histogram of isolated open times departs substantially from the best-fitting distribution, indicating that the distribution of isolated openings is a poor representation of the kinetics of closure of a single channel. 2° The similarity between the two functions (the solid curves of Fig. 2A,B) obtained by fitting Eq. (9) indi~9F. Guhary and F. Sachs, J. PhysioL (London) 363, 119 (1985). 20 M. B. Jackson, £ Physiol. (London) 397, 555 (1988).
[51 ]
lOOO
STATIONARYSINGLE-CHANNELANALYSIS A
741
B
N 100
I0 [
2
l 4
i 6 t
(msec)
I
/
i
8
I0
12
i
I
I
I
I
I
L
2
4
6
8
lO
12
d
t (msec)
FIG. 2. Comparison of cumulative distributions of lifetimes of isolated openings (doUed curve) with the best-fitting single-channel open time probability density function (smooth curve), obtained by curve fitting with Eq. (8) (from Ref. 20). In (A) the level of channel activity was low, and there were almost never simultaneous channel openings. In (B) the level of channel activity was high; most channel openings were interrupted by other channel openings, and isolated openings were relatively rare. The data were obtained from cultured mouse muscle using 2 # M carbachol. The difference in channel activity between the two patches reflects primarily differences in receptor density. Note that whereas the isolated open time distributions are very similar. The parameters for the best fits in (A) were #t = 6.0 m see-I, #2 = 0.24 msec-t, x~ = 0.26 for open times, with a frequency of channel opening of 0.029 msec- L In (B) the parameters were #~ = 5.5 msec-~, g2 = 0.021 msec-t, and xt = 0.30, with a frequency of channel opening of 0.21 msec- ~.
cates that with both sparse and crowded channel data, the same kinetics operates at the single-channel level. When fitting to these equations using likelihood maximization, it is important to incorporate the out-of-range event correction in an analogous manner to that described for Eq. (4). The integrations are lengthy but straightforward. Equation (9) can be used to analyze open times in very crowded data sets, provided that at least several hundred isolated openings are observed. When the frequency of channel opening is so high that the slowest time constant of the observed dosed time distribution is shorter than the slowest time constant of the single-channel open time distribution, then the parameters determined by fitting Eq. (9) will have large standard errors. 2° It will then be necessary to carry out such fits on many data sets to obtain accurate estimates of parameters.
742
DATA STORAGE AND ANALYSIS
[51]
Burst Analysis
Many channels exhibit "bursting behavior," in which openings occur in clusters. The openings within such a burst are separated by short closed time intervals. A great deal of jargon has arisen around this phenomenon. Terms such as "flickers, .... gaps," and "nachschlags" refer to various types of events in a "bursty" channel record. In the case of chemically activated channels, the bursts of openings are generally thought to be repeated openings of the same liganded receptor-channel complex. This mechanistic interpretation has provided a framework for a useful simplified approach to the kinetic analysis. Treating ligand binding and channel gating as separate steps of a kinetic scheme gives the following model: R + nA .-~- RA,,~'~R*A,~
(10)
where A represents the ligand, R the receptor, and R* indicates that the channel is open. A gap or flicker is a sojourn in state RA, between two openings. The mean burst duration, tb, mean gap duration, tg, and mean number of gaps per burst, n8, can be derived for this model and solved for the kinetic parameters to give2! a = (n 8 + 1)/(t b
-
-
nets)
k_~ = I/[2ts(n s + 1)1 fl = 2 k _ l n s
(1 1)
It then remains to estimate tb, ns, and ts from the data. The mean gap duration is best determined as the time constant of the fast component of the best-fitting closed time distribution. By inspecting a closed time distribution, one can define an operational cutoff time and assign all shorter closed intervals as gaps. One can then reconstruct bursts and determine burst duration and number of gaps. Although the assignment of all closures shorter than a cutoff as gaps may seem arbitrary, when there are two time constants of a closed time distribution differing by more than a factor of 10, then the number of incorrect assignments becomes insignificant. This type of burst analysis has been used in the study of both acetylcholine~7and ;,-aminobutyric acid (GABA)22 receptor channel data. Burst analysis is particularly useful with data showing frequently occurring subconductance states, or with which idealization is difficult owing to a low signal-to-noise ratio. By including subeonductance states as part of a burst, one can avoid the problem of how to idealize these events, which would require treating the subconductances as separate states. 2~ D. Colquhoun and A. G. Hawkes, Philos. Trans. R. Soc. London, Ser. B300, 1 (1982). 22 j. Bormann and D. E. Clapham, Proc. Natl. Acad. Sci. U.S.A. 82, 2168 (1985).
[51 ]
STATIONARY SINGLE-CHANNEL ANALYSIS
743
Burst analysis has the advantage of being simple and the disadvantage of being too simple. Kinetic models very different from Eq. (10) can generate bursting behavior. Equation (1 1) may not be valid for these models. If there are other pathways leading away from state RAn, such as a desensitization step, 23 then what is determined as k_ 1 will actually be the sum of k_ t and the rates of the other steps leading away from RA n. If there are multiple open states, but only one state that can gate, then burst analysis will still work. If there are three or more components in the closed time distribution, then it is very difficult to decide where to place the cutoff in closed times to classify gaps. If there is an open channel blocking process, then tg will no longer be a function solely of k_l and a. If experiments have established the validity of the model embodied in Eq. (10), then burst analysis can be carried out on the system under a wide variety of conditions to see how the rate constants a, t , and k_~ are affected by experimental manipulations. Correlation Analysis Correlations are evident in much single-channel data. Bursts of long openings may be separated by very brief closures, whereas brief openings are often separated by long closed times. To detect these correlations in event lifetimes, one employs the linear correlation coefficient or covariance function, which are, respectively (t~tb - t~tb)/aaab and t#b - t~tb. In general, it is easy to write a program that will accept as input an idealized data set and compute these functions for pairs of events related in any way. The estimation of the statistical significance of a linear correlation coefficient is a standard part of linear regression analysis. Alternatively, events can be classified as long and short with reference to a cutoff time, and a two-by-two contingency test can be used to evaluate the significance of the correlation. 24 The simple demonstration of statistically significant correlations allows one to make an important statement about the gating mechanism: there is more than one gating pathway. This can most easily be established by determining the correlation coefficient of adjacent open and closed times, or for successive openings separated by a single closed time. ~ For short closed times, there is a high probability that the two successive open times are of the same channel. Thus, focusing on the open intervals separated by relatively brief closures improves the detection of significant correlations. 23 j. Dudel, C. Franke, and H. Hatt, Biophys. J. 57, 533 0990). 24 M. B. Jackson, B. S. Wong, C. E. Morris, H. Lecar, and C. N. Christian, Biophys. J. 42, 109 0983).
744
D A T A S T O R A G E A N D ANALYSIS
A
[51]
1.o .8 .6 .4 .2 0
B
16 12 8
4
m 0 ~ range of closed 0.075- 0.125- 0 . 5 2 5 - 1.725time intervals 0.125 0.:525 1.725 21
@
//
• -
2141
>41
42
507
(msec) # poirs
128
89
69
150
FIO. 3. Single-channel currents through nicotinic receptor channels in cultured embryonic mouse muscle were activated by the agonist suberyldicholine and analyzed for correlations in successive open times (from Ref. 25). Pairs of openings were separated on the basis of the intervening closure lifetime (abscissa). The probability that a particular opening arose from the rapid or slow exponential component was calculated for a sum of two exponentials fitted to all open times. This probability, denoted as PF, is plotted in (A) versus the closed time between pairs of openings. Openings adjacent to long closed times have higher values of PF, meaning they are more likely to be brief. Openings adjacent to brief closed times have lower values of PF, meaning they are more likely to be long. This indicates an inverse correlation between open times and the adjacent closed time. (B) Correlation between successive open times was tested with a two-by-two contingency test. Z2 reflects the deviation from random mixing of long duration and short duration consecutive open times. A high Z2 means that pairs of consecutive openings are more often similar in open time; members of each exponential component arc segregated. A low)~ suggests that long and short open times are mixed randomly. Pairs of openings separated by dosed times between 0.125 and 1.725 msec are highly correlated in open time. This indicates that there arc at least two distinct gating pathways. Pairs of openings separated by long open times are uncorrelated, as expected for the gating of independent channels. Pairs of openings separated by very brief closures are uncorrelated, suggesting that most of these brief closures gate only to a long duration open state. This is consistent with the analysis shown in (A).
[5 l ]
STATIONARY SINGLE-CHANNEL ANALYSIS
745
A measure of the correlation for successive openings of acetylcholine receptor channels is displayed for different values of the intervening dosed time~ (Fig. 325). Significant correlations for intermediate closed times (Fig. 3B) show that there is more than one gating pathway. The same result can be inferred from Fig. 3A, which shows that long duration openings are more likely to be adjacent to short duration closures. Figure 3A thus shows that one of the gating pathways is between a short lifetime dosed state and a long lifetime open state. Similar results have been obtained for a Ca 2÷activated K + channel. ~ Figure 3B suggests that closed states with very short durations have access to only one gating pathway but that dosed states with intermediate durations can gate by two pathways to either long duration or short duration open states. With long duration dosed times there are no significant correlations between the durations of the two openings, as expected, since such openings are likely to be openings of different channels. When large numbers of events are available (typically more than 10,000), the condition on the lifetime of a neighboring event can be exploited much more effectively. The quantitative dependence of the lifetime distribution can be examined. This method has been termed adjacent state analysis, and it has been used to discriminate between very complex models in the gating of a C1- channel from rat skeletal muscle.27 As noted above, the covariance as a function of n is a sum of exponenrials, where n is the number of events separating two events for which the covariance is being computed. The number of exponential terms in the covariance function is one less than the minimum number of gating pathways. 3,4 Although specifying the number of dosed states, open states, and gating pathways is often deemed substantial progress toward understanding the gating mechanism, it generally does not identify a unique model. Correlation analysis can be pushed further to discriminate between models. The most general way to obtain additional information starts with a derivation of the two-dimensional probability density function for adjacent open and dosed times for each candidate model with the requisite number of open states, closed states, and gating pathways. 3 For some of the candidate models, there may be no difference in the form of the two-dimensional probability density function, and discrimination between these models cannot be achieved solely by the analysis of stationary single-channel data) ,4 The two-dimensional probability density functions can be curve fitted, using likelihood maximization, to the relevant list of pairs of dwell 25 M. B. Jackson, Adv. Neurol. 44, 171 (1986). 26 O. B. McManus, A. L. Blatz, and K. L. Magleby, Nature (London) 317, 625 (1985). 27 A. L. Blatz and K. L. Magleby, J. Physiol. (London) 410, 561 (1989).
746
DATA STORAGE AND ANALYSIS
[52]
times, prepared from the idealized data record. 3,2s The fits provided by the different models can then be compared and ranked according to a goodness-of-fit criterion. 28 These methods have been developed in detail and applied to the glutamate receptor channel of locust muscle. 29
Acknowledgments I thank Larry Truss¢ll for critical comments on the manuscript. 2s F. G. Ball and M. S. P. Sansom, Proc. R. Soc. London B 236, 385 (1989). 29 S. E. Bates, M. S. P. Sansom, F. G. Ball, R. L. Ramsey, and P. N. R. Usherwood, Biophys. J. 58, 219 (1990).
[52]
Analysis of Nonstationary
Single-Channel
Currents
By F. J. SIGWORTH and J. Z n o u Introduction A nonstationary process is a random process whose behavior changes with time. In this chapter we consider the analysis of channel activity as a particular kind of nonstationary process. A channel is prepared in a certain state or configuration of states, for example, by holding the membrane at a negative potential. Then a sudden perturbation is applied, for instance, a step depolarization, causing the channel to relax from this configuration, passing through various states on its way to a new extuilibdum configuration; the time course of this channel activity is what is analyzed. Some kinds of channels must be studied in pulsed experiments of this kind. Sodium channels and A-current potassium channels inactivate very strongly, so that one has little hope of observing much channel activity unless one gives pulsed depolarizations. Similarly, some neurotransmitter receptors desensitize rapidly, making pulsed application of agonist necessary. Analyzing the results of such pulsed experiments involves extra difficulty compared to continuous recordings: there are artifacts from the stimulus to deal with, and special consideration has to be given to dwell times immediately after the stimulus. However a great deal of extra information can be gained from such data. For example, if the channel starts in a particular closed state So, the time before the channel first opens (the first latency) can be analyzed to give a description of the kinetic steps between So and the open states of the channel. METHODS IN ENZYMOLOGY,VOL. 207
Cop~d~t © 1992by AcademicPress,Inc. Allr~htsof reproductionin any formr~c~'rved.
