%\VignetteIndexEntry{R packages: LaTeX vignettes}
%\VignetteEngine{R.rsp::tex}

\documentclass[journal, onecolumn, draftclsnofoot]{IEEEtran}

\usepackage{xr}
\usepackage{amsthm}
\usepackage[mathscr]{eucal}
\usepackage{amsfonts}
\usepackage[cmex10]{amsmath}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage{caption}
\usepackage{subcaption}
\usepackage{pstricks}
\usepackage{hyperref}
\usepackage{booktabs}
\usepackage{colortbl}
\usepackage{comment}
\usepackage{tabularx}

\usepackage{xcolor}

\usepackage[authoryear,round]{natbib}


\newcolumntype{N}{>{\centering\arraybackslash}m{.85in}}

%\newrgbcolor{myblue}{0.2 0.42 1}
\newrgbcolor{myblue}{0 0 1}
\newrgbcolor{myred}{0.466667 0 0}
\newrgbcolor{mygreen}{0 0.466667 0}
\newrgbcolor{mylightred}{1 0 0}
\newrgbcolor{mylightblue}{0 0.8 1}

\newcommand{\ttblue}[1]{\textcolor{myblue}{#1}}
\newcommand{\ttorange}[1]{\textcolor{myorange}{#1}}
\newcommand{\ttred}[1]{\textcolor{myred}{#1}}
\newcommand{\ttlred}[1]{\textcolor{mylightred}{#1}}
\newcommand{\ttlblue}[1]{\textcolor{mylightblue}{#1}}
\newcommand{\ttgreen}[1]{\textcolor{mygreen}{#1}}


\usepackage{tabularx,ragged2e,booktabs,caption}
\newcolumntype{C}[1]{>{\Centering}m{#1}}
\renewcommand\tabularxcolumn[1]{C{#1}}

\newcommand\solidrule[1][0.5cm]{\rule[0.5ex]{#1}{.8pt}}
\newcommand\dashedrule{\mbox{%
  \solidrule[1mm]\hspace{1mm}\solidrule[1mm]\hspace{1mm}\solidrule[1mm]}}

\newcommand\dottedrule{\mbox{%
  \solidrule[0.5mm]\hspace{1mm}\solidrule[1mm]\hspace{1mm}\solidrule[0.5mm]\hspace{1mm}\solidrule[1mm]}}

\newcommand\dotsrule{\mbox{%
  \solidrule[0.5mm]\hspace{.45mm}\solidrule[0.5mm]\hspace{.45mm}\solidrule[0.5mm]}}

\usepackage{setspace}
\doublespacing

\hyphenation{op-tical net-works semi-conduc-tor}

\newcommand{\mb}[1]{\mathbb{#1}}
\newcommand{\mc}[1]{\mathcal{#1}}
\newcommand{\msc}[1]{\mathscr{#1}}
\newcommand{\bs}[1]{\boldsymbol{#1}}
\newcommand{\mf}[1]{\mathfrak{#1}}
\newcommand{\fbm}[1]{\bs{\msc{#1}}}
\newcommand{\fmb}[1]{\bs{\mb{#1}}}
\newcommand{\fmc}[1]{\bs{\mc{#1}}}
\newcommand{\msf}[1]{\mathsf{#1}}


\newcommand{\tn}[1]{\textnormal{#1}}
\newcommand{\inner}[2]{\langle #1,#2\rangle}
\newcommand{\Bf}[1]{{\bf #1}}
\newcommand{\B}[1]{{\bf #1}}
\newcommand{\ind}{{1\hspace{-2.5pt}\tn{l}}}


\newcommand{\ic}{\bs{i}}
\newcommand{\weak}[1]{\stackrel{#1}{\rightsquigarrow}}
\newcommand{\snorm}[1]{|\hspace{-0.99mm}|\hspace{-0.99mm}| #1 |\hspace{-0.99mm}|\hspace{-0.99mm}|}
\newcommand{\normi}{\left\|}
\newcommand{\normd}{\right\|}
\newcommand{\tv}[2]{\| #1 - #2 \|_{\tt{TV}}}
\newcommand{\Fnorm}[1]{\normi #1 \normd_{\tt{F}}}
% ----------------------------------------------------------------------------------
\DeclareMathOperator*{\argmin}{argmin}

% ----------------------------------------------------------------------------------
\newcommand{\e}{\B{e}}
\newcommand{\R}{\mb{R}}
\newcommand{\Z}{\mb{Z}}
\newcommand{\N}{\mb{N}}
\newcommand{\E}{\mb{E}}
\newcommand{\Q}{\mb{Q}}
\newcommand{\T}{\mb{T}}
% \newcommand{\I}{\mbol{I}}
\newcommand{\p}{\mb{P}}
\newcommand{\Sb}{\mb{S}}
\newcommand{\Sm}{\mc{S}}
\newcommand{\V}{\mb{V}}
\newcommand{\C}{\mb{C}}
% \newcommand{\F}{\mbol{F}}
\newcommand{\prob}{\B{P}}
\newcommand{\probQ}{\B{Q}}
% \newcommand{\snorm}[1]{\pmb{|}#1\pmb{|}}
% ----------------------------------------------------------------------------------
\newcommand{\wh}[1]{\widehat{#1}}
\newcommand{\wt}[1]{\widetilde{#1}}
\newcommand{\br}[1]{\breve{#1}}
\newcommand{\vt}{\vartheta}
% ----------------------------------------------------------------------------------
\newcommand{\ME}{\msf{E}}
\newcommand{\VAR}{\msf{VAR}}
\newcommand{\SD}{\msf{SD}}
\newcommand{\MSE}{\msf{MSE}}
\newcommand{\COV}{\msf{Cov}}
\newcommand{\BIAS}{\msf{BIAS}}



\newcommand{\todo}[1]{\textcolor{red}{TODO: #1}}

\newcommand{\kernel}{F}
\newcommand{\sr}{f_{\bs{s}}}
\newcommand{\cutoff}{f_{\bs{c}}}

\newcommand{\JULES}{\operatorname{JULES}}
\newcommand{\HSMUCE}{\operatorname{HSMUCE}}
\newcommand{\HMM}{\operatorname{HMM}}
\newcommand{\JSMURF}{\operatorname{JSMURF}}
\newcommand{\TRANSIT}{\operatorname{TRANSIT}}
\newcommand{\MDL}{\operatorname{MDL}}
\newcommand{\HILDE}{\operatorname{HILDE}}
\newcommand{\SMUCE}{\operatorname{SMUCE}}
\newcommand{\valmu}{c}
\newcommand{\median}{\operatorname{median}}

\usepackage{siunitx}
\DeclareSIUnit{\Molar}{\textsc{m}}
\sisetup{range-phrase = --}

\usepackage{algorithm2e}
\usepackage{tikz}
\usetikzlibrary{shapes,arrows,backgrounds,fit,calc,positioning}

\numberwithin{equation}{section}
\newcounter{satze}
\numberwithin{satze}{section}
\theoremstyle{plain} 
\newtheorem{Proposition}[satze]{Proposition}
\newtheorem{Lemma}[satze]{Lemma}
\newtheorem{Corollary}[satze]{Corollary}
\newtheorem{Theorem}[satze]{Theorem}


\makeatletter
\def\maxwidth{ %
  \ifdim\Gin@nat@width>\linewidth
    \linewidth
  \else
    \Gin@nat@width
  \fi
}
\makeatother

\definecolor{fgcolor}{rgb}{0.345, 0.345, 0.345}
\newcommand{\hlnum}[1]{\textcolor[rgb]{0.686,0.059,0.569}{#1}}%
\newcommand{\hlstr}[1]{\textcolor[rgb]{0.192,0.494,0.8}{#1}}%
\newcommand{\hlcom}[1]{\textcolor[rgb]{0.678,0.584,0.686}{\textit{#1}}}%
\newcommand{\hlopt}[1]{\textcolor[rgb]{0,0,0}{#1}}%
\newcommand{\hlstd}[1]{\textcolor[rgb]{0.345,0.345,0.345}{#1}}%
\newcommand{\hlkwa}[1]{\textcolor[rgb]{0.161,0.373,0.58}{\textbf{#1}}}%
\newcommand{\hlkwb}[1]{\textcolor[rgb]{0.69,0.353,0.396}{#1}}%
\newcommand{\hlkwc}[1]{\textcolor[rgb]{0.333,0.667,0.333}{#1}}%
\newcommand{\hlkwd}[1]{\textcolor[rgb]{0.737,0.353,0.396}{\textbf{#1}}}%
\let\hlipl\hlkwb

\usepackage{framed}
\makeatletter
\newenvironment{kframe}{%
 \def\at@end@of@kframe{}%
 \ifinner\ifhmode%
  \def\at@end@of@kframe{\end{minipage}}%
  \begin{minipage}{\columnwidth}%
 \fi\fi%
 \def\FrameCommand##1{\hskip\@totalleftmargin \hskip-\fboxsep
 \colorbox{shadecolor}{##1}\hskip-\fboxsep
     % There is no \\@totalrightmargin, so:
     \hskip-\linewidth \hskip-\@totalleftmargin \hskip\columnwidth}%
 \MakeFramed {\advance\hsize-\width
   \@totalleftmargin\z@ \linewidth\hsize
   \@setminipage}}%
 {\par\unskip\endMakeFramed%
 \at@end@of@kframe}
\makeatother

\definecolor{shadecolor}{rgb}{.97, .97, .97}
\definecolor{messagecolor}{rgb}{0, 0, 0}
\definecolor{warningcolor}{rgb}{1, 0, 1}
\definecolor{errorcolor}{rgb}{1, 0, 0}
\newenvironment{knitrout}{}{} % an empty environment to be redefined in TeX

\usepackage{alltt}
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}


\begin{document}

\title{Analysis of Patchclamp Recordings:\\ Model-Free Multiscale Methods and Software}

\author{Florian~Pein, Benjamin~Eltzner and Axel~Munk% 
\thanks{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.}%
\thanks{Benjamin~Eltzner and Axel Munk are with the Institute for Mathematical Stochastics, Georg-August-University of Goettingen, Goldschmidtstr. 7, 37077 G\"ottingen, Germany.}%
\thanks{Axel Munk is also with the Max Planck Institute for Biophysical Chemistry, Am Fassberg 11, 37077 G\"ottingen, Germany, and with the Felix Bernstein Institute for Mathematical Statistics in the Biosciences, Goldschmidtstr. 7, 37077 G\"ottingen, Germany.}%
\thanks{Financially support of DFG (CRC803, project Z02) over the past twelve years is gratefully acknowledge. We are grateful to our collaborators Annika~Bartsch, Ulf~Diederichsen, Manuel~Diehn, Thomas~Hotz, Ingo~P.~Mey, Tatjana~Polupanow, Ole~M.~Sch\"utte, Ivo~Siekmann, Hannes~Sieling, Claudia~Steinem, Inder~Tecuapetla-G\'omez and Laura~Yineth~Jula~Vanegas. Our special thank goes to Florian~Ebmeier, Mariyam~Khan and Stanislav~Syekirin for their work on the graphical user interface. F. Pein was also supported by EPSRC (Statscale programme). A. Munk was also supported by DFG (Cluster of excellence 2067 MBExC Multiscale Bioimaging: From Molecular Machines to Networks of Excitable Cells) and the Volkswagen Foundation (FBMS).}}%

\maketitle

\begin{abstract}
Analysis of patchclamp recordings is often a challenging issue. We give practical guidance how such recordings can be analyzed using the model-free multiscale idealization methodology $\JSMURF$, $\JULES$ and $\HILDE$. We provide an operational manual how to use the accompanying software available as an R-package and as a graphical user interface. This includes selection of the right approach and tuning of parameters. We also discuss advantages and disadvantages of model-free approaches in comparison to hidden Markov model approaches and explain how they complement each other.
\end{abstract}

\begin{IEEEkeywords}
Deconvolution, flickering, fully-automatic, hidden Markov models, homogeneous and heterogeneous noise, ion channel recordings, lowpass filtering, open channel noise, PorB, subconductance states
\end{IEEEkeywords}

\IEEEpeerreviewmaketitle

\section{Introduction}\label{sec:introduction}
The patchclamp technique has been and still is a fundamental tool for the quantitative analysis of electrophysiological processes of transmembrane proteins, in particular of ion channels, \citep{Neher.Sakmann.76, Sakmann.Neher.95}. A detailed understanding of the dynamics of transmembrane proteins and their manifold interactions with their surrounding is of high importance in medicine and biochemistry, for instance for the development of new drugs \citep{kass2005channelopathies, overington2006many}. However, most electrophysiologists will agree that conducting patchclamp experiments but also the analysis of their recordings is a challenging issue, and the latter is far from being a routine data analysis in general \citep{sivilotti2016praise}. In this work we provide practical guidance on how to analyze such recordings. We focus mainly on model-free multiscale idealizations (explained below), which we have developed over the last decade.

\paragraph*{Patchclamp recordings}
The patchclamp technique allows to measure the conductance of a channel (i.e., the recorded current divided by the applied voltage) over time. An example is given in Figure \ref{fig:PorBData}. It shows a recording of the outer membrane porin PorB from \textit{Neisseria meningitidis}, a pathogenic bacterium in the human nose and throat region \citep{virji2009pathogenic}. PorB is a trimeric porin and the second most abundant protein in the outer membrane of \textit{Neisseria meningitidis}. The added antibiotic ampicillin blocks the ion flow for short periods of time which allows to draw conclusions about the transport of antibiotics into the cell, which is relevant for the understanding of antibiotic resistances. For further details see \citep{Bartschetal18}. In addition, we will also use a PorB dataset without ampicillin (Figure \ref{fig:PorBhetData}) and a Gramicidin A dataset (Figure \ref{fig:GramicidinData}) throughout this work as illustrating examples\footnote{All shown measurements were taken at the Steinem lab, Institute of Organic and Biomolecular Chemistry, Georg-August-University of G{\"o}ttingen.}.

\begin{figure}[!htb]
\centering
\includegraphics[width = 0.99\textwidth]{figures/PorBwithAmpData.pdf}
\caption{From seconds to microseconds: Patchclamp recording (grey points) displayed at the level of seconds (top panel), of milliseconds (middle panel) and of microseconds (bottom panels). Data points result from a representative conductance recording of PorB wild type with $1$ mM ampicillin by the patchclamp technique using black lipid membranes at $\SI{80}{\milli\volt}$. Data points are explicitly displayed instead of a line plot, which provides an accurate representation at fine resolution levels.}
\label{fig:PorBData}
\end{figure} 

\paragraph*{Idealization}
Important dynamics such as the number of conductance levels, their values and how long each level persists can be examined provided the conductance recordings (data points) are properly idealized \citep{Colquhoun.87, Sakmann.Neher.95}, i.e., the conductance trace over time (the underlying signal) is accurately reconstructed (estimated, denoised).

An idealization can either be obtained model-free\footnote{Strictly speaking the terminology 'model-free' is misleading as also model-free approaches require an underlying model to be valid. However, we use this terminology mainly for historical reasons. The precise (non-parametric) models underlying our methodology are reviewed in Section \ref{sec:models}. We only assume that the underlying conductance is piecewise constant, but make no further assumptions about the gating dynamics.}, i.e., without prior assumptions about the gating dynamics, or in a model based way by assuming an underlying statistical (parametric) model with a few parameters for the gating dynamics. For the latter most commonly hidden Markov models (HMMs) are used, see \citep{Ball.Rice.92} for an early reference, where parameters correspond to states, transition probabilities and noise characteristics.

\paragraph*{Filtering}
The noise before filtering is often assumed to be Gaussian white noise. However, patchclamp recordings are usually lowpass filtered, integrated in the hardware of the measurement device, to stay in the transmission range of the amplifier. Such filtering introduces colored noise and smooths the underlying conductance, see Figure \ref{fig:convolution} in Section \ref{sec:existingMethods}. Ignoring filtering typically results in the detection of false positives (additional wrong events). This is illustrated in Figure \ref{fig:PorBSMUCE}, where we used $\SMUCE$ \citep{frick2014asymptotic}, a multiscale method that does not include filtering but is otherwise similar in spirit to our model-free idealization approaches, to be explained later. Filtering especially affects short temporal scales (at and below the magnitude of the filter length, say) and is therefore particularly relevant to the analysis of short events, also called \textit{flickering}.

\begin{figure}[!htb]
\centering
\includegraphics[width = 0.99\textwidth]{figures/PorBwithAmpSMUCE.pdf}
\caption{Idealization (red) of the observations in Figure \ref{fig:PorBData} by $\SMUCE$ \citep{frick2014asymptotic} displayed on three different temporal scales. In the lower panels we also show the convolution of the idealization with the lowpass filter (blue). $\SMUCE$ is a multiscale method similar to our model-free idealization approaches, but does \textit{not} take into account filtering. Hence, it overestimates the number of conductance changes massively, since it misinterprets local variations due to colored noise as events and splits abrupt conductance changes in multiple steps since the convolution of the conductance with the lowpass filter is ignored.}
\label{fig:PorBSMUCE}
\end{figure} 

\paragraph*{Flickering and subgating}
Flickering typically has its own dynamics and can result from various molecular processes like conformational changes of the protein \citep{Grosse.etal.14} or by the passage of larger molecules blocking the ions' pathway through the protein \citep{Singh.etal.12}. An example for the latter is the PorB analysis in Figures \ref{fig:PorBData} and \ref{fig:PorBHILDE}. A second potential challenge in the analysis are \textit{subconductance} states \citep{fox1987ion}, meaning that two or more conductance levels are close to each other, as illustrated in Figures \ref{fig:GramicidinData} and \ref{fig:GramicidinJSMURF}.

\begin{figure}[!htb]
\centering
\includegraphics[width = 0.99\textwidth]{figures/PorBwithAmpHILDE.pdf}
\caption{Idealization (red) of the observations in Figure \ref{fig:PorBData} by $\HILDE$ \citep{HILDE} displayed on three different temporal scales. Lower panels: Convolution of the idealization with the lowpass filter (blue). Events are well idealized down to microseconds.}
\label{fig:PorBHILDE}
\end{figure}

\begin{figure}[!htb]
\centering
\includegraphics[width = 0.99\textwidth]{figures/GramicidinData.pdf}
\caption{From seconds to milliseconds: Patchclamp recording (grey points) displayed at the level of seconds (top panel and middle panel) and of milliseconds (bottom panels). Data points result from a representative conductance recording of an acylated gramicidin A derivative by the patch clamp technique using solvent-free bilayers at $\SI{100}{\milli\volt}$. Small conductance changes, which most likely result from subconductance states, occur frequently.}
\label{fig:GramicidinData}
\end{figure} 

\begin{figure}[!htb]
\centering
\includegraphics[width = 0.99\textwidth]{figures/GramicidinJSMURF.pdf}
\caption{Idealization (red) of the observations in Figure \ref{fig:GramicidinData} by $\JSMURF$ \citep{JSMURF} displayed on three different temporal scales. It idealizes the observations well and finds in particular a larger number of gating events with a small conductance change (subconductance states).}
\label{fig:GramicidinJSMURF}
\end{figure} 

\paragraph*{Model-free idealizations}
In this paper, we discuss mainly our model-free idealization methods $\JSMURF$ \citep{JSMURF}, $\JULES$ \citep{JULES} and $\HILDE$ \citep{HILDE}, which are primarily designed as versatile tools to analyze patchclamp recordings in a multiscale fashion, for instance to deal well with \textit{subconductance} states and \textit{flickering}. Due to their multiscale character they act on various temporal scales simultaneously and hence are able to idealize events of different lengths well in a single step. Moreover, all parameters, i.e., locations of conductance changes and conductance levels, are obtained by (local) deconvolution, hence they take into account lowpass filtering explicitly. Furthermore, all three approaches control the overestimation of the number of conductance changes. More precisely, the probability to detect at least one false positive is bounded approximately by the error level $\alpha$, a tuning parameter. A more detailed review of these and further model-free idealization approaches is given in Section \ref{sec:existingMethods}.

All three methods can be used when homogeneous noise is assumed (i.e. the error variability does not change over time, see Section \ref{sec:models} for details), but $\JSMURF$ and $\HILDE$ in addition allow for heterogeneous noise. The latter means that different parts of the data, for instance different states, may have different noise levels, as for instance caused by open channel noise. Moreover, $\JSMURF$ requires that events are slightly longer, while $\JULES$ or $\HILDE$ are able to deal with flickering (short events) at the possible expense of longer computational time. Figure \ref{fig:methodSelection} summarizes for which datasets which one of them is most suitable and Section \ref{sec:guidanceMethod} explains those choices in full detail.

\begin{figure}[!htb]
\centering
\includegraphics[width = 0.75\textwidth]{figures/methodSelection.pdf}
\caption{Table for the selection of the right model-free idealization method. For more details see Section \ref{sec:guidanceMethod}.}
\label{fig:methodSelection}
\end{figure} 

\paragraph*{Software}
$\JSMURF$, $\JULES$ and $\HILDE$ are implemented as \texttt{R} functions in the package \textit{clampSeg} \citep{clampSeg}. In Section \ref{sec:softwareUse} use of those methods is demonstrated. They can be combined with the packages \textit{readABF} \citep{readABF} to load recordings, and \textit{lowpassFilter} \citep{lowpassFilter} for certain data processing steps around filtering such as computing the convolution of an idealization with the kernel of a lowpass filter.

Alternatively, a graphical user interface, available at \url{https://github.com/FlorianPein/clampSegGUI} together with detailed manuals on how to install it and on how to use it, allows access without requiring any \texttt{R} or other programming knowledge. The idealizations can be visualized in the interface but also saved as \textit{csv} files and hence post-processed by any other program.

\paragraph*{Interplay between model-free idealizations and hidden Markov models}
There is wide agreement that, except in few counterexamples \citep{Fulinski.98, Mercik.Weron.01, Goychuk.etal.05, Shelley.etal.10}, the gating dynamics of ion and many other channels are usually Markovian. Hence, hidden Markov model (HMM) based approaches are widely used to analyze partchclamp recordings. However, from our own experience we stress that the hidden part is usually the critical part of the assumption of a (homogeneous) HMM. This is because of artifacts, which are for instance caused by the electronics, external vibrations or small holes in the membrane, or because of additional high frequency $f^2$ (violet) and long tailed $1/f$ (pink) noise components, see for instance \citep{Neher.Sakmann.76, venkataramanan1998identification, levis1993use}. Hence, HMM based analyses often rely on intensive preprocessing or on more complicated models: For instance \citep{venkataramanan1998identification} assumed an HMM that allows additional colored noise, and \citep{diehn2019maximum} provided modifications to incorporate inhomogeneous errors. Moreover, lowpass filtering often requires further, computationally demanding extensions, see for instance \citep{venkataramanan1998identification, deGunst.etal.01, Diehn17, almanjahie2019moving}.

In contrast, model-free idealizations do not assume a specific (parametric) model for the gating dynamics and immediately provide an idealization without such an assumption. Moreover, they usually act rather locally on the data. Hence, they are typically more robust to artifacts and often simpler to apply. Contrarily, HMM based approaches have (potentially) a finer time resolution and provide more concise results. A more detailed review of HMM based analyses, their advantages and disadvantages in comparison to model-free approaches and their interplay is given in Section \ref{sec:HMM}.

Model-free idealizations allow a flexible analysis of the number of conductance states, their values and which transitions are possible. To fit a Markov model, one has to cluster the conductance values, e.g., by fitting a Gaussian mixture distribution and assigning each value to the nearest mean value. The outcome will then have only a small number of conductance levels. This can be used to select and verify a Markov model and to estimate its parameters, which often requires taking into account missing of short events. Further details and tools which can be used for those steps are described in Section \ref{sec:analysis}. Finally, model-free idealizations can be used to assist HMM based approaches in any of their analysis steps, e.g., they can be used to remove artifacts, to select and verify a specific Markov model, to provide starting values for iterative procedures such as the Baum-Welch algorithm and to verify estimated parameters and the provided idealization.
 
All in all, model-free and HMM approaches have different strengths and weaknesses and hence should less be seen as competing approaches, but rather as tools that benefit from and complement each other. In fact, as an indication of a proper data analysis it can be checked whether the results of model-free and HMM based analyses are in compliance.

\paragraph*{Organization of this work}
In Section \ref{sec:mf-idealization} we give detailed instructions how to use our methods to obtain model-free idealizations. In Section \ref{sec:models} we provide details of the statistical models underlying the presented model-free idealization methodology. In Section \ref{sec:idealizations} we review in more detail existing approaches for the analysis of patchclamp recordings. We start in Section \ref{sec:HMM} with a review of HMM based methodology, their advantages and disadvantages in comparison to model-free approaches and the interplay of model-free idealization with them. Afterwards, in Section \ref{sec:existingMethods} we review some existing model-free idealization methods, with a particular focus on our approaches $\JSMURF$, $\JULES$ and $\HILDE$. This is complemented by a brief summary of simulation results. Finally, in Section \ref{sec:analysis} we discuss how (model-free) idealizations can be used to analyze patchclamp recordings. The paper concludes with a discussion in Section \ref{sec:discussion}, in which we highlight open research questions.


\section{Model-Free Idealizations}\label{sec:mf-idealization}
This section provides a comprehensive guide on how to use our methods $\JSMURF$, $\JULES$ and $\HILDE$, which have different strengths and weaknesses depending on certain structural features of the measured data. We explain in detail how to use our software (Section \ref{sec:softwareUse}) and provide guidance for which method is preferable in which situations (Section \ref{sec:guidanceMethod}).

\subsection{Using our software}\label{sec:softwareUse}
For this section we use \texttt{R} code\footnote{https://www.r-project.org/}. We start by describing how recordings can be loaded, how the lowpass filter can be specified and how our methods can be called. To this end, we require the R packages \textit{readABF} \citep{readABF}, \textit{lowpassFilter} \citep{lowpassFilter} and \textit{clampSeg} \citep{clampSeg}. All three packages are available on CRAN\footnote{https://cran.r-project.org/}. For users who are not familiar with \texttt{R}, we also provide a graphical user interface\footnote{https://github.com/FlorianPein/clampSegGUI}, which contains detailed manuals on how to install it and on how to use it.

\paragraph*{Loading the data}
Patchclamp recordings are typically stored as \textit{abf} files. The \textit{readABF} package allows to read such files in \texttt{R}.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{library}\hlstd{(readABF)}
\hlcom{# path and filename have to adjusted}
\hlstd{data} \hlkwb{<-} \hlkwd{readABF}\hlstd{(}\hlstr{"../data/1 M KCl 10 mM HEPES pH 75_0268.abf"}\hlstd{)}
\hlcom{# convert it to a data.frame with two columns:}
\hlcom{# time and conductance (current divided by voltage)}
\hlstd{data} \hlkwb{<-} \hlkwd{as.data.frame}\hlstd{(data,} \hlkwc{type} \hlstd{=} \hlstr{"one"}\hlstd{,} \hlkwc{channels} \hlstd{=} \hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,} \hlnum{2}\hlstd{),} \hlkwc{unit} \hlstd{=} \hlstr{"nS"}\hlstd{)}
\hlstd{time} \hlkwb{<-} \hlstd{data}\hlopt{$}\hlstd{`Time [s]`}
\hlstd{data} \hlkwb{<-} \hlstd{data}\hlopt{$}\hlstd{`Data [nS]`}
\end{alltt}
\end{kframe}
\end{knitrout}

Our idealization methods but also any other approach should not be used as a black-box. We strongly recommend to start with an empirical and visual data analysis to gain understanding of the datasets and major features that can be used to direct further analysis. In alignment with the underlying multiscale philosophy of our methods we recommend always to plot on various temporal scales. See Figure \ref{fig:PorBData} for an example of such plots on three different scales ranging from a minute to milliseconds. Moreover, histograms of the raw data (point amplitude histograms), see for instance Figure \ref{fig:pointHistogram} and for code the paragraph 'Interpreting, plotting and verification of the output' below, are helpful visual cues. As detailed below, this can already help to decide whether the noise is homogeneous or heterogeneous and whether short events occur in the dataset. Moreover, we recommend to identify potential artifacts that might disturb analysis and interpretation. However, we have found that model-free idealizations are usually quite robust to artifacts. Hence, our default suggestion is first to apply the idealization methods on the unmodified dataset and to decide later whether artifacts require a more careful analysis.

\paragraph*{Lowpass filter}
Our methodology requires to specify correctly the lowpass filter in the measurement device. The type, often a Bessel filter with an even number of poles, should be specified in the hardware documentation. The sampling rate and cut-off frequency can typically be varied by the user. In our example (Figure \ref{fig:PorBData}), the recordings were sampled at $50\,000$ Hz and lowpass filtered by a $4$-pole Bessel filter with normalized cutoff frequency of $0.1$ ($5\,000$ Hz cutoff frequency in time domain). 

For simplification and since the error is negligible we truncate the kernel of the lowpass filter after $m$ data points, for sufficiently large $m$. As a working rule, we choose $m$ such that the autocorrelation function of the untruncated analogue lowpass filter is below $10^{-3}$ afterwards, which leads for instance to $m = 11$ in the example above.

This is implemented in the function \textit{lowpassFilter} in the package \textit{lowpassFilter} \citep{lowpassFilter} (currently only Bessel filters are supported). The following code creates the filter object.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{library}\hlstd{(lowpassFilter)}
\hlstd{filter} \hlkwb{<-} \hlkwd{lowpassFilter}\hlstd{(}\hlkwc{type} \hlstd{=} \hlstr{"bessel"}\hlstd{,}
                        \hlkwc{param} \hlstd{=} \hlkwd{list}\hlstd{(}\hlkwc{pole} \hlstd{=} \hlnum{4L}\hlstd{,} \hlkwc{cutoff} \hlstd{=} \hlnum{0.1}\hlstd{),} \hlkwc{sr} \hlstd{=} \hlnum{5e4}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}

We strongly recommend to verify that the filter is correctly specified by zooming into single events and checking whether the obtained idealization convolved with the lowpass filter fits the observations well. This is detailed in paragraph 'Interpreting, verification and storing of the output' below, where we discuss in more generality how to assess the quality of an obtained idealization. Additionally, one can compare the auto-correlation resulting from the filter, \textit{filter\$acf}, with the estimated auto-correlation of the recordings. To this end, one can either apply standard time series estimators, as for instance offered by the \textit{acf} function in \texttt{R}, to long segments without conductance changes or use the robust difference based estimators of \citep{Munk.Tecuapetla.15} on the raw data, available in the \texttt{R}-package \textit{dbacf}\footnote{http://www.stochastik.math.uni-goettingen.de/dbacf}.

\paragraph*{Obtaining an idealization}
$\JSMURF$, $\JULES$ and $\HILDE$ are available in the package \textit{clampSeg}. All three functions can be applied when homogeneous noise is assumed, but only $\JSMURF$ and $\HILDE$ allow for heterogeneous noise. The following code illustrates how to call those functions depending on whether the noise is homogeneous or heterogeneous. See Section \ref{sec:guidanceMethod} for guidance which method and which noise option should be chosen to idealize a given measurement.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{library}\hlstd{(clampSeg)}
\hlcom{# JSMURF assuming homogeneous noise}
\hlstd{fit1} \hlkwb{<-} \hlkwd{jsmurf}\hlstd{(data,} \hlkwc{filter} \hlstd{= filter,} \hlkwc{family} \hlstd{=} \hlstr{"jsmurfPS"}\hlstd{,}
               \hlkwc{alpha} \hlstd{=} \hlnum{0.05}\hlstd{,} \hlkwc{r} \hlstd{=} \hlnum{1e4}\hlstd{)}
\hlcom{# JSMURF allowing heterogeneous noise}
\hlstd{fit2} \hlkwb{<-} \hlkwd{jsmurf}\hlstd{(data,} \hlkwc{filter} \hlstd{= filter,} \hlkwc{family} \hlstd{=} \hlstr{"hjsmurf"}\hlstd{,}
               \hlkwc{alpha} \hlstd{=} \hlnum{0.05}\hlstd{,} \hlkwc{r} \hlstd{=} \hlnum{1e4}\hlstd{)}
\hlcom{# JULES assuming homogeneous noise}
\hlstd{fit3} \hlkwb{<-} \hlkwd{jules}\hlstd{(data,} \hlkwc{filter} \hlstd{= filter,} \hlkwc{alpha} \hlstd{=} \hlnum{0.05}\hlstd{,} \hlkwc{r} \hlstd{=} \hlnum{1e4}\hlstd{)}
\hlcom{# HILDE assuming homogeneous noise}
\hlstd{fit4} \hlkwb{<-} \hlkwd{hilde}\hlstd{(data,} \hlkwc{filter} \hlstd{= filter,} \hlkwc{family} \hlstd{=} \hlstr{"jsmurfPS"}\hlstd{,} \hlkwc{method} \hlstd{=} \hlstr{"LR"}\hlstd{,}
              \hlkwc{alpha1} \hlstd{=} \hlnum{0.01}\hlstd{,} \hlkwc{alpha2} \hlstd{=} \hlnum{0.04}\hlstd{,} \hlkwc{r} \hlstd{=} \hlnum{1e3}\hlstd{,} \hlkwc{lengths} \hlstd{=} \hlnum{1}\hlopt{:}\hlnum{20}\hlstd{)}
\hlcom{# HILDE allowing heterogeneous noise}
\hlstd{fit5} \hlkwb{<-} \hlkwd{hilde}\hlstd{(data,} \hlkwc{filter} \hlstd{= filter,} \hlkwc{family} \hlstd{=} \hlstr{"hjsmurf"}\hlstd{,} \hlkwc{method} \hlstd{=} \hlstr{"2Param"}\hlstd{,}
              \hlkwc{alpha1} \hlstd{=} \hlnum{0.01}\hlstd{,} \hlkwc{alpha2} \hlstd{=} \hlnum{0.04}\hlstd{,} \hlkwc{r} \hlstd{=} \hlnum{1e3}\hlstd{,} \hlkwc{lengths} \hlstd{=} \hlnum{1}\hlopt{:}\hlnum{65}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}

The followings paragraphs discuss run time, required Monte-Carlo simulations, the output of the approaches and how to proceed with it. Further, we explain a potentially occurring warning and how to choose tuning parameters, e.g., \textit{r} and \textit{alpha}.

\paragraph*{Run time and Monte-Carlo simulations}
The run time of all approaches depends on the size of the dataset, but also on the number of detected events. The primary reason are Monte-Carlo simulations which are required to obtain critical values that balance the probabilities of detection of true events and of false positives. Monte-Carlo simulations depend on the number of data points and on the lowpass filter. Hence, a new Monte-Carlo simulation is required when new values for those parameters occur or when more repetitions are requested. Depending on the number of data points and the total number of repetitions $r$, Monte-Carlo simulations may take long, even up to several hours. Hence, we store and load their results such that they have to be performed only once. They are fully automatically stored in the workspace and on the disk of the local computing machine, for more details see the documentation of the function \textit{getCritVal} in the package \textit{clampSeg} \citep{clampSeg}. To keep track of the progress of a Monte-Carlo simulation one can set the argument \textit{messages} to a positive integer value $m$ to print a message every $m$ repetitions.

While a larger number of repetitions increases the run time of the simulations, it also reduces statistical errors in the computation of the critical values. For a final analysis we recommend to use the default values, $10\,000$ for $\JSMURF$ and $\JULES$ and $1\,000$ for $\HILDE$. For a quick analysis, for instance to decide whether further measurements or analyses are required, few hundreds up to $1\,000$ repetitions usually suffice.

Additionally, also the main computation of the idealization can take some time, usually between few seconds and few minutes, depending on the used idealization method, on the size of the dataset and on the number of detected events. Usually, the run time increases with the complexity of the idealization approach, $\JSMURF$ is the fastest and $\HILDE$ the slowest. A situation which is computationally particularly demanding is displayed in Figure \ref{fig:PorBwithAmpJSMURF} (see below). $\JSMURF$ detects almost no events. Due to internals in the dynamic programming algorithm, this causes a considerably long run time, in this example of roughly half an hour. We stress that $\HILDE$ uses $\JSMURF$ as an initial step and hence also $\HILDE$ is slow in such a situation, though it detects many events as it is able to resolve events on smaller temporal scales at and below the magnitude of the filter length.

\paragraph*{Interpreting, plotting and verification of the output}
All shown idealization methods return an object of the classes \textit{stepblock} and \textit{localDeconvolution}. We omit the exact structure of it (and refer to the man files of the called functions), but demonstrate important ways how to proceed with such an object. First of all, the idealization can be plotted using standard functions in \texttt{R}. Further, the convolution of the idealization with the kernel of the lowpass filter can be computed using the function \textit{getConvolution} in the package \textit{lowpassFilter}. The following code demonstrates how to do so. It provides the lower left panel in Figure \ref{fig:PorBHILDE}.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{fit} \hlkwb{<-} \hlstd{fit4}

\hlkwd{par}\hlstd{(}\hlkwc{mar} \hlstd{=} \hlkwd{c}\hlstd{(}\hlnum{4.5}\hlstd{,} \hlnum{4.5}\hlstd{,} \hlnum{0.5}\hlstd{,} \hlnum{0.5}\hlstd{))}
\hlkwd{plot}\hlstd{(time, data,} \hlkwc{pch} \hlstd{=} \hlnum{16}\hlstd{,} \hlkwc{col} \hlstd{=} \hlstr{"grey30"}\hlstd{,} \hlkwc{ylim} \hlstd{=} \hlkwd{c}\hlstd{(}\hlnum{2.5}\hlstd{,} \hlnum{4.2}\hlstd{),}
     \hlkwc{xlim} \hlstd{=} \hlkwd{c}\hlstd{(}\hlnum{0.0460798}\hlstd{,} \hlnum{0.0463768}\hlstd{),} \hlkwc{ylab} \hlstd{=} \hlstr{"Conductance in nS"}\hlstd{,}
     \hlkwc{xlab} \hlstd{=} \hlstr{"Time in s"}\hlstd{,} \hlkwc{cex.lab} \hlstd{=} \hlnum{1.8}\hlstd{,} \hlkwc{cex.axis} \hlstd{=} \hlnum{1.8}\hlstd{)}

\hlstd{ind} \hlkwb{<-} \hlkwd{seq}\hlstd{(}\hlnum{0.0460}\hlstd{,} \hlnum{0.0464}\hlstd{,} \hlnum{1e-6}\hlstd{)}
\hlstd{convolvedSignal} \hlkwb{<-} \hlstd{lowpassFilter}\hlopt{::}\hlkwd{getConvolution}\hlstd{(ind, fit, filter)}
\hlkwd{lines}\hlstd{(ind, convolvedSignal,} \hlkwc{col} \hlstd{=} \hlstr{"blue"}\hlstd{,} \hlkwc{lwd} \hlstd{=} \hlnum{3}\hlstd{)}

\hlkwd{lines}\hlstd{(fit,} \hlkwc{col} \hlstd{=} \hlstr{"#FF0000"}\hlstd{,} \hlkwc{lwd} \hlstd{=} \hlnum{3}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}

In Figure \ref{fig:PorBHILDE} we found that the convolution fits the recorded observations well, which is a confirmation for our idealization, but also for a correct specification of the model, in particular of the underlying lowpass filter. We always recommend such a graphical inspection to evaluate the quality of the idealization. If the idealization is not sufficiently good, one might modify tuning parameters (see the paragraph below), try a different idealization method (see Section \ref{sec:guidanceMethod}), remove artifacts or seek to improve the quality of the recordings.

Obtaining a model-free idealization is usually only one step in a data analysis. In Section \ref{sec:analysis} we discuss typical follow up steps. The idealized conductance values and the start and end times of the segments are given in \textit{fit\$values}, \textit{fit\$leftEnd} and \textit{fit\$rightEnd}, respectively. For instance, the following code creates histograms such as the ones in Figure \ref{fig:histograms}. We use the half sample mode \citep{robertson1974iterative}, implemented in the \textit{R}-package \textit{modeest}, to determine the underlying conductance levels, see the paragraph 'Analysis of the conductance levels' in Section \ref{sec:analysis} for further discussion.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlcom{# point amplitude histogram (histogram of the raw data)}
\hlkwd{hist}\hlstd{(data[data} \hlopt{>=} \hlnum{2.5} \hlopt{&} \hlstd{data} \hlopt{<=} \hlnum{4.5}\hlstd{],} \hlkwc{breaks} \hlstd{=} \hlkwd{seq}\hlstd{(}\hlnum{2.5}\hlstd{,} \hlnum{4.5}\hlstd{,} \hlnum{0.05}\hlstd{),}
     \hlkwc{main} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{xlab} \hlstd{=} \hlstr{"Conductance in nS"}\hlstd{)}
\hlkwd{abline}\hlstd{(}\hlkwc{v} \hlstd{= modeest}\hlopt{::}\hlkwd{mlv}\hlstd{(data[data} \hlopt{>=} \hlnum{2.5} \hlopt{&} \hlstd{data} \hlopt{<=} \hlnum{4.5}\hlstd{],}
                        \hlkwc{method} \hlstd{=} \hlstr{"hsm"}\hlstd{),} \hlkwc{col} \hlstd{=} \hlstr{"red"}\hlstd{,} \hlkwc{lwd} \hlstd{=} \hlnum{2}\hlstd{)}

\hlcom{# event histogram (histogram of the idealized conductance levels)}
\hlkwd{hist}\hlstd{(fit}\hlopt{$}\hlstd{value[fit}\hlopt{$}\hlstd{value} \hlopt{>=} \hlnum{2.5} \hlopt{&} \hlstd{fit}\hlopt{$}\hlstd{value} \hlopt{<=} \hlnum{4.5}\hlstd{],}
     \hlkwc{breaks} \hlstd{=} \hlkwd{seq}\hlstd{(}\hlnum{2.5}\hlstd{,} \hlnum{4.5}\hlstd{,} \hlnum{0.05}\hlstd{),}
     \hlkwc{main} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{xlab} \hlstd{=} \hlstr{"Codunctance in nS"}\hlstd{)}
\hlkwd{abline}\hlstd{(}\hlkwc{v} \hlstd{= modeest}\hlopt{::}\hlkwd{mlv}\hlstd{(fit}\hlopt{$}\hlstd{value[fit}\hlopt{$}\hlstd{value} \hlopt{>=} \hlnum{2.5} \hlopt{&} \hlstd{fit}\hlopt{$}\hlstd{value} \hlopt{<=} \hlnum{4.5}\hlstd{],}
                        \hlkwc{method} \hlstd{=} \hlstr{"hsm"}\hlstd{),} \hlkwc{col} \hlstd{=} \hlstr{"red"}\hlstd{,} \hlkwc{lwd} \hlstd{=} \hlnum{2}\hlstd{)}
\hlkwd{abline}\hlstd{(}\hlkwc{v} \hlstd{= modeest}\hlopt{::}\hlkwd{mlv}\hlstd{(fit}\hlopt{$}\hlstd{value[fit}\hlopt{$}\hlstd{value} \hlopt{>=} \hlnum{2.5} \hlopt{&} \hlstd{fit}\hlopt{$}\hlstd{value} \hlopt{<=} \hlnum{3.5}\hlstd{],}
                        \hlkwc{method} \hlstd{=} \hlstr{"hsm"}\hlstd{),} \hlkwc{col} \hlstd{=} \hlstr{"red"}\hlstd{,} \hlkwc{lwd} \hlstd{=} \hlnum{2}\hlstd{)}

\hlcom{# amplitude histogram (difference between idealized conductance levels)}
\hlstd{amp} \hlkwb{<-} \hlkwd{diff}\hlstd{(fit}\hlopt{$}\hlstd{value[fit}\hlopt{$}\hlstd{right} \hlopt{-} \hlstd{fit}\hlopt{$}\hlstd{left} \hlopt{<=} \hlnum{20} \hlopt{/} \hlstd{filter}\hlopt{$}\hlstd{sr} \hlopt{&}
                        \hlstd{fit}\hlopt{$}\hlstd{right} \hlopt{-} \hlstd{fit}\hlopt{$}\hlstd{left} \hlopt{>=} \hlnum{4} \hlopt{/} \hlstd{filter}\hlopt{$}\hlstd{sr])}

\hlkwd{hist}\hlstd{(amp[amp} \hlopt{<} \hlnum{2} \hlopt{&} \hlstd{amp} \hlopt{>} \hlnum{0}\hlstd{],} \hlkwc{breaks} \hlstd{=} \hlkwd{seq}\hlstd{(}\hlnum{0}\hlstd{,} \hlnum{2}\hlstd{,} \hlnum{0.05}\hlstd{),}
     \hlkwc{main} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{xlab} \hlstd{=} \hlstr{"Amplitude in nS"}\hlstd{)}
\hlkwd{abline}\hlstd{(}\hlkwc{v} \hlstd{= modeest}\hlopt{::}\hlkwd{mlv}\hlstd{(amp[amp} \hlopt{<} \hlnum{2} \hlopt{&} \hlstd{amp} \hlopt{>} \hlnum{0.7}\hlstd{],}
                        \hlkwc{method} \hlstd{=} \hlstr{"hsm"}\hlstd{),} \hlkwc{col} \hlstd{=} \hlstr{"red"}\hlstd{,} \hlkwc{lwd} \hlstd{=} \hlnum{2}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}


\paragraph*{Warning}
Users may experience a \textit{warning} saying \textit{``at least one segment could not be deconvolved since two successive short segments occurred''}. This is caused by the fact that the deconvolution approach incorporated in our methods can only deal with single changes or with isolated peaks (two changes in quick succession but separated by few more observations from other events). Obtaining a deconvolution for three or more  changes in quick succession is complicated and time-consuming and hence we decided to ignore such events when applying a deconvolution, but to mark them in \textit{attr(fit, "noDeconvolution")}. For a further analysis we usually recommend to ignore such events as they might even indicate artifacts. This can be done by setting all marked values to \textit{NA}.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{fit}\hlopt{$}\hlstd{value[}\hlkwd{attr}\hlstd{(fit,} \hlstr{"noDeconvolution"}\hlstd{)]} \hlkwb{<-} \hlnum{NA}
\hlstd{fit}\hlopt{$}\hlstd{rightEnd[}\hlkwd{attr}\hlstd{(fit,} \hlstr{"noDeconvolution"}\hlstd{)} \hlopt{-} \hlnum{1}\hlstd{]} \hlkwb{<-} \hlnum{NA}
\hlstd{fit}\hlopt{$}\hlstd{rightEnd[}\hlkwd{attr}\hlstd{(fit,} \hlstr{"noDeconvolution"}\hlstd{)]} \hlkwb{<-} \hlnum{NA}
\hlstd{fit}\hlopt{$}\hlstd{leftEnd[}\hlkwd{attr}\hlstd{(fit,} \hlstr{"noDeconvolution"}\hlstd{)]} \hlkwb{<-} \hlnum{NA}
\hlkwa{if} \hlstd{(}\hlkwd{length}\hlstd{(}\hlkwd{attr}\hlstd{(fit,} \hlstr{"noDeconvolution"}\hlstd{))} \hlopt{>} \hlnum{0}\hlstd{) \{}
  \hlkwa{if} \hlstd{(}\hlkwd{attr}\hlstd{(fit,} \hlstr{"noDeconvolution"}\hlstd{)[}
    \hlkwd{length}\hlstd{(}\hlkwd{attr}\hlstd{(fit,} \hlstr{"noDeconvolution"}\hlstd{))]} \hlopt{==} \hlkwd{length}\hlstd{(fit}\hlopt{$}\hlstd{leftEnd)) \{}
    \hlstd{fit}\hlopt{$}\hlstd{leftEnd[}\hlkwd{attr}\hlstd{(fit,} \hlstr{"noDeconvolution"}\hlstd{)[}
      \hlopt{-}\hlkwd{length}\hlstd{(}\hlkwd{attr}\hlstd{(fit,} \hlstr{"noDeconvolution"}\hlstd{))]} \hlopt{+} \hlnum{1}\hlstd{]} \hlkwb{<-} \hlnum{NA}
  \hlstd{\}} \hlkwa{else} \hlstd{\{}
    \hlstd{fit}\hlopt{$}\hlstd{leftEnd[}\hlkwd{attr}\hlstd{(fit,} \hlstr{"noDeconvolution"}\hlstd{)} \hlopt{+} \hlnum{1}\hlstd{]} \hlkwb{<-} \hlnum{NA}
  \hlstd{\}}
\hlstd{\}}
\end{alltt}
\end{kframe}
\end{knitrout}

%$

If too many segments are marked and they appear to be important for the given dataset, we cannot recommend to use our approaches, in this situation of \textit{extreme / high flickering} a better alternative might be approaches based on conductance distribution fitting, for further details see our review in Section \ref{sec:HMM}.

\paragraph*{Storing of the output}
To allow proceeding in a different program, one can store the idealization for instance in a \textit{csv} file as demonstrated by the following code. Note that we also remove the first and last segment, since their true start and end, respectively, cannot be identified by the data.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{fit} \hlkwb{<-} \hlkwd{data.frame}\hlstd{(}\hlkwc{left} \hlstd{= fit}\hlopt{$}\hlstd{leftEnd[}\hlopt{-}\hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,} \hlkwd{length}\hlstd{(fit}\hlopt{$}\hlstd{leftEnd))],}
                  \hlkwc{right} \hlstd{= fit}\hlopt{$}\hlstd{rightEnd[}\hlopt{-}\hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,} \hlkwd{length}\hlstd{(fit}\hlopt{$}\hlstd{rightEnd))],}
                  \hlkwc{value} \hlstd{= fit}\hlopt{$}\hlstd{value[}\hlopt{-}\hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,} \hlkwd{length}\hlstd{(fit}\hlopt{$}\hlstd{value))])}
\hlkwd{write.csv}\hlstd{(fit,} \hlkwc{file} \hlstd{=} \hlstr{"fit.csv"}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}


\paragraph*{Tuning parameters}
All three methods have multiple parameters which can be tuned to adapt to particular needs. Nonetheless, it is advisable to leave them unchanged unless specific reasons exists. All parameters are described in the man files of the called functions and in the referenced papers. Hence, in the following we will only give a brief overview about the most important ones. Further details are also provided in the review of our idealization approaches in Section \ref{sec:existingMethods}.

The choice of the number of repetitions of the Monte-Carlo simulations, the argument \textit{r}, was already discussed above in paragraph 'Run time and Monte-Carlo simulations'. The parameters \textit{alpha}, \textit{alpha1} and \textit{alpha2} are error levels $\alpha, \alpha_1, \alpha_2$ that bound approximately the probability of detecting one or more false positives (under the idealized scenario that the observations follow exactly the assumed model). As a default choice we suggest $\alpha = 0.05$. Larger values increase the chance to detect true events, but also to detect more false positives. One may use larger $\alpha$ values to 'screen' if important events are difficult to detect.

For $\HILDE$ the error level $\alpha := \alpha_1 + \alpha_2$ is split between the multiscale criterion of $\JSMURF$ (error level $\alpha_1$) and the local tests (error level $\alpha_2$). As default values, we suggest $\alpha_1=0.01$ and $\alpha_2=0.04$, since the focus of $\HILDE$ is typically on detecting short events primarily, while events on larger scales are often easier to detect. More weight can be put on $\alpha_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, for instance because of subconductance states. $\HILDE$ requires to specify the largest scale $l_{\max}$ (in \texttt{R}: lengths $= 1:l_{\max}$), this value should be chosen such that all events on larger scales are reliably detected by $\JSMURF$. If required, this can be tested by applying $\JSMURF$ or by Monte-Carlo simulations.

\subsection{Choosing the right method}\label{sec:guidanceMethod}
A guide which method to use is given in Figure \ref{fig:methodSelection}. The two main criteria are whether the noise is homogeneous or heterogeneous and whether short events are present and relevant. Recall that $\JSMURF$, $\JULES$ and $\HILDE$ are all suitable when one assumes homogeneous noise, but only $\JULES$ and $\HILDE$ allow for heterogeneous noise. Moreover, $\JULES$ and $\HILDE$ are designed to deal with short events, while $\JSMURF$ requires that events are slightly longer. Because of run time and precision, we generally recommend to use the simplest approach that is suitable for a dataset. Unless the dataset demands otherwise, we recommend $\JSMURF$ over $\JULES$ over $\HILDE$ and a homogeneous over a heterogeneous noise setting.

\paragraph*{Visual inspection}
\emph{Homogeneous noise} means that the noise distribution is the same at all times and for all conductance levels, otherwise the noise is called \emph{heterogeneous}. Heterogeneous noise is often clearly visible by naked eye, as in Figure \ref{fig:PorBhetData} where the noise level is higher for the higher conductance level, whence one should use either $\JSMURF$ or $\HILDE$ with heterogeneous noise setting. In most cases, if heterogeneous noise is not clearly visible, approaches that assume homogeneous noise are suitable.

A \emph{short event} is defined by two conductance changes in quick succession, e.g., a channel opening only very briefly before closing again. $\JSMURF$ should be used if it is not expected to miss relevant short events. Which events are too short depends not only on the absolute length, but also on the magnitude of the conductance change, noise levels, filtering and tuning parameters. At least, events shorter than filter length will certainly be missed by $\JSMURF$. Figure~\ref{fig:PorBData} shows such an example, whence one should use either $\JULES$ or $\HILDE$. 

\paragraph*{Empirical comparison}
If visual inspection is not sufficient, we suggest the following empirical procedure. The user should apply all potentially suitable methods to a small excerpt of the data and decide which leads to the best idealization. In general, if the idealizations are similar, the simpler approach should be preferred.

To illustrate the procedure for short events, consider Figure~\ref{fig:PorBHILDE} where we see that $\HILDE$ detects a large number of short events. In comparison, we see in Figure~\ref{fig:PorBwithAmpJSMURF} that $\JSMURF$ is not able to detect those events and hence is unsuitable for this dataset. In this case, $\HILDE$ appears to be more suitable. Contrarily, Figure~\ref{fig:GramicidinJSMURF} demonstrates that $\JSMURF$ is very suitable to idealize the Gramicidin dataset, where no short events occur, but events with small conductance changes, while the more complex method $\JULES$ struggles to detect all of them, as seen in Figure~\ref{fig:GramicidinJULES}.

To illustrate the procedure for heterogeneous noise, we idealized the observations in Figures~\ref{fig:PorBhetData}, which have visibly heterogeneous noise, with $\HILDE$, which is designed to deal with heterogeneous noise. Results are displayed in Figure~\ref{fig:PorBhetHILDE}. For comparison, an idealization by $\JULES$, assuming homogeneous noise, are displayed in Figure~\ref{fig:PorBhetJULES}. We see that $\JULES$ detects many additional events in the open state, which has higher noise level and while it is able to detect the short events, the fit is visibly worse than the fit by $\HILDE$.

To decide whether the noise is heterogeneous, we recommend to more advanced users also the following systematic approach: If longer segments without gating events are present, one can use them to estimate the noise level. Alternatively, one can idealize the data with $\JSMURF$ or $\HILDE$ with heterogeneous noise setting and use the idealization to determine noise levels as detailed in \citep[Section VI-C]{HILDE}.

Finally, if homogeneous noise is assumed and short events are relevant, we usually recommend to use $\JULES$ instead of $\HILDE$ as it is simpler and faster. Only if events are very short, such as in Figure \ref{fig:PorBHILDE}, $\HILDE$ should be used as it detects such events more likely.

%Good scientific practice demands reporting all results and reasons why which methods are chosen. One should never perform the whole analysis with two or more different approaches and only report the one that confirms the own hypotheses (performing multiple, but reporting and discussing all is of course absolutely fine).

\begin{figure}[!htb]
\centering
\includegraphics[width = 0.99\textwidth]{figures/PorBhetData.pdf}
\caption{From seconds to microseconds: Patchclamp recording (grey points) displayed at the 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 $\SI{20}{\milli\volt}$. The observations in the open state have visibly a larger noise level than the ones in closed state.}
\label{fig:PorBhetData}
\end{figure}
 

\begin{figure}[!htb]
\centering
\includegraphics[width = 0.99\textwidth]{figures/PorBhetHILDE.pdf}
\caption{Idealization (red) of the observations in Figure \ref{fig:PorBhetData} by $\HILDE$ \citep{HILDE} displayed on three different temporal scales. Lower panels: Convolution of the idealization with the lowpass filter (blue). Events are well idealized down to very short temporal resolutions.}
\label{fig:PorBhetHILDE}
\end{figure} 

\begin{figure}[!htb]
\centering
\includegraphics[width = 0.99\textwidth]{figures/PorBhetJULES.pdf}
\caption{Idealization (red) of the observations in Figure \ref{fig:PorBhetData} by $\JULES$ \citep{JULES} displayed on three different temporal scales. Lower panels: Convolution of the idealization with the lowpass filter (blue). $\JULES$ detects short events, but finds many small events, which are most likely false positives, at areas of high conductance and high variance (see for instance the idealization of the observations around \SI{0.36}{\nano\siemens} in the middle panel). These detections prevent $\JULES$ from performing a deconvolution at those positions (see for instance the lower left panel) and make the idealization unreliable.}
\label{fig:PorBhetJULES}
\end{figure}

\begin{figure}[!htb]
\centering
\includegraphics[width = 0.99\textwidth]{figures/PorBwithAmpJSMURF.pdf}
\caption{Idealization (red) of the observations in Figure \ref{fig:PorBhetData} by $\JSMURF$ \citep{JSMURF} displayed on three different temporal scales. Almost no events are detected, since all events are around or below the magnitude of the filter length.}
\label{fig:PorBwithAmpJSMURF}
\end{figure} 

\begin{figure}[!htb]
\centering
\includegraphics[width = 0.99\textwidth]{figures/GramicidinJULES.pdf}
\caption{Idealization (red) of the observations in Figure \ref{fig:GramicidinData} by $\JULES$ \citep{JULES} displayed on three different temporal scales. It misses many of the small conductance changes, since the fact that it also looks for very short events (which are not present in this dataset) slightly decreases detection power for long events compared to $\JSMURF$.}
\label{fig:GramicidinJULES}
\end{figure} 


\section{Models}\label{sec:models}
In this section we explain the statistical models underlying our methodology. For more details see \citep{JSMURF, JULES, HILDE}.

We assume that the recorded data $Y_1,\ldots,Y_n$ (the measured conductance at time points $t_i=i/\sr,\ i = 1,\ldots,n$, equidistantly sampled at rate $\sr$) result from a conductance $f$ perturbed by a centered Gaussian white noise process $\eta$. The noise is scaled by the noise level $\sigma$. Furthermore, conductance and noise are convolved with an analogue lowpass filter, with (truncated) kernel $\kernel_m$. Hence, after digitization at sampling rate $\sr = n / \tau_{\operatorname{end}}$, we obtain
\begin{equation}\label{eq:model}
Y_i = \big(\kernel_m \ast (f + \sigma \eta)\big)\left(i/\sr\right) = (\kernel_m \ast f)(i/\sr) + \epsilon_i,\quad i=1,\ldots,n,
\end{equation}
with $\ast$ the convolution operator. Here, $n$ denotes the total number of data points (typically several hundred thousands up to few millions). Hence, the resulting errors $\epsilon_1,\ldots,\epsilon_n$ are Gaussian and centered, but correlated (colored noise).\\
The conductance $f$ is assumed to be piecewise constant with potentially many different (unknown) segments of (unknown) length and size. The noise can either be homogeneous, i.e., the noise level $\sigma$ does not vary over time, or heterogeneous. In the latter case, we assume the noise level $\sigma$ to be an unknown piecewise constant function with potential jumps at the locations where the conductance changes, since changes of the noise level also depend on gating events\footnote{Strictly speaking, this models only \textit{heteroscedasticity}, one special form of heterogeneous noise that is for instance caused by open channel noise. But we expect our methods also to be robust to other forms of heterogeneous noise, confer the simulation results in \citep{HILDE}.}. More precisely, we model the conductance $f$ and the noise level $\sigma$ by
\begin{equation}\label{eq:signal}
\begin{array}{lll}
f(t) = \sum_{j=0}^{K}\,\valmu_j\,\ind_{[\tau_j, \tau_{j+1})}(t) & \text{ and } \sigma \equiv \sigma_0\in \R, & \text{if homogeneous noise is assumed},\\
f(t) = \sum_{j=0}^{K}\,\valmu_j\,\ind_{[\tau_j, \tau_{j+1})}(t) & \text{ and } \sigma(t) = \sum_{k=0}^{K}\,s_k\,\ind_{[\tau_k, \tau_{k+1})}(t), & \text{if heterogeneous noise is assumed},
\end{array} 
\end{equation}
where $t$ denotes physical time. The (unknown) conductance levels are denoted as $\valmu_0, \ldots, \valmu_K$, the (unknown) noise levels as $s_0,\ldots,s_K>0$, the (unknown) number of gating event as $K$ and the (unknown) locations of the gating events as $-\infty =: \tau_0 < \tau_1 < \cdots < \tau_K < \tau_{K+1}:=\tau_{\operatorname{end}}$. We stress that the class of signals in \eqref{eq:signal} is very flexible as potentially any arbitrary number of gating events at arbitrary conductance levels and arbitrary noise levels can be imposed, see Figure \ref{fig:PorBHILDE} for an example.

\section{Review of tools to analyze patchclamp recordings}\label{sec:idealizations}
In this section we give a review about methods for the analysis of patchclamp recordings. We start in Section \ref{sec:HMM} with a HMM-based analysis and discuss also their interplay with model-free idealizations as well as their advantages and disadvantages in comparison to model-free approaches. The analysis by and the interplay between the different approaches is also illustrated in Figure \ref{fig:interplay}. Secondly, we review existing model-free idealization methods in Section \ref{sec:existingMethods}. Finally, we discuss in Section \ref{sec:analysis} how idealizations can be used to analyze patchclamp recordings. Given the large amount of different methodology, we are by far not able to give a full review. The following only intends to summarize major ideas in order to help the reader to put $\JSMURF$, $\JULES$ and $\HILDE$ in the right context.

\subsection{Hidden Markov Models}\label{sec:HMM}
\paragraph*{HMM based analysis}
We limit our discussion mostly to homogeneous HMMs, which means that the parameters, which describe state transition properties and noise distribution, are constant in time. Inhomogeneous HMMs, see for instance \citep{diehn2019maximum}, are rarely used, as they are computationally more challenging and theoretical guarantees for parameter estimates are much harder to prove. As already discussed in the introduction, the assumption of a homogeneous Markov chain underlying the gating dynamics is almost always appropriate, but the assumption of a homogeneous error distribution to obtain a homogeneous HMM is more critical, since, e.g. because of artifacts, often intensive data cleaning or more complicated models are required. We stress that the quality of a HMM based analysis crucially depends on the stringent modeling assumption given by a HMM.

Obtaining an idealization by a HMM proceeds in several steps (illustrated in the right hand side of Figure \ref{fig:interplay}): First, a specific hidden Markov model has to be selected and ideally verified. This includes to find a Markov model for the gating dynamics, e.g., to fix the number of states and which transitions are possible. Note, that often multiple Markov states are required for one conductance level, e.g., to accommodate different noise levels or dwell times. Though data-driven model-selection tools are available, see e.g. \citep{gassiat2000likelihood, gassiat2003optimal, celeux2008selecting, chambaz2009minimum, lehericy2019consistent} and the references therein, this is often done manually by an empirical data analysis or by repeating the steps below until results are satisfying, which can be time-consuming and introduces subjectivity.

\begin{figure}[!htb]
  \centering
  \includegraphics[width = 0.99\textwidth]{figures/Interplay.pdf}
  \caption{Illustration of the interplay between HMM and model-free approaches.}
  \label{fig:interplay}
\end{figure}

As soon as a specific HMM is selected, parameters of the Markov model can either be estimated by the Baum-Welch algorithm, see \citep{Venkataramanan.etal.00, Qin.etal.00}, by Bayesian approaches, in particular MCMC sampling, see \citep{deGunst.etal.01, Siekmann.etal.11}, or by approaches based on the conductance (current) distribution, see \citep{yellen1984ionic, heinemann1991open, schroeder2015resolve} and the references therein.

Finally, an idealization can be obtained by the Viterbi algorithm \citep{viterbi1967error} or by Bayesian methods, in particular particle filtering, see \citep{fearnhead2018particle} and the references therein. Recently, a deep neural network approach has been proposed \citep{celik2020deep}, which skips the parameter estimation step and directly obtains an idealization. This approach can be seen as a hybrid method in between parametric and model-free approaches. It does not require a specific HMM to obtain an idealization but training in advance is required, which was done by assuming classes of hidden Markov models with hyperparameters.

Once an idealization is obtained, it can be used in reverse to estimate the parameters of the Markov model. We postpone details to Section \ref{sec:analysis}, since one proceeds as for model-free idealization. %However, one has to note that an idealization that is obtained using a wrong model or wrong parameters is not a reliable source to determine the model and the parameters. 
Using a layered HMM on simulated filtered signals \citep[Section IV-D]{HILDE} as well as in real data applications \citep[Section V]{HILDE}, \citep{Bartschetal18} we observed that the thus estimated parameters were significantly better than the parameters obtained directly by a Baum-Welch algorithm, most likely because of the applied missed event correction.

\paragraph*{Interplay}
As we will demonstrate in Section \ref{sec:analysis}, model-free idealizations allow a standalone analysis of patchclamp recordings. Moreover, they can assist a HMM based analysis in various forms (illustrated in Figure \ref{fig:interplay}): Model-free idealization can help to identify and remove artifacts, we used for instance $\JULES$ in \citep{Bartschetal18} to assists a HMM based analysis in that way. They can be used to determine the number of conductance levels (paragraph 'Analysis of the conductance levels' in Section \ref{sec:analysis}) and help select and verify a specific Markov model (paragraph 'Selection and verification of a Markov model' in Section \ref{sec:analysis}). Furthermore, most HMM based parameter estimation approaches are iterative procedures which require starting values. Those are particularly crucial when the procedure converges to a local optimum only. Such starting values can be provided by previously obtained values using model-free idealizations. Finally, model-free idealizations and the resulting parameter estimates using a missed event correction can be used to verify HMM based idealization and parameter estimates, and vice versa. This is particularly valuable as they have different strengths and weaknesses as outlined in the following paragraph.

We also note that the local deconvolution approach used in our model-free idealization methods, see \citep{JULES}, can be used to improve HMM based idealizations, obtained for instance by a Viterbi algorithm, as our approach not only takes into account explicitly the filtering but is also time-continuous. It only relies on a prior fit that fixes the number of conductance changes and their rough locations. It can be called by the function \textit{deconvolveLocally} in the package \textit{clampSeg}.

\paragraph*{HMM versus model-free idealization: Compared and contrasted}
In general, HMM based approaches achieve a higher temporal resolution of gating dynamics, because of their stronger assumptions. Hence, parameter estimates might be more accurate as they rely on more detected events. Moreover, HMMs allow for immediate parameter estimation and interpretation, which is often the main goal of an analysis. And since the HMM state space is fixed in advance, the idealization immediately assigns every time point to one of the states. In contrast, model-free idealizations often have to be postprocessed, (e.g., by clustering or thresholding) to identify discrete states because conductance levels are determined freely.

On the other hand, there are several disadvantages, some of which are closely entangled with the advantages. As discussed above, the need for often extensive preprocessing adds subjectivity and also more potential sources of data analysis errors. In contrast, in such situations model-free methods may right away provide a reasonable idealization as they can potentially handle inhomogeneity in a more flexible way, in particular those which act locally on the dataset. The state space and a model for the noise must be fixed in advance, thereby strongly limiting the possible results. Model selection always has a subjective component and can lead to a flawed idealization, for example by inadvertently modeling two states with similar but subtly different conductance or noise levels as only one state, or by prescribing an unsuitable noise model which can lead to detection of spurious state changes. Within a HMM framework one can only incompletely determine whether the data is compatible with the underlying model assumptions. Hence, despite the above described advantages, at least in simulations and real data examples in \citep{HILDE, Bartschetal18} we observed that parameter estimates based on an idealization (either obtained by model-free approaches or by the Viterbi algorithm) appear to be more accurate than direct estimates by the Baum-Welch algorithm. Gating dynamics are time-continuous processes, but for simplification many HMM approaches underlie a time-discrete Markov chain as an approximation. A time-discrete approximation is also implied by most model-free approaches as they allow gating events only at the sampling points. An exception is the local deconvolution approach used in \citep{JULES} and \citep{HILDE}.

Some of the subjectivity and other problems in HMM modeling can be mitigated by conducting a model-free idealization to inform preprocessing and model selection. In summary, HMM-based and model-free approaches can (and should) be used to complement each other to verify each other's results. Artifacts and missed events might be reasons for some differences, but otherwise results should be similar.

\subsection{Review of existing model-free idealization approaches}\label{sec:existingMethods}
Many analyses are still performed by visual inspection, often with manually chosen event times or in a semi-automatic way, for instance by amplitude thresholding \citep{Colquhoun.87, Sakmann.Neher.95}, as e.g. offered by pCLAMP 10 software (Molecular Devices), or by the semi-automatic \textit{SCAN} software \citep{colquhoun1995fitting} which allows \textit{time-course fitting}. Hence, those approaches are typically time-consuming and subjective. Moreover, approaches which are based on additional filtering (often by lowpass Gauss filters) aggravate detection of small events. A first approach for a fully automatic idealization was slope thresholding \citep{Basseville.Benveniste.83}, for instance $\TRANSIT$ \citep{VanDongen.96}. Recently, \citet{gnanasambandam2017unsupervised} proposed idealizations based on the minimal description length ($\MDL$). All of them (except the semi-automatic \textit{SCAN} software) ignore lowpass filtering and hence may have difficulties to idealize events correctly on small temporal scales. Furthermore, if events are present on multiple scales (recall Figures \ref{fig:PorBData}, \ref{fig:GramicidinData}, \ref{fig:PorBhetData}), uniscale thresholding procedures will usually fail. 

As mentioned in the introduction, $\JSMURF$ \citep{JSMURF}, $\JULES$ \citep{JULES} and $\HILDE$ \citep{HILDE} are multiscale procedures combined with local deconvolution and hence take into account both issues. Consequently, they provide usually more accurate results as demonstrated in simulations and real data applications. As described in the following they mostly differ in how they take into account the filter when detecting events and hence whether they are suitable to detect short events, but also whether they incorporate the possibility to allow for heterogeneous noise. To understand the methodology better, it is illustrative to plot the convolution of a single gating event and single peaks with the kernel of a lowpass filter, see Figure \ref{fig:convolution}. We stress that for the short event displayed in Figure \ref{fig:convolutionpeak} the filtered signal does not reach the lower conductance level of the original signal. This is generally the case for peaks shorter than the filter length $m/\sr$. Hence, if such short events are present, deconvolution techniques are indispensable to idealize those conductance levels correctly. 

\begin{figure}[htb]
  \begin{subfigure}[b]{0.33\linewidth}
    \centering
    \includegraphics[width=0.9\linewidth]{figures/Jump.png} 
    \caption{Single gating event.} 
    \label{fig:convolutionjump} 
    \vspace{4ex}
  \end{subfigure}%% 
  \begin{subfigure}[b]{0.33\linewidth}
    \centering
    \includegraphics[width=0.9\linewidth]{figures/Peak.png} 
    \caption{Short peak.} 
    \label{fig:convolutionpeak} 
    \vspace{4ex}
  \end{subfigure}%%
    \begin{subfigure}[b]{0.33\linewidth}
    \centering
    \includegraphics[width=0.9\linewidth]{figures/Long.png} 
    \caption{Long peak.} 
    \label{fig:convolutionlong} 
    \vspace{4ex}
  \end{subfigure}  
  \vspace{-.5cm}
\caption{\footnotesize Signals (black line) containing a single gating event, a short peak and a longer peak and their convolutions (blue line) with a four-pole lowpass Bessel filter with normalized cutoff frequency of $0.1$ and sampling rate $10^4$. Vertical red lines indicate the event time plus the filter length $m / \sr$.}
  \label{fig:convolution} 
\end{figure}

\paragraph*{$\JSMURF$}
The \textbf{J}ump-\textbf{S}egmentation by \textbf{MU}lti\textbf{R}esolution \textbf{F}ilter, $\JSMURF$, from \citet{JSMURF}, combines a multiscale criterion with rigorous error control to reliably detect events on various temporal scales simultaneously. More precisely, it takes into account all scales above the filter length and for each of those intervals it ignores the first $m$ data-points. As illustrated in Figure \ref{fig:convolution} only during these $m$ point long transitions the convolution is not matching the conductance $f$. It provides the following strict error control. The probability that the idealization contains at least one false positive event (an event that is not contained in the true conductance $f$) is bounded by the error level $\alpha$. The original work \citep{JSMURF} assumed homogeneous noise, \citep{HILDE} proposed an extension to heterogeneous noise.

\paragraph*{$\JULES$}
The \textbf{JU}mp \textbf{L}ocal d\textbf{E}convolution \textbf{S}egmentation filter, $\JULES$, from \citet{JULES}, applies a multiscale criterion to all temporal scales and combines it with a postfilter step to remove incremental steps as for instance occurring in Figure \ref{fig:PorBSMUCE}. Finally, a local deconvolution approach is proposed to idealize short events well. The error level $\alpha$ bounds the probability of detecting a false positive approximately. All in all, $\JULES$ is particularly designed for homogeneous noise and short events.

\paragraph*{$\HILDE$}
\textbf{H}eterogeneous \textbf{I}dealization by \textbf{L}ocal testing and \textbf{DE}convolution \citep{HILDE} obtains idealizations in three steps: It applies $\JSMURF$ to detect events on large temporal scales, afterwards it tests locally for additional short events. Those tests explicitly take filtering into account. The final idealization is once again obtained by local deconvolution. Local tests are performed on scales up to length $l_{\max}$. The error level $\alpha := \alpha_1 + \alpha_2$ is split between the multiscale criterion of $\JSMURF$ (error level $\alpha_1$) and the local tests (error level $\alpha_2$). False positives occur again only with probability approximately $\alpha$.

%$\HILDE$ was particularly designed to deal with heterogeneous noise but can also be used when homogeneous noise is assumed. It is, to the best of our knowledge, together with the extension of $\JSMURF$ which was proposed in \citep{HILDE} as well, the first model-free idealization method that admits heterogeneous noise explicitly and hence deals efficiently with it. Additionally, it is the first method that takes into account filtering explicitly when detecting (short) events. Therefore, $\HILDE$ is particularly suited to idealize events that are short in time, both when the noise is homogeneous or heterogeneous.

\paragraph*{Simulation results}
In the following we give a brief qualitative summary about the simulation results in \citep{JSMURF, Pein17, JULES, HILDE}. Generally speaking, such computer simulations are a systematic but also computation intensive way to determine precisely how long an event has to be such that an idealization method is able to reliably detect it. However, we stress that all quantitative results depend on the signal to noise ratio, the filter and on tuning parameters.

We found that $\JSMURF$ reliably idealizes events of medium or large length (usually an event has to be at least few times the filter length) even when the conductance change is small, confer \citep{JSMURF}. This is essential to idealize \textit{subgating} events. In comparison, $\JULES$ and $\HILDE$ are able to reliably idealize much shorter peaks, if they are isolated. Isolated means that two events have to be separated by at least three times the filter length if homogeneous noise is assumed but at least five times the filter length if heterogeneous noise is assumed (for a filter truncated after $m=11$ sampling points). Moreover, for a good idealization events have to be usually only few sampling points long but can be shorter than the filter length. Hence, those two approaches are suitable to idealize \textit{flickering}. $\HILDE$ allows events to be a bit shorter than $\JULES$.

$\JSMURF$ is usually the fastest of our three approaches and an idealization of several hundred thousands up to a few million data points take often only seconds (when Monte-Carlo simulations are already performed). In comparison, an idealization of the same dataset with $\JULES$ may last around a minute and with $\HILDE$ few minutes. All run times are measured on a standard laptop and increase typically linearly in the number of data points. A notable exception are situations in which $\JSMURF$ detects almost no change-points, then the run time increases quadratically in the number of observations. For instance the idealization in Figure \ref{fig:PorBwithAmpJSMURF} took roughly half an hour. Since $\HILDE$ uses $\JSMURF$ as a first step, its run time is similarly slow.

\subsection{Analysis of patchclamp recordings}\label{sec:analysis}
In this section we provide a step-by-step guide on how to analyze patchclamp recordings using model-free idealizations. In addition we describe their interplay with Markov model based analyses, see also the introduction and in particular Figure \ref{fig:interplay} for an illustration. Of course, any analysis depends on the specific datasets and its goals. Hence, the following steps should be seen as more of general guidance that has to be interpreted flexibly. We also stress that it contains time consuming verification steps which might not be necessary in every analysis.

\paragraph*{Analysis of the conductance levels}
For this step, we assume that the underlying protein attains only a finite number of conformations and hence that only a finite number of conductance levels occur. We aim to determine this number, the values of the conductance levels and possible transitions between those levels. This can be done in various ways and we will only sketch important ideas. Event histograms (histograms of the idealized conductance levels) and amplitude histograms (histogram of the differences between consecutive segments in the idealization) should be used as a visualization of the underlying conductance levels, see Figure \ref{fig:histograms}.

The idealized conductance levels form a mixture distribution, typically a Gaussian mixture, around the true conductance levels, where randomness results from measurement and idealization errors (and hence the peaks are narrower if those are better performed). Modes correspond to the true conductance levels. An example can be found in Figure \ref{fig:histograms}. One can use simple approaches based on a Gaussian assumption to estimate modes. We obtained good results using the half sample mode \citep{robertson1974iterative}, because it is quite robust against outliers. In more difficult cases, where peaks cannot be identified that clearly, more involved statistical methodology to estimate the components of a mixture distribution has to be used, for an overview see \citep{mclachlan2004finite} and the references therein, or the accuracy of the measurement or idealization has to be increased.

\begin{figure}[htb]
  \begin{subfigure}[b]{0.33\linewidth}
    \centering
    \includegraphics[width=0.9\linewidth]{figures/HistogramRaw_PorBamp.png} 
    \caption{Point amplitude histogram.} 
    \label{fig:pointHistogram} 
    \vspace{4ex}
  \end{subfigure}%% 
  \begin{subfigure}[b]{0.33\linewidth}
    \centering
    \includegraphics[width=0.9\linewidth]{figures/HistogramFit_PorBamp.png} 
    \caption{Event histogram.} 
    \label{fig:eventHistogram} 
    \vspace{4ex}
  \end{subfigure}%%
    \begin{subfigure}[b]{0.33\linewidth}
    \centering
    \includegraphics[width=0.9\linewidth]{figures/Amplitude_PorBamp.png} 
    \caption{Amplitudes.} 
    \label{fig:amplitudes} 
    \vspace{4ex}
  \end{subfigure}  
  \vspace{-.5cm}
\caption{\footnotesize Histograms of the PorB measurement with ampicillin in Figure \ref{fig:PorBData}. Histograms are based on the visualized and on ten additional traces. Code to obtain such histograms was explained in the paragraph 'Interpreting, plotting and verification of the output' in Section \ref{sec:softwareUse}. In the point amplitude histogram we found one dominant conductance level of $3.9664$ nS (estimated by the half sample mode). Smaller conductance levels are not visual, since they are too short and smoothed by the lowpass filter (in total $2\,476$ data points are between $2.5$ nS and $3.5$ nS). The event histogram confirms this conductance level. Note that the peak is much narrower in the event than in the point amplitude histogram. This is usually the case and improves identification of conductance levels. Moreover, because of the deconvolution step in our idealization and since dwell times are not represented in the event histogram, we were also able to identify a second conductance level of $2.7956$ nS, i.e., the amplitude (difference, blockage effect of the ampicillin) is $1.1708$ nS. The amplitude histogram confirms this finding with a pronounced mode at $1.1662$ nS. A simple mean is roughly the same with $1.1546$ nS. The amplitude histogram but also the event histogram shows further events. Those events could be matched to processes unrelated to the interaction of ProB and ampicillin, and hence should be ignored for the ampicillin influence, confer \citep{Bartschetal18}.}
  \label{fig:histograms} 
\end{figure}

Afterwards, one often wants to map each idealized conductance level to its mixture component, for instance by defining thresholds based on such histograms. Note that the idealization methods are often sensitive enough to detect baseline fluctuations and fluctuations due to pink noise as events. As a result, often several events in a row are mapped to the same conductance level. In other words, this process can also merge segments and thus remove spurious events. Moreover, it is often a good idea to remove segments that are far from any estimated conductance level from the subsequent analysis, since they typically result from artifacts.

\paragraph*{Selection and verification of a Markov model}
As discussed before, a time-continuous Markov model is a common assumption to analyze patchclamp recordings. Since model-free idealizations are obtained without any prior assumption on the gating dynamics, they can be used to determine and verify a Markov model. In order to avoid statistical dependency, a careful analysis involves splitting the measurements and using the first part to select a Markov model and the second part to verify the model. Recall that a Markov model has two key properties: dwell times (how long a channel stays in one Markov state) are independent of each other and are exponentially distributed. Since it is often simpler, one might aim to verify uncorrelated, instead of independent, dwell times, though lack of correlation does not imply independence. When checking whether dwell times follow a Markov model, one has to take into account that short events might be missed. Nonetheless, at least in simple Markov models with only few states one readily can check for uncorrelated and exponentially distributed dwell times, for an example see \citep[Figure S4 in the supplement]{Bartschetal18}.

\paragraph*{Parameter estimation}
Once a specific Markov model is assumed, one has to estimate its parameters. To this end, it is essential to take into account missed events. Missing events shorter than a certain resolution limit are widely discussed in the literature. The exact distribution is calculated by \citet{hawkes.etal.90}, an estimator called $\operatorname{MIL}$ of the Q-matrix is suggested by \citet{qin1996estimating} and integrated in the $\operatorname{QuB}$ software package \citep{nicolai2013solving}, the exact maximum likelihood estimator for the Q-matrix for two conductance levels is obtained by \citet{colquhoun1996joint} and recently a Bayesian approach was proposed by \citet{epstein2016bayesian}. In \citep{JULES, HILDE, Bartschetal18} we applied simpler approximations, which worked well since the measurements could be modeled well by Markov models with only two or three states.

\paragraph*{Verification by using hidden Markov approaches}
This step was already discussed in Section \ref{sec:HMM}. The previous analysis using model-free idealizations can be an essential help to perform an analysis using HMM based approaches. HMM based approaches are however potentially able to achieve better temporal resolution. Hence, both approaches should be used to verify each other's results, both in terms of parameter estimation and of idealizations.

\section{Discussion}\label{sec:discussion}
We gave detailed guidance on how to obtain model-free idealizations using $\JSMURF$, $\JULES$ and $\HILDE$ and on how to use those idealizations together with HMM based approaches to analyze patchclamp recordings. We believe that this provides a rather comprehensive toolkit for the analysis of many patchclamp recordings.

A notable exception are experiments with varying conductance. Such experiments are interesting, since not only the present value of voltage affects the channel, some channels are also affected by the present rate of voltage change. This includes channels that show no gating when the voltage is constant, but can be activated by a varying voltage. For other channels different dynamics are observed when the voltage changes. One example is the protein channel Tim23 which tends to close when larger voltage levels are applied constantly \citep{denkert2017cation}. Moreover, experiments with a constant voltage only allow to examine the gating dynamics at few voltage levels (or require large experimental effort), while with varying voltage the dynamics can be analyzed for a whole range of voltages by a single experiment. Brief ideas were discussed in \citep[Section 6.2.1]{Pein17} and \citep{diehn2019maximum}. 

%Another interesting research idea could be the development of model-free idealization approaches that allow only a finite number of conductance levels. They would be an intermediate approach between model-free and HMM and hence share some benefits and disadvantages with both of them. In comparison to traditional model-free approaches, they will benefit from the additional assumptions by having a better sensitivity to detect small events (either short in time or with a small conductance change) and by providing more concise results. But they also share with HMM based approaches the disadvantage that they are heavily affected by artifacts. Detection of changes with the restriction of a finite level set (alphabet) was studied in \citep{behr2018multiscale}, but their application to detect copy-number aberrations from genetic sequencing data does not require filtering and hence an application to patchclamp recordings is not straightforward.

Though model-free approaches are in general more robust to artifacts than HMM based approaches, confer \citep{JULES, HILDE} who demonstrated for $\JULES$ and $\HILDE$ certain robustness to model violations, there is need for improved methodology (either model-free or HMM based) with a larger focus on robustness.

%Finally, in this paper we focused primarily on obtaining good idealizations, and hence could only discuss typical analysis steps applied afterwards. It would be of interest to study those steps more systematically and to aim to automatize them at least to some extent. We also think that good software for those steps is lacking in general. This includes the estimation of the number of conductance states, a missed event analysis to estimate parameters of a Markov model, but also model verification.


\bibliographystyle{apalike}
%\bibliographystyle{acm}
%\bibliographystyle{IEEEtran}
\bibliography{Literature}

\end{document}