A new approach for processing of piecewise-constant signals is proposed. It is based on modeling the observed data as a sum of a random signal and noise. The random signal has a Gibbs distribution, and the noise is Gaussian. A MAP criterion in derived for joint estimation of the number of signal levels and reconstruction of signal. The criterion comprises of three terms, one corresponding to the likelihood of the data and two to penalties. One penalty term penalizes for unnecessary transitions, and the other, for unnecessary levels. The method has been tested on synthesized data and applied to single ion channel recording

Djuric, P.M.

Fwu, Jong-Kae

Jovanovic, S.

Lynn, K.

UNT Digital Library

ol\l F- 77605--/0 8- - 1 BNL-62775 4 Y I ON THE PROCESSING OF PIECEWISE-CONSTANT SIGNALS BY HIERARCHICAL MODELS WITH APPLICATION TO SINGLE ION CHANNEL CURRENTS Petar M. Djuri6,* Jong-Kae Fm,* Slobodan Jouanovit,** and Kelvin Lynn'B-k * Department of Electrical Engineering State University of New York Stony Brook, NY 11794-2350 ** Physics Departmen Brookhaven National Lab0 Upton, NY 11973 djuric@sbee.sunysb.edu jovanovic@bnlarm.bnl.gov fwu@sbee.sunysb.edu kgl@aps.org ABSTRACT easily interpreted. The method is implemented itera- tively in stages by using a treestructure (TS) initial- ization scheme [5] and the Iterated Conditional Modes (ICM) method [Z]. In the first stage, it is assumed that in the third three, and so on. After the completion of each stage, the best sequence is selected for that stage stages. If based on the derived criterion the sequence from the latest stage is better, it is kept, otherwise it is removed. The procedure ends with the processing of of levels. We are interested to apply our method to patch clamping recordings. One of the main objectives of patch clamping is the study of ion permeation mecha- nisms in biological membranes. Patch clamping allows for the isolation of small patches of membranes and in- volves measurement of ion channel currents. The ion channels are large proteins embedded in the membranes of all living cells. These macromolecules form pores across the cell membrane, and in certain conforma- tions, they allow the flow of ions into or out of the cell, thereby controlling the electrical function of the cell. The measured currents are noisy piecewise constant sig- nals which reflect the gating kinetics of the individual channels. Previous attempts to process such signals include the half amplitude analysis [8], mean-variance histograms [8], the stategram [7], various forms of the Hinkley detector [3], and Bayesian methods based on Markov chain prior distributions [4]. For an overview of statistical analysis methods in the study of ion chan- nels, see [l]. In this paper, we provide the derivation of our prc+ cedure, interpret the results, and show the necessary steps for its implementation. We also test the method on synthesized signals and apply it to single ion channel currents. A new approach for processing of piecewise-constant signals is propoged. It is based on modeling the ob- served data as a sum of a random signal and noise. The Gaussian. A MAP criterion is derived for joint estima- tion of the number of signal levels and reconstruction one corresponding to the likelihood of the data and two to penalties. One penalty term penalizes for unneces- sary transitions, and the other, for unnecessary levels. applied to single ion channel recordings. 1. INTRODUCTION Piecewise constant signals arise frequently in many ar- eas of science and engineering. They are characterized by several constant levels and are commonly corrupted by unknown noise. In many applications, the number of levels and their associated values are not known. The signals themselves also change levels randomly. The main processing task is that of detecting the number of levels and reconstructing the noiseless signal. There are a variety of methods for addressing this problem. In this paper, we propose a procedure based on hierarchical models. The observed data are repre- sented as a sum of two random processes, one corre- sponding to the signal and the other to the noise. The unobservable (nod-) signal is modeled by a Gibbs distribution and the noise by a Gaussian distribution. We apply a maximum a posteriori (MAP) criterion to obtain the optimal number of levels and the estimate of the associated signal. The resulting criterion is a penalized likelihood function with terms that can be This work WM supported by the National Science Foundation random signa1 h a  a Gibbs distribution, and the noise is si@a1 Only one level, in the second two levels, of the signal. The criterion comprises of three terms, and to the best sequence from the Previous The method has been tested on synthesized data and the hypothesis associated with the maximum number i under Award No. MIP-9506743 and, in ar t  b U.S. DOE under contract no. DE-ACd2-Y6Cii00016. DISTRIBUTION OF THIS DQCU 2. PROBLEM STATEMENT AND TERMINOLOGY Let x k  be a signal which represents a random sequence of N samples, with k being the number of possible signal values (leveb), and let Xk be a realization of Xk. The set of various levels is denoted by s k  = {SI, $2,. - Q } ,  and the collection of all possible signal sequencies x k  by the set x k ,  where Xk = {Xk = [Zk[l] zk[2] - * *  zk[N]IT : Zk[d E &, 1 5 i 5 N}. The Signal remains constant with a value equal to 8j E &, i = 1,2,. .C, for some time, which is referred to as the dwell time, and then moves to a different level. It then remains there for another dwell time, then changes lev- els again, etc. The signal Xk is unobservable. Instead, we observe the noise corrupted version y ,  which is ob- tained as a sum of x k  and a white noise vector w, Le., y = X k + W .  (1) The number of signal levels k, the signal levels & , and the noise parameters are unknown. Also unknown are the dwell times of the signal. Based on the model (1) and the assumptions 0 the signal x k  has a finite number of well defined levels k, where k 5 K, for a known K, 0 the signal is described by a Gibbs distribution, 0 the noise samples are independent and zero mean Gaussian with variances that depend on the sig- nal levels, and e the introduced random processes are time revers- ible, the objective is to determine from y the number of levels h, k E (1,2, - - - , K}, and estimate the signal xk, which includes the signal levels Si, i = 1,2, * e ,  C, and the associated dwell times. One of the key assumptions is the Gibbs distribu- tion of the random signal Xk. If xk is one outcome from the set .&, its probability is given by [2], [6], where Z t  is a normalizing constant, and V(xk)  is an energy function. With the assumption (2), we adopt a neighborhood system N, which is a collection of subsets of {1,2,.-.,N}, that is N = {Ni : i E (1,2,...,N}}, where Ni denotes the neighbors of i and satisfies the conditions (a) i Nj and (b) i E Nj iff j E Ni. Thus, Xk is a one-dimensional Markov random field with re- spect t o N .  Finally, we introduce the sequence of k different la- bels l k ,  where It[d € { 1,2, . - - , k}. This sequence is associated with the signal sequence xk via I (3) where a(-) is a function that uniquely maps the label sequence It to the signal Xk according to Z k [ q  = s([~J, i ~ { 1 , 2 , . . . , N }  . Forexample,ifIk=[I 1 2 1 3  '..IT, x k  = (81 81 sa SI 83 . -  -IT. In quantifying the probabil- ity p(xk), we choose energy functions that allow us to write d x k )  = dlk).  (4) 3. ESTIMATION CRITERION We want to apply the maximum a posteriori probability (MAP) criterion for estimating k. The marginalized MAP estimate is given by ( 5 )  where p(k1y) is the posterior probability mass function of k given the observed data y- This criterion will be very difficult to implement because its solution is based on where f(yIlt,&) is the density of the data given the label sequence l k  and the signal and noise parameters 6k with prior density f(&), 8 k  is the parameter space of 81, and p ( k )  is the prior probability of k levels. Note that the summation in (6) has an extremely large num- ber of terms. Another MAP estimator jointly provides the MAP estimates of k and Xk, and it is expressed by This criterion can readily be evaluated by using one of several existing iterative techniques. Finally, we introduce a third criterion, which we find more appropriate than (7) for reasons to be ex- plained below. The form of the criterion is or .- .-  Clearly, the criterion (9) selects the joint MAP esti- mates of E and the label sequence l k .  Its main dif- ference from (7) is that it integrates out the unknown signal and noise parameters. Of course, once k and 1 k  are estimated, it is trivial to determine the signal levels si and obtain the estimate of Xk. The evaluation of (9) is easier than that of (6) but is still computationally intractable. One difficulty is the evaluation of the normdizing constants zk needed due to p(1k). To alleviate this problem, we substitute p(lt) by the pseudolikelihood (21 where 4 denotes the parameters of the Gibbs distri- bution, and &[q is the set of neighboring samples of Another problem in evaluating (9) is the computa- tion of the integral JB, f (Yllk,  &)f(8k)d&.  However, by Taylor expending f(yllk, 6,) around the maximum likelihood estimates 8 k ,  and using asymptotic assump tions, we can approximate it by l k  [i]. where it is the estimated label sequence, and n j  are the total number of samples whose label is It[-] = j. With these approximations and the assumptions that the additive noise is zero mean Gaussian with level dependent variance and p(E) = 1/K, the MAP crite- rion results in (12) where 6; is the estimated variance of the samples la- beled as a(.] = j. The interpretation of the three terms in (12) is straightforward. The first term is the likeli- hood which decreases with increasing k. The second and third terms are penalties for introducing signal transitions from one level to another and additional levels, respectively, and they grow with the number of transitions and k. Note that penalization for the pa- rameters of the Gibbs distribution is not necessary be- cause we assume the same parameters are present for all the models used in (12). 4. IMPLEMENTATION OF THE PROPOSED CRITERION The criterion (12) can be implemented by applying the iterated conditional modes (ICM) algorithm [l]. For good performance, the ICM requires relatively good initial conditions. To obtain them Fe propose to use the tree-structure (TS) method introfluced in [5], which is a completely data driven scheme. The overall (Ts ICM) procedure consists of the following steps. First, assume that k = 1, estimate the only level 61, and evaluate the criterion function. Set h = 2, and as initial conditions use iy’ = i 1  - e and i?’ = 61 + e, where e is some small number, and 61 is the result for k = 1. Apply the ICM method, estimate 61 and 62, and evaluate the criterion. Set k = 3 and use as initial conditions iy’ = 81 - e ,  if) = 61 + E, and ic) = 82, with 61 and 62 being the results for k = 2, and continue with the ICM and the evaluation of the criterion. As another possible set of initial conditions use 6(,0) = i l l  6‘2”’ = 02-6, and 6:’ = &+e, and again apply the ICM and evaluate the criterion. As a final sequence of three levels, 13, choose the one that has the smaller criterion value. Next set k = 4, and continue along the same lines until the testing of k = K leyels is completed. The solution is the label sequence l k  that yields the smallest criterion value. Once i k  is determined, the signal x i  is easily estimated. 5. SIMULATION RESULTS Our method was tested on synthesized data and ap- plied to real patch clamp recordings. The synthesized data records had 2000 samples, three different levels, and 6 level changes. The signal-to-noise ratio (SNR) was defined by minIA/ai, where A is the difference between the levels of two adjacent segments and u the standard deviation of the noise, which was the same for every segment. The SNR was varied between 1 and 5 in steps of 1. For each SNR, there were 100 trials. The probability p(/k[~1&[8~, 4) in (10) was defined by where 4 = p, and um[d denotes the number of neigh- bors of i having the label m. The number of neighbors in the experiment was equal to four, that is, two neigh- bors on each side. The results are shown in Table 1. For SNR’s of 2 and higher, the method had very good performance. Figure 1 at the top shows a typical real- ization with SNR=2 and at  the bottom, the histogram of the detected level changes in 100 trials. The peaks of the histogram are at the correct locations of the signal transitions. Figure 2 at the top displays a real patch clamp recording, which is quite complex and has several con- ductance levels. We applied our procedure to these data to determine the number of levels and estimate kk. The maximum number of hypothesized levels was 15. The procedure found 10 levels and estimated the noise-frke signal shown at the bottom of Figure 2. Table 1: The entries represent the number of times k levels were estimated in 100 trials. l e o r  . 1 6. REFERENCES [l] F. G. Ball and J. A. Rice, “Stochastic models for ion channels: Introduction and bibliography,” Mathe- matical Bioscience, vol. 112, pp. 189-206, 1992. [2] J. Besag, “On the statistical analysis of dirty pic- tures,” Journal of tbe Royal Statistical Society B, vol. 48, no. 3, pp. 259-302,1986. [3] S. Draber and R. Schultze, UDetection ofjumps in single-channel data containing subconductance lev- els,” Biophysical Journal, vol. 67, pp. 1404-1413, 1994. [4] D. R. F’redkin and J. A Rice, “Bayesian restoration of single channel patch-clamp recordings,” Biomet- r i a ,  vol. 48, pp. 427-448, 1992. [5] J.-K. Fwu and P. M. DjuriC, “Unsupervised vec- tor image segmentation by a tree structure-ICM al- gorithm,” submitted to the IEEE “kansactions on Medical Imaging. [6] S. Geman and D. Geman, “Stochastic relaxation, Gibbs distributions, and the Bayesian reconstruc- tion of images,” IEEE “kansactions on Pattern Analysis and Machine Intelligence, vol. PAMI-6, no. [7] S. Jovanovic, “A system for on-line analysis of sin- gle ion channel currents,” Master Thesis, State Uni- versity of New York, 1994. [8] J. B. Patlak, “Measuring kineties of complex single ion channel data using mean-variance histograms,” Biophysical Journal, vol. 65, pp. 29-42, 1993. 6, pp. 721- 741, 1984. Figure 1: Top: A realization with 3 levels and 6 transi- tions for SNR=Z. Bottom: Histogram of the estimated transitions in 100 trials. 1 800 1 1400 1200 1000 600 e00 400 200 0 1800 1400 1200 1000 800 e00 400 200 0 0 0 1 0 0  ZOO 300 400 SO0 600 700 a00 BOO IO00 Figure 2: Top: Red data with unknown number of levels. Bottom: Estimated signal with 10 levels. i DISCLAIMER This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsi- bility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Refer- ence herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recom- mendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof. 

On the processing pf piecewise-constant signals by hierarchical models with application to single ion channel currents

https://digital.library.unt.edu/ark:/67531/metadc665224/m2/1/high_res_d/224109.pdf

On the processing pf piecewise-constant signals by hierarchical models with application to single ion channel currents

Abstract

Similar works

Full text

Available Versions

UNT Digital Library