Heterogeneous Idealization of Ion Channel Recordings – Open Channel Noise

We propose a new model-free segmentation method for idealizing ion channel recordings. This method is designed to deal with heterogeneity of measurement errors. This in particular applies to open channel noise which, in general, is particularly difficult to cope with for model-free approaches. Our methodology is able to deal with lowpass filtered data which provides a further computational challenge. To this end we propose a multiresolution testing approach, combined with local deconvolution to resolve the lowpass filter. Simulations and statistical theory confirm that the proposed idealization recovers the underlying signal very accurately at presence of heterogeneous noise, even when events are shorter than the filter length. The method is compared to existing approaches in computer experiments and on real data. We find that it is the only one which allows to identify openings of the PorB porine at two different temporal scales. An implementation is available as an R package.


I. INTRODUCTION
The voltage patch clamp technique is a major tool to quantify the electrophysiological dynamics of ion channels in the cell membrane (Neher and Sakmann, 1976;Sakmann and Neher, 1995).It allows to record the conductance trace Florian Pein is with the Statistical Laboratory of the Department of Pure Mathematics and Mathematical Statistics (DPMMS) at the University of Cambridge, Wilberforce Road, Cambridge, CB3 0WB, United Kingdom.Annika Bartsch and Claudia Steinem are with the Institute of Organic and Biomolecular Chemistry, Georg-August University of Goettingen, Tammannstr.2, 37077 Göttingen, Germany.
Axel Munk is with the Institute for Mathematical Stochastics, Georg-August-University of Goettingen, Goldschmidtstr.7, 37077 Göttingen, Germany.Axel Munk is also with the Max Planck Institute for Biophysical Chemistry, Am Fassberg 11, 37077 Göttingen, Germany, and with the Felix Bernstein Institute for Mathematical Statistics in the Biosciences, Goldschmidtstr.7, 37077 Göttingen, Germany.Support of DFG (CRC803, projects Z02 and A01, Cluster of excellence 2067 MBExC Multiscale Bioimaging: From Molecular Machines to Networks of Excitable Cells), the Volkswagen Foundation (FBMS) and of EPSRC (EP/N031938/1 -Statscale programme) is gratefully acknowledged.We thank Timo Aspelmeier, Manuel Diehn, Benjamin Eltzner, Ingo P. Mey, Ole M. Schütte, Ivo Siekmann, Inder Tecuapetla-Gómez and Frank Werner for fruitful discussions.arXiv:2008.02658v1[stat.ME] 6 Aug 2020 (i.e., the recorded current trace divided by the applied voltage) of a single ion channel in time, which is for instance important in medical research for the development of new drugs (Kass, 2005;Overington et al., 2006).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;Pein et al., 2018).To obtain such an idealization an extensive amount of methodology is available nowadays, a selective review is given below.
Open channel noise: In this paper, we focus on recordings that are affected by open channel noise, i.e., have larger noise on segments with a larger conductance.The name open channel noise refers to the fact that a larger conductance results from an open pore.The additional noise when the channel is open can for instance be explained by current interruptions lasting approximately 1 microsecond (Sigworth, 1985(Sigworth, , 1986;;Sigworth et al., 1987;Heinemann and Sigworth, 1988, 1990, 1991).We analyze these recordings in a 'model-free' manner, i.e., without assuming a hidden Markov or related model, as a time series which is obtained by equidistant sampling from the convolution of a piecewise constant signal contaminated by white noise with the kernel of a lowpass filter.The white noise is scaled by an unknown piecewise constant standard deviation function to allow variance heterogeneity caused by open channel noise.We stress, that this modeling is rather general and also allows to deal with heterogeneity of measurement errors, not necessarily due to open channel noise.It will be explained in full detail in Section II.
Data: the outer membrane porin PorB: Figure 1 shows exemplarily a conductance trace of the outer membrane porin PorB from Neisseria meningitidis, a pathogenic bacterium in the human nose and throat region (Virji, 2009).
PorB is a trimeric porin and the second most abundant protein in the outer membrane of Neisseria meningitidis.
It is for instance relevant for the transport of antibiotics into the cell and hence of current interest to understand antibiotic resistance better.Recordings are obtained by the patch clamp technique using solvent-free bilayers.In Figure 1 it is clearly visible that the two conductance levels around 0.04 nS and 0.36 nS are affected from open channel noise as the variance of observations around 0.36 nS is much larger than of the ones around 0.04 nS.Such recordings were manually analyzed in (Bartsch et al., 2019) (Figure 1 and its explanation).They have been a major motivation for our work as they show distinct heterogeneous noise but also short event times, which we could not tackle satisfactorily by existing idealization methods, but also not by a manual analysis.In fact, in (Bartsch et al., 2019) only the conductance levels were investigated but not the full gating dynamics, since events on short time scales could not be idealized.
Methods for open channel noise: Idealization methodology can be divided into so called model-free methods (Colquhoun, 1987;VanDongen, 1996;Hotz et al., 2013;Gnanasambandam et al., 2017;Pein et al., 2018) which do not rely on a specific model for the gating dynamics, to methodology based on hidden Markov models (HMM) (Ball and Rice, 1992;Venkataramanan et al., 2000;Qin et al., 2000;de Gunst et al., 2001;Siekmann et al., 2011; Fig. 1: From seconds to microseconds: Ion channel recordings (grey points) displayed at a level of seconds (top panel), of milliseconds (middle panel) and of microseconds (bottom panels).Data points result from a representative conductance recording of PorB by the patch clamp technique using solvent-free bilayers at 20 mV.Diehn et al., 2019) and to current distribution fitting (Yellen, 1984;Heinemann and Sigworth, 1991;Schroeder, 2015;Hartel et al., 2019).The latter often assume a hidden Markov model as well but focus on parameter estimation directly.An idealization can be obtained by the Viterbi algorithmus (Viterbi, 1967) as soon as the parameters are determined.
Most HMM methods can deal with heterogeneous noise.Moreover, they allow to extrapolate information from larger (observable) to smaller (not observable) time scales and hence can provide a good idealization on small temporal scales.However, they rely heavily on the correctness of the assumed model assumptions.Up to few exceptions, see (Fuliński et al., 1998;Goychuk et al., 2005;Mercik and Weron, 2001;Shelley et al., 2010), a Markov model is a reasonable assumption for the underlying ion channel dynamics.However, artifacts in the data observed, for instance base line fluctuations, occur frequently in ion channel recordings and require elaborate data cleaning before a HMM can be fitted.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 electronics or from building vibrations (although damped).The PorB measurements display in Figure 1 show several artifacts of this type (see for instance the waviness of the observations or the conductance increase around 10.1 s, which severely hinders straightforward fitting by a HMM: We tried to fit this data set with in total four different hidden Markov model approaches.We achieved the best results when we assumed three states, but with the assumption that two states (with small conductivity) share the same expectation and variance.More details, also on parameter choices, are given in Section X in the supplement.The obtained idealization is shown in Figure 2. It fits long events well, but misses very short events, see for instance the lower left panel.Fitting such events well requires to take into account the filtering which is computationally very demanding for HMMs.We will discuss such an approach in Section X in the supplement as well.In summary, in addition, to the low robustness against artifacts, the choice of a specific Markov model, especially the determination of the number of states, can be a demanding task and often involves subjective choices by the analyst.
Contrary, model-free approaches can deal way more flexible with artifacts as they act rather locally on the time series without the assumption of a global model.Hence they are more robust than HMMs to model violations.
Therefore, they complement HMMs well, e.g., as a preprocessing step.For example, model-free methods can be used to select or verify a specific Markov model, in particular to determine the number of states and possible transitions, as they explore and potentially remove artifacts in a model-free manner.See also (Pein et al., 2018) for a more extensive discussion of further aspects of the different approaches.
To the best of our knowledge, all existing model-free approaches assume (implicitly or explicitly) homogeneous noise and hence produce unreliable results when open channel noise is present.Among the first methods which fall into this category is TRANSIT (VanDongen, 1996).An idealization by this approach, details of its limitations in our setting and further discussions can be found in Section X in the supplement.In Figure 3 we display JULES (Pein et al., 2018), a novel multiscale approach that also falls into this type of methods.It detects many small events on segments with a larger conductance and variance, but none on ones with a smaller conductance and variance.
These additional events are most likely artifacts caused by open channel noise.Indeed, in Section IV-D we found that the rates of a simulated hidden Markov model with parameters similar to them underlying the observations in Figure 1 could not be recovered when we used JULES to idealize the underlying signal.This effect is even more Fig.2: Idealization (red) of the data in Figure 1 assuming a HMM displayed on three different temporal scales.The model consists of three states, whereby two states are assumed to have the same expectation and variance.In the lower panels we also show the convolution of the idealization with the lowpass filter (blue).It fits most part of the data well, but misses short events, for instance the event displayed in lower left panel.
severe when the variance heterogeneity is larger.
Recently, there has been made some progress to adjust for heterogeneous noise in the context of model-free methods.However, they are not dedicated to idealize ion channel recordings, which means in particular that they do not incorporate lowpass filtering.Obviously, ignoring the filtering will deteriorate results.For illustration purposes Fig. 3: Idealization (red) of the data in Figure 1 by JULES displayed on three different temporal scales.In the lower panels we also show the convolution of the idealization with the lowpass filter (blue).It detects short events, but finds many small events (which are most likely false positives) at parts of high conductance and high variance (see for instance the idealization of the observations around 0.36 nS in the middle panel).These detections hinder the decovolution (see for instance the lower left panel) and make the idealization unreliable.
we display the heterogeneous multiscale approach HSMUCE (Pein et al., 2017) in Figure 4. We found that it provides reasonable results on larger temporal scales.However, due to filtering HSMUCE misses shorter events, see for instance the missed peaks around 10.2 s (lower left panel), 10.4 s or 10.7 s.Also for this type of methods we provide further examples and discussions in Section X in the supplement.Fig. 4: Idealization (red) of the data in Figure 1 by HSMUCE displayed on three different temporal scales.In the lower panels we also show the convolution of the idealization with the lowpass filter (blue).It detects events on larger temporal scales well (middle panel, lower right panel), but misses short events, see for instance at around 10.2 s (lower left panel), 10.4 s or 10.7 s.
The occurrence of short events is often called flickering.Missing them does not only potentially disturb the analysis of the general channel behavior, the analysis and hence the idealization of flickering events is also of its own interest in many applications, since flickering has often its own dynamics and can result from different molecular processes.Typical examples are conformational changes of the ion channel (Grosse et al., 2014) or the passage of larger molecules blocking the ions pathway through the channel (Raj Singh et al., 2012;Bartsch et al., 2019).Hence, one main goal of this paper will be to idealize and detect such events as well.
To this end, we introduce in Section II a statistical model which resembles all features (open channel noise, events on a large range of scales, filtering) of such complex data as in the previous example.In summary, we then ask for a model-free idealization method that adapts automatically to heterogeneous noise, hence is able to detect and idealize events on a large range of relevant scales accurately, but in particular also events shorter than the filter length.Furthermore, we aim to provide theoretical justification for the detected events (controlling false positives) and for a computationally efficient method to deal with large data sets.

HILDE:
To address these tasks, we propose in this paper a new method called Heterogeneous Idealization by Local testing and DEconvolution, HILDE.This method combines multiscale regression and deconvolution as it takes into account the convolution of the signal with the lowpass filter explicitly for detecting events that are short in time.Before we will explain our (quite involved) methodology further, we discuss firstly the general challenges: A major difficulty for any such method due to the presence of heterogeneous noise is to distinguish between small jumps in the signal and random fluctuations caused by the noise of unknown level.However, simultaneously estimating the signal and the noise level locally is notoriously difficult in general (Pein et al., 2017) and further hampered in our situation since the unknown signal and noise are both smoothed by the filter and hence deconvolution is required when shorter temporal scales are considered.We solve this by means of a multiresolution approach in combination with a local deconvolution to idealize events on all relevant temporal scales accurately.Whereas statistical multiresolution idealization that ignores the deconvolution can be computed efficiently by dynamic programming, see for instance (Hotz et al., 2013;Frick et al., 2014;Pein et al., 2017), combining multiresolution procedures with deconvolution is algorithmically difficult, since due to the coupling of all observations in the idealization, dynamic programming is not applicable without further ado.We will overcome this burden by focusing firstly on larger temporal scales and then improving the idealization on smaller temporal scales.More precisely, HILDE consists of the following three steps: a) detection of long events, b) detection of short events and c) parameter estimation by deconvolution.A summary about all three steps is given in Algorithm 1 (see Section III).a) Detection of long events: We will obtain an idealization by multiresolution regression that covers all important features on larger temporal scales (for the data set analyzed here large means events of length 6.5 ms at least, i.e., ≥ 65 sampling points).This step is discussed in Section III-A and technical details are given in Section VII in the supplement.b) Detection of short events: Our data set contains several short events that will be missed by the previous step, see for instance in Figure 1 at around 10.2 s (lower left panel), 10.5 s, 10.7 s and 10.8 s.To detect such events, we test locally whether additional events on smaller temporal scales have to be incorporated.This is impaired by the lowpass filter and the resulting convolution has to be taken into account explicitly.To this end, we assume that signal and noise left and right of the interval on which we test are given by the idealization from the previous multiresolution step.These tests are detailed in Section III-B, while technical details are postponed to Section VIII in the supplement.
Steps a) and b) determine the number of events and their rough locations.The final idealization in Figure 5 (see e.g. the lower left panel) confirms that step b) is indeed able to detect short events (up to 0.2 ms, corresponding to only two subsequent observations).c) Parameter estimation by deconvolution: Finally, the precise locations of the events and the conductance levels have to be obtained.This will done in an additional deconvolution step, as the recordings are filtered.To this end, we use the local deconvolution approach from Pein et al. (2018) with minor modifications.This step is discussed in Section III-C and technical details are explained in Section IX in the supplement.
Figure 5 shows the final idealization by HILDE of the observations in Figure 1.Despite distinct heterogeneous noise, the idealization covers all main features on all relevant scales, in particular also short events, while at the same time it does not include systematically additional artificial changes.The zooms into single peaks (lower panels) show that HILDE fits the observations well down to a scale of microseconds, which is also a confirmation of our approach, including the modeling.We stress that HILDE is not only robust against heterogeneous noise but has typically also a larger detection power for event detection than JULES (even when the noise is homogeneous), since it takes into account the convolution explicitly for detection.This is discussed in more detail in Section VI-B, where we also outline a version of HILDE that assumes homogeneous noise to improve detection power even further if the homogeneous noise assumption is justified.
While the first and the third step are mostly useful modifications of existing methodologies, we want to stress that this is not true for the second step.To the best of our knowledge, no other model-free ion channel idealization method is able to take the convolution explicitly into account when detecting events.As discussed before, this is however indispensable to detect short events when filtering and heterogeneous noise are present.
Implementation and run time: Each step of HILDE can be computed separately.Hence, HILDE can be applied and modified in modular fashion.This allows for instance to skip the second step if a data set contains only longer events, hence saving computation time.Another usage might be to modify the local tests in the second step, for instance to increase the detection power in a data set with small conductance changes but large difference in the noise levels, without modifying the first or third step.Such modifications are discussed in Section VI-A.
The first multiresolution regression step can be computed by a pruned dynamic program.The computation of the local tests in the second step is straightforward and the deconvolution in the third step can be computed by an iterative grid search.These steps are detailed in Section III-D and summarized in Algorithm 1.An implementation Fig. 5: Idealization (purple) of the data in Figure 1 by HILDE displayed on three different temporal scales.In lower panels we also show the convolution of the idealization with the lowpass filter (orange).HILDE idealizes events on all relevant temporal scales well. is available by the function hilde in the R package clampSeg accompanying this paper.The package is available on request and has been submitted parallel to CRAN (Pein et al., 2019b).
The worst case computational complexity is quadratic in the number of observations, but in most ion channel recordings conductance changes occur frequently which reduces the complexity to linear in the number of observations.For instance the 600 000 observations in Figure 1 can be idealized in a few minutes on a standard laptop.
A detailed discussion of the computational complexity is given in Section III-D.
Simulations: In Section IV we investigate the performance of HILDE in Monte-Carlo simulations which resemble the characteristics of the data in the application in Section V. Based on this we confirm that HILDE works very well for data sets like the one shown in Figure 1.In more detail, it can detect events which last 0.2 ms, corresponding to only two subsequent observations and being less than one fifth of the filter length long, with probability almost one.Furthermore, all parameters (conductance levels and the locations of the changes) are estimated very accurately, see Section IV-B for more details.Moreover, two subsequent events can be separated reliably as soon as the distance between them is larger than five times the filter length, see Section IV-C.In Section IV-D we simulate data from a hidden Markov model.Our method is not assuming a HMM, but such a model is still illustrative to simulate as it is a standard assumption for the analysis of ion channel recordings.We will also see in Section V-C that a Markov model is reasonable for the PorB recordings.We find in Section IV-D that HILDE recovers all parameters with high precision.Those parameters are chosen similar to those which we have estimated in Section V-C.Finally, we investigate robustness issues against f 2 -and 1/f -noise in Section IV-E.We omit most of the time a systematic comparison with other approaches, since, as discussed in the introduction before, to the best of our knowledge all existing approaches assume a more restrictive model which hinders a fair comparison.However, we include JULES, HSMUCE and a HMM based approach in the simulations in Section IV-D to illustrate the shortcomings (and benefits) of these approaches further.
Application to PorB recordings: Our analysis of single channel recordings of PorB in Section V confirms all major results from (Bartsch et al., 2019) about this data set.Moreover, a novel finding of our analysis is that the dwell times do not fit a single exponential distribution, but suggests that two different regimes for the dwell times are underlying: very short openings of 2.31 ms estimated average duration and longer openings of 51.62 ms estimated average duration.To best of our knowledge, fast and slow gating at the same time was not observed for PorB before, but for another porine OmpG (Grosse et al., 2014).We stress that all results obtained by HILDE could be confirmed by at least one other approach.However, none of the other methods were able to reproduce all results obtained by HILDE.
In summary, in this work we proposed with HILDE the first fully automatic model-free method for the analysis of ion channel recordings affected from open channel noise, i.e., to the best of our knowledge no other existing methodology is able to estimate a piecewise constant function in a model-free manner when filtering and heterogeneous noise are present at the same time.Simulations confirm that HILDE deals efficiently with heterogeneous noise and filtered data at the same time and idealizes events on various time scales efficiently.More precisely, to obtain a good idealization events have to be only at least two subsequent observations long but separated from each other by at least five times the filter length (at signal and noise ratio and filtering as in the analyzed data).This allowed us to obtain novel findings for the PorB channel, e.g., that it can have shorter and longer opening processes at the same time.

