Fully-Automatic Multiresolution Idealization for Filtered Ion Channel Recordings: Flickering Event Detection

TL;DR

JULES is a new model-free segmentation method that detects flickering ion channel events at small scales, providing precise dwell times and amplitudes, and can enhance hidden Markov analysis.

Contribution

The paper introduces JULES, a novel multiresolution and deconvolution-based segmentation method for ion channel recordings, capable of detecting small-scale flickering events.

Findings

JULES accurately detects flickering events below filter frequency.

It provides precise dwell times and amplitude levels.

It enhances hidden Markov analysis for ion channel data.

Abstract

We propose a new model-free segmentation method, JULES, which combines recent statistical multiresolution techniques with local deconvolution for idealization of ion channel recordings. The multiresolution criterion takes into account scales down to the sampling rate enabling the detection of flickering events, i.e., events on small temporal scales, even below the filter frequency. For such small scales the deconvolution step allows for a precise determination of dwell times and, in particular, of amplitude levels, a task which is not possible with common thresholding methods. This is confirmed theoretically and in a comprehensive simulation study. In addition, JULES can be applied as a preprocessing method for a refined hidden Markov analysis. Our new methodolodgy allows us to show that gramicidin A flickering events have the same amplitude as the slow gating events. JULES is available as an R function jules in the package clampSeg.

Tables

Table 1. Table 1 . Performance of JULES JULES \operatorname{JULES} , TRANSIT TRANSIT \operatorname{TRANSIT} and JSMURF JSMURF \operatorname{JSMURF} in idealizing a signal with an isolated peak having levels l 0 = l 2 = 40 subscript 𝑙 0 subscript 𝑙 2 40 l_{0}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=l_{2}}=40 , l 1 = 20 subscript 𝑙 1 20 l_{1}=20 , change-points τ 1 = 0.2 subscript 𝜏 1 0.2 \tau_{1}=0.2 and τ 2 = τ 1 + ℓ / f 𝒔 subscript 𝜏 2 subscript 𝜏 1 ℓ subscript 𝑓 𝒔 \tau_{2}=\tau_{1}+\ell/f_{\boldsymbol{s}} , ℓ = 2 , 3 , 5 ℓ 2 3 5 \ell=2,3,5 . The standard deviation of the error used in each simulation is σ 0 = 1.4 subscript 𝜎 0 1.4 \sigma_{0}=1.4 . Results are based on 10 000 10000 10\,000 pseudo samples.

Method	Length ($\ell$)	Correctly identified ($\%$)	Detected ($\%$)	False positive (Mean)
$\mathrm{JULES}$	2	65.17	65.17	0.0290
$\mathrm{TRANSIT}$	2	0.01	96.02	19.9692
$\mathrm{JSMURF}$	2	0.00	0.00	0.2471
$\mathrm{JULES}$	3	99.82	99.82	0.0004
$\mathrm{TRANSIT}$	3	0.02	74.65	19.8013
$\mathrm{JSMURF}$	3	0.00	0.00	0.1366
$\mathrm{JULES}$	5	100.00	100.00	0.0000
$\mathrm{TRANSIT}$	5	0.00	98.20	19.8484
$\mathrm{JSMURF}$	5	0.00	0.00	0.0000
$\mathrm{MDL}$	5	0.02	100.00	96.8718

Table 2. Table 2 . Performance of JULES JULES \operatorname{JULES} and TRANSIT TRANSIT \operatorname{TRANSIT} in idealizing a signal with an isolated peak having levels l 0 = l 2 = 40 subscript 𝑙 0 subscript 𝑙 2 40 l_{0}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=l_{2}}=40 , l 1 = 20 subscript 𝑙 1 20 l_{1}=20 , change-points τ 1 = 0.2 subscript 𝜏 1 0.2 \tau_{1}=0.2 and τ 2 = τ 1 + ℓ / f 𝒔 subscript 𝜏 2 subscript 𝜏 1 ℓ subscript 𝑓 𝒔 \tau_{2}=\tau_{1}+\ell/f_{\boldsymbol{s}} , ℓ = 2 , 3 , 5 ℓ 2 3 5 \ell=2,3,5 . The standard deviation of the error used in each simulation is σ 0 = 1.4 subscript 𝜎 0 1.4 \sigma_{0}=1.4 . Results are based on 10 000 10000 10\,000 pseudo samples. Results are given as multiples of the sampling rate f 𝒔 = 10 4 subscript 𝑓 𝒔 superscript 10 4 f_{\boldsymbol{s}}=10^{4} .

Method	Length ($\ell$)	$f_{𝒔}^2\mathrm{𝖬𝖲𝖤}\left({\hat{\tau }}_1\right)$	$f_{𝒔}\mathrm{𝖡𝖨𝖠𝖲}\left({\hat{\tau }}_1\right)$	$f_{𝒔}\mathrm{𝖲𝖣}\left({\hat{\tau }}_1\right)$	$f_{𝒔}^2\mathrm{𝖬𝖲𝖤}\left({\hat{\tau }}_2\right)$	$f_{𝒔}\mathrm{𝖡𝖨𝖠𝖲}\left({\hat{\tau }}_2\right)$	$f_{𝒔}\mathrm{𝖲𝖣}\left({\hat{\tau }}_2\right)$
$\mathrm{JULES}$	2	0.4022	-0.1047	0.6255	0.2677	0.0587	0.5141
$\mathrm{TRANSIT}$	2	9.3427	2.6218	1.5713	24.1103	4.7920	1.0709
$\mathrm{JULES}$	3	0.1170	0.0044	0.3420	0.1087	-0.0012	0.3297
$\mathrm{TRANSIT}$	3	13.0761	2.9386	2.1074	24.2032	4.7115	1.4162
$\mathrm{JULES}$	5	0.0670	-0.0026	0.2588	0.0669	0.0025	0.2587
$\mathrm{TRANSIT}$	5	15.4710	3.5622	1.6679	18.6559	4.0973	1.3670

Table 3. Table 3 . Performance of JULES JULES \operatorname{JULES} and TRANSIT TRANSIT \operatorname{TRANSIT} in idealizing a signal with an isolated peak having levels l 0 = l 2 = 40 subscript 𝑙 0 subscript 𝑙 2 40 l_{0}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=l_{2}}=40 , l 1 = 20 subscript 𝑙 1 20 l_{1}=20 , change-points τ 1 = 0.2 subscript 𝜏 1 0.2 \tau_{1}=0.2 and τ 2 = τ 1 + ℓ / f 𝒔 subscript 𝜏 2 subscript 𝜏 1 ℓ subscript 𝑓 𝒔 \tau_{2}=\tau_{1}+\ell/f_{\boldsymbol{s}} , ℓ = 2 , 3 , 5 ℓ 2 3 5 \ell=2,3,5 . The standard deviation of the error used in each simulation is σ 0 = 1.4 subscript 𝜎 0 1.4 \sigma_{0}=1.4 . Results are based on 10 000 10000 10\,000 pseudo samples. Note that JULES JULES \operatorname{JULES} overestimates the amplitude in very few cases (less than 0.2 % percent 0.2 0.2\% ) heavily while in all other simulations the estimation is very well. To illustrate this issue we also reported trimmed versions of all error quantities, where we only considered values between 0 0 and 40 40 40 .

Method	Length ($\ell$)	$\mathrm{𝖬𝖲𝖤}\left({\hat{l}}_1\right)$	$\mathrm{𝖡𝖨𝖠𝖲}\left({\hat{l}}_1\right)$	$\mathrm{𝖲𝖣}\left({\hat{l}}_1\right)$	${\mathrm{𝖬𝖲𝖤}}_{\mathrm{Trim}}\left({\hat{l}}_1\right)$	${\mathrm{𝖡𝖨𝖠𝖲}}_{\mathrm{Trim}}\left({\hat{l}}_1\right)$	${\mathrm{𝖲𝖣}}_{\mathrm{Trim}}\left({\hat{l}}_1\right)$
$\mathrm{JULES}$	2	222978.6771	-83.6452	464.7745	32.5343	-0.6010	5.6726
$\mathrm{TRANSIT}$	2	112.7163	10.1733	3.0367	99.3660	9.7529	2.0609
$\mathrm{JULES}$	3	552.8490	-1.0170	23.4919	12.3488	-0.7219	3.4393
$\mathrm{TRANSIT}$	3	69.6075	7.0327	4.4891	47.2984	6.2367	2.8988
$\mathrm{JULES}$	5	2.7763	-0.1081	1.6628	2.7763	-0.1081	1.6628
$\mathrm{TRANSIT}$	5	26.8554	2.0811	4.7462	11.5351	1.4210	3.0849