II. MODELING
We assume that the recordings result from equidistant sampling from the convolution of an unknown piecewise constant signal corrupted by Gaussian white noise with the (known) kernel of a lowpass filter.We stress however that our methodology can be extended to an unknown filter by using the methodology of (Tecuapetla-Gómez and Munk, 2017).To incorporate heterogeneity, the white noise is scaled by an unknown piecewise constant function to allow a larger variance on segments on which the conductance is larger.We only allow potential variance changes when the conductance changes, since variance changes also depend on gating events of the channel.More precisely, we model the conductivity and the standard deviation by piecewise constant signals f and σ, where t denotes physical time.The (unknown) conductance levels are denoted as c 0 , . . ., c K , the (unknown) standard deviations as s 0 , . . ., s K > 0, the (unknown) number of changes as K and the (unknown) locations as −∞ =: A and zero otherwise.The signals are extended to τ 0 = −∞ to define the convolution correctly but we will see at the end of this section that only a very short time period before recordings started, i.e., before t = 0, will be relevant.We assume c k = c k+1 to define the number of changes unambiguously, i.e., to obtain an identifiable model.But we allow s k = s k+1 , i.e., the standard deviation does not have to change between different events and in particular homogeneous noise is still part of the model (σ ≡ s 0 ).We stress that the class of signals in (II.1) is very flexible as potentially any arbitrary number of changes at arbitrary conductance levels and arbitrary standard deviations can be imposed, see Figure 5 for an example.
We assume further that the recorded data points Y 1 , . . ., Y n (the measured conductivity at time points t i = i/f s , i = 1, . . ., n, equidistantly sampled at rate f s ) result from convolving the signal f perturbed by Gaussian white noise η scaled by the standard deviation function σ with an analogue lowpass filter, with (truncated) kernel F m , and digitization at sampling rate f s = n/τ end , i.e., with * the convolution operator.Here, n denotes the total number of data points (typically several hundred thousands up to few millions).Like in (Hotz et al., 2013;Pein et al., 2018) we truncate (and rescale) the kernel of the lowpass filter and the covariance function at m/f s to simplify our model.This is implemented in the R function lowpassFilter (Pein et al., 2019b).As a working rule, we choose m such that the autocorrelation function of the untruncated analogue lowpass filter is below 10 −3 afterwards.For the later analyzed PorB traces, which are filtered by a 4-pole lowpass Bessel filter with 1 kHz cut-off frequency and sampled at 10 kHz, this choice leads to m = 11.Hence, the resulting errors 1 , . . ., n are Gaussian and centered, E[ i ] = 0, and have covariance Note that we have in (II.3) an unknown non-stationary covariance structure.However, the covariance can be decomposed into a known stationary autocorrelation given by the lowpass filter and an unknown non-stationary variance, which is modeled by a piecewise constant function that shares its change-points with the mean function.
An analytic expression of A m (t, l) is implemented in the R function lowpassFilter.Hence, (II.3) can be computed exactly and efficiently.
The major aim will be now to idealize (reconstruct) the unknown signal f taking into account the convolution, the heterogeneous noise given by (II.3) and the specific structure of f in (II.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 certain errors levels α = α 1 + α 2 , the maximal scale on which local tests are performed l max and two filter specific parameters have to be selected in advance, see Section III-E.

III. METHODOLOGY: HILDE
In this section we detail the three steps of our Heterogeneous Idealization by Local testing and DEconvolution (HILDE) approach.A summary of these steps is given in the Meta-algorithm 1.
Data Y 1 , . . ., Y n , error levels α = α 1 + α 2 , maximal scale l max , regularization parameter γ 2 , filter with kernel F m truncated at m/f s Detection of long events (> l max ): Multiresolution regression at error level α 1 Detection of short events: (≤ l max ): Local tests that take into account the convolution explixitly, at joint level α 2

Parameter estimation:
Local deconvolution, regularized by γ 2 Idealization f , i.e., all event times and conductance levels Algorithm 1: Steps of HILDE.

A. Detection of long events
To detect events on larger temporal scales, we use a modification of the Heterogeneous Simulataneous MUltiscale Change-point Estimator, HSMUCE from (Pein et al., 2017), which is a multiresolution procedure that is robust against heterogeneous noise.To avoid false positives due to the filter, we omit on each interval the first m observations and do not test on very short intervals.Since we truncated the filter, the signal and the convolution of the signal with the lowpass filter differ only at the beginning of each segment.More precisely, if the signal is constant on an interval [i/f s , j/f s ] with conductance level c ij and the first m observations Y i , . . ., Y i+m−1 are ignored, all other observations Y i+m , . . ., Y j have constant expectation equal to the conductance level c ij .Hence, we take into account only intervals longer than m/f s and ignore the first m observations of each interval.This leads to an estimator that detects change-points at presence of heterogeneous noise and filtering, i.e., when the heterogeneous ion channel model from Section II is assumed, while at the same time the probability to overestimate the number of events is controlled, i.e., a false positive is only added with probability at most equal to the tuning parameter α 1 , see Theorem VII.1 in Section VII in the supplement.A detailed definition of this estimator is given in Section VII in the supplement.Note that it does not take into account the convolution explicitly but still has good detection properties if events are long enough, but almost no detection power on small scales.Simulations (not displayed) show that for our data set events with of length at least 6.5 ms, corresponding to 65 sampling points, are detected reliably.
In the following two sections we will present a refinement of this idealization to detect and idealize events on smaller time scales, too, which proves to be relevant for our data example.Note that in this section and in the next section (a refinement will be provided in the local deconvolution step) we restrict all changes to the grid on which the observations are given, in other words, we assume that f s τ i are integers.

B. Detection of short events
To detect short events, we test on all intervals containing l = 1, . . ., l max (to be defined later) observations whether the previous idealization is the underlying signal or whether the inclusion of an additional event on the considered interval is significantly better.More precisely, let [τ L , τ R ] = [i/f s , j/f s ] be the interval on which we test.And assume for the moment that τ is the only change in [(i − m + 1)/f s , (j + m − 1)/f s ] with conductance levels c L before and c R afterwards.Note that this also includes the scenario of no change in by setting c L = c R .Then, we decide whether an additional event on [τ L , τ R ] is required by testing the hypothesis The form of these hypotheses allows us to construct a test statistic that takes into account the convolution explicitly.
Moreover, information provided by potential variance changes can be used as well.We provide details of the corresponding test in the paragraph 'Local testing' in Section VIII in the supplement.All choices there are motivated by a trade-off between a good detection power for events in the measurements in Section V, see Figure 1, and a reasonable computational complexity.
If a hypothesis is rejected, we replace the single change-point at τ by a short peak.Temporary locations will be placed at τ L = i/f s and τ R = j/sr, but exact locations and the conductance level c will be obtained in the upcoming deconvolution step.However, note that usually one event in the data causes rejections of multiple tests.
Therefore, we only consider the event with the largest test statistic among all rejections on intervals that intersect or adjoin each other.More details are provided in the paragraph 'Multiple dependent rejections' in Section VIII in the supplement.

C. Parameter estimation by local deconvolution
The final idealization is obtained by local deconvolution as described in Section 3.2 of (Pein et al., 2018) with two adjustments, which will be discussed in Section IX in the supplement.This means in particular that we still use the likelihood function of observations with homogeneous noise, although heterogeneous noise is assumed.
Simulations show, see Section IV, that this works reasonably well for the recordings we analyze in Section V.
Alternatives for recordings with more pronounced noise heterogeneity are discussed in Section VI-A.

D. Computation and run time
The multiresolution regression step in Section III-A can be computed by a pruned dynamic program as described in Section A.1 in the supplement of (Pein et al., 2017).For related ideas, see also (Killick et al., 2012;Frick et al., 2014;Li et al., 2016;Maidstone et al., 2017) and the references given there.The implementation of the local tests in Section III-B is straightforward.The local deconvolution in Section III-C can be computed by an iterative grid search as described in Section 3.2 of (Pein et al., 2018).An implementation of HILDE is available by the R function hilde in the package clampSeg.The package is available on request and has been submitted parallel to CRAN (Pein et al., 2019b).All run time critical parts are written in C++ and are interfaced by the R code.
The worse case computation complexity of the dynamic program is quadratic in the number of observations n.
However, in most ion channel recordings conductance changes occur frequently which reduces the complexity to be linear O(n), see Section A.3 in the supplement of Pein et al. (2017).The local tests in Section III-B are of complexity O(l 2 max n), since for each of the l max scales 1, . . ., l max roughly n tests have to be performed and the complexity to compute a single test is at most of order O(l max ).The computation time of the local deconvolution is dominated by the iterative grid search to deconvolve a single event.The deconvolution of a single event is constant in the number of observations, since the number of involved observations and the grid sizes do not increase.Moreover, the number of involved observations is small and the covariance matrix is a band matrix, with band size equal to m, which allows fast computation.Hence, the complexity of the deconvolution increases linearly in the number of events which increases for ion channel recordings typically linearly in the number of observations.In summary, for a typical channel trace the complexity to compute HILDE increases only linearly in the number of observations.This is confirmed by a run time of less than five minutes for idealizing the 600 000 observations in Figure 1 on a Dell Latitude E6530 with Intel(R) Core(TM) i5-3340M CPU 2.70GHz processor.Similar run times are obtained for the traces generated in Section IV-D.Thus, the theoretical considerations as well as the empirical run times confirm that HILDE can be computed efficiently, which is important since large data sets have to be analyzed.

E. Parameter choices
HILDE can be tuned by the parameters α 1 , α 2 , l max and γ 2 , see Algorithm 1 and the referenced sections for a definition.The probability to overestimate the number of conductance changes is approximately controlled by the sum of the error levels α = α 1 + α 2 .Hence, if such an overestimation control is desired, α should be chosen small.
As a default choice we suggest α = 0.05.Increasing α yields to a larger detection power (at the price of including more false positives).Hence, one may 'screen' for different α if important events are difficult to detect.The levels α 1 and α 2 allocates the power between the multiresolution test for detecting events on large scales (> l max ) and the local tests to detect events on small scales (≤ l max ).We have chosen α 2 = 0.04 and α 1 = 0.01 in our data analysis, since our focus was on detecting short events primarily, while events on larger scales were easier to detect.
More weight can be put on α 1 if either short events are of less interest or if long events are difficult to detect as well, e.g.since they have a smaller jump size than the short events.The latter is often called subgating and was for instance studied in (Hotz et al., 2013).The tuning parameter l max , the largest scale on which local tests are performed to find short events, should be chosen such that all events on larger scales are detected by the previous multiresolution test.This can for instance be determined by Monte-Carlo simulations.In our setting, we choose l max = 65, since simulations (not displayed) showed that the multiresolution step in Section III-A is able to detect events which contain more than 65 observations with probability almost one.The correlation matrix is regularized with parameter γ 2 = 1, further details can be found in Section III B in (Pein et al., 2018).And, as mentioned before, we truncate the kernel and autocorrelation function of the filter at m = 11 as the autocorrelation function is below 10 −3 afterwards.All of these choices are the default parameters of the function hilde and are used in the simulations in Section IV and in the real data application in Section V.

IV. SIMULATIONS
In this section we examine the performance of HILDE in Monte-Carlo simulations.Since to our best knowledge no other model-free method is known that takes into account heterogeneous noise and filtering explicitly, it is difficult to compare HILDE with other methods.Most similar in spirit are JULES (Pein et al., 2018), HSMUCE and an HMM based approach (Diehn, 2017).These have been included in a simulation in Section IV-D for purpose of comparison.The simulation study consists of four parts.First of all, we investigate the detection and idealization of isolated peaks.Secondly, we identify the minimal distance at which HILDE is able to separate two consecutive peaks.Thirdly, although HILDE does not rely on a hidden Markov model assumption, we examine its ability to recover the parameters of a Markov model, since a hidden Markov model is a common assumption for ion channel recordings.Finally, we investigate its robustness against violations of the model in Section II, in particular against additional f 2 and 1/f noise.

A. Data generation
We generate all signals and observations accordingly to the heterogeneous ion channel model we described in Section II and such that they are in line with the measured data we analyze in Section V.This means in particular that amplitudes, dwell times and noise levels of the generated observations are chosen such that they are similar to those of the analyzed datasets.We also simulate a 4-pole Bessel filter with 1 kHz cut-off frequency and sample the observation at 10 kHz.
The expectation of the observations, given by the convolution of the signal with the truncated kernel F m of the lowpass Bessel filter, can be computed explicitly.For the errors we oversample by a factor of 100, i.e., we generate 100 times as many independent Gaussian observations, discretize the filter accordingly, compute a discrete convolution and rescale the observations such that they have the desired standard deviation.

B. Isolated peak
In this simulation with 4 000 observations we examine the detection and idealization of a single isolated peak.
More precisely, in accordance with the model in Section II and with the estimated values in Section V for the observations in Figure 1, we choose conductance levels c 0 = c 2 = 0, c 1 = 0.32, variances s 2 0 = s 2 2 = 6.1 • 10 −5 and varying variance s 2 1 ∈ {2 • 10 −4 , 5 • 10 −4 , 10 −3 , 2 • 10 −3 , 5 • 10 −3 } to examine the influence of different noise levels.Note that s 2 1 = 10 −3 is roughly the noise level in the measurements in Section V.Moreover, we simulate changes at τ 1 = 2000/f s and τ 2 = (2000 + )/f s , c.f. (II.1), and are interested in how well HILDE detects the peak and idealizes the locations τ 1 and τ 2 and the level l 1 as a function of , the length (relative to the sampling rate f s ) of the peak.For = 5, Figure 6 shows an example of the simulated data as well as the idealizations by HILDE and their convolutions with the Bessel filter in a neighborhood of the peak.Tables I-III summarize our results based on 10 000 repetitions for = 2, 3, 5.
To this end, we count how often 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 |τ j −τ 1 | < m/f s and |τ j+1 −τ 2 | < m/f s as a peak is shifted at most m/f s by the filter.If only one change but not a peak is within these boundaries 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 locations and the level we only consider cases where the peak is detected and report the mean square error, the bias and the standard deviation.In most scenarios, HILDE has a good detection power and detects almost no false positives, see Table I.Only for a five times larger variance than in the real data and when l ≤ 3 few events are missed.In Tables II and III we found that idealization of the locations τ 1 and τ 2 and the conductance value c 1 works well for variances similar to the real data, but has some issues when the variance of the peak is larger, in particular in the scenario of a five times larger variance.For such observations it might be desirable to take into account the heterogeneous noise in the deconvolution step, see Section VI-A for more details.For smaller variances the results for estimating the locations are better when the peak is longer, but for larger variances results are even worse when the peak is longer.noise and the length of the peak.More precisely, it has changes at τ 1 = 0.2 and τ 2 = τ 1 + /fs, = 2, 3, 5, conductance levels c 0 = c 2 = 0, c 1 = 0.32, variances s 2 0 = s 2 2 = 6.1 • 10 −5 and varying variance s 2 1 ∈ {2 • 10 −4 , 5 • 10 −4 , 10 −3 , 2 • 10 −3 , 5 • 10 −3 }. Results are based on 10 000 pseudo samples.An example, s 2 1 = 10 −3 and = 5, is given in Figure 6.An explanation might be two effects with opposite influences.The conductance change provide more information when the peak is longer, but then also the overall variance of the observations is larger which reduces estimation accuracy.Estimation of the level c 1 is always more accurate when the peak is longer.It seems that here the first effect dominates.

Setting
All in all, these simulations confirm that HILDE performs very well for observations comparable to them in Section V.

C. Separation of two consecutive peaks
To examine how well HILDE separates two consecutive peaks we perform the same simulations as in Section 4.3 in (Pein et al., 2018), since results are identical for homogeneous and heterogeneous noise as separation depends Results are based on 10 000 pseudo samples and are given as multiples of the sampling rate fs = 10 4 .An example, s 2 1 = 10 −3 and = 5, is given in Figure 6.

Setting
Length ( )  channel noise and the length of the peak.More precisely, it has changes at τ 1 = 0.2 and τ 2 = τ 1 + /fs, = 2, 3, 5, conductance levels c 0 = c 2 = 0, c 1 = 0.32, variances s 2 0 = s 2 2 = 6.1 • 10 −5 and varying variance s 2 1 ∈ {2 • 10 −4 , 5 • 10 −4 , 10 −3 , 2 • 10 −3 , 5 • 10 −3 }. Results are based on 10 000 pseudo samples.An example, s 2 1 = 10 −3 and = 5, is given in Figure 6. on the method and distance between the peaks but not on the noise level.More precisely, we consider a signal f with changes at τ 1 = 2 000/f s , τ 2 = τ 1 + 5/f s , τ 3 = τ 2 + d and τ 4 = τ 3 + 5/f s , , with τ 0 = 0 and τ end = 4 000/f s and levels l 0 = l 2 = l 4 = 40 pS and l 1 = l 3 = 20 pS.Hence d is the distance between the two peaks.We distinguish between perfect separation, i.e., the detection step of HILDE identifies the two peaks (4 changes) and the local deconvolution yields idealizations for the four levels (illustrated in Figure 7c).Secondly, separation fails in the detection step, i.e., the multiresolution reconstruction recognizes only 2 changes and identifies one peak whose level can be further deconvolved (illustrated in Figure 7a).Finally, separation fails in the deconvolution step, i.e., HILDE identifies two peaks but the distance is so small that the deconvolution step cannot separate them, in other words, no long segment is in between (illustrated in Figure 7b).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, . . ., 70}.We found that the two peaks are detected if d > 6, but separation in the detection step and hence an appropriate idealization requires d ≥ 62.Hence, in Section V events have to be separated by more than 6.2 ms to be idealized appropriately.In comparison, we found that events are on average separated by 1/(2.67 + 4.50) s ≈ 0.14 s which shows that this limitation is not an issue for the analyzed PorB recordings.