Table 4. Table 4 . Robustness of JULES JULES \operatorname{JULES} against additional noise components in idealizing a signal with an isolated peak having levels l 0 = l 2 = 40 subscript 𝑙 0 subscript 𝑙 2 40 l_{0}=l_{2}=40 , l 1 = 20 subscript 𝑙 1 20 l_{1}=20 , change-points τ 1 = 0.2 subscript 𝜏 1 0.2 \tau_{1}=0.2 and τ 2 = τ 1 + ℓ / f 𝒔 subscript 𝜏 2 subscript 𝜏 1 ℓ subscript 𝑓 𝒔 \tau_{2}=\tau_{1}+\ell/f_{\boldsymbol{s}} , ℓ = 2 , 3 , 5 ℓ 2 3 5 \ell=2,3,5 . The standard deviation of the error used in each simulation is σ 0 = 1.4 subscript 𝜎 0 1.4 \sigma_{0}=1.4 . Results are based on 10 000 10000 10\,000 pseudo samples.

Method	Length ($\ell$)	Correctly identified ($\%$)	Detected ($\%$)	False positive (Mean)
White noise	2	64.42	64.42	0.0333
$f^2$ noise	2	65.56	65.56	0.0219
$1/f$ noise	2	78.43	99.93	0.2306
Heterogeneous noise	2	59.05	59.05	0.0214
White noise	3	99.82	99.82	0.0005
$f^2$ noise	3	99.88	99.88	0.0000
$1/f$ noise	3	78.02	100.00	0.2379
Heterogeneous noise	3	95.23	95.23	0.0023
White noise	5	100.00	100.00	0.0000
$f^2$ noise	5	100.00	100.00	0.0000
$1/f$ noise	5	73.60	100.00	0.2911
Heterogeneous noise	5	99.95	99.95	0.0000

Equations

$$f(t)=\sum _{j=0}^K{l}_j{1\text{l}}_{[{\tau }_j,{\tau }_{j+1})}(t),$$

$$Y_{i}=\big{(}F\ast(f+\sigma_{0}\eta)\big{)}\left(i/f_{\boldsymbol{s}}\right)=(F\ast f)(i/f_{\boldsymbol{s}})+\epsilon_{i},\quad i=1,\ldots,n,$$

$$Y_i=(F_m\ast f)(i/f_{\mathbf{s}})+{ϵ}_i,i=1,\dots ,n,$$