D. Hidden Markov model
In this section we simulate data from a three state hidden Markov model.Since hidden Markov models are often assumed for ion channel recordings, it is instructive to investigate the methods in such a scenario.We simulate observations that resemble the PorB data we analyze in Section V.More precisely, we have expectations 0 nS, 0 nS and 0.32 nS as well as standard deviations 0.0078 nS, 0.0078 nS and 0.0316 nS, i.e. the variances are 6.1 • 10 −5 (nS) 2 and 10 −3 (nS) 2 .The dwell times in the first, second and third state are exponentially distributed with rates 20 Hz, 400 Hz and 7 Hz, respectively.The process always jumps from the first or second state to the third state, i.e., no transitions between the first and second state are allowed.And it jumps from the third state with probability 2/3 to the first state and with probability 1/3 to the second state.We generate five time series with 600 000 observations, each.Each trace looks similar to the observations in Figure 1 and hence we refrain from showing an example.
We analyze these data sets with HILDE and for purpose of comparison with JULES (Pein et al., 2018), HSMUCE Fig. 8: Results for HILDE assuming heterogeneous noise in idealizing two consecutive peaks separated by distance d.Its frequencies for no separation in the detection step (green), for successful detection, but no separation in the deconvolution step (red) and for successful detection and deconvolution (blue).Results are based on 10 000 simulations for each value of d. (Pein et al., 2017) and an HMM based approach which assumes the true three state model, i.e. three states, whereby two have the same expectations and variances and no transitions are allowed between them.We used µ = (0, 0, 0.32) T , s = (0.0102, 0.0102, 0.0321) T , P = .
We will discuss the estimated transition matrix later in comparison with the other approaches and when we also discuss the results using the Viterbi algorithm.The estimated expectations and standard deviations are accurate.For the other approaches we show in Figure 9 histograms of the estimated amplitudes of all events with an amplitude between 0.2 nS and 0.5 nS.We found in Figure 9 that all approaches estimate the amplitude accurately.The estimated amplitudes of HSMUCE are skewed, but the final estimation is still decent.
We continue with an analysis of the dwell times.To this end, we consider from now on all events with estimated conductance level between −0.05 nS and 0.05 nS as a closed event and between 0.15 nS and 2 nS as an open event, while all other events are considered as artifacts and are ignored.Figure 10 shows histograms of the dwell times in the closed state for various approaches.We rescaled all lines such that the area under them are standardized to one to make them comparable to the histograms.
We see that with the exception of HSMUCE (it misses short events) none on the histograms look exponentially distributed, since we have a mixture of short and long events.Hence, in Figures 11 and 12 we will analyze short and long events separately.To this end, we say an event is short if its dwell time is between 0.1 ms and 5 ms and long if its dwell time is between 20 ms and 200 ms.To estimate the rates, we apply a missed event correction like in (Pein et al., 2018).Fig. 11: Histograms of the dwell times in the closed state with a length between 0.1 ms and 5 ms to analyze short events, together with the true exponential distribution (black line) and exponential fits (red) that are corrected for missed events.We rescaled all lines such that the area under them are standardized to one to make them comparable to the histograms.This explains why the true distribution does not look always the same.The blue line in the HMM plot indicates an exponential fit when the same missed event correction is used that is applied to the results obtained by HILDE.
We found from Figures 11 and 12 that HILDE recovers in both cases the exponential distribution very well and estimates both rates of 400 Hz and 20 Hz with 332.59 Hz and 19.2629 Hz accurately.In comparison, JULES is not able to deconvolve all events due to the detection of additional spurious events, compare Figure 3.The rate for events, together with the true exponential distribution (black line) and exponential fits (red) that are corrected for missed events.We rescaled all lines such that the area under them are standardized to one to make them comparable to the histograms.This explains why the true distribution does not look always the same.
the short events is with 520.16 Hz still accurately, but the rate for the long events is with 9.0682 Hz significantly underestimated.Notably the dwell times are still (almost) exponentially distributed.HSMUCE misses short events, in total it has detected only 18 short events.Hence, a rate for the short events cannot be estimated.The rate for the long events is with 6.0182 Hz underestimated as well.The hidden Markov approach estimated with 347.44 Hz and 19.3949 Hz both rates accurately.However, since this approach misses very short events, for the rate for the short events we had to apply a stricter missed event correction that takes into account only events with a length of at least 0.75 ms.Hence, at least in the used form the hidden Markov approach is less favorable to analyze short events (since its corrected estimate is based on less event and hence will have a larger variance).This is remarkably, since the idealization on very short temporal scales is considered to be a strength of hidden Markov approaches.
Finally, the estimated exit probabilities by the Baum-Welch algorithmus of 0.002671 and 0.000644 corresponds to estimated rates of 26.71 Hz and 6.44 Hz which is much worse than the rates estimated using the idealizations obtained by the Viterbi algorithm.
We are now analyzing how often closing events occur.To this end, we analyze the dwell times in the open state or in other words the distance between two closing events.Moreover, we analyze the proportions of short and long events.Therefore, we divide the number of detected events by the estimated probability that such an event is detected assuming an exponential distribution for the dwell times.
We found from Figure 13 that HILDE, HSMUCE and the hidden Markov approach recover the exponential distribution very well and estimate the rate of 7 Hz with 6.4007 Hz, 6.4919 Hz and 6.2389 Hz accurately.Only JULES underestimates the rate because of previously explained reasons with 4.3362 Hz a bit.HILDE, JULES and HMM estimated with 0.3608, 0.2982 and 0.4082, respectively, the proportion of short events decently, recall that the truth is 1/3.This number could not be determined using HSMUCE, since it misses almost all short events.
Once again, the Baum-Welch provides with 27.06 Hz and 0.0329 much worse results.and exponential fits (red) that are corrected for missed events.We rescaled all lines such that the area under them are standardized to one to make them comparable to the histograms.This explains why the true distribution does not look always the same.
All in all, we found that HILDE was indeed able to recover all parameters very well.All other model-free idealization methods had at least one massive problem.The hidden Markov approach might be usable, but requires a more restrictive missed event correction and is also more complicated to apply.One should also keep in mind that we used the true parametric model class as prior knowledge.

E. Robustness
The model we proposed in Section II is a good assumption for ion channel recording at presence of open channel noise.However, in some patch clamp recordings additional high frequency f 2 (violet) and long tailed 1/f (pink) noise components have been observed, for a more detailed discussion see (Neher and Sakmann, 1976;Venkataramanan et al., 1998b;Levis and Rae, 1993) and the references therein.Thereto, in this section we examine how robust HILDE is against such noise components.To this end, we revisit the simulation setting from Section IV-B with s 1 = 10 −3 only.
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.
We assume that the pink noise is globally present.More precisely, we reduce the previously present noise by a factor of 1/2 and add pink noise which is scaled such that its standard deviation is equal to 1/2 √ 6.1 • 10 −5 (half of the standard deviation in the background in Section IV-B).For the high frequency violet noise we consider the setting that the new noise component is state-dependent as well.In other words, we generated errors from such a moving average process and convolved them with the kernel of the lowpass filter instead of assuming white noise errors.
We found in Table IV that HILDE is very robust against the additional f 2 but effected by 1/f noise.At presence of the latter noise, the standard deviation estimation on the long segments is wrong which causes the detection of false positives in roughly a quarter of the cases.Note, that false positives are caused by the underestimated standard deviation but also by the long range dependency itself.However, the false positives have a small amplitude and therefore do not influence the analysis severely or can be removed by postfiltering.Parameter estimation (not displayed) is slightly effected by 1/f noise (estimation of the change-point locations is slightly worse, but estimation of the size of the change is even improved), but not affected by presence of f 2 noise.