$$\mathsf{Cov}\big{[}\epsilon_{i},\epsilon_{i+j}\big{]}=\left\{\begin{array}[]{cl}\sigma_{0}^{2}(F\ast F)(j/f_{\boldsymbol{s}}),&\text{for }|j|=0,\ldots,m,\\ 0,&\text{for }|j|>m.\end{array}\right.$$

$$Y_i=g(i/f_{\mathbf{s}})+{ϵ}_i,i=1,\dots ,n.$$

$${\hat{\sigma }}_0=\frac{\mathrm{IQR}\left(Y_{m+1}-Y_1,\dots ,Y_n-Y_{n-m}\right)}{2{\Phi }^{-1}(0.75)\sqrt{2(F\ast F)(0)}},$$

$$M(\tilde{g}):=\max_{\begin{subarray}{c}1\leq i\leq j\leq n\\ \tilde{g}([i/f_{\boldsymbol{s}},j/f_{\boldsymbol{s}}])\text{ const.}\end{subarray}}\,\left\{|S_{i,j}|+p_{i,j}\right\}\\ :=\max_{\begin{subarray}{c}1\leq i\leq j\leq n\\ \tilde{g}([i/f_{\boldsymbol{s}},j/f_{\boldsymbol{s}}])\text{ const.}\end{subarray}}\,\left\{\left|\frac{1}{\sigma_{i,j}}\,\sum_{k=i}^{j}\,\big{(}Y_{k}-\tilde{g}(k/f_{\boldsymbol{s}})\big{)}\right|+p_{i,j}\right\}\leq q,$$

$$\hat{g}(t):=\sum_{j=0}^{\hat{K}}\,\hat{l}_{j}\,{1\textnormal{l}}_{[\hat{\tau}_{j},\hat{\tau}_{j+1})}(t):=\operatorname*{argmin}_{\tilde{g}\in G_{\hat{K}},M(\tilde{g})\leq q}\,\sum_{i=1}^{n}\big{(}Y_{k}-\tilde{g}(k/f_{\boldsymbol{s}})\big{)}^{2}.$$

$$(Y_{i,j}-{\hat{\mu }}_{i,j}{)}^T{\Sigma }_{i,j}^{-1}(Y_{i,j}-{\hat{\mu }}_{i,j}),$$

$$\hat{\mu}_{i,j}=\big{(}(K\ast\hat{f})((f_{\boldsymbol{s}}\hat{\tau}_{i}-m+1)/n),\ldots,(K\ast\hat{f})((f_{\boldsymbol{s}}\hat{\tau}_{j}+m-1)/n)\big{)}^{T}$$

Full text

1. Introduction

Ion channels are fundamental to many vital human functions like such as breathing and heartbeating as they allow for a regulative transport of ions through the otherwise ion impermeable lipid membrane along electrochemical gradients. The channels can switch between different conformations (gating) resulting in different conductance levels (states) controlled by external stimuli such as voltage, ligand binding or mechanical stress Chung et al. (2007). Understanding of ion channel gating has a long history and plays an important role for many purposes, e.g., for developing drugs Kass (2005); Overington et al. (2006). A major tool to analyze the gating behavior is the patch clamp technique which allows to record the conductance trace (i.e., the recorded current trace divided by the applied voltage) of a single ion channel Neher and Sakmann (1976); Sakmann and Neher (1995).

To stay in the transmission range of the amplifier the recordings are convolved by a lowpass filter (typically a 4, 6 or 8 pole Bessel filter) and digitized afterwards Sakmann and Neher (1995). This attenuates high frequent noise components, e.g., caused by shot noise, for an example see Figure 1. Therefore, lowpass filters are in general integrated in the hardware of the technical measurement device, since they are indispensable for noise reduction and visualization. Important channel characteristics such as amplitudes and dwell times can be obtained provided the conductance changes of the traces are idealized (underlying signal is reconstructed) from these recordings (data points) Colquhoun (1987); Sakmann and Neher (1995); Hotz et al. (2013).

Besides subjective and time-consuming analysis by visual inspection, often with manually chosen event times or in a semi-automatic way, as e.g. offered by pCLAMP 10 software (Molecular Devices), there is a demand for fully automatic methods, which provide idealizations on an objective basis, preferably with statistical guarantees for the recovered events. These automatic methods can be roughly divided into two groups.

HMM based methods

Firstly, methods that assume an underlying statistical model with a few parameters describing the gating dynamics are utilized. Most common are hidden Markov models (HMMs), see Ball and Rice (1992), where parameters correspond to transition probabilities and levels. These parameters can either be estimated by the Baum-Welch algorithm, see Venkataramanan et al. (2000); Qin et al. (2000), or by MCMC sampling, see de Gunst et al. (2001); Siekmann et al. (2011), or by approaches based on the conductance (current) distribution, see Yellen (1984); Heinemann and Sigworth (1991); Schroeder (2015) and the references therein. Although, they do not provide an immediate idealization of the underlying signal, the most likely path (idealization) can be obtained from these models, e.g., by the Viterbi algorithm Viterbi (1967) and related methods. In addition, and relevant to this paper, parameter estimation can often be corrected for missing short events, see Hawkes et al. (1990); Qin et al. (1996); Colquhoun et al. (1996); Epstein et al. (2016).

Model-free methods

Secondly, and more in J. Tukey’s spirit of letting the data speak for itself, non-parametric, often called model-free, methods that do not rely on a specific (parametric) model are used. Examples of such methods are amplitude thresholding Colquhoun (1987); Sakmann and Neher (1995), slope thresholding Basseville and Benveniste (1983); VanDongen (1996), idealization based on the minimal description length ($\operatorname{MDL}$) Gnanasambandam et al. (2017) and the jump segmentation multiresolution filter $\operatorname{JSMURF}$ Hotz et al. (2013).

The methodology suggested in this paper belongs to this second class of models. We stress, however, that it can also be seen as a companion for HMM based methods, e.g., to determine the number of states and for prefiltering purposes. This will be specified in the following.

Methods - compared and contrasted

Usually, gating ($1\text{\,}\mathrm{ns}100\text{\,}\mathrm{ns}$) is much faster than the sampling rates ($1\text{\,}\mathrm{kHz}100\text{\,}\mathrm{kHz}$) of the recordings and, hence, channel recordings have the appearance of abrupt random changes between states Hamill et al. (1981). Approaches that assume a HMM, e.g., conductance distribution fitting, allow to resolve the gating dynamics below the measured temporal scales, potentially. Nevertheless, often a rigorous HMM based analysis is hindered by several issues. To specify the Markov model for the gating dynamics it is required to fix the number of states. Although some model selection approaches exist, this is not an easy task, in general. Further, in rare cases the Markovian property itself is questionable, cf. Fuliński et al. (1998); Mercik and Weron (2001); Goychuk et al. (2005); Shelley et al. (2010) or, more commonly, the implicit homogeneity assumption underlying a hidden Markov model is violated, e.g., due to artifacts providing an inhomogeneous error structure. Hence, in the presence of artifacts like for instance base line fluctuations (cf. Figure 1) the assumption of few conductance levels that occur repeatedly requires elaborate data cleaning before an HMM can be fitted. In contrast, in such situations model-free methods may provide right away a reasonable idealization as they potentially can handle inhomogeneity in a more flexible way. Additionally, the computational complexity of many model-free methods increases almost linearly in the number of data points, an important issue, since often large data sets have to be analyzed. Finally, model-free methods, in addition to model selection approaches, can be used to select or verify a specific Markov model, in particular to determine the number of states and possible transitions, and to explore and potentially remove artifacts. Furthermore, fitting the conductance distribution requires often an initial guess for the parameters which can be provided by model-free methods, as well. Finally, if a Markov model is assumed, missing event corrections can also be combined with model-free methods to improve the parameter estimation based on their idealizations, see Section 4.4.

Consequently, hidden Markov models, fitting the conductance distribution, and model-free methods complement each other well and should not be seen as competing approaches. While model-free methods ease idealization for e.g. nonhomogeneous data, approaches based on a Markov model allows to extrapolate information from larger (observable) to smaller (not observable) time scales.

While the literature on hidden Markov models and fitting the conductance distribution is comprehensive, we see a demand for better working model-free approaches and its theoretical justification. Many approaches are based on additional filtering (often by lowpass Gauss filters) and manually chosen thresholds which are difficult to determine in an objective way. It is well known that filtering threatens to miss events in the idealization Colquhoun (1987); Sakmann and Neher (1995); Hotz et al. (2013) as the filter length predetermines the scales on which events can be detected. Further, such methods often ignore the lowpass filter of the amplifier. But ignoring this filtering threatens to include false positives which falsify the analysis as well.

Hence, providing rigorous statistical guarantees for the idealized traces is lacking to model-free methods, in general, and filtering hinders this additionally. The multiresolution method $\operatorname{JSMURF}$ Hotz et al. (2013) offers a model-free idealization which comes with statistical guarantees for its occurrences of conductance changes at a range of time scales. Roughly speaking, at a given error level $\alpha$ it guarantees detected conductance changes to be significant (controls the false positives). However, this method is limited to achieve accuracy at the filter resolution and beyond. As a consequence, it fails to detect flickering events (events shorter than the filter length), see Figure 3 for an illustration and the simulations in Section 4.

This is a crucial issue, as flickering typically has its own dynamics and can result from different molecular processes like conformational changes of the ion channel Grosse et al. (2014) or by the passage of larger molecules blocking the ion’s pathway through the channel Raj Singh et al. (2012). Thus, and to the best of our knowledge, no model-free method is able to detect flickering events reliably, i.e., at a given statistical accuracy, while at the same time detecting events at larger scales. Moreover, even if an event is detected, precise idealization of its amplitude and its dwell times requires new tools, since the convolution of the signal with the filter smoothes the peak, see Figures 2 and 5 for illustration and Colquhoun (1987); Sakmann and Neher (1995). An exception is the semi-automatic SCAN software Colquhoun and Sigworth (1995) which allows time-course fitting. This means to idealize events by least squares fitting based on an approximation of the Bessel filter kernel by a Gauss kernel, but interventions by the experimenter are required. Moreover, a conductance value can only be idealized by least squares fitting if the event is long enough, otherwise it is guessed by previously obtained values.

Therefore, to overcome this lack of methodology, in this paper we will present a model-free method for detection and precise idealization of event times and amplitudes of flickering events, while at the same hand all relevant events at larger time scales will be recovered as well.

JULES

We first introduce in Section 2 an extremely flexible and simple model for the recordings, which does not rely on a Markov assumption. Then, we propose based on this model a new JUmp Local dEconvolution Segmentation filter, $\operatorname{JULES}$, that consists of a (model-free) detection step by a multiresolution constraint (Section 3.1) and an estimation (idealization) step for the signal by local deconvolution (Section 3.2). The level of resolution (encoded in the multiresolution constraint) of $\operatorname{JULES}$ will be of the magnitude of the sampling rate and allows us even detection of events below the filter frequency. To this end, the final idealization of the time-continuous jump locations and in particular of the amplitudes is achieved by local deconvolution for which statistical guarantees can be provided as well.

Implementation

These steps are complex and are summarized in Algorithm 1. $\operatorname{JULES}$ is implemented as an R R Core Team (2016) function jules in the package clampSeg Pein et al. (2017b). The detection step is computed via a pruned dynamic program implemented in C++ to speed up computations. The deconvolution step is a maximum likelihood estimator computed by an iterative grid search algorithm, see Section 3.2.

Evaluation

The performance of $\operatorname{JULES}$ is assessed in several simulation studies (Section 4) and on gramicidin A (gA) traces (Section 5). The simulations indeed confirm that $\operatorname{JULES}$ is able to detect and idealize events that are short in time with very high precision. More precisely, we found that $\operatorname{JULES}$ detects (for the signal-to-noise ratio and filter setup of the later analysed gramicidin A traces) with high probability peaks with length at least one quarter of the filter length, see Figure 2 for an illustration and Figure 6 for a systematic analysis, and separates two consecutive peaks with high probability once the distance between them is at least three times the filter length, see Figure 8. Although not designed under a Hidden Markov assumption, Hidden Markov model simulations with two flickering states show that $\operatorname{JULES}$ is able to separate them if the amplitude of the two states differs at least by $3\text{\,}\mathrm{pS}$ (which is roughly two times the noise level and one sixth of the amplitude). In comparison, $\operatorname{JSMURF}$ misses all peaks below the filter length, see Figure 3 for illustration, whereas $\operatorname{TRANSIT}$ VanDongen (1996) includes many additional wrong detections, illustrated in Figure 4. Moreover, its idealization shifts peaks to the right and estimates too small amplitudes for shorter lengths when the signal is smoothed by the filter, see Figure 5(b) and Table 3.

JULES can be nicely combined with HMMs: If in addition the gating dynamics are assumed to follow a Markov model, even shorter dwell times can be analyzed by using missed event corrections. More precisely, for the simulated signal to noise ratio and the used filter setup, dwell times slightly below the sampling rate ($10$ kHz) are analyzable by $\operatorname{JULES}$, since then still enough events are long enough to be detected. As discussed before, methods that directly use the Markov model, in particular conductance distribution fitting, might be able to resolve faster gating, however at the price of for instance a much smaller robustness against model violations.

Data analysis

Gramicidn A is a small peptide produced by the soil bacterium Bacillus brevis showing antimicrobial activity O’Connell et al. (1990); Andersen et al. (2007). It forms an ion channel by interleaflet dimerization making its putative gating mechanism sensible to properties of the surrounding lipid membrane O’Connell et al. (1990); Andersen et al. (2005, 2007). Apart from the classical slow gating by dimerization also short closing events of the channel have been observed which could be linked to the physical and chemical properties of the lipid membrane; especially its hydrophobic thickness Armstrong and Cukierman (2002); Ring (1986).

Idealization of the gramicidin A traces supports the finding of the simulations with high evidence, see Figures 2-4. Zooming into single peaks demonstrates that $\operatorname{JULES}$ idealizes the observation extremely well (see Figure 2 lower panels), which also supports our underlying model. As all simulations (see above) haven been performed for models which are resembling our gramicidin A data sets very well, this establishes $\operatorname{JULES}$ as a reliable tool to identify and idealize flickering events in filtered ion channel recordings. We found that the flickering events are full-sized, i.e., have the same amplitude than the slow gating events, see Figure 12. In the near future we will continue with an analysis of further transmembrane proteins. Currently we examine for instance whether certain mutations influence the transport of antibiotics in a bacterial porin and analyze the gating dynamics of the translocase channel Tim23.

2. Modeling

As argued above (see also Hotz et al. (2013)) the conductivity will be modeled as a piecewise constant signal

[TABLE]

where $t$ denotes physical time. The (unknown) levels are denoted as $l_{0},\ldots,l_{K}$, its (unknown) number of conductivity changes as $K$ and the (unknown) changes between states as {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-\infty}=:\tau_{0}<\tau_{1}<\cdots<\tau_{K}<\tau_{K+1}:=\tau_{\operatorname{end}}. Note, this implies that we only consider conductivity changes after recording started at time [math] and assume for simplicity that the signal before recordings started is constant and equal to the first function value. When truncating the filter we will see that only a very short time period before recordings started will be relevant, details are given at the end of this section. For an analysis of how long the channel dwells in the same conductance level we exclude always the first and last segment. We stress that the class of signals in (1) is very flexible as potentially an arbitrary number of states and changes at arbitrary levels can be imposed, see again Figure 2 for an example. We assume further that the recorded data points $Y_{1},\ldots,Y_{n}$ (the measured conductivity at time points $t_{i}=i/f_{\boldsymbol{s}},\ i=1,\ldots,n$, equidistantly sampled at rate $f_{\boldsymbol{s}}$) result from convolving the signal $f$ perturbed by gaussian white noise $\eta$ with an analogue lowpass filter, with kernel $F$, and digitization at sampling rate $f_{\boldsymbol{s}}=n/\tau_{\operatorname{end}}$, i.e.,

[TABLE]

with noise level $\sigma_{0}>0$ and $\ast$ the convolution operator. Here $n$ denotes the total number of data (typically several hundred thousands up to few millions). Hence, the resulting errors $\epsilon_{1},\ldots,\epsilon_{n}$ are gaussian and centered, $\mathbb{E}[\epsilon_{i}]=0$, and have covariance $\mathsf{Cov}[\epsilon_{i},\epsilon_{i+j}]=\sigma_{0}^{2}(F\ast F)(j/f_{\boldsymbol{s}})$. For example the gA traces in our subsequent analysis (see Figure 2 and Section 5) are filtered by a 4-pole lowpass Bessel filter with $1\text{\,}\mathrm{kHz}$ cut-off frequency and sampled at $10\text{\,}\mathrm{kHz}$. The major aim will be now to idealize (reconstruct) the unknown signal $f$ taking into account the convolution and the specific structure of $f$ in (1). This will be done fully automatically and with statistically error control. By fully automatic we mean that no user action is required during the idealization process, only the error level $\alpha$ and two filter specific parameters have to be selected in advance, see Section 4.1. This is no restriction when the filter is known (as it is the case in most applications) and when it is unknown it can be pre-estimated from the data Tecuapetla-Gómez and Munk (2016).

Typically used lowpass filters have infinite impulse responses but a quickly decaying kernel $F$, for instance the kernels of Bessel filters decay exponentially fast Proakis and Manolakis (1996). Thus, we truncate the kernel $F$ at some point, say m{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}/f_{\boldsymbol{s}}}, which provides further simplification of our model. The trucated kernel is denoted by $F_{m}$ and rescaled such that $\int_{-\infty}^{\infty}{F_{m}}=1$, see Hotz et al. (2013). This is implemented in the R function lowpassFilter Pein et al. (2017b). We also truncate the autocorrelation function of the noise, resulting in an $m$-dependent process. All in all, model (2) reduces to