A. Measurements
We analyze single channel recordings of PorB from Neisseria meningitidis (recall the last paragraph in the introduction).In the following we analyze six traces, each of them is one minute long and consists of 600 000 observations.An example is shown in Figure 1, which shows distinct heterogeneous noise.
Measurements were performed on solvent-free planar bilayers using the Port-a-Patch (Nanion Technologies).Giant

B. Idealization
Idealizations are obtained by HILDE with parameter choices as in Section III-E.Moreover, an illustrative comparison with other approaches was discussed in the introduction (recall Figures 19-4).In Figure 1 we see that the channel switches frequently between two conductance levels, roughly between 0.04 nS and 0.36 nS, the variance is roughly 6.1 • 10 −5 (nS) 2 in the closed state and 10 −3 (nS) 2 in the open state.Moreover, several artifacts seem to be present, see for instance the fluctuating conductance in the open state in the first ten seconds.We stress that such artifacts heavily disturb any idealization that assumes a HMM, confer Figure 19.Contrarily, the modelfree idealization by HILDE (Figure 5) recovers all visible features on small as well as on large temporal scales accurately.In particular, the zooms into single peaks (Figure 5, lower panels) shows that HILDE fits the observations well which is also a confirmation of our model.Since PorB forms three pores, four different conductance levels are possible.However, in this measurement we see only two different conductance levels.Such a cooperative opening and closing was observed before, see for instance (Song et al., 1998).