[TABLE]

where $\epsilon_{1},\ldots,\epsilon_{n}$ have covariance

[TABLE]

Note that we truncate the autocorrelation function itself as this leads to more accurate values for lags close to $m$ instead of convolving the truncated kernel. As a working rule, we choose $m$ such that the autocorrelation function {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(F\ast F)(\cdot)/(F\ast F)(0)} of the analogue lowpass filter is below $10^{-3}$. For gA traces, this choice leads to $m=11$ and is confirmed visually by comparing the observations with the convolution of the kernel with the idealized signal, see Figures 2 and 5.

3. Methodology: JULES

3.1. Detection step: multiresolution reconstruction

In this section we present our novel multiresolution methodology for model (3) that detects all significant conductance changes in the recordings and allows us, at the same time, to gain precise statistical error control of the reconstructed quantities. One major challenge is to take temporal scales below the filter frequency into account. We begin by disregarding for the moment the underlying convolution of the signal with the lowpass filter and developing a reconstruction method for the model

[TABLE]

Here $g$ is also a function of type (1), which may have different changes and levels. In other words, with respect to model (3), in a first step we aim for a piecewise constant approximation of $F_{m}\ast f$ which will then be further refined at the end of this section via postfiltering to a reconstruction of $f$. These two steps fix the number of conductance changes and their rough location. Afterwards, the final idealization $\operatorname{JULES}$ of $f$ will be obtained by local deconvolution in Section 3.2. A summary of these steps is given in Algorithm 1.

Observe that with (5) we intend to model an unknown piecewise constant function which is perturbed by an $m$-dependent (colored) noise. Its correlation structure, however, is given by the filter and therefore explicitly known, see (4). Alternatively, it could be preestimated as in Tecuapetla-Gómez and Munk (2016). We preestimate the noise level $\sigma_{0}$ by

[TABLE]

used in Hotz et al. (2013) and implemented in the R function sdrobnorm Pein et al. (2017a). At this point (a refinement will be provided in the local deconvolution step) we restrict all change-point locations to the grid on which the observations are given, in other words, we assume that $f_{\boldsymbol{s}}\tau_{i}$ are integers.

Multiresolution statistics

Suppose that we are given a piecewise constant candidate $\tilde{g}$ from model space (1) and we have to decide whether or not $\tilde{g}$ is a good fit (reconstruction) for the unknown signal $g$. To this end, we test simultaneously on all scales (resolution levels) whether $\tilde{g}$ fits the data appropriately. More precisely, $\tilde{g}$ is a suitable candidate if

[TABLE]

where $p_{i,j}=\sqrt{2\log(e(j-i+1)/n)}$ is a penalty term which is used to set different scales on equal footing, $\sigma_{i,j}:=\sqrt{\mathsf{VAR}(\epsilon_{i}+\cdots+\epsilon_{j})}=\sigma_{0}\,\sqrt{(j-i+1)(F\ast F)(0)+2\,\sum_{k=1}^{m}{(j-i+1-k)_{+}(F\ast F)(k/f_{\boldsymbol{s}})}}$ and $q=q(\alpha)$ is a critical value that can be obtained in a universal manner by Monte Carlo simulations such that (7) is a level $\alpha$-test. Here $\alpha$ is fixed in advance (by the experimenter) to control the number of fluctuations of the solution $\tilde{g}$. For details how to generate observations from model (3) to simulate $q$, see Section 4. However, to keep it computationally feasible we only evaluate the maximum in (7) over the system of all intervals that contain a dyadic number of observations, i.e., the maximum in (7) is only taken over all $1\leq i\leq j\leq n$ such that $j-i+1=2^{l}$ for some $l\in\mathbb{N}_{0}$. This reduces the complexity of the system from $n(n+1)/2\in\mathcal{O}(n^{2})$ to $\lfloor\log_{2}(n)\rfloor(2n-\lfloor\log_{2}(n)\rfloor+1)/2\in\mathcal{O}(n\log(n))$ intervals.

In order to provide a reconstruction for $g$ we proceed as follows. Let $G_{k}$ denote the class of all candidates in (1) with $K=k$ changes, $k\geq 0$, and let $\hat{K}$ denote the minimal $k$ for which there exits a $\tilde{g}\in G_{k}$ such that $M(\tilde{g})\leq q$, cf. (7). Then, the reconstruction for $g$ is given by

[TABLE]

Such an estimator has many favorable properties, for instance it detects changes at the optimal detection rate and allows for confidence statements. For more details to the theory behind multiresolution tests as in (7) and the choice of the penalty see e.g. Siegmund and Yakir (2000); Dümbgen and Spokoiny (2001); Dümbgen and Walther (2008) and for estimates of type (8) see Boysen et al. (2009b); Frick et al. (2014a) and the references therein.

Dynamic Programming

The estimator in (8) can be computed by dynamic programming as described in section 3 of Frick et al. (2014a). For related ideas, see also Killick et al. (2012); Li et al. (2016); Maidstone et al. (2017) and the references given there. However, its worst case computation time is quadratic in the number of observations. In order to speed up computations we included several pruning steps in the backward as well as in the forward step of the dynamic program, for more details see section A.1 in the supplement of Pein et al. (2016). In the situation of routinely changes (which is typical for patch clamp recordings) the runtime of such a program is almost linear in the number of observations such that typical traces (several hundred thousands up to several million data points) can be processed in less than a minute on a standard laptop. An implementation of $\operatorname{JULES}$ (including the local deconvolution step, see next section) can be found in the CRAN-R package clampSeg Pein et al. (2017b).

Postfiltering

Our preliminary reconstruction disregarded that also the signal is convolved with the filter and hence we obtained a piecewise constant approximation on the convolution $F_{m}\ast f$ instead of a reconstruction of the signal $f$ itself. However, they only differ at the beginning of each segment, i.e., when $f$ changes abruptly at $\tau$ from $l_{0}$ to $l_{1}$, then $F_{m}\ast f$ changes on $(\tau,\tau+m/f_{\boldsymbol{s}})$ from $l_{0}$ to $l_{1}$ and is constant afterwards (as long as no other change occur). As a consequence the reconstruction $\hat{g}$ in (8) potentially changes in incremental steps, see Li et al. (2017) for a mathematically rigorous statement. To correct for this we merge a segment with all subsequent segments for which the distance between the two starting points is less than $m/f_{\boldsymbol{s}}$ and if all changes are in the same direction. More precisely, we set the starting point of the merged segments to the starting point of the first segment and its level to the level of the last segment.