C. Analysis of flickering dynamics
We now use the obtained idealizations to analyze the gating dynamics in a similar fashion as the simulated data in Section IV-D.We will focus in this section on HILDE, but we will compare it in Section XI in the supplement with analyses based on JULES, HSMUCE and HMM.We say a channel opens (a gating event from the lower conductance level to the higher conductance level) if the idealized level is between 0.25 nS and 2 nS and the previous level is between 0 nS and 0.1 nS.To study the amplitude, we consider the conductance difference of all such events.
Figure 14 shows a histogram of the so obtained amplitudes between 0.2 nS and 0.5 nS.All other events are either closing events or are considered as artifacts.Such artifacts can for instance be base line fluctuations as discussed in the introduction.We stress that an analysis of the closing events leads to very similar results.
The histogram in Figure 14 shows only one mode.Hence, all events have the same amplitude up to measurements and idealization errors.This means especially that also the flickering events are full-sized.An amplitude of 0.3194 nS is estimated by the half sample mode (Robertson and Cryer, 1974), computed in R by using the modeest package.
Note that other mode estimators or Gaussian mean estimation lead to similar results.This amplitude coincides with the one obtained by a manual analysis using the pClamp 10.2 software package (Axon Instruments), see (Bartsch et al., 2019).
We now analyze the dwell time in the open state and how frequently the channel opens.We take into account events with an amplitude between 0.2 nS and 0.5 nS and with a dwell time between 0.1 ms and 200 ms, since shorter events cannot be detected reliably and longer events are rare and often interrupted by artifacts.Histograms of the dwell time in the open state are shown in Figure 15 together with an exponential fit using a missed event correction like in (Pein et al., 2018).(a) All events between 0.1 ms and 200 ms.
(b) Short events between 0.1 ms and 5 ms.
(c) Long events between 20 ms and 200 ms.
Fig. 15: Histograms of the dwell times in the open state of all opening events with amplitude between 0.2 nS and 0.5 nS together with exponential fits using a missed event correction (red line).
Interestingly, the dwell times do not fit a single exponential distribution, but when we split the events in short (shorter than 5 ms) and long (longer than 20 ms) ones, both fit exponential distributions very well, with an estimated average duration of 51.62 ms and 2.31 ms, respectively.Note, that these estimations are approximations, since an exponential distribution with a large / small rate generates with a small probability a long / short event, but since the average dwell times are very different this error is negligible.To best of our knowledge, fast and slow gating at the same time was not observed for PorB before.However, Grosse et al. (2014) showed that the loop within the pore structure of OmpG leads to fast flickering (fast time constant).If the loop is removed, there is still gating observed but less frequent (slower time constant).Even though this is not the same protein, in PorB we have a loop L3 which is also localized in the pore and forms an α-helix in its center, which constricts the pore to its narrowest point.Hence, our findings support that similar dynamics might occur for PorB as well.
We are now analyzing the distance between two opening events.This is shown in Figure 16.Moreover, we analyze how many of the openings are short or long.Once again we apply a correction for missed events.
Fig. 16: Histograms of the distances between two opening events with amplitude between 0.2 nS and 0.5 nS together with exponential fits using a missed event correction (red line).
The distance between two events seem to be exponentially distributed and the estimated rate is 5.75 Hz.We found that 39.08% of all opening events were short events.Moreover, we found in Section XI in the supplement that all results obtained by HILDE could be confirmed by at least one other approach, but none of the other methods was able to reproduce all results obtained by HILDE.

VI. DISCUSSION AND OUTLOOK
In this paper we proposed a new model-free idealization method for ion channel recordings, called HILDE.
In comparison to existing approaches, HILDE provides still reasonable idealizations under heterogeneous noise, for instance caused by open channel noise.Moreover, it detects and idealizes flickering events reliable, is fullyautomatic and can be computed efficiently.It offers great flexibility in adapting to the needs of a specific data analysis by modifying the error probabilities α = α 1 + α 2 and the scale l max that distinguishes short and long events.Its precise idealization is confirmed by simulations and a real data application to PorB recordings.We found that these recordings contain opening events of significantly different length.
We stress that HILDE is modular, i.e., single components like the choice of the test statistics and functionals to optimize can be changed without further modifications.This can be used to adapt HILDE to specific challenges in the measurements.We will discuss several such possibilities in the following.Some of them are implemented in the clampSeg package and just require to choose different parameters, for others few lines of code have to be modified.

A. Alternative approaches
A different underlying interval set can be used for the multiresolution test.The set of all intervals of dyadic length provides in general a good compromise between detection power and computation time.But, if a larger detection power is required, the set of all intervals can be used at the price of a larger computational complexity.The other way around, if faster computation is demanded, a smaller interval set, for instance the dyadic partition like in (Pein et al., 2017), can be used.This might be particularly beneficial in situations in which the multiresolution test detects almost no events which results in a large computation time.Interesting alternatives are also the approaches in (Chan and Walther, 2013;Kovács et al., 2020) which require only a slightly larger computational effort than the use of all intervals of dyadic length but detects change-points in a certain sense statistically optimally.A different way to increase the detection power is to use likelihood ratio tests, again at computational expenses.We found in simulations (not displayed) that the likelihood ratio test statistic is slightly more powerful on small scales, but much slower to compute.However, a slightly worse detection power on small scales should not be a big concern, since a refinement by local tests will be done in the next step.Also for detecting events on small scales by local tests, see Section III-B, different statistics can be used to increase the detection power.For instance the likelihood ratio test or maximum likelihood estimators for the parameters (c, s 2 ) can be considered.However, they are computationally very demanding, since the likelihood function involves the inverse and the determinant of the covariance matrix given by (VIII.5).
Finally, our deconvolution approach assumes still homogeneous noise which we found in simulations works still well at presence of open channel noise, see Tables II and III.Taking into account the heterogeneous noise might be beneficial, in particular if the noise level differences are large, but difficult, maybe even impossible, since avoiding an ill-conditioned matrix by regularization and keeping the variance levels might be impossible to achieve at the same time.

B. Homogeneous noise
We designed HILDE particularly to deal with heterogeneous noise.However, taking into account the convolution explicitly when detecting changes is also beneficial if the noise is homogeneous, i.e., a constant variance is assumed.
In this situation, its detection power can be further improved by small modifications that utilize the assumption of a constant variance, they are explained in Section XII in the supplement.We found that HILDE has a better detection power than JULES (Pein et al., 2018), but at the price of worse separation properties and a larger computation time.More precisely, in the simulations in Section IV-B in (Pein et al., 2018) we found that JULES is able to detect an isolated peak of length 3.1/f s with probability almost one.In comparison, HILDE requires only 2.3/f s if homogeneous noise is assumed and 2.8/f s if heterogeneous noise is assumed (see Section 3.9.3. in (Pein, 2017)).
Remarkably, the detection power of HILDE is even larger than the one of JULES if HILDE does not use the assumption of homogeneous noise which illustrates how much detection power is lost by not taking into account the convolution.

C. Idealizing the variance
Our focus was on idealizing the conductance while the unknown variance was considered as a nuisance parameter.
However, as a byproduct HILDE can easily be extended to an idealization of the variance which offers for instance a residual analysis of the noise to validate a given model.To this end, we use HILDE to estimate the change-point location of the conductance and assume that these are the change-points of the variance as well.Note that the model of Section II allows the variance to stay constant at such a location but precludes further variance changes.
With the definitions from before (see Section III), if a segment is long, the square of the estimator in (VIII.1)can be used.Afterwards, the variance on short segments can be estimated by the estimator in (VIII.8).The resulting function will be an idealization of the variance.
SUPPLEMENT TO HETEROGENEOUS IDEALIZATION OF ION CHANNEL RECORDINGS -OPEN CHANNEL NOISE VII.LARGE SCALES Following the ideas of HSMUCE (Pein et al., 2017) we propose the idealization f by . ., i}.In other words, f is the maximum likelihood estimator restricted to all solutions of the optimization problem Thereby, F k is the set of all candidate signals in (II.1) with K = k changes, k ≥ 0. K is the estimated number of changes defined by the minimal number k for which there exits an f ∈ F k such that M ( f ) ≤ 0. Finally, M ( f ) denotes the multiresolution test statistic This tests simultaneously on all scales (resolution levels) whether f fits the data well.If it does not, the local test statistic T i,j will be larger than the scale dependent critical value q |j−i+1| , exact definitions are given below, and f will not be considered as a potential idealization.Such a method has many favorable properties, for more details see (Frick et al., 2014;Pein et al., 2017Pein et al., , 2018)).Most importantly, the number of false positives is controlled if the scale dependent critical values are defined appropriately (one possible choice is outlined below), see Theorem VII.1.
We use as test statistic the statistic in (1.5) in (Pein et al., 2017) without taking into account the first m observations (JSMURF principle), i.e., with fij the conductance level of f on the interval [i/n, j/n], Y i+m,j = (j − i + 1 − m) −1 j l=i+m Y l and ŝi+m,j = (j − i − m) −1 j l=i+m (Y l − Y i+m,j ) 2 .This statistic estimates the variance locally and is large if the mean of the observations differ significantly from the conductance level fij .The scale dependent critical values q m+1 , . . ., q n are obtained in a universal manner by Monte Carlo simulations as described in Section 2 of (Pein et al., 2017) such that (VII.3) is a level α 1 -test and different scales are balanced by weights.Here, we use the default choice of uniform weights.The error level α 1 ∈ (0, 1) has to be fixed in advance (by the experimenter) to control the number of false positives of the estimate f as stated in the following Theorem VII.1.See Section III-E for a discussion how to choose α 1 .
Theorem VII.1.Assume the heterogeneous ion channel model from Section II and let f be defined as in (VII.1) with error level α 1 .Then, for any f and σ in (II.1) the probability that f overestimates the number of changes is bounded by α 1 , i.e.
To keep the method computationally feasible, we only evaluate the maximum in (VII.3) over the system of all intervals that contain a dyadic number of observations, i.e., the maximum in (VII.3) is only taken over all 1 ≤ i ≤ j ≤ n such that j − i + 1 = 2 l for some l ∈ N 0 .This reduces the complexity of the system from O(n 2 ) to O(n log(n)) intervals.This speeds up the pruned dynamic program and the simulations in Section IV show that this interval set is large enough to allow a good performance.More details and a discussion of the run time are given in Section III-D.Note that (Hotz et al., 2013) performed tests on a slightly different interval set.
They required that j − i + 1 − m is a dyadic number.Although depending on the true signal, their choice leads in general to a better detection power on scales slightly larger than the filter length m/f s , but its computation lasts much longer.And missed short events can be detected by the upcoming local tests.The very long computation time was one major criticism in (Gnanasambandam et al., 2017).To be fair, both implementations differ in other points, too, and a fast implementation of their interval set might be possible as well, but our approach was easier to integrate in the dynamic programming framework of the stepR package (Pein et al., 2019a).Finally, we remark that the restricted maximum likelihood estimator ignores the convolution and hence the locations of the detected changes are typically a little bit shifted to the right which will be corrected in the upcoming deconvolution step in Section III-C.Different to this, (Hotz et al., 2013) suggested to move all locations by a constant factor t 0 , only depending on the filter, to the left, but our deconvolution step will (usually) be more precise.

VIII. SMALL SCALES
The upcoming three paragraphs describe precisely how the local tests are performed to detect short events that are missed in the idealization obtained in Section III-A.If no change is contained, we test a constant signal against the alternative of an additional event on [i/f s , j/f s ] with an arbitrary conductance level.If one change is contained, we test a signal with one change, exact details are discussed below, against the alternative of an additional event on [i/f s , j/f s ] with an arbitrary conductance level.
If the test rejects, the single change is replaced by two and the exact locations and the conductance level between these two changes are obtained in the upcoming deconvolution step in Section III-C.In the rare situation that two or more changes are present no local test is performed to save computation time, since the parameters of more than two changes can anyway not be estimated in the upcoming deconvolution step.Moreover, we only test on intervals with start and end point at the observation grid.Both limitations can be narrowed as discussed in Section VI-A, at the price of a larger computation time.But, we found that our choices are sufficient for the data we analyze.
We now describe how we obtain the parameters τ, c L , c R , s L and s R in the hypotheses in (III.1) and alternatives in (III.2).To this end, note that changes in the idealization from Section III-A are typically slightly shifted to the right, since the convolution is ignored and the lowpass filter acts only in the past.Hence, if we simply obtain the parameters from the previous idealization, many hypotheses will be wrongly rejected, even if the true underlying signal has only one change in [(i − m + 1)/f s , (j + m − 1)/f s ].To correct for this, we reestimate the locations of all isolated changes by deconvolution.This is performed locally as described in the upcoming Section III-C.
Since we assume for testing that all changes are on the observation grid, we perform the deconvolution only at the observation grid, without any refinement at finer scales.This includes a reestimation of the conductance levels on long segments by medians.In other words, as the hypothesis we assume the signal which will be obtained by deconvolution if no test rejects, up to refinements using finer grids.At the same time, estimation of the conductance levels on long segments by the median guarantees that they are not too badly estimated even if few short peaks are missed.
On long segments, in addition to the expectation, the standard deviation s is estimated by using the same observations as used for estimating the expectation.Here, Φ −1 denotes the quantile of the standard normal distribution.
Finally, we recommend to choose l max such that events on all larger scales are already detected by the previous idealization (or have such a small jump size that they are also not detectable by the tests in this step).We found in simulations (not displayed) that l max = 65 is a suitable default choice.

Local testing:
In this paragraph we propose a test that provides a good trade-off between detection power for events in the measurements in Section V and computational complexity.We start with estimating the unknown parameters c and s under the alternative.
To estimate c we use the least squares estimator where with F m the antiderivative (step function) of the truncated filter kernel.Moreover, it follows from (II.3) that under the alternative (III.2) the covariance is given by Hence, for estimating the variance s 2 we use the weighted estimator (VIII.9) Note that the random variable of which we take the expectation in (VIII.9)can be written as a quadratic form ), where Y i,j = (Y i , . . ., Y j ) t and all entries of the matrix C are non-negative and depend only on v l and w l,0 , l = i + 1, . . ., j + m − 1.This combined with (VIII.5)confirms the proposed structure in (VIII.9)follows and allows to computed A and B(s L , s R ) explicitly.
Note that the estimator ĉ is unbiased, while for ŝ this would be true without the projection of negative values to zero in (VIII.8),which however reduces the mean square error.
Using these estimators, under the alternative the observation Y l has estimated expectation ĉ1,l := v l ĉ + c LR,l and estimated variance ŝ2 1,l := w l,0 ŝ2 + s 2 LR,l,0 .Under the null hypothesis the observation Y l has expectation Finally, using these estimators we propose the test statistic (VIII.10) We are aware that this test statistic and its underlying estimators might be improvable with respect to efficiency of the estimators and the power of the resulting test for its various alternatives, for a more detailed discussion and potential alternatives see Section VI.But, as stressed before, we aimed for a test that has at least a good power for the recordings in Section V and can be computed efficiently.This will be confirmed by the simulations in Section IV.
Moreover, note that this test uses information provided by potential standard deviation changes as well.This is different to the multiresolution test we used in Section III-A to detect events on large scales and to HSMUCE.
Note that, different to similar ideas in these settings, it can be computed efficiently, since only testing is required and not regression based on these tests.This is another gain of the three step procedure we propose in this paper.
The test problem is also of a different type than the one in (Enikeeva et al., 2018), since we allow the standard deviation to be constant (s = s L or s = s R ) when the conductance changes.
Critical values: It remains to choose critical values that balance the different tests appropriately.To this end, we obtain again scale depend critical values by using the approach from Section 2 in (Pein et al., 2017).We apply it with significance level α 2 and equal weights β 1 , . . ., β lmax = 1/l max .By this we aim to control the overall probability of detecting an false positive by α := α 1 + α 2 .While we showed in Theorem VII.1 such a control for the multiresolution procedure in Section III-A, we are not able to prove such a bound for the local tests as well, since the previous idealization might not be exactly the true signal up to events on shorter temporal scales and hence the observations are not generated exactly according to the hypotheses (III.1) and alternatives (III.2).Moreover, to speed up the required Monte-Carlo simulations we use the following simplification.When computing the test statics in the Monte-Carlo simulations, we ignore the previous idealization step and assume instead a constant signal.Since the idealization by JSMURF leads with probability at least 1 − α 1 to a constant idealization, this error is negligible.
All in all, we found in simulations that the local tests keep the error level α 2 well.
Multiple dependent rejections: Usually one event in the data causes rejections of multiple tests.Hence, we only add the event that corresponds to the rejection with the largest test statistic among all rejections on intervals that intersect or adjoin each other.More precisely, two rejections on intervals [i 1 /f s , j 1 /f s ] and [i 2 /f s , j 2 /f s ] are only considered as two separated events if the intervals are disjoint and (w.l.o.g.let j 1 < i 2 ) there exists an l ∈ {j 1 + 1, . . ., i 2 − 1} such that all tests on intervals containing l/f s accept the hypothesis.The choice to consider the rejection with the largest test statistic is a natural choice for all tests on intervals of the same length, since they share the same distribution (under their respective null hypotheses and alternatives).For tests on intervals of different lengths this is not exactly true.Nonetheless, we found that considering the rejection with the largest test statistic works very well in practice.That is because usually the test statistics are much larger when their alternative is true than when their null hypothesis is true, which outweighs the (slightly) different distributions (under their respective null hypotheses and alternatives).Also note that a slight missestimation of a location does not have a noticeable effect, since the final estimation of them is obtained in the upcoming deconvolution step.

IX. LOCAL DECONVOLUTION
In this section we describe two minor modifications on the local deconvolution approach of (Pein et al., 2018) which we made to adapt to the heterogeneous noise setting.As summarized at the end of the introduction the local deconvolution approach of (Pein et al., 2018) requires that two short events are separated by at least one long event.Parameters on long events can be estimated without deconvolution.In our setting we have to estimate mean and standard deviation instead of only the mean as in (Pein et al., 2018).Hence, our first modification is that we require at least 25 instead of ten observations in the definition of a long segment to guarantee a reasonable well parameter estimation.Secondly, we adapt the choice of the grids.For the idealization in Section III-A and for the detection step of JULES (Pein et al., 2018)  estimated precisely, but also not systematically shifted to one side.The reestimation of the conductance levels on long segments is adapted in the same way.Everything else is performed in the same way as explained in Section 3.2 of (Pein et al., 2018).