3.2. Estimation step: local deconvolution

In this parameter estimation step we obtain the final idealization by determining the precise changes $\tau_{1},\ldots,\tau_{\hat{K}}$ and levels $l_{0},\ldots,l_{\hat{K}}$ by local deconvolution. This step is crucial for idealizing (reconstructing) flickering events with high accuracy. To this end, we split the constant segments in short and long ones. Basically, long segments are those in which we can determine its level without taking the lowpass filter explicitly into account, details are given below. Pre-estimating the level of long segments enables to deconvolve the signal locally, i.e., to compute the (regularized) maximum likelihood estimator of the remaining parameters (change-point locations and levels) by an iterative grid search, see below for algorithmic details.

Long segments

We will define a segment $[\hat{\tau}_{i},\hat{\tau}_{i+1})$ as long if we can determine its level $\hat{l}_{i}$ without taking the lowpass filter explicitly into account. If the signal $f$ changes at $\tau_{i}$ its convolution $F_{m}\ast f$ will change continuously on $(\tau_{i},\tau_{i}+m/f_{\boldsymbol{s}})$ and will be constant on $[\tau_{i}+m/f_{\boldsymbol{s}},\tau_{i+1}]$. Hence, the detection step will (if the change is detected) find a change in the interval $[\tau_{i},\tau_{i}+m/f_{\boldsymbol{s}}]$. The other way around, a detected change $\hat{\tau}_{i}$ implies a change in $\hat{f}$ in the interval $[\hat{\tau}_{i}-m/f_{\boldsymbol{s}},\hat{\tau}_{i}]$ and, hence, a constant convolution $F_{m}\ast\hat{f}$ on $[\hat{\tau}_{i}+m/f_{\boldsymbol{s}},\hat{\tau}_{i+1}-m/f_{\boldsymbol{s}}]$. If the observations in this interval are enough to determine the level, say ten or more, by the median, i.e., $\hat{l}_{i}=\operatorname{median}(Y_{f_{\boldsymbol{s}}\hat{\tau}_{i}+m},\ldots,Y_{f_{\boldsymbol{s}}\hat{\tau}_{i+1}-m})$, we define the segment as long.

Local deconvolution

Knowing these levels enables us to perform the deconvolution locally, i.e., only few observations have to be taken into account and only few parameters have to be estimated. More precisely, let $[\hat{\tau}_{i-1},\hat{\tau}_{i})$ and $[\hat{\tau}_{j},\hat{\tau}_{j+1})$, $i\leq j$, be two consecutive long segments (no long segment inside, but potentially some short segments). Then, we aim to determine the precise locations of the changes $\hat{\tau}_{i},\ldots,\hat{\tau}_{j}$ and the precise levels $\hat{l}_{i},\ldots,\hat{l}_{j-1}$. Note that the levels $\hat{l}_{i-1}$ and $\hat{l}_{j}$ are already determined. And recall that a detected change-point $\hat{\tau}$ implies a change in $[\hat{\tau}-m/f_{\boldsymbol{s}},\hat{\tau}]$ and that the level $\hat{l}_{i}$ affects the convolution on $(\hat{\tau}_{i},\hat{\tau}_{i+1}+m/f_{\boldsymbol{s}})$. Therefore, we only have to take into account the observations $Y_{f_{\boldsymbol{s}}\hat{\tau}_{i}-m+1},\ldots,Y_{f_{\boldsymbol{s}}\hat{\tau}_{j}+m-1}$ and maximize its likelihood in the unknown parameters which is equivalent to minimize

[TABLE]

with $Y_{i,j}:=(Y_{f_{\boldsymbol{s}}\hat{\tau}_{i}-m+1},\ldots,Y_{f_{\boldsymbol{s}}\hat{\tau}_{j}+m-1})^{T}$ the vector of the observations,

[TABLE]

the vector of the expected conductance levels $K\ast\hat{f}$ for a candidate signal $\hat{f}$ and $\Sigma_{i,j}$ the (known) correlation matrix of the observation vector, see the explanations after (3). For the signal $\hat{f}$ the number of change-points and the function values $\hat{l}_{i-1}$ and $\hat{l}_{j}$ are fixed by the prior reconstruction and the change-point locations are restricted to the intervals $[\hat{\tau}_{k}-m/f_{\boldsymbol{s}},\hat{\tau}_{k}]$, $k=i,\ldots,j$. The minimization is performed by an iterative grid search, see the next but one paragraph.

Regularization

For the Bessel filter used in the later analyzed gA traces we have to regularize the correlation matrix $\Sigma_{i,j}$ in (9), since the matrix is ill-conditioned, i.e., has a condition number around $10^{-4}$. This characteristic is typical for most analogue lowpass filters. Hence, small errors in the model (for instance in ion channel recordings often base line fluctuations occur) will be amplified by inversion of $\Sigma_{i,j}$ and may result in a large missestimation. A standard approach to cope with this is Tikhonov / $L_{2}$ regularization, i.e., we replace $\Sigma_{i,j}$ by $\Sigma_{i,j}+\gamma^{2}I$, with $I$ the identity matrix and $\gamma^{2}\geq 0$ a regularization parameter, see e.g. (Engl et al., 2000, chapter 5). Such a regularization allows quick computation of (9) by the Cholesky decomposition (which can be stored) and by solving a triangular system of equations. Note that the correlation structure is determined by the filter and hence the parameter $\gamma^{2}$ can be chosen universally for all recordings on the same system. For our data we choose the regularization parameter such that the convolution fits well to the recordings, which is nicely confirmed by visual inspection as illustrated in Figure 2. This precalibration step has to be done only on a small data excerpt in advance. We found that the results are robust in $\gamma^{2}$, and simply setting $\gamma^{2}=1$ led to satisfactory results in all cases. For automatic choices, e.g. by cross-validation, and for other regularization, e.g. truncated SVD, see Vogel (2002) and references therein.

Iterative grid search

In this paragraph we describe how we compute the (local) maximum likelihood estimator by an iterative grid search. In general, computing the estimator, i.e., performing the deconvolution, is difficult due to the non-convexity of the optimization problem to minimize (9). Grid search means that we fix for each change $\hat{\tau}_{i},\ldots,\hat{\tau}_{j}$ a set of possible locations, the grid, and compute (9) for all combinations and take the solution with the minimal value. Note that for given changes close formulas for the optimal levels $\hat{l}_{i},\ldots,\hat{l}_{j-1}$ exist which allows fast computation of (9). In theory such an optimization could be done for arbitrarily fine grids. However, to keep it computationally feasible, we start with the observations grid, i.e., $\{\hat{\tau}_{k}-m/f_{\boldsymbol{s}},\ldots,\hat{\tau}_{k}\}$ for the change $\hat{\tau}_{k},k=i,\ldots,j-1$. Recall that a detected change $\hat{\tau}_{i}$ implies a change in $f$ in the interval $[\hat{\tau}_{i}-m/f_{\boldsymbol{s}},\hat{\tau}_{i}]$. Afterwards, a refinement can be done by repeating the grid search iteratively with finer and finer grids around the changes found in the previous step. We found that this refinement for determining the changes with arbitrarily precision works very well in practice. We found it sufficient to iterate such a refinement twice, each with a ten times finer grid between the neighboring grid points of the candidate for the location of the change.

This approach is computationally feasible for a jump between two long segments, i.e., only the jump location is unknown, and for a peak, i.e., a short segment between two long segments, here the unknowns are the level of the short segment and the two jump locations of the short segment between two long segments. See Figures 2 and 5 for illustration. In rare situations more than one short segment between two long segments occur. The same approach could then be done for $k$ consecutive short segments by looking for $k+1$ changes between two long segments and determining all $k$ local levels. However, such an optimization is time consuming and not necessary for our data. Hence, we abandon and simply drop these segments from the analysis. For the idealization we keep the reconstruction from the detection step. This leads to potentially miss of events if they are too close together. Figure 8 shows that for the analyzed gramicidin A traces this only occurs if the distance between two events is less than $3.2\text{\,}\mathrm{ms}$ (roughly $3$ times the filter length). This is much smaller than estimated average distance between two events of $0.3\text{\,}\mathrm{s}$ and hence the whole effect seems indeed negligible.

4. Simulations

The purpose of this section is threefold. We will assess first the ability of $\operatorname{JULES}$ in detecting and idealizing an isolated peak. Then, we will identify the minimal distance at which $\operatorname{JULES}$ is able to separate two consecutive peaks. Finally, although $\operatorname{JULES}$ does not rely on a Hidden Markov model assumption we will examine the ability of $\operatorname{JULES}$ to recover the parameter of a Hidden Markov model with one open state and two closed flickering states, in particular how well the two flickering states can be separated. We will generate all signals and observations accordingly to (1) and (3) such that they are in line with the data we analyze in Section 5. We also choose amplitudes, dwell time and noise level of the generated observations such that they are similar to those of the datasets of our applications and we also simulate a 4-pole Bessel filter with $1\text{\,}\mathrm{kHz}$ cut-off frequency and sample the observation at $10\text{\,}\mathrm{kHz}$. For data generation we use the Durbin-Levinson algorithm (Brockwell and Davis, 2006, proposition 5.2.1) to compute the coefficients of the moving average process corresponding to the covariance of the errors in (3) and compute the convolution of the simulated signal with the truncated kernel $F_{m}$ exactly. We also examine briefly robustness issues against $f^{2}$ and $1/f$ noise.

For comparison, we also apply the slope thresholding method $\operatorname{TRANSIT}$ VanDongen (1996) and the jump segmentation multiresolution filter $\operatorname{JSMURF}$ Hotz et al. (2013). To this end, we used the R functions transit (with default parameters) and jsmurf (also with the default interval system of dyadic lengths and $\alpha=0.05$) in the package stepR Pein et al. (2017a). We also briefly investigate the performance of $\operatorname{MDL}$ Gnanasambandam et al. (2017) for which we used the Matlab code provided by the authors.

4.1. Parameter choices

The following specifications will be used in this and in the upcoming section where we analyze the gA traces. For the detection step (multiresolution reconstruction) of $\operatorname{JULES}$ we aim to control the probability to overestimate the number of change-points at significance level $\alpha$ and hence use a conservative choice of $\alpha=0.05$, resulting in a quantile of $q=1.4539$. In the estimation step (local deconvolution) of $\operatorname{JULES}$ we use a regularized correlation matrix with regularization parameter $\gamma^{2}=1$, see Section 3.2. And, as mentioned before, we truncate the kernel and autocorrelation function of the filter at $m=11$ which corresponds to the fact that the autocorrelation function is below $10^{-1}$ afterwards. These choices are used as default by the R function jules in the package clampSeg Pein et al. (2017b).

4.2. Isolated peak

In this simulation with $4\,000$ observations we examine the detection and idealization of a single isolated peak. We consider a situation that is comparable to the analyzed data in Section 5. More precisely, we simulate from model (3) with a signal with levels $l_{0}=l_{2}=40$ and $l_{1}=20$ and changes at $\tau_{1}=2000/f_{\boldsymbol{s}}$ and $\tau_{2}=(2000+\ell)/f_{\boldsymbol{s}}$, c.f. (1) and Figure 5. The standard deviation of the errors $\epsilon_{i}$ in (3) used in each simulation is $\sigma_{0}=1.4$. We are interested in the performance of the methods in detecting the peak and estimating the change-point locations $\tau_{1}$ and $\tau_{2}$ and the level $l_{1}$ as a function of $\ell$, the length (relative to the sampling rate $f_{\boldsymbol{s}}$) of the peak.

Tables 1-3 summarize our results based on $10\,000$ repetitions for $\ell=2,3,5$. For $\ell=5$, Figure 5 shows an example of the simulated data as well as the idealizations by $\operatorname{JULES}$, $\operatorname{TRANSIT}$ and $\operatorname{JSMURF}$ and their convolutions with the Bessel filter in a neighborhood of the peak. Finally, Figure 6 shows the relative amount of true detections among all detections for $\ell\in\{0,0.1,\ldots,5\}$. To this end, we count in how many simulations the signal is correctly identified, i.e., only the peak and no other change is detected. More precisely, we define the peak as detected if there exists a $j$ such that $|\hat{\tau}_{j}-\tau_{1}|<m/f_{\boldsymbol{s}}$ and $|\hat{\tau}_{j+1}-\tau_{2}|<m/f_{\boldsymbol{s}}$ as a peak is shifted at most $m/f_{\boldsymbol{s}}$ by the filter. If only one change but not a peak is within these boundaries (see the zooms in Figure 3 for an illustration) we do not count it as a true detection, but also not as a false positive, whereas all other changes are counted as false positives. For the estimated change-point locations and the level we only consider cases where the peak is detected and report the mean squared error, the bias and the standard deviation. We also report trimmed versions for the estimated level, where we compute these error quantities only for all values between [math] and $40$.

We found that $\operatorname{JULES}$ detection power increases with the length $\ell$ in range between $0.8$ to $2.4$ times the sampling rate and is (almost) one for larger lengths while at the same time it provides almost no false positives, resulting in a relative amount of true detections of almost one for lengths larger than one quarter of the filter lengths. In comparison, $\operatorname{TRANSIT}$ has a slightly larger detection power for very small lengths, but detects much more false positives, resulting in almost no correctly identified signals. Moreover, its detection power fluctuates for larger lengths in an uncontrollable way. $\operatorname{JSMURF}$ is not able to detect such small filtered peaks as it does not take into account the corresponding scales, but also does not detect false positives. We also briefly investigated $\operatorname{MDL}$ Gnanasambandam et al. (2017), see Table 1. To decrease the amount of false positives the authors suggested to assume a minimal length for the events. Although this might be problematic for real data applications with short events, we used for the simulations with $l=5$ the assumption of that length as prior information. But even then the number of events is heavily overestimated. More precisely, on average $98.87$ change-points are discovered instead of the two change-points of the peak. Note, that the simulated observations are filtered by a 4-pole lowpass Bessel filter with normalized cutoff frequency of $0.1$ and as noted by Gnanasambandam et al. (2017) their approach works better when recordings are unfiltered or only slightly filtered.

$\operatorname{JULES}$ idealizes almost all detected peaks with high accuracy which increases with the length of the peak. Only in very rare cases (less than $0.2\%$) and only if the peak is short the amplitude is heavily overestimated, resulting in large mean squared error, bias and standard deviation, as the corresponding trimmed values are much smaller. We still found that the convolved idealization fits the observations always very well, see Figure 5(a) for an illustration. We will also see in Section 4.4 that $\operatorname{JULES}$ is nevertheless able to identify two levels close to each other. In comparison, Figure 5(b) and the biases in Tables 2 and 3 show that $\operatorname{TRANSIT}$ miss-estimates the peak systematically, in particular the amplitude is underestimated if the peak is smoothed by the filter. The major difficulty of the devonvolution problem is that convolutions of signals with larger amplitude but smaller length can look very similar to the convolution of a signal with smaller amplitude but larger length. In fact, it is possible to show that the change-points locations $\tau_{i}$ can by no method be estimated better than at $1/\sqrt{n}$ rate (instead of the sampling rate $1/n$ without convolution), see Boysen et al. (2009a) for a similar setting and also Frick et al. (2014b); Goldenshluger et al. (2006) for further theoretical results of estimation a peak from filtered data. Note that the filter is rather short ranged with a filter length of eleven but the signal to noise ratio is rather small with $20/1.4\approx 14.29$, since gA has a small conductance in comparison to other proteins, and the lengths are very short with two to five observations.

Similar results were obtained in additional simulation studies (not included) with different amplitude and noise level. Also a simulation (not displayed) where we shift the two changes by $0.1/f_{\boldsymbol{s}},0.25/f_{\boldsymbol{s}},0.5/f_{\boldsymbol{s}},0.75/f_{\boldsymbol{s}}$ and $0.9/f_{\boldsymbol{s}}$ leads to almost identical results for $\operatorname{JULES}$. This confirms the ability of $\operatorname{JULES}$ to identify correctly the location (and amplitude) of changes between the sampling points.

4.3. Separation of two consecutive peaks

In this simulation study we are interested in the performance of $\operatorname{JULES}$ in separating two consecutive peaks as a function of $d$, the distance between them. Separation is necessary in the detection step, i.e., two peaks have to be detected, but also in the deconvolution step as we are only able to deconvolve both peaks individually if they are isolated. To this end, we consider a signal $f$ with change-points $\tau_{1}=2\,000/f_{\boldsymbol{s}}$, $\tau_{2}=\tau_{1}+5/f_{\boldsymbol{s}}$, $\tau_{3}=\tau_{2}+d$ and $\tau_{4}=\tau_{3}+5/f_{\boldsymbol{s}}$, with $\tau_{0}=0$ and $\tau_{\operatorname{end}}=4\,000/f_{\boldsymbol{s}}$ and levels $l_{0}=l_{2}=l_{4}=$40\text{,}\mathrm{pS} and $l_{1}=l_{3}=$20\text{\,}\mathrm{pS}. Here we fixed the length of the consecutive peaks at $5/f_{\boldsymbol{s}}$, since want to focus on separation only and we found in Section 4.2 that peaks of this size can be detected with (almost) probability one. In this section we focus only on $\operatorname{JULES}$ as separation signifies an issue for this method.