X. IDEALIZATIONS BY EXISTING APPROACHES
In this section we discuss in more detail than in the introduction the idealization of the observations in Figure 1 by various approaches.We begin with approach that assume a hidden Markov model.We observed two different conductance levels and hence started with two states.We used the following starting values   In the lower panels we also show the convolution of the idealization with the lowpass filter (blue).It detects very short events but finds a huge number of (most likely false positive) events at parts of low conductance (see for instance the lower right panel).
scales, but is not able to detect shorter events.If we increase α, this effect reduces, but at the precise of additional false positives, see Figure 22, where we used α = 0.99.This effect is even more pronounced when we use CBS (Venkatraman et al., 2004), see Figure 23, a method that like HSMUCE takes into account heterogeneous noise, but ignores the filtering.However, in comparison to HSMUCE it puts less emphasize on avoiding false positives, In the lower panels we also show the convolution of the idealization with the lowpass filter (blue).It detects short events, but finds many small events (which are most likely false positives) at parts of high conductance and high variance (see for instance the idealization of the observations around 0.36 nS in the middle panel).These detections hinder the decovolution (see for instance the lower left panel) and make the idealization unreliable.

XII. HOMOGENEOUS NOISE
This section details how HILDE can be adapted to the assumption of homogeneous noise.A constant variance can be preestimated as in ( 6) in (Pein et al., 2017).The first multiresolution step to detect events on larger scales can be performed by JSMURF as defined in (Hotz et al., 2013).And for the local tests to detect events on smaller    Fig. 26: Histograms of the dwell times of opening events with amplitude between 0.2 nS and 0.5 nS together with exponential fits using a missed event correction (red line).We consider short dwell times between 0.1 ms and 5 ms.
III.2) with c ∈ R arbitrary.The same structure is assumed for standard deviation functions σ 0 and σ 1 with values s L , s and s R .The precise hypotheses and alternatives, i.e., the values for τ, c L , c R , s L and s R , are determined by the previous idealization step.If more than one change is contained in [(i − m + 1)/f s , (j + m − 1)/f s ], no local test will be performed on this interval.The reasoning behind this and how to obtain τ, c L , c R , s L and s R exactly are explained in the paragraph 'Obtaining the hypotheses and alternatives' in Section VIII in the supplement.All tests are performed at simultaneous error level α 2 > 0.
for the Baum-Welch algorithm.Those standard deviations were determined by taking the empirical standard deviation of all observations below and above 0.2, respectively.Idealizations are obtained by using a Viterbi algorithm.The Baum-Welch algorithm estimated the following parameters µ = (0, 0, 0.3177) T , s = (0.0079, 0.0079, 0.0374) T , P =

Fig. 9 :
Fig. 9: Histograms of the estimated amplitudes using various idealization approaches.The black line (for HILDE hidden by the red line) indicates the true amplitude of 0.32 and the red line the estimated amplitudes of 0.3205, 0.3161 and 0.3148 by the half sample mode using the idealizations from HILDE, JULES and HSMUCE, respectively.

Fig. 10 :
Fig.10: Histograms of the dwell times in the closed state with exponential fits (red).We rescaled all lines such that the area under them are standardized to one to make them comparable to the histograms.

Fig. 12 :
Fig.12: Histograms of the dwell times in the closed state with a length between 20 ms and 200 ms to analyze long events, together with the true exponential distribution (black line) and exponential fits (red) that are corrected for missed events.We rescaled all lines such that the area under them are standardized to one to make them comparable to the histograms.This explains why the true distribution does not look always the same.

Fig. 13 :
Fig. 13: Histograms of the dwell times in the open state, together with the true exponential distribution (black line)and exponential fits (red) that are corrected for missed events.We rescaled all lines such that the area under them are standardized to one to make them comparable to the histograms.This explains why the true distribution does not look always the same.
unilamellar vesicles (GUVs) composed of 1,2-diphytanoyl-sn-glycero-3-phosphocholine (DPhPC)/cholesterol (9:1) were prepared by electroformation (AC, U = 3 V, peak-to-peak, f = 5 Hz, t = 2 h) in the presence of 1 M sucrose at 20 • C. Spreading of a GUV in 1 M KCl, 10 mM HEPES, pH 7.5 on an aperture (d = 1-5 µm) in a borosilicate chip by applying 10-40 mbar negative pressure resulted in a solvent-free membrane with a resistance in the GΩ range.Once the membrane with a GΩ seal was formed, varying amounts of a PorB stock solution (2.2 µM in 200 mM NaCl, 20 mM Tris, 0.1% (w/w) LDAO, pH 7.5) were added to the buffer solution (50 µL) at an applied DC potential of +40 mV.Current traces were recorded at a sampling rate of 10 kHz and filtered with a low-pass four-pole Bessel filter of 1 kHz using an Axopatch 200B amplifier (Axon Instruments).For digitalization, an A/D converter (Digidata 1322; Axon Instruments) was used.

Fig. 14 :
Fig. 14: Histograms of the amplitudes between 0.2 nS and 0.5 nS.The vertical red line indicates the estimated amplitude of 0.3194 nS by the half sample mode.
Obtaining the hypotheses and alternatives: In this paragraph we give details for the construction of the hypothesis (III.1) and alternative (III.2).Assuming the model from Section II, we see that the expectation of the observations Y i , . . ., Y j is determined by the signal on the interval [(i − m)/f s , j/f s ].The other way around, information about the underlying signal on an interval [i/f s , j/f s ] is provided by the observations Y i+1 , . . ., Y j+m−1 and the signal on [(i − m + 1)/f s , (j + m − 1)/f s ] effects the expectation of these observations.Hence, for a local test on an interval [i/f s , j/f s ] we distinguish few scenarios depending on how many changes the previous idealization has in [(i − m + 1)/f s , (j + m − 1)/f s ].
the locations of the estimated changes are shifted to the right, since the convolution was ignored.Hence, we use for a change τ detected in Section III-A (and not replaced by two detected changes by the local tests) still the grid {τ − m/f s , τ − (m + 1)/f s , . . ., τ }.However, for a change τ detected by the local tests in Section III-B we use instead {τ − m/2 /f s , . . ., τ + m/2 /f s }, since the locations are not mean, standard deviation and and transition matrix, respectively.The standard deviations were determined by taking the standard deviation of all observations below and above 0.2, respectively.

Fig. 17 :
Fig.17: Idealization (red) of the data in Figure1assuming a HMM with two states displayed on three different temporal scales.In the lower panels we also show the convolution of the idealization with the lowpass filter (blue).It fits most part of the data well, but misses short events, for instance the event displayed in lower left panel.

Fig. 19 :
Fig.19: Idealization (red) of the data in Figure1assuming a HMM with filtering(Diehn, 2017) displayed on three different temporal scales.In the lower panels we also show the convolution of the idealization with the lowpass filter (blue).It detects very short events but finds a huge number of (most likely false positive) events at parts of low conductance (see for instance the lower right panel).

Fig. 20 :
Fig.20: Idealization (red) of the data in Figure1by TRANSIT displayed on three different temporal scales.In the lower panels we also show the convolution of the idealization with the lowpass filter (blue).It detects short events, but finds many small events (which are most likely false positives) at parts of high conductance and high variance (see for instance the idealization of the observations around 0.36 nS in the middle panel).These detections hinder the decovolution (see for instance the lower left panel) and make the idealization unreliable.

Fig. 21 :
Fig. 21: Idealization by TRANSIT (red) of the signal underlying the observations in Figure 1 and its convolution with the lowpass filter (darkred).It detects a huge amount of false positives at segments with high conductance and high variance.

Fig. 22 :
Fig.22: Idealization (red) of the data in Figure1by HSMUCE using a large significance of α = 0.99 displayed on three different temporal scales.In the lower panels we also show the convolution of the idealization with the lowpass filter (blue).It still misses short events (for instance the one displayed in the lower left panel), but also includes several false positives.

Fig. 23 :
Fig. 23: Idealization (red) of the data in Figure 1 by CBS displayed on three different temporal scales.In the lower panels we also show the convolution of the idealization with the lowpass filter (blue).It detects a huge amount of false positives at segments with high conductance and high variance.

Fig. 24 :
Fig.24: Histograms of the estimated amplitudes using various idealization approaches.The the red line indicates the estimated amplitudes of 0.3194 nS, 0.3182 nS and 0.3216 nS by the half sample mode using the idealizations from HILDE, JULES and HSMUCE, respectively.

Fig. 25 :
Fig.25: Histograms of the dwell times of opening events with amplitude between 0.2 nS and 0.5 nS together with exponential fits using a missed event correction (red line).We consider all dwell times between 0.1 ms and 200 ms.

Fig. 27 :
Fig.27: Histograms of the dwell times of opening events with amplitude between 0.2 nS and 0.5 nS together with exponential fits using a missed event correction (red line).We consider long dwell times between 20 ms and 200 ms.

Fig. 28 :
Fig. 28: Histograms of the dwell times in the closed state, together with exponential fits (red) that are corrected for missed events.

TABLE I :
Performance of HILDE in idealizing a signal with an isolated peak in different settings.They differ in the amount of open channel

TABLE II :
Performance of HILDE in idealizing a signal with an isolated peak in different settings.They differ in the amount of open channel noise and the length of the peak.More precisely, it has changes at τ 1 = 0.2 and τ 2 = τ 1 + /fs, = 2, 3, 5, conductance levels

TABLE III :
Performance of HILDE in idealizing a signal with an isolated peak in different settings.They differ in the amount of open