We identified three outcomes for $\operatorname{JULES}$, all illustrated in Figure 7. First of all, perfect separation, i.e., the multiresolution reconstruction of $\operatorname{JULES}$ identifies the two peaks (4 change-points) and the local deconvolution yields idealizations for the four levels. Secondly, separation fails in the detection step, i.e., the multiresolution reconstruction recognizes only 2 change-points and identifies one peak whose level can be further deconvolved. Finally, separation fails in the deconvolution step, i.e., $\operatorname{JULES}$ identifies two peaks but the distance is so small that the deconvolution method cannot separate them, in other words, no long segment is in between.

Figure 8 shows the frequency at which each scenario occurred as a function of $d$, the distance between the two peaks, in $10\,000$ simulations for each value of $d=\{1,2,\ldots,43\}$. We found that the detection of both peaks fails potentially when $d\leq 12$, but the frequency decreases rapidly with $d$. For $d\in[4,31]$ detection is often possible, but not the separation in the deconvolution, whereas for $d>33$ this scenario is no longer observed. Finally, separation of the two peaks is possible with high probability as soon as the distance between them is at least $32/f_{\boldsymbol{s}}$, roughly three times the filter lengths. This is equal to the minimal distance of a long segment $10$ plus the two shifts on the left and right side of the segment by $m=11$ to take the filter into account. This corresponds to $3.2\text{\,}\mathrm{ms}$ for the gA traces, whereas the estimated average distance is with 1\text{/}3.28\text{\,}\mathrm{s}$\approx$0.3\text{\,}\mathrm{s} (see Section 5) much larger. Moreover, we found that the distance required to detect two separated peaks is only roughly three times the length that is required to detect a single, isolated peak (confer Figure 6). This shows that $\operatorname{JULES}$ separates peaks well with respect to detection.

4.4. Hidden Markov model

In this section we simulate data from a three state Hidden Markov model. None of the considered methods require a Hidden Markov assumption, but as such a model is often assumed for ion channel recordings, it is instructive to investigate these methods in such a scenario. More precisely, we simulate one open state with $40\text{\,}\mathrm{pS}$ and two closed states with $20\text{\,}\mathrm{pS}$ and 20+\Delta$$\text{\,}\mathrm{pS}, i.e., the amplitudes of the two flickering states ($20\text{\,}\mathrm{pS}$ and 20-\Delta$$\text{\,}\mathrm{pS}) differ by $\Delta$. The dwell time in the open state is exponentially distributed with rate $2.5$ and the channel switches to both closed states equally likely. For both closed states the channel reopens quickly with rate $800$ ($0.8\text{\,}\mathrm{kHz}$). We generate five traces with $600\,000$ observations, each.

In the following we analyze these data in a similar fashion as we will do it for the gramicidin A traces in Section 5. In accordance with the definition of a short segment in Section 3.2 we define a closing event as a flickering event if its dwell time is smaller than or equal to $2.6\text{\,}\mathrm{ms}$, but results are qualitative the same (amplitude histograms are even quantitative very similar) if we increase or decrease this threshold within a reasonable range. We also exclude events with a dwell time smaller than $0.24\text{\,}\mathrm{ms}$ as we cannot detect such events reliably and they would disturb the analysis. We do not include results for $\operatorname{JSMURF}$, since it is not able to detect such small flickering events, c.f. Figure 6. Figure 9 shows histograms of the idealized amplitudes of flickering events with a small bin width of $0.5$ (visually chosen) to already see small indications of two different levels.

In Table 3 we found that $\operatorname{JULES}$ level estimate has a large mean squared error, mainly due to heavy overestimations in rare cases. Nevertheless, we see in this simulation that $\operatorname{JULES}$ idealizes and separates the two amplitudes very well. For $\Delta=2$ a tiny indication of two different amplitudes can be seen, for $\Delta=3$ the two peaks are distinct and for larger differences the two different states can be clearly detected. In comparison, $\operatorname{TRANSIT}$ finds some smaller amplitudes and hence the separation is slightly less clear, but still possible for $\Delta\geq 3$. We already saw in Figure 5 and Table 3 that $\operatorname{TRANSIT}$ estimates a too small amplitude if a short peak is smoothed by the filter.

Figure 10 shows the dwell times in the closed state and an exponential fit. More precisely, we only consider events with a dwell time such that we can detect the events reliably, i.e., events with dwell time between $0.24$ and $2.6\text{\,}\mathrm{ms}$ and estimate the rate by maximizing the likelihood of these events assuming an exponential distribution for the dwell times. Note that the results for the dwell time and also for the later analyzed distances are similar for all considered $\Delta$. We show results exemplary for $\Delta=2$.

Apart from the fact that both methods miss extremely short events ($<0.3\text{\,}\mathrm{ms}$) which coincides with Figure 6, $\operatorname{JULES}$ confirms an exponential distribution and estimates with $0.9$ roughly the right parameter. $\operatorname{TRANSIT}$ detects additional spurious events with lengths between $0.4$ and $0.5\text{\,}\mathrm{ms}$ and hence overestimates the rate with $1.08$ slightly more. A similar behavior could be observed in Table 2, too.

Figure 11 shows the distance between two flickering events which coincides with the dwell time in the open state if other closing events are considered as spurious events. We include events with a distance between $0.032$ and $1\text{\,}\mathrm{s}$ as Figure 8 shows that $\operatorname{JULES}$ is not able to separate peaks with a smaller distance. However, such a fit is not enough as we miss events. Hence, a precise recovery of the rates requires to recalculate the distribution by taking into account missed events, see for instance Hawkes et al. (1990); Colquhoun et al. (1996); Qin et al. (1996); Epstein et al. (2016). In this case simply dividing the estimated rate by the probability that an event has dwell time between $0.24$ and $2.6\text{\,}\mathrm{ms}$ assuming an exponential distribution with the fitted rate appears to be appropriate.

For both methods the distribution of the distances seems to be exponential, however, due to the missing events the rate $2.5$ cannot be recovered by a simple exponential fit as it is the case for the dwell times. For $\operatorname{JULES}$ we find a rate of $1.86$ and with $\operatorname{TRANSIT}$ of $1.88$. Including the missed events correction leads with $2.69$ and $2.72$, respectively, to satisfactory results.

4.5. Robustness

While filtered white noise is already a good approximation for the noise occurring in the data analyzed here, see Section 5, additional high frequent $f^{2}$ (violet) and long tailed $1/f$ (pink) noise components have been observed in patch clamp recordings as well, for a more detailed discussion see Neher and Sakmann (1976); Venkataramanan et al. (1998b); Levis and Rae (1993) and the references therein. Moreover, heterogeneous noise can be caused by open channel noise Venkataramanan et al. (1998b), i.e., the noise level can be larger when the conductance is larger. Therefore, we focus in this section on analyzing how robust $\operatorname{JULES}$ is against such noise. To simplify the data generation we use a discrete convolution for the filter. For the violet noise we use as suggested by Venkataramanan et al. (1998a) a moving average process with coefficients $0.8$ and $-0.6$. For the pink noise we use the algorithm available on https://github.com/Stenzel/newshadeofpink. In both cases we scale the noise such that each single observation has standard deviation $1.4$ (to make it comparable to the simulations in the previous sections) and such that half of it results from the white noise and half of it from the additional noise component. For the heterogeneous noise case we simulate standard deviation $1.4$ for the background noise and $2.8$ for the events by oversampling the recordings by a factor of $100$.

From Table 4 we conclude that $\operatorname{JULES}$ is very robust against the additional $f^{2}$ and heterogeneous noise but influenced by $1/f$ noise. At presence of the latter noise, the global standard deviation is underestimated by (6) which leads to a larger detection power, but at the price of additional false positives. Note, that false positives are caused by the underestimated standard deviation but also by the long range dependency itself. However, these false positives have a small amplitude and therefore do not influence the analysis severely or can be removed by postfiltering. At presence of heterogeneous noise the detection power is slightly decreased due to the larger variance. However, we stress that $\operatorname{JULES}$ is only robust against heterogeneous noise as long as only short events occur. If we would have also some longer events the non-stationarity caused by open channel noise would be a severe issue. Parameter estimation (not displayed) is slightly worse at presence of $f^{2}$ or heterogeneous noise, but not affected by presence of $1/f$ noise.

5. Data analysis

5.1. Measurements

Gramicidin A was obtained from Sigma-Aldrich (Schnelldorf, Germany) and used without further purification. 1,2-Diphytanoyl-sn-glycero-3-phosphocholine (Avanti Polar Lipids, Alabaster, USA) and cholesterol (Sigma) were dissolved in chloroform and mixed in a molar ratio of 9:1. Giant-unilamellar vesicles (GUVs) containing the peptide (peptide to lipid ratio of 1:10000 up to 1:1000) were obtained by electroformation Angelova and Dimitrov (1986). Voltage clamp experiments of solvent-free lipid bilayers were performed using the Port-a-Patch (Nanion Technologies, Munich, Germany) in buffer ($1\text{\,}\mathrm{\textsc{m}}$ KCl, $10\text{\,}\mathrm{m\textsc{m}}$ HEPES, pH 7.4). Electrically insulating membranes with resistances $>1\text{\,}\mathrm{G\SIUnitSymbolOhm}$ were obtained by spreading GUVs on the aperture of a glass chip applying a slight negative pressure. A DC potential of $+$100\text{,}\mathrm{mV}$$ was applied and data was recorded at a sampling rate of $10\text{\,}\mathrm{kHz}$ using a $1\text{\,}\mathrm{kHz}$ 4-pole Bessel filter.

5.2. Idealization

Idealizations are obtained by $\operatorname{JULES}$, $\operatorname{TRANSIT}$ and $\operatorname{JSMURF}$ with parameter choices as in Section 4.1. The recorded observations in Figure 1 show gating events between two states on various time scales, but also several noise effects like outliers or varying conductivity, compare for instance the conductivity from $10$ to $25$ seconds and from $55$ to $60$ seconds. Note that such effects raise substantial difficulties for methods that assumes a Hidden Markov model or similar models. Contrary, $\operatorname{JULES}$ (Figure 2) provides a reasonable idealization covering all major features of the data and some smaller fluctuations. In comparison, $\operatorname{TRANSIT}$ (Figure 4) is able to detect such events too but at the same hand detects far too many false positives. $\operatorname{JSMURF}$ (Figure 3) works well on larger time scales but misses flickering events below the filter length. The zooms into single peaks (Figures 2-4 lower panels) demonstrates that $\operatorname{JULES}$ fits the observation well, whereas $\operatorname{TRANSIT}$ estimates a too small amplitude due to the smoothing of the filter and $\operatorname{JSMURF}$ misses the peaks.

5.3. Analysis of flickering dynamics

In this section we analyze the flickering dynamics of gramicidin A based on the obtained idealizations in a similar fashion as the simulated data in Section 4.4. Note that flickering does not occur in all data traces, for this analysis we focus on five traces, each with $600\,000$ data points, which show significant flickering. As before we define a closing event as a flickering event if its dwell time is smaller than or equal to $2.6\text{\,}\mathrm{ms}$. For comparison we also study briefly the amplitude of the slow gating (closing events between segments with dwell time longer than $10\text{\,}\mathrm{ms}$). We focus on events with an amplitude larger than $\geq 10\text{\,}\mathrm{pS}$ as events with a smaller amplitude are more likely base line fluctuations than channel gating. Base line fluctuations are for instance caused by small defects in the membrane, which is unavoidable in the recordings. There might be also periodic oscillations, resulting from the electronic or from building vibrations (although damped). Figure 12 compare the idealized amplitudes of flickering events between $10$ and $30\text{\,}\mathrm{pS}$ ($307$ events) with the amplitudes of slow gating events ($44$ events) of the same magnitude by kernel density estimates. The latter are performed by the R function bkde in the package KernSmooth Wand and Jones (1994) with bandwidth $2$ (visually chosen).

We see no distinct peak as in the simulations in Figure 9. Note that this also true for histograms instead of kernel density fits, but kernel density fits allow a better comparison of flickering and the gating on larger time scales. Hence, either multiple levels are underlying or more likely additional errors occur. Recall for instance the conductivity fluctuations in Figure 1. More importantly, we found that the flickering events are full-sized (have the same amplitude as the slow gating events). In comparison, $\operatorname{TRANSIT}$ suggests a smaller amplitude which can be explained by the fact that $\operatorname{TRANSIT}$ ignores the smoothing effect by the filter. We already saw this effect in the simulations in Section 4, see Figures 5 and 9 as well as Table 3. We do not show results for $\operatorname{JSMURF}$, since only very few flickering events have been detected, also this was observed before in Table 1 and Figure 6 in the simulations.

For the same reason we exclude $\operatorname{JSMURF}$ from the following analysis of the dwell times and distances between two flickering events. Figure 13 shows the dwell times in the closed state and Figure 14 shows the distance between two flickering events. The latter coincides with the dwell time in the open state if other events, in particular events with smaller amplitude, are considered as noise and not as gating.

Although the lipid system is totally different, our results are of a similar order than in Ring (1986); Armstrong and Cukierman (2002). Flickering events occur on average roughly every second and are only around a millisecond long. We miss extremely short events ($<0.24\text{\,}\mathrm{ms}$) in accordance with the simulations, see Figures 6 and 10. Up to this, we can confirm an exponential distribution with rate $0.73$ (in $\text{\,}\mathrm{ms}$) for the dwell times. In comparison, $\operatorname{TRANSIT}$ detects additional short events which, however, seem to be artifacts, see Figure 4 and Table 1 in the simulations, resulting in a rate of $1.46$.

From Figure 14 we draw that too many short distances are observed for a good exponential fit. Since such an effect was not observed in the simulations, cf. Figure 11, we speculate that this is caused by artifacts. Using also here a missed event correction leads to an estimated rate of $3.28\text{\,}\mathrm{Hz}$ based on the idealization by $\operatorname{JULES}$. In contrast, when we use $\operatorname{TRANSIT}$ we estimate a higher rate of $6.5\text{\,}\mathrm{Hz}$ due to its many additional findings which are likely to be wrong, confer the number of false positives in Table 1. We stress that we cannot exclude that a mixture of two (or more) exponential distributions is underlying which could be for instance caused by multiple underlying biochemical processes provoking flickering events, but a more detailed analysis would require more events or less artifacts.

6. Discussion and Outlook

We proposed a new method $\operatorname{JULES}$ for idealizing ion channel recordings, particularly suited for detecting and idealizing very short flickering events, while at the same time it idealizes events on intermediate and large scales. Detection and idealization of flickering is a difficult task due to pre-lowpass filtering incorporated in the amplifier. We confirmed its ability to detect and idealize such peaks in simulations and in an application on gramicidin A traces. One limitation is that two consecutive peaks have to be separated by enough observations (roughly three times the filter length in the given setting, see Figure 8) to enable a quick local deconvolution as in the paper. This can be refined at computational expenses. For the given data this does not seem to be a significant limitation as the estimated average distance is with $0.3\text{\,}\mathrm{s}$ much larger than the needed distance of $3.3\text{\,}\mathrm{ms}$.

Another improvement could be done in the detection step. The detection power can be increased by using the interval system of all intervals instead of the system with all intervals containing a dyadic number of observations. However, this increases the computational complexity from $\mathcal{O}(n\log(n))$ to $\mathcal{O}(n^{2})$. In the same spirit, the multiscale statistic $M$ in (7) could be based on local likelihood ratio test statistics instead on partial sums. However, the likelihood ratio test statistic is computationally more expensive as it involves inverting the corresponding covariance matrices on each interval considered in (7) and an empirical study (not displayed) has shown that reconstructions based on local likelihood ratio tests does only slightly improve detection power.

The current version of $\operatorname{JULES}$ requires an explicit knowledge of the filter which is usually known. The clampSeg package has an implementation for the most commonly used Bessel filter (with an arbitrary number of poles and cutoff frequency). If the filter is not known or if no implementation is available, the covariance could be estimated purely data driven by difference based methods Tecuapetla-Gómez and Munk (2016) from the recordings as well. Then, from the estimated covariance a discrete kernel for the filter could be computed and used for the deconvolution instead of the time-continuous kernel.

While we assume a non-stationary signal, we require a stationary noise. We found in Section 4.5 that $\operatorname{JULES}$ is rather robust against heterogeneous noise as long as all events are very short. In general, $\operatorname{JULES}$ assumes a global standard deviation and will be influenced by changes of the standard deviation during the recording, as e.g. caused by open channel noise. We are currently working on a follow up paper which will be able to deal with heterogeneous noise to overcome this limitation. We expect such a method also to be robust against changes of the correlation structure.

Besides these possible improvements, we will continue with using $\operatorname{JULES}$ to analyze further recordings with high biological relevance, for instance whether certain mutations influence the transport of antibiotics in a bacterial porin and to analyze the gating dynamics of the translocase channel Tim23.

Acknowledgment

Financial support of DFG (CRC803, projects Z02 and A01) and of the Volkswagen Foundation (FBMS) is gratefully acknowledged. O.M.S. thanks the GGNB PBCS for a Ph.D. fellowship. We thank Timo Aspelmeier, Annika Bartsch, Niels Denkert, Manuel Diehn, Thomas Hotz, Housen Li, Michael Meinecke, Ingo P. Mey, Erwin Neher, Dominic Schuhmacher, Ivo Siekmann and Frank Werner for fruitful discussions. We also thank two anonymous referees for their helpful comments and suggestions which helped us to improve this paper.