Here we treat an exemplary ligand-gated channel with two ligand binding steps and one open-closed isomerization described by an HMM (see Figure 1a). For this model, confocal patch-clamp fluorometry (cPCF) data were simulated: time courses of ligand binding and channel current upon concentration jumps were generated (see Appendix 5 and Materials and methods section). Idealized example data with added white noise are shown in Figure 1b–d. We added realistic instrumental noise to the simulated data (see Appendix 5). A qualitative description of the statistical problem that needs to be addressed when modeling time series data such as the simulated is outlined in Box. 1.

Figure 1

Download asset Open asset

Figure 1—source data 1

The example data is provided.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig1-data1-v3.zip

Download elife-62714-fig1-data1-v3.zip

Here and in the Materials and methods section, we derive the mathematical tools to account correctly for the stochastic Markov dynamics of single molecules in the fluctuations of macroscopic signals. The KF is a Bayesian filter (see Materials and methods), that is a continuous state HMM with a multivariate normal transition probability Ghahramani, 1997 (Figure 2a). We define the hidden ensemble state vector

which counts the number of channels in each state ${\displaystyle \mathbf{s}}$ (see Methods). To make use of the KF, we assume the following general form of the dynamic model: The evolution of ${\displaystyle \mathbf{n}(t)}$ is determined by a linear model that is parametrized by the state evolution matrix ${\displaystyle \mathrm{T}}$

where ∼ means sampled from and ${\displaystyle \mathcal{N}(\cdot |\mathit{\mu },\mathbf{\Sigma })}$ is a shorthand for the multivariate normal distribution, with the mean μ and the variance-covariance matrix ${\displaystyle \mathbf{\Sigma }}$. The state evolution matrix (transition matrix) is related to the rate matrix ${\displaystyle \mathbf{K}}$ by the matrix exponential ${\displaystyle \mathbf{T}=\exp (\mathbf{K}\mathrm{\Delta }t)}$. The mean of the hidden state evolves according to the equation ${\displaystyle \mathbb{E}[{\mathbf{n}}_{t+1}|{\mathbf{n}}_t]=\mathbf{T}{\mathbf{n}}_t}$. It is perturbed by normally distributed white process noise ${\displaystyle \mathit{\omega }}$ with the following properties: The mean value of the noise fulfills ${\displaystyle \mathbb{E}[{\mathit{\omega }}_t]=0}$ and the variance-covariance matrix of the noise is ${\displaystyle \mathrm{cov}[{\mathit{\omega }}_t,{\mathit{\omega }}_t]=\mathbf{Q}(\mathbf{T},{\mathbf{n}}_t)}$ (see Materials and methods Equation 38d, Ball, 2016). In short, Equation 3 defines a Gaussian Markov process.

Figure 2

Download asset Open asset

The observations ${\displaystyle {\mathbf{y}}_t}$ depend linearly on the hidden state ${\displaystyle {\mathbf{n}}_t}$. The linear map is determined by an observation matrix ${\displaystyle \mathbf{H}}$.

The noise of the measurement setup (Appendix 5 and Equation 43) is modeled as a random perturbation of the mean observation vector. The noise fulfills ${\displaystyle \mathbb{E}[\mathit{\nu }]=0}$ and ${\displaystyle \mathrm{cov}[{\mathit{\nu }}_t,{\mathit{\nu }}_t]={\mathbf{\Sigma }}_t}$. Equation 4 defines the state-conditioned observation distribution ${\displaystyle \mathbb{O}}$ (Figure 2a). The set of all measurements up to time ${\displaystyle t}$ is defined by ${\displaystyle {\mathcal{Y}}_t=\left\{y_1,\dots ,y_t\right\}}$. If the system strictly obeys Equation 3 and Equation 4 then the KF is optimal in the sense that it is the minimum variance filter of that system Anderson and Moore, 2012. If the distributions of ν and ω are not normal, the KF is still the minimum variance filter in the class of all linear filters but there might be better non-linear filters. In case of colored noise ν and ω the filtering equations (see Materials and methods) can be reformulated by state augmentation or measurement-time-difference approach techniques Chang, 2014. For each element in a sequence of hidden states ${\displaystyle \{{\mathbf{n}}_t:0<t<T\}}$ and for a fixed set of parameters ${\displaystyle \mathit{\theta }}$, an algorithm based on a Bayesian filter (Figure 2a), explicitly exploits the conditional dependencies of the assumed Markov process. A Bayesian filter recursively predicts prior distributions for the next ${\displaystyle {\mathbf{n}}_t}$

given what is known about ${\displaystyle {\mathbf{n}}_{t-1}}$ due to ${\displaystyle y_{t-1}}$. The KF as a special Bayesian filter assumes that the transition probability is multivariate normal according to Equation 3

Note, that Equation 6 is a central approximation of the KF. While the exact transition distribution of an ensemble of ion channels is the generalized-multinomial distribution (Methods Equation 32), the quality of normal approximations to multinomial Milescu et al., 2005 or generalized-multinomial Moffatt, 2007 distributions depends on the number of ion channels $N_{\mathrm{ch}}$ in the patch and on the position of the probability vector in the simplex space. The difference between the log-likelihoods of the true generalized-multinomial dynamics and Equation 6 type approximation scales as ${\displaystyle 1/N_{\mathrm{c}\mathrm{h}}}$ Moffatt, 2007. As a rule of thumb one should be careful with both algorithms for time traces with ${\displaystyle N_{\mathrm{c}\mathrm{h}}\in \left[{10}^1,{10}^2\right]}$. Below or even inside this interval there are more qualified concepts such as the forward algorithm or even particle filters (Golightly and Wilkinson, 2011; Gillespie and Golightly, 2012) which avoid the normal approximation.

Each prediction of ${\displaystyle {\mathbf{n}}_t}$ (Equation 6) is followed by a correction step,

that allows to incorporate the current data point into the estimate, based on the Bayesian theorem (Chen, 2003). Additionally, the KF assumes (Anderson and Moore, 2012; Moffatt, 2007) a multivariate normal observation distribution

If the initial prior distribution is multivariate normal then due to the mathematical properties of the normal distributions the prior and posterior $\mathbb{P}(\cdot )$ in Equation 8 become multivariate normal Chen, 2003 for each time step. In this case, one can derive algebraic equations for the prediction (Materials and methods Equation 37, Equation 38d) and correction (Materials and methods Equation 58 and Equation 58) of the mean and covariance. The algebraic equations originate from the fact that a normal prior is the conjugated prior for the mean value of a normal likelihood. Due to the recursiveness of its equations, the KF has a time complexity that is linear in the number of data points, allowing a fast algorithm. The denominator of Equation 8 is the normal distributed marginal likelihood ${\displaystyle \mathbb{L}({\mathbf{y}}_t|{\mathcal{Y}}_{t-1},\mathit{\theta })}$ for each data point, which constructs by

a product marginal likelihood of normal distributions of the whole time trace ${\displaystyle {\mathcal{Y}}_{\mathrm{T}}=\{{\mathbf{y}}_1,\dots ,{\mathbf{y}}_{N_{\mathrm{T}}}\}}$ of length $N_{\mathrm{T}}$ for the KF. For the derivation of ${\displaystyle {\mathbf{P}}_t}$ and ${\displaystyle {\mathbf{\Sigma }}_t}$ see Materials and methods (Equation 38d) and Equation 43. ${\displaystyle {\mathbf{P}}_t}$ is the covariance of the prior distribution over ${\displaystyle \mathbf{n}(t)}$ before the KF took ${\displaystyle \mathbf{y}(t)}$ into account. The likelihood for the data allows to ascribe a probability to the parameters ${\displaystyle \mathit{\theta }}$, given the observed data (Methods Equation 20). An illustration for the operation of the KF on the observation space is given in Figure 2b. The predicted mean signal ${\displaystyle \mathbf{H}\mathbb{E}[\mathbf{n}(t)]}$ corresponds to binding degree ${\displaystyle B(t)=\frac{\mathbf{H}\mathbb{E}[\mathbf{n}(t){]}_1}{N_{\mathrm{c}\mathrm{h}}}}$ and open probability ${\displaystyle P_O(t)=\frac{\mathbf{H}\mathbb{E}[\mathbf{n}(t){]}_2}{N_{\mathrm{c}\mathrm{h}}}}$. These values are plotted as vector trajectories.

The standard KF (Moffatt, 2007; Anderson and Moore, 2012; Chen, 2003) has additive constant noise ${\displaystyle {\mathrm{\Sigma }}_t=const}$ in the observation model. Thus, in this case a constant variance term ${\displaystyle \mathrm{\Sigma }}$ is added, in Equation 9 to the aleatory variance ${\displaystyle \mathbf{H}{\mathbf{P}}_t{\mathbf{H}}^{\mathrm{\top }}}$ which, as mentioned above, originates (Equation 38d) from the the fact that we do not know the true system state ${\displaystyle \mathbf{n}(t)}$. For signals with Poisson-distributed photon counting or open-channel noise, we need to generalize the noise model to account for additional white-noise fluctuations with ${\displaystyle \mathbf{n}(t)}$-dependent variance. For instance, in single-channel currents additional noise is often observed whose variance is referred to by ${\displaystyle {\sigma }_{\mathrm{o}\mathrm{p}}^2}$. In macroscopic currents this additional noise can be modeled by a term ${\displaystyle {\sigma }_{\mathrm{o}\mathrm{p}}^2{n}_4(t)}$, causing state-dependency of our noise model.

The second noise term ${\nu }_{\mathrm{op}}$ is defined in terms of the first two moments ${\displaystyle \mathbb{E}({\nu }_{\mathrm{o}\mathrm{p}})=0}$ and ${\displaystyle \mathrm{v}\mathrm{a}\mathrm{r}({\nu }_{\mathrm{o}\mathrm{p}})=\mathbb{E}({\nu }_{\mathrm{o}\mathrm{p}}^2)={\sigma }_{\mathrm{o}\mathrm{p}}^2{n}_4(t)}$. To the best of our knowledge such a state-dependent noise makes the integration of the denominator of Equation 8 (which is also the incremental likelihood) intractable

This is because the state distribution ${\displaystyle \mathcal{N}(\mathbf{n}|\overline{\mathbf{n}}(t),\mathbf{P}(t))}$ as the prior also influences the variance parameter of the likelihood which means that the conjugacy property is lost. While a normal distribution is the conjugated prior of the mean of a normal likelihood, it is not the conjugated prior for the variance. However, by applying the theorem of total variance decomposition Equation 46a we deduce a normal approximation to Equation 8 and to the related problem of Poisson-distributed noise in fluorescence Equation 57, Equation 55a data. By computing the mean and the variance or covariance matrix of the signal, we can reformulate the noise model to fit the form of the traditional KF framework. Note, that the derived equations for the covariance matrix are still exact for the more general noise model. Mean and covariance just do not form a set of sufficient statistics anymore.

Our derivation is not limited to ligand-gated ion channels. For example, when investigating voltage-gated channels, the corresponding noise model can be easily adapted. This holds also when using the P/n protocol for which the noise model resembles that of the additional variance in the fluorescence signal. The additional variance is induced because the mean signal from the ligands swimming in the bulk (Materials amd methods Equation 43 Appendix 5) is eliminated by subtracting scaled mean reference signal which itself has an error. This manipulation adds additional variance to the resulting signal comparable to P/n protocol. Other experimental challenges, as for example series resistance compensation promoting oscillatory behavior of the amplifier, deserve certainly advanced treatment. Nevertheless, for voltage-clamp experiments with a rate equation approach it also becomes clear (Lei et al., 2020) that modeling of the actual experimental limitations, including series resistance, membrane and pipette capacitance, voltage offsets, imperfect compensations by the amplifier, and leak currents are necessary for consistent kinetic scheme inference.

The Bayesian posterior distribution

encodes all information from model assumptions and experimental data used during model training (see Materials and methods). A full Bayesian inference is usually not an optimization (finding the global maximum or mode of the posterior or likelihood) but calculates all sorts of quantities derived from the posterior distribution such as mean values of any function $f$ including the mean value or covariance matrix of the parameters themselves or even the likelihood of the data.

Besides the covariance matrix of the parameter to express parameter uncertainty, the posterior allows to calculate a credibility volume. The smallest volume ${\displaystyle V_P}$ that encloses a probability mass ${\displaystyle P}$ of

is called the Highest Density Credibility Volume/Interval (HDCV/HDCI). Those credibility volumes should not be confused with confidence volumes although under certain conditions they can become equivalent. Given that our model sufficiently captures the true process, the true values ${\displaystyle {\mathit{\theta }}_{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}}}$ will be inside that volume with a probability ${\displaystyle P}$. Unfortunately, typically there is no analytical solution to Equation 12 . However, it can be solved numerically with Monte Carlo techniques, enabling to calculate all quantities related to Equation 13 and Equation 14 . Our algorithm uses automatic differentiation of the statistical model to sample from the posterior (Appendix 1—figure 1a) via Hamiltonian Monte Carlo (HMC) (Betancourt, 2017), see Appendix 7 , as provided by the Stan software (Hoffman and Gelman, 2014; Gelman et al., 2015).

We compare the posterior distribution (Figure 3) of our algorithm against Bayesian versions of the deterministic (Milescu et al., 2005) and stochastic (Moffatt, 2007) algorithms, which we consider as the gold standard algorithms for macroscopic patch-clamp data. Simulated currents of a patch with ${\displaystyle N_{ch}=5\cdot {10}^3}$ are shown in (Figure 3d). The resulting posteriors (Figure 3a) show that both former algorithms are further away from the true parameter values with their maxima or mean values (Figure 3a). E.g., the relative error of the maximum of the posterior are ${\displaystyle \mathrm{\Delta }{k}_{21}\approx 200\mathrm{\%}}$ for Milescu et al., 2005 and ${\displaystyle \mathrm{\Delta }{k}_{32}\approx 240\mathrm{\%}}$ for Moffatt, 2007 . The four other parameters including the three equilibrium constants behave less problematic as judged by their relative error. Additionally, if one does not only judge the performance by the relative distance of maximum (or some other significant point) of the posterior but considers the spread of the posterior as well, it becomes apparent, that the marginal posterior of both former algorithms fail to cover the true values within at least the reasonable parts of their tails. Accordingly, for maximum likelihood inferences the true value would be far outside the estimated confidence interval. For the RE approach only the marginal posterior of ${\displaystyle {\tilde{K}}_{21}}$ is nicely centered over the true values and the marginal of ${\displaystyle {\tilde{K}}_{32}}$ could be considered to cover within a reasonable part of the distribution the true value. Uncertainty quantification is investigated in more detail further down (Figures 4—9). Note that in Figure 3a, parameter unidentifiability by heavy tails/ multiple maxima of the posterior distribution or (anti-) correlation is easily visible as non axial symmetric patterns.

Figure 3

Download asset Open asset

Figure 3—source data 1

The data folder includes all 15 sets of time traces.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig3-data1-v3.zip

Download elife-62714-fig3-data1-v3.zip

Figure 4

Download asset Open asset

Figure 4—source data 1

6µMol.

Each of the source data folders contains for a specific ligand concentration the time traces of cPCF data for all $N_{\mathrm{ch}}$.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig4-data1-v3.zip

Download elife-62714-fig4-data1-v3.zip

Figure 4—source data 2

64µMol.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig4-data2-v3.zip

Download elife-62714-fig4-data2-v3.zip

Figure 4—source data 3

32µMol.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig4-data3-v3.zip

Download elife-62714-fig4-data3-v3.zip

Figure 4—source data 4

8µMol.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig4-data4-v3.zip

Download elife-62714-fig4-data4-v3.zip

Figure 4—source data 5

4µMol.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig4-data5-v3.zip

Download elife-62714-fig4-data5-v3.zip

Figure 4—source data 6

1µMol.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig4-data6-v3.zip

Download elife-62714-fig4-data6-v3.zip

Figure 4—source data 7

2µMol.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig4-data7-v3.zip

Download elife-62714-fig4-data7-v3.zip

Figure 4—source data 8

025µMol.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig4-data8-v3.zip

Download elife-62714-fig4-data8-v3.zip

Figure 4—source data 9

05µMol.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig4-data9-v3.zip

Download elife-62714-fig4-data9-v3.zip

Figure 4—source data 10

00625µMol.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig4-data10-v3.zip

Download elife-62714-fig4-data10-v3.zip

Figure 4—source data 11

003125µMol.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig4-data11-v3.zip

Download elife-62714-fig4-data11-v3.zip

Figure 4—source data 12

0125µMol.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig4-data12-v3.zip

Download elife-62714-fig4-data12-v3.zip

Figure 5

Download asset Open asset

Figure 6

Download asset Open asset

Figure 7

Download asset Open asset

Figure 7—source data 1

The folder of the five-state model includes 15 sets of time traces in the interval $N_{\mathrm{ch}}\in [{10}^3,5\cdot {10}^5]$.

Each of them has 10 ligand concentrations. The number in the file name reports the amount of ion channels in the patch.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig7-data1-v3.zip

Download elife-62714-fig7-data1-v3.zip

Figure 7—source data 2

The folder of the six-state model includes 14 sets of time traces in the interval $N_{\mathrm{ch}}\in [{10}^3,{10}^6]$.

Each of them has 10 ligand concentrations. The number in the file name reports the amount of ion channels in the patch.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig7-data2-v3.zip

Download elife-62714-fig7-data2-v3.zip

Figure 8

Download asset Open asset

Figure 9

Download asset Open asset

To assess the location of the posterior conditioned on $N_{\mathrm{ch}}$, we select the median vector $\mathit{\theta }$ of the marginal posteriors and calculate its Euclidean distance to the true values by:

This defines a single value to judge the overall accuracy of the posterior. Varying ${\sigma }_{\text{op}}/i$ reveals the range of the validity (Figure 3b) of the algorithm (red) from Moffatt, 2007 . While both stochastic approaches are nearly equivalent for low open-channel noise, the RE (blue) performs consistently poorer. It may seem surprising that even for ${\displaystyle {\sigma }_{op}/i<0.01}$ the two stochastic algorithms start to produce different results. But considering the scaling (Materials and methods Equation 46a) of the total open-channel noise (top axis) from currents of an ensemble patch $\propto {(N_{\mathrm{ch}}{P}_{\mathrm{open},\max }0.5)}^{0.5}{\sigma }_{\mathrm{open}}$ one sees that if $\propto {(N_{\mathrm{ch}}{P}_{\mathrm{open},\max }0.5)}^{0.5}{\sigma }_{\mathrm{open}}$ approaches σ the traditional KF suffers from ignoring state dependent noise contributions. The lower scale changes with experiments (e.g. $N_{\mathrm{ch}}$ and ${\sigma }_{op}$). In contrast, the upper scale is largely independent of the particular measurements. The two different normalizations indicate an experimental intuition: “ Why should I consider the extra noise from the open state of the single channel if only ${\sigma }_{\mathrm{op}}/i={\sigma }_{\mathrm{op}}/\sigma \approx 0.01$” is misleading. The small advantage of our algorithm for small ${\sigma }_{\mathrm{op}}/i$ over Moffatt, 2007 is due to the fact that we could apply an informative prior in the formulation of the inference problem on ${\sigma }_{\exp }\sim \mathrm{normal}({\sigma }_{\exp ,\mathrm{true}}^2,{\sigma }_{\exp ,\mathrm{true}}^2\cdot 0.01)$ by taking advantage of our generalization (Equation 46a) Bayesian filter. Further, Figure 3b indicates the importance that the functional form of the likelihood is flexible enough to capture the second order statistics of the noise of the data sufficiently.

For an increasing data quality, which in our benchmark is an increasing $N_{\mathrm{ch}}$ per trace, we show (Figure 3c) that the deterministic RE and our Bayesian filter are consistent estimators, that is they converge in distribution to the true parameter values with their posterior maxima or median for increasing data quality. The scaling of the RE approach (blue) and our Bayesian filter (green) vs. $N_{\mathrm{ch}}$ shows that for large $N_{\mathrm{ch}}$ both algorithms seem to have a constant error ratio relative to each other. They are both well described by $\mathrm{error}(N_{\mathrm{ch}})\propto a/\sqrt{N_{\mathrm{ch}}}$ with an error ratio computed from the fit of 4.4. Thus, although our statistical model is singular (meaning that the fisher information matrix is singular Watanabe, 2007), its asymptotic learning behaviour is similar to a regular model (Figure 4c) which, however, means that the euclidean error from both algorithms stays different also for large $N_{\mathrm{ch}}$. For data with ${\displaystyle N_{\mathrm{c}\mathrm{h}}<2\cdot {10}^3}$ the samples from the posterior typically indicate that the posterior is improper which is defined as

We consider this as the case of unidentified parameters. This data-driven definition is in so far different from structural and practical identifiability definitions (Middendorf and Aldrich, 2017a; Middendorf and Aldrich, 2017b) as the two latter cases are not distinguished. Still the practical consequence of structural or practical unidentifiability, which is usually an improper posterior, is captured. Cases of structural or practical unidentifiability which lead to a confined region of constant posterior density will be considered identified as the posterior is still normalizable thus the uncertainty quantification will still be correct, even when this finding is not sufficient to answer the research question at hand.

For the simulated time traces with an optimistically high signal-to-noise assumption, the posterior of the KF (from hereon KF denotes our Bayesian Filter) and a RE (Milescu et al., 2005) approach are compared for cPCF data (Figure 4a–d). For a brief introduction of the RE approach, see Appendix 8 . The failure to analyze PC data with moderate open-channel noise (Moffatt, 2007; Figure 3a) disqualifies the classical KF with its constant noise variance also as a useful algorithm for fluorescence data, because here the Poisson distribution of the signal generates an even stronger state dependency of the signal variance.

By “high signal-to-noise assumption” , we refer to an experimental situation with a standard deviation of the current recordings ${\sigma }_{\mathrm{ex}}/i=0.5$, a low additional ${\sigma }_{\mathrm{op}}/i=0.05$, and a high mean photon rate per bound ligand and frame ${\lambda }_{\mathrm{b}}=5$. Additionally, we assume vanishing fluorescence background noise generated by the ligands in the bulk. The benefit of the high signal-to-noise is that the limitations of the two different approximations to the stochastic process of binding and gating can be investigated without running into the risk of being compensated or obscured by the noise from the experimental setup. For these experimental settings (Figure 4a), we calculate the Euclidean distance of the median (Equation 15) for different $N_{\mathrm{ch}}$. For ${\displaystyle N_{\mathrm{c}\mathrm{h}}<500}$ (gray shaded area in Figure 4a), the Euclidean error of both algorithms is roughly the same. On the single parameter level (Figure 4b), this can be seen as an onset of correlated deviations from the true value for both algorithms. Each marginal posterior has for each $N_{\mathrm{ch}}$ a similar deviation in magnitude and direction. That is in particular true for ${\tilde{k}}_{32}$ and ${\tilde{K}}_{32}$ which dominate Equation 15 . In spite of the correlation in direction of the errors of ${\tilde{k}}_{21}$ and ${\tilde{K}}_{21}$ their magnitude is still smaller for the KF. In summary, this indicates that in this regime the approximations to the involved multinomial distributions fail in a similar manner for both algorithms. That implies that treating the autocorrelation of the gating and binding becomes similar important compared to the error induced by normal approximations (which are used by the KF and the RE approach). For larger $N_{\mathrm{ch}}$, the Euclidean error of the RE is on average 1.6 times larger than the corresponding error of the posterior mode of the KF, which we deduce by fitting the function $\mathrm{error}(N_{\mathrm{ch}})=\frac{a}{\sqrt{N_{\mathrm{ch}}}}$. On the one hand, both algorithms are better in approaching the true values than with patch-clamp data alone. On the other hand, the smaller error ratio means, that adding a second observable constrains the posterior, such that much of the overfitting is prevented for the RE approach. By overfitting, we define the adaptation of any inference algorithm to the specific details of the used data set due to experimental and intrinsic noise which is aggravated if too complex kinetic schemes are used. Similarly, (Milescu et al., 2005) showed that the over fitting tendency of the RE can be reduced if the autocorrelation of the data is eliminated. The dotted green curve derives from PC data. The Euclidean error is roughly an order of magnitude larger for ${\displaystyle N_{\mathrm{c}\mathrm{h}}>2000}$. Thus, in this regime the cPCF data set is equivalent to 10² fold more time traces or 10² more $N_{\mathrm{ch}}$ in a similar PC data set. For ${\displaystyle N_{\mathrm{c}\mathrm{h}}<2000}$ only cPCF establishes parameter identifiability (given a data set of 10 ligand concentrations and no other prior information). In Figure 4b(1-6), we demonstrate the 0.95-HDCI (Equation 14) of all parameters and their modes vs. $N_{\mathrm{ch}}$. Even though the Bayesian filter and the RE approach are both consistent estimators, the RE approach covers the true values with its 0.95-HDCI only occasionally. The modeling assumption of the RE approach of treating each data point as if it does not come from a Markov process but from an individual draw from a multinomial distribution with deterministically evolving mean and variance makes the parameter estimates overly confident (Figure 4b(1-6)) . A likely explanation can be found by analyzing the extreme case where data points are sampled at high frequency relative to the time scales of the channel dynamics. The RE approach treats each data point as a new draw from Equation 67 while in reality the ion channel ensemble had no time to evolve into a new state. In contrast, the KF updates its information about the ensemble state after incorporating the current data point and then predicts from this updated information the generalised multinomial distribution of the next data point. For ${\displaystyle N_{\mathrm{c}\mathrm{h}}>200}$, the marginal posterior of the KF usually contains the true value. Nevertheless, one might depict a bias in both algorithms, in particular (Figure 4b 2,4) for ${\tilde{k}}_{32}$ and ${\tilde{K}}_{32}$ for ${\displaystyle N_{\text{ch}}<2\cdot {10}^3}$, similar to the findings of Moffatt, 2007 . A proper investigation of bias can be found in Figure 11 and 12 and in the Appendix. Notably, with the more realistic higher experimental noise level, in those tests the bias is hardly or not all detectable (consider the unfiltered or infinitely fast integrated data). A plausible explanation is that the bias only occurs (Figure 4 2,4) because the data are that perfect that the discrete nature of the ensemble dynamics is almost visually detectable, thus deviating from to the modeling assumption of multi-variate normal distributions.

To investigate the six one-dimensional 0.95-HDCIs simultaneously, we declare the analysis of a data set as successful if all 0.95-HDCIs include the true values. Otherwise we define it as a failure. This enables to determine the true probability at which the probability mass of the KF and the RE approach covers the true values in a binomial setting. The left blue vertical line in Figure 4c indicates $p={0.95}^6\approx 0.735$ which is the lower limit and which would be the true success probability for an ideal model whose six 0.95-HDCIs are drawn from $y\sim \mathrm{binomial}(0.95,6)$. This is the probability of getting 6 successes in 6 trials. The right blue vertical line equals $p=0.95$, signifying the upper limit obtained by treating the six $0.95-$ HDCIs as being drawn from $y\sim \mathrm{binomial}(0.95,1)$ each, which is a rather loose approximation. All marginal distributions are computed from the same high-dimensional posterior which is formed by one data set for each trial. Thus, the six $0.95-$ HDCIs $y\sim \mathrm{binomial}(0.95,1)$ must have success rates between those two extremes if the algorithm creates an accurate posterior. We next combine the binomial likelihood with the conjugated beta prior (Hines et al., 2014) for mathematical convenience. On this occasion, for the sake of the argument, $\mathrm{beta}(1,1)$ seems sufficient. A $\mathrm{beta}(1,1)$ prior is a uniform prior on the open interval $(0,1)$. The estimated true success rate of the RE approach (blue) is $\approx 0.15$ and therefore far away from the success probability an algorithm should have when it is based on an exact likelihood of the data. In contrast, the posterior (green) of the true success probability of the KF resides with a large probability mass between the lower and upper limit of the success probability of an optimal algorithm (given the correct kinetic scheme). As both algorithms use the same prior distribution, the different performance is not induced by the prior.

Exploiting six one-dimensional posterior distributions does not necessarily answer whether the posterior is accurate in 6 dimensions but we can refine the used binomial setting. In Figure 4d $\mathbb{P}({\tilde{k}}_{32},{\tilde{K}}_{43})$, we see that 2-D marginal distributions can, due to their additional degree of freedom, twist around the true value without covering it with HDCV (Equation 14) of reasonable size while simultaneously the two 1–D marginal distribution do cover it with a reasonable HDCI. In general, the KF posterior distribution has its mode much closer to the true value for various parameter combinations and it seems that the posterior is approximately multivariate normal. Further, we recognize that the probability mass of the reasonably sized HDCV of the KF posterior includes the true values whereas the HDCV from the RE does not. In 6 dimensions we lack visual representations of the posterior. Since we showed that both algorithms are consistent for a given identifiable model, we are looking for a way to ask whether the posterior is accurate (has the posterior distribution the right shape). We can answer that question by asking, how much probability mass around the mode (or around multiple modes) needs to be counted to construct a HDCV Equation 14 which includes the true values. Then we can ask for $N_{\mathrm{set}}$ data sets how often did we find the true values inside a volume $V(P)$ of a specific probability mass $P$ of the posterior distribution

An algorithm which estimates the parameters of the true process should fulfill this property simultaneously to being consistent. Otherwise credibility volumes or confidence volumes are meaningless. Noteworthy, that this is a empirical test of how sufficient the Bayesian filter and the RE approach hold frequentist coverage property of their HDCVs (Rubin and Schenker, 1986). We explain (Appendix 8) in detail how to quantify the overall shape and $n$-dimensional posterior and comment on its geometrical meaning. One way is to use an analytical approximation via the cumulative Chi-squared distribution (Figure 5a and b), The other way is to count the probability mass of $n$-dimensional histogram bins starting with the highest value until the first bin includes the true values (Figure 5b).

Knowing how much volume/probability mass is needed to include the true rate matrix allows us to test whether all HDCVs constructed from the two probability distributions match the binomial distributions of the ideal model. For each data set and for each HDCV of a fixed probability mass, there are two possible outcomes: The true rate matrix is inside or outside of that volume. For a chosen HDCV with a fixed probability volume, as indicated by a gray space in Figure 5b , we count how many times the true matrix is included in the volume of that probability mass for each trail in a fixed amount of trials. Since the success is binomially distributed, we plot the expected mean of a perfect model $E[y]=N_{\mathrm{trials}}{P}_{\mathrm{true}}$ and binomial quantiles and compare them with the success rate found in our test runs (Figure 5c) for both algorithms with both methods to determine the posterior shape. The posterior of the KF distributes the probability mass in a consistent manner such that each volume includes the true rate matrix within the quantile range. In contrast, the RE approach fails for all data sets for all HDCVs (from 0 – 0.95 probability mass) and does not include the true values in one single case. Note, that all the binomial trials for each HDCV are made from the same set of data sets which explains the correlated deviation from the mean. For lower but realistic signal to noise ratios, where the fit quality decreases, for example by producing larger errors/wider posterior distributions (Figure 6a), the statistics of the HDCV from the RE approach improve but are still outperformed by the KF. In particular, in our tested case of realistic experimental noise we never find the true values within a 0.65-HDCV if the data are analyzed with a RE approach. Even for the highest noise level (Figure 6b), the probability mass of the KF posterior needed to include the true rate matrix remains almost always smaller then the posterior mass of the RE approach. That means that the posterior mass of the KF is much closer to the true value distributed than the posterior mass of the RE. With the KF we find the true rate matrix for one data set in small volume ${\displaystyle P<0.05}$ around the mode. To achieve the same with the RE approach we need at least a probability mass of 0.3.

In the inset of Figure 6a and b we do not observe a trend, thus no indication that the RE approach has a better performance for large values $N_{\mathrm{ch}}$ in this regard. This challenges the common argument that the RE approach should be equivalent to the KF for large $N_{\mathrm{ch}}$ because the ratio of mean signal vs. the intrinsic binding and gating noise is so large. Thus, including the autocorrelation into the analysis is important even for unrealistic large $N_{\mathrm{ch}}$. One possible explanation is model a signal-to-noise ratio which scales $\propto N_{\mathrm{ch}}$. From the multinomial distribution both algorithms inherit mean signals which scale $\propto N_{\mathrm{ch}}$ and variances which scale in the terms dominating for large $N_{\mathrm{ch}}$ similarly with $\propto N_{\mathrm{ch}}$. Thus, identical to the real signal, both algorithms model the scaling of the signal-to-noise ratio $\propto \sqrt{N}$. It is plausible, that both algorithms remain sensitive for the occurrence of autocorrelation of the noise even for largest signal-to-noise ratios. In Figure 5c we compare the Euclidean error of the pessimistic high white noise case with an over-optimistic low noise case. We see, that when increasing $N_{\mathrm{ch}}$ there is a regime where the Euclidean error increases faster than ${\sqrt{N_{\mathrm{ch}}}}^{-1}$ which we indicate with a coarse approximate fit $\propto N_{\mathrm{ch}}^{-1}$. In that regime two effects happen simultaneously. First, the mean and the intrinsic fluctuations of the signal become more and more dominant over the experimental noise. Second, the standard deviation of intrinsic fluctuations becomes smaller relative to the mean signal. We speculate, that this produces together a learning rate which is faster than the usual asymptotic learning rate ${\sqrt{N_{\mathrm{ch}}}}^{-1}$ of a regular model but relaxes asymptotically towards ${\sqrt{N_{\mathrm{ch}}}}^{-1}$.

We have seen in Figure 3c and Figure 4a that the RE and the KF algorithm are consistent estimators, while their error ratio (Figure 7a) seems to have no trend to approach 1 with increasing $N_{\mathrm{ch}}$. Adding a second observable increases parameter accuracy and adds identifiability for both algorithms since less aspects of the dynamics need to be statistically inferred (Figure 4a). Furthermore, the second observable takes away much of the tendency (compare Figure 4b 1 – 6 with 8) of the RE approach to overinterpret (overfit) which leads to a shrinking of the error ratio $5.6\pm 1.4$ for PC data to smaller values for cPCF data (Figure 7a) (red) which are on average still bigger than one, while the Euclidean error is reduced (Figure 4a). If we then keep the amount and quality of the PC/cPCF data but increase the complexity of the model which produced the data (Figure 7b and d) from a four-state to a five-state model (see kinetic schemes above Figure 7a–c), we see that for cPCF data the error ratio stays roughly the same (difference between Figure 7a and b). For PC data instead both algorithms deliver an unidentified k₂₁ for $N_{\mathrm{ch}}\leqq 2\cdot {10}^3$ (defined as an improper posterior). For larger $N_{\mathrm{ch}}$ the KF always identifies all parameters while the RE fails at $N_{\mathrm{ch}}\in \{7000,2000,75000\}$ to identify k₅₄. Thus, the KF reduces the risk of unidentified parameters. To calculate the mean error ratio, we exclude the values were some of the parameters are unidentified in total that still amounts to $6.8\pm 2.7$ thus the advantage of the KF (given all parameters are identified) might increase with model complexity for PC data. The lowest Euclidean error for this kinetic scheme has cPCF data analyzed with the KF. (Figure 7d). A 6-state-1-open-states model with cPCF data has again an error ratio of the the usual scale (Figure 7c). As expected, the Euclidean error continuously increases with model complexity (Figure 7d and e). For PC data of the 6-state-1-open-states model even the likelihood of the KF is that weak (Figure 7e) that it delivers unidentified parameters even for and we can detect heavy tailed distributions up until . Using RE on PC data alone does not lead to parameter identification, thus no error ratio can be calculated.

Consistent with our findings, fluorescence data itself, should lower the advantage of the KF compared to PC data simply by signal-to-noise arguments. The stochastic aspect of the ligand binding is usually more dominated by the noise of Photon counting and background noise than the stochastic gating is dominated in current data by experimental noise. In terms of uncertainty quantification the advantage of the KF with cPCF varies with the model complexity (see, Appendix 9).

Besides analyzing what causes the changes in the Euclidean error (Figure 7a and b) at the single parameter, we now investigate whether the posterior is a proper representation of uncertainty. Thus, we look back at the HDCIs. The HDCIs of the 4-state-1-open-state (Figure 8) of the PC data from Figure 3 reveal an exacerbated over-confidence problem of the RE approach (blue) compared to cPCF-data (Figures 4b1—6). This, underlines our conclusion of Figures 5 and 6 that the Bayesian posterior sampled by the RE approach is misshaped. As a consequence a confidence volume derived from the curvature at the ML estimate of the RE algorithm understates parameter uncertainty. A possible way for ML methods to derive correct uncertainty quantification is by using bootstrapping data methods (Joshi et al., 2006). Furthermore, the error ratios of each single parameter from its true value in the last column ${\tilde{k}}_{43}$${\tilde{K}}_{43}$ strongly increased their magnitudes (insets Figure 8). Even error ratios of $5\cdot {10}^2$ are possible. Note, that the way we defined Equation 15 suppresses the influence of the smaller parameter errors in the overall error ratio. Thus the advantage (error ratio) of the KF over RE approach for a single parameter can be much larger or lower compared to the error ratio derived from the Euclidean error if the respective parameter is contributing less to the Euclidean error. The posterior of the KF (green) seems to be unbiased after the transition into the regime ${\displaystyle N_{\mathrm{c}\mathrm{h}}>2\cdot {10}^3}$ where all parameters are identified. Similarly, for the RE algorithm there is no obvious bias in the inference. If we use the RE algorithm and change from the four-state to the five-state model (PC data from Figure 7b), bias occurs (Figure 9) in many inferred parameters, even for the highest $N_{\mathrm{ch}}$ investigated. Milescu et al., 2005 showed that one or the reason of the biased inference of the RE approach is its ignorance of autocorrelation of the intrinsic noise. We add here that the bias problem clearly aggravates with an increased model complexity. It is even present in unrealistically large patches which in principle could be generated by summing up 10² time traces with $N_{\mathrm{ch}}={10}^3$. In contrast, the KF algorithm reveals that its parameter inference is either unbiased or at least much less biased in the displayed $N_{\mathrm{ch}}$ regime. Furthermore, for both algorithms the position of the HDCI relative to the true value is for some parameters highly correlated, which corresponds to the correlation between optima of the ML method of Milescu et al., 2005 ; Moffatt, 2007.

As a side note, unbiased parameter estimates are a highly desirable feature of an inference algorithm. For example, with a bias in the inference, repeated experiments do not lead to the true value if the arithmetic mean of the parameter inferences is taken. With bias even bootstrapping methods fail to produce reliable uncertainty quantification. Due to the variation of the data the k₅₄ parameter is either identified in some neighbourhood of the true value or complete unidentified (Figure 9), if the RE algorithm is used. The unidentified k₅₄ occurs even at high-quality data such as $N_{\mathrm{ch}}=7.5\cdot {10}^4$. Only because of the nonphysical prior (Figure 9) of k₅₄ induced by the limits of the sampling box of the sampling algorithm the posterior appears to be proper but is in fact either unidentified or or more than two orders of magnitude away from the true value. For the same data using the KF did not result in any unidentified parameters. Note, that comparable inference pathologies such as multimodal distributions of inferred parameter were also reported for the maximum likelihood RE algorithm for low quality PC data or too simple stimulation protocols (Milescu et al., 2005).

In conclusion, the two different perspectives on parameter uncertainty: On the one hand distributions of ML estimates due to the random data (Milescu et al., 2005 ; Moffatt, 2007) and the Bayesian posterior distribution loose their tightly linked (and necessary) connection if the RE algorithm is used. Thus, KF robustifies also ML inferences of the rate matrix. Our findings are consistent with the findings for gene regulatory networks (Gillespie and Golightly, 2012) which show that RE approaches deliver a too narrow posterior in contrast to stochastic approximations which deliver an acceptable posterior compared to the true posterior (defined by a particle filter algorithm). On the data side of the inference problem adding cPCF data eliminates the bias, reduces the variance of the position of the HDCI and eliminates unidentified parameters (Appendix 9—figures 1 and 2) for both investigated algorithms. This advantage increases with model-complexity.

For the five-state and six-state model, we applied microscopic-reversibility (Colquhoun et al., 2004). We enforced it by hierarchical prior distribution (Materials and methods Equation 60) whose parameters can be chosen such that they allow only arbitrarily small violations of microscopic-reversibility. But the prior distribution can also be used to enforce some softer regularization around microscopic-reversibility. Thus, we can transfer the usually strictly applied algebraic constraint (Salari et al., 2018) of microscopic-reversibility to a constraint with scalable softness. In that way we can model the lack of information if microscopic-reversibility is exactly fulfilled (Colquhoun et al., 2004) by the given ion channel instead of enforcing the strict constraint upon the model.

In the Bayesian framework, the likelihood of the data and the prior generate the posterior. Thus, the performance of both algorithms can be influenced by appropriate prior distributions. We used a uniform prior over the rate matrix which is not optimal. Note, that uniform priors are widely used by several reasons. They appear to be unbiased, and are assumed to be a ‘no prior’ option (which they are not). This is true for location parameters like mean values. In contrast, for other parameters, such as scaling parameters like rates or variances, a uniform prior adds bias to the inference towards faster rates (Zwickl and Holder, 2004). We suspect, that for the PC data even in the simplest model discussed here the lower data quality limit below which we detected unidentified parameters (improper posteriors) is caused by the uniform prior. This lower limit for the KF also increases with the complexity of the model from ${\displaystyle N_{\mathrm{c}\mathrm{h}}<2\cdot {10}^3}$ for the foue-state model till $N_{\mathrm{ch}}\leqq 2\cdot {10}^4$ for 6-state-1-open-state model. Note, that it is hardly possible to fit the 6-state-1-open-state model with the RE approach for the same amount of PC data. We observe cPCF data eases this problem because the likelihood becomes more concentrated for all parameters. The likelihood dominates the uniform prior. Nevertheless, for most parts of the paper we used a uniform prior over the rates and equilibrium constants to be comparable with the usual default method: a plain ML which influences our results in data regimes in which the data is not strong enough to dominate the bias from the uniform prior. Thus, both algorithms perform better with smarter informative or at least unbiased prior choices for the rate matrix.

In principle, to rule out an influence of the prior, unbiased priors should be used for the rates. The standard concept for unbiased least informative priors is to construct a Jeffreys prior Jeffreys, 1946 for the rate matrix which is, however, beyond the scope of the paper.

To evaluate the advantage of cPCF data Biskup et al., 2007 with respect to PC data only (Figure 10), we compare different types of ligands: Idealized ligands with brightness ${\lambda }_{\text{b}}$, emitting light only when bound to the channels, ‘real’ ligands which also produce background fluorescence when diffusing in the bath solution (Appendix 5) and current data alone. For datasets including fluorescence, the increased precision for the dissociation rate of the first ligand, $k_{2,1}$, is that strong that the variance of the posterior $\mathbb{P}(k_{2,1},k_{3,2})$ nearly vanishes in the combined plot with the current data (nearly all probability mass is concentrated in a single point in Figure 10a). The effect on the error of the equilibrium constants $K_i$ is less strong. Additionally, the bias is reduced and even the estimation of $N_{\mathrm{ch}}$ is improved. The brighter the ligands are, the more the posterior of the rates decorrelates, in particular $\mathbb{P}(k_{2,1},k_{3,2})$ (Figure 10a). All median estimates of nine different cPCF data sets (Figure 10b) differ by less than a factor 1.1 from the true parameter except $k_{3,2}$, which does not profit as much from the fluorescence data as $k_{2,1}$ (Figure 10c). The 95th percentiles, l₉₅ of $\mathbb{P}(k_{2,1})$ and $\mathbb{P}(K_1)$ follow $l_{95}\sim 1/\sqrt{{\lambda }_{\mathrm{b}}}$. Thus, with increasing magnitude of ligand brightness λ, the estimation of $k_{2,1}$ becomes increasingly better compared to that of $k_{3,2}$ (Figure 10c). The posterior of the binding and unbinding rates of the first ligand contracts with increasing ${\lambda }_{\mathrm{b}}$. The l₉₅ percentiles of other parameters exhibit a weaker dependency on the brightness ($l_{95}\sim {\lambda }^{-0.1}$). For ${\lambda }_{\mathrm{b}}=0.01$ photons per bound ligand and frame, which corresponds to a maximum mean signal of 20 photons per frame, the normal approximation to the Poisson noise hardly captures the asymmetry of photon counting noise included in the time traces. Nevertheless, l₉₅ decreases about ten times when cPCF data are used (Figure 10c). The estimated variance of

with the mean predicted signal $\mathbf{H}E[\mathbf{n}(t_i)]$, for PC or cPCF data is ${\sigma }^2(r_i)\approx 1$ (Figure 10d) which means that the modeling predicts the stochastic process correctly up to the variance of the signal. Note that the mean value and covariance of the signal and the state form sufficient statistics of the process, since all involved distributions are approximately multivariate normal. The fat tails and skewness of $\mathbb{P}(k_{21})$ and $\mathbb{P}(k_{12})$ arises because the true model is too flexible for current data without further prior information. The KF allows to determine the variance (Figure 10e) of the open-channel current noise for ${\sigma }_{\mathrm{op}}=0.1i$. Adding fluorescence data has roughly the same effect on the estimation of ${\sigma }_{\mathrm{op}}$ like using five times more ion channels to estimate ${\sigma }_{\text{op}}^2$.

Figure 10

Download asset Open asset

Figure 10—source data 1

Five different sets of time traces of panel b.

All instrumental noise is already added.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig10-data1-v3.zip

Download elife-62714-fig10-data1-v3.zip

Figure 10—source data 2

Eighteen sets of time traces of panel c.

The number in the file name indicates the brightness.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig10-data2-v3.zip

Download elife-62714-fig10-data2-v3.zip

On the one side, every analog signal to be digitized needs analog filtering for antialiasing according to the Nyquist theorem. On the other side, every analog filter does not only suppress unwanted white noise but also distorts the dynamics (Figure 11a) of the signal of interest (Silberberg and Magleby, 1993). Therefore, (Qin et al., 2000) recommend to avoid analog filtering as much as possible in single-channel analysis and let the HMM analyze the data in the rawest available form, even with simultaneous drift correction (Sgouralis and Pressé, 2017a). One can also expect that analog filtering of a macroscopic signal is harmful for the inference of the KF and the RE approach. For the CCCO model considered herein we investigated the mean behavior (accuracy and precision) of the posterior of both algorithms with seven data sets (simulated at 100 kHz to mimic an analog signal). A digital fourth-order Bessel filter (Virtanen et al., 2020) was then applied. The maximum analysing frequency f_ana of the KF used is $100-400$ Hz to be comparable to cPCF setups. The slower frequency at which the Bayesian filter analyzes the data is necessary because the applied Bessel filter has caused additional time correlations in the originally white noise of the signal. Thus, an all-data-points fit would immediately violate the white noise assumption of Equation 4 which we restore by analyzing at a much lower frequency. We then let the time scales of the induced time correlations become larger and larger by decreasing f_cut. Physically, the absolute cut-off frequency f_cut is irrelevant; what matters is the magnitude of f_cut relative to f_ana and to the eigenvalues ${\alpha }_i$ of the ensemble (see, Appendix 3), since the eigenvalues determine the time evolution of the mean ensemble state, the autocorrelation, and Fourier spectrum of the fluctuations around the equilibrium distribution (Colquhoun et al., 1997b). The eigenvalues depend on the ligand concentration such that for a four-state model for each ligand concentration there are three relevant time scales $-1/{\alpha }_i$ (where $i=2,3,4$) plus the equilibrium solution which satisfies ${\alpha }_1=0$. For 10 different time series $3\cdot 10+3$ the outcome is to have different values of ${\alpha }_i$.Each eigenvalue is the inverse of the time constant of an exponential decay (see, Appendix 3). For this reason, we normalize in the following (Figure 11) the cut-off frequencies by ${\alpha }_2$ at the highest ligand concentration. We analyze the arithmetic mean from 7 different data sets of the median of the posterior of the rate matrix. The mean Euclidean error of the median (Figure 11b) and a series of quantiles demonstrate that overall the error of the mean median of the posterior KF (green) is smaller than that obtained by the RE. For unfiltered data, the accuracy of the mean median of the KF is increased by $\approx 1.6$. Based on the Euclidean error both algorithms benefit slightly from careful analog filtering for $f_{\mathrm{cut}}/{\alpha }_2\ge 1$ while the offset remains rather constant. A strong negative effect of analog filtering starts for both algorithms around $f_{\mathrm{cut}}\approx 1\mathrm{kHZ}$. This is induced by $f_{\mathrm{cut}}\to f_{\mathrm{ana}}$ (see, Appendix 10). In contrast, based on the level of each individual parameter of the rate matrix (Figure 11c 1–6) the bias induced by analog filtering immediately starts with $f_{\mathrm{cut}}=70\mathrm{kHz}$ (Figure 11c 1–3). Note, that visual inspection of the signal (Figure 11a) does not reveal signal distortions $f_{\mathrm{cut}}\ge 10\mathrm{kHz}$ though they are detected by both algorithms. For unfiltered data, the maximum of the posterior for the RE approach is a biased estimate $E[{\theta }_{\mathrm{ME}}]\ne {\theta }_{\mathrm{true}}$ for at least the parameters ${\tilde{k}}_{21},{\tilde{K}}_{21},{\tilde{K}}_{32}$ of the true value ${\theta }_{\mathrm{true}}$, which is explained (Milescu et al., 2005) by the fact that RE approaches ignore the autocorrelation of the intrinsic noise. Additionally, the data indicate that for ${\tilde{K}}_{43}$ the maximum of the posterior is even for the KF a biased estimate which we interpret as limitations induced by the fact that the mean vector and covariance-matrix do not constitute sufficient statistics as soon as Poisson distributed photon counting or open-channel noise blurs the signal. For the RE approach, the additional bias induced by the analog filter on the mean maximum of all parameters of the posterior starts with $f_{\mathrm{cut}}\approx 70$ kHz or, in other words, at the fastest time scale in the whole data set. The total bias in the estimate is reduced for k₂₁ with the additional bias from the analog filtering but increased for k₃₂ which for the Euclidean error leads at first to a small increase in accuracy. The KF is more robust towards analog filtering, as the results alter less with f_cut (given a reasonable f_cut), and less biased for unfiltered data in the estimates of these parameters. On the one hand, the Euclidean error shrinks for ${\displaystyle f_{\mathrm{c}\mathrm{u}\mathrm{t}}>10}$ kHz (Figure 11b). On the other hand, on the single-parameter level (Figure 11c 1–6), the parameter estimates pick up bias due the analog filtering even for high filter frequencies, in particular for the RE approach. Only for k₄₃ the KF is more biased than the RE approach.

Figure 11

Download asset Open asset

Figure 11—source data 1

Representative data set of cPCF data whose current has been filtered with decreasing f_cut.

In the name of each file is f_cut encoded in units of 100 kHz. The time axis is identical for all time traces and saved in the file Time.txt. The folder is similar to six other provided source data folders which have in total been used in this figure.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig11-data1-v3.zip

Download elife-62714-fig11-data1-v3.zip

Figure 11—source data 2

Representative data set of cPCF data whose current has been filtered with decreasing f_cut.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig11-data2-v3.zip

Download elife-62714-fig11-data2-v3.zip

Figure 11—source data 3

Representative data set of cPCF data whose current has been filtered with decreasing f_cut.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig11-data3-v3.zip

Download elife-62714-fig11-data3-v3.zip

Figure 11—source data 4

Representative data set of cPCF data whose current has been filtered with decreasing f_cut.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig11-data4-v3.zip

Download elife-62714-fig11-data4-v3.zip

Figure 11—source data 5

Representative data set of cPCF data whose current has been filtered with decreasing f_cut.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig11-data5-v3.zip

Download elife-62714-fig11-data5-v3.zip

Figure 11—source data 6

Representative data set of cPCF data whose current has been filtered with decreasing f_cut.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig11-data6-v3.zip

Download elife-62714-fig11-data6-v3.zip

The KF is the unique minimal variance Bayesian filter for a linear Gaussian process (Anderson and Moore, 2012) which means, given that the assumptions of the KF are fulfilled by the true process of interest, the KF is mathematically proven the best model-based filter to apply. Consequently, analog filtering does not provide an advantage unless it removes specific high frequency external noise sources (colored noise). We demonstrate (Appendix 10) this for PC data and varied f_cut and f_ana. On the downside, increasing f_ana makes the results of both algorithms more fragile if $f_{\mathrm{cut}}\gg f_{\mathrm{ana}}$ does not hold. Thus, the critical edge in Figure 11b is indeed induced by f_cut approaching f_ana. This suggests that the white noise assumption of both algorithms is violated. On the upside, if $f_{\mathrm{cut}}\gg f_{\mathrm{ana}}$ is given, the KF with an order of magnitude higher f_ana has a reduced bias of up to 20% for $f_{\mathrm{cut}}\to \mathrm{\infty }$ for individual parameters compared to the KF with lower f_ana. Additionally, a higher f_ana reduces the variance. To reduce the bias of parameter estimates to a minimum, the experimental design offers two remedies, either doing cPCF experiments with additional discussed advantages or using the KF at a high f_ana with even much higher f_cut.

By theoretical grounds a further argument for doing less analog filtering is that this benchmark analyzes data of a finite state Markov process, which is a coarse proxy for the true process. In reality, relaxation of a protein is a high-dimensional continuous-state Markov process with infinitely many relaxation time scales (eigenvalues) (Frauenfelder et al., 1991) which, however, might be grouped in slower experimentally accessible and non-accessible faster time scales (Noé et al., 2013). With larger data sets of higher quality from better experiments, the faster time scales might become accessible if they are not distorted by analog filtering. In conclusion, deciding on a specific kinetic scheme and inferring its parameters means finding a model which accommodates in the best way to the set of observed eigenvalues. Analog filtering hampers the RE, KF or HMM forward-backward algorithm (Qin et al., 2000) to correctly describe the faster time scales.

So far, we idealized the fluorescence data integration time as being instantaneously relative to the time scales of ensemble dynamics. In real experiments, the fluorescence signal of cPCF data has orders of magnitude longer minimal integration time $T_{\mathrm{int}}$ (time to record all voxels of a frame) or maximal integration frequency $f_{\mathrm{int}}=1/T_{\mathrm{int}}$, than the possible sampling frequency of current recordings. We mimic the finite integration time

by summing with a sliding window the 100 kHz signal including the white noise to obtain data at an effectively lower sampling frequency (Figure 12a). Additionally we set the Bessel filter for the current data to $f_{\mathrm{cut}}/{\alpha }_2=4.59$ or $f_{\mathrm{cut}}=90$ kHz. The fastest used analysing frequency is $f_{\mathrm{ana}}=500\mathrm{Hz}$. We scale mean photo brightness ${\lambda }_{\mathrm{b}}$ and background noise down such that the signal-to-noise ratio of the lower integration frequency data is the same as of the high-frequency data ${\lambda }_{\mathrm{b}}/T_{\mathrm{int}}=const$ . We do that in order to separate the bias from the finite integration time from other effects such as a better signal to noise ratios for each integrated point. Note that we only analyzed the plot until $f_{\mathrm{int}}=f_{\mathrm{ana}}$. Both algorithms incur very similar bias due to the finite integration time (Figure 12b). The KF (green) is more precise for high integration frequencies $f_{\mathrm{cut}}/{\alpha }_2$ until $f_{\mathrm{cut}}/{\alpha }_2\approx 0.08$ then the RE approach becomes more robust. Similar to Bessel-filtered current data (Figure 11b) on the single parameter level the systematic deviations start early for example $f_{\mathrm{int}}=10$ kHz for K₂₁ (Figure 12c4). Possibly the systematic deviations start (Figure 12c2) already at $f_{\mathrm{int}}=50$ kHz for k₃₂. The sudden increase of the Euclidean error (Figure 12b) of the mean median at $f_{\mathrm{cut}}/{\alpha }_2\approx 0.2$ occurs in this case not due to f_int approaching f_cut but due to $f_{\mathrm{int}}⪅{\alpha }_{1,2}$ for many ligand concentrations. To show this we plot the results of the fitting of five different data sets without including the highest 4 ligand concentrations (red) which means the largest eigenvalues are much smaller (Figure 12b,C1-6). Additionally, we keep $f_{\mathrm{int}}=const$. Although fluctuations of the posterior medians are higher, the KF becomes more robust. Note, that the fastest eigenvalues of these reduced data sets are indicated by the blue bars (Figure 12b and c4). Based on the Euclidean error (Figures 11a and 12a) the robustness of both algorithms against the cut-off frequency is compared with the robustness against the integration frequency found to be about an order of magnitude higher. That is related to a specific detail of the model used: the binding reaction, corresponds to the fastest time scales of the overall dynamic (difference between Figure 1b and c), which is exposed by the fluorescence signal. Thus, kinetic analysis of any data should make sure that the corresponding frequency of the most dominant timescales of the time series are much slower than the respective f_intf_cut independently of the investigated algorithms.

Figure 12

Download asset Open asset

Figure 12—source data 1

The data folder contains one data set of cPCF data.

The fluorescence has been integrated with decreasing f_int. The number in the file name encodes the amount of $100kHz$ data points which have been summed up to mimic the integration time. The time axis is identical for all time traces and saved in the file Time.txt. The folder is similar to six other provided source data folders which have in total been used in this figure.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig12-data1-v3.zip

Download elife-62714-fig12-data1-v3.zip

Figure 12—source data 2

The data folder contains one data set of cPCF data.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig12-data2-v3.zip

Download elife-62714-fig12-data2-v3.zip

Figure 12—source data 3

The data folder contains one data set of cPCF data.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig12-data3-v3.zip

Download elife-62714-fig12-data3-v3.zip

Figure 12—source data 4

The data folder contains one data set of cPCF data.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig12-data4-v3.zip

Download elife-62714-fig12-data4-v3.zip

Figure 12—source data 5

The data folder contains one data set of cPCF data.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig12-data5-v3.zip

Download elife-62714-fig12-data5-v3.zip

Figure 12—source data 6

The data folder contains one data set of cPCF data.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig12-data6-v3.zip

Download elife-62714-fig12-data6-v3.zip

Figure 12—source data 7

The data folder contains one data set of cPCF data.

https://cdn.elifesciences.org/articles/62714/elife-62714-fig12-data7-v3.zip

Download elife-62714-fig12-data7-v3.zip

We generalized the filter equations (Methods Equation 37, 38d, 57, 58 and 59) of the KF for analyzing the gating and binding dynamics of ligand-gated ion channels with a realistic signal-generating model for isolated patch-clamp (PC) and confocal patch-clamp fluorometry (cPCF) data including open-channel noise, photon-counting noise and background noise. Any other type of linear kinetic scheme (e.g. for voltage-dependent channels) and signal can be applied as long as the characteristics of the signal are sufficiently described by normal distributions. Our approach is derived by approximating the chemical master equation of a first order chemical reaction network (which ion channel experiments usually are) which is exact up to the second statistical moment. For first-order chemical reaction networks, the linear noise approximation (Wallace et al., 2012) are exact up to the second moment too (Grima, 2015). Thus, we can conclude that our Bayesian filter uses a time integrated version of the linear noise approximation. To our understanding of Wallace et al., 2012 our approach is thus equivalent to approaches based on the chemical Langevin or Fokker-Planck equations (Gillespie, 2002). Consequently, this also makes the considerations of the quality of the chemical Langevin equation as an approximation (Gillespie, 2000) of the chemical master equation valid for our approach. Compared to previous attempts Moffatt, 2007, this mathematical generalization is necessary (Figure 3b) in order to use Bayesian filters on macroscopic PC or cPCF data. With our algorithm, we demonstrate (Figures 3c and 7) that the common assumption that for large ensembles of ion channels simpler deterministic modeling by RE approaches is on par with stochastic modeling, such as a KF, is wrong in terms of Euclidean error and uncertainty quantification (Figures 5a–c ,–6a–b).

Enriching the data by fluorescence-based ligand binding reveals two regimes. In one regime, the two-dimensional data increase the accuracy of the parameter estimates up to $\approx 10$-fold (Figure 4a and c). In the other regime of lower channel expression, enriching the data by the second observable, makes non-identified parameters to identified parameters. The second observable in cPCF data decreases the overfitting tendency (Figure 4a, b and d) of the RE approach on the true process. Thus, in this regard the advantage of the KF becomes smaller. However, by exploiting Bayesian HDCV we gain a second perspective: We show for various signal-to-noise ratios (Figures 5a–c ,–6a–b) that the posterior sampled by a RE approach never covers the true values within a reasonable HDCV. Thus, the central feature of Bayesian statistics, exact uncertainty quantification by having the full posterior, is meaningless in combination with an RE approach (considering the type of data and set of signal-to-noise ratios that we tested). This even holds true for very pessimistic signal-to-noise assumptions Figure 6b. If HDCVs based on an RE approach cannot be trusted, the same applies to confidence volumes based on the curvature of the likelihood. This is not the case for the KF which delivers properly shaped posteriors (Figures 6a–c ,–5a–c). Increasing the model complexity, at unchanged PC data quality (Figure 7) shows that the RE approach displays unidentified rates even for large ion channel ensembles while our approach identified all parameters for the same data. We also investigated the robustness of both algorithms against the cut-off frequency of a Bessel filter (Figure 11) and showed the overall superior robustness of the KF against errors of analog filtering compared to the RE approach. Analog filtering has its limitations due to distorting the higher frequencies of the Fourier spectrum of the signal. Thus, one should let the KF sample as fast as possible, with a cut-off frequency of at least one order of magnitude higher than the sampling frequency of the KF.

Similar to the Bessel filter, the KF is more robust than the RE approach against errors due to the finite integration time. Nevertheless, it is crucial for both algorithms (Figure 12), that the intrinsic time scales (1/eigenvalues) of the process to be analyzed are slower than the integration time of the data points. Otherwise the accuracy of the inference deteriorates.

Altogether, we demonstrated the performance of the generalized Kalman filter on ion channel data for inference of kinetic schemes. Nevertheless, our approach can approximate any other stochastic system and signal distributions of linear (pseudo-first-order) kinetics (Sorenson and Alspach, 1971). Prospective extensions of the Bayesian filter, for example by Bayesian Gaussian sum filters or similar numerically brute force concepts such as particle filters (Golightly and Wilkinson, 2011; Gillespie and Golightly, 2012), can overcome modeling errors at low ion channel numbers or low photon fluxes.

Hereinafter, we derive the equations for our Bayesian filter for time series analysis of hidden linear chemical reaction networks (kinetic schemes). A detailed description of the experimental noise is provided in the Appendix 5.

The following conventions are generally used: Bold symbols are used for multi-dimensional objects such as vectors or matrices. Calligraphic letters are used for (some) vectorial time series and double-strike letters are used for probabilities and probability densities. Within the Bayesian paradigm (Hines, 2015; Ball, 2016), each unknown quantity, including model parameters $\mathit{\theta }$ and time series of occupancies of hidden states $N_T={\{\mathbf{n}(t_i)\}}_{i=1}^T$, are treated as random variables conditioned on observed time series data ${\displaystyle {\mathcal{Y}}_T={\mathbf{y}(t_i)}_{i=1}^T}$. The prior $\mathbb{P}(\mathit{\theta })={\prod }_j^{N_{\mathrm{par}}}\mathbb{P}({\theta }_j)$ or posterior distribution ${\displaystyle \mathbb{P}(\mathit{\theta }|{\mathcal{Y}}_T)}$ encodes the available information about the parameter values before and after analysing the data, respectively. According to the Bayesian theorem, the posterior distribution

is a probability distribution of a parameter set $\mathit{\theta }$ conditioned on ${\displaystyle {\mathcal{Y}}_T}$. The likelihood ${\displaystyle \mathbb{L}({\mathcal{Y}}_T|\mathit{\theta })}$ encodes the distribution of the data by modelling the intrinsic fluctuations of the protein as well as noise coming from the experimental devices. The prior provides either assumptions before measuring data or what has been learnt from previous experiments about $\mathit{\theta }$. The normalization constant

ensures that the posterior is a normalized distribution. The KF is a special class of models in the family of Bayesian filters (Ghahramani, 1997), which is a generalisation of the classical KF. Due to its linear time evolution (Equation 1), the KF is particularly useful for modeling time series data of ensembles dynamics of first order chemical networks. It delivers a set of recursive algebraic equations (Materials and methods Equation 28 and Equation 32) for each time point, which allows to express the prior ${\displaystyle \mathbb{P}(\mathbf{n}(t)|{\mathcal{Y}}_{t-1})}$ and (after incorporating $\mathbf{y}(t)$) the posterior ${\displaystyle \mathbb{P}(\mathbf{n}(t)|{\mathcal{Y}}_t)}$ occupancies of hidden states $\mathbf{n}(t)$ for all $t$ given a set of parameters $\mathit{\theta }$. This means the KF solves the filtering problem (inference of $N_T$) by explicitly modeling the time evolution of $\mathbf{n}(t)$ by multivariate normal distributions. This allows us to replace ${\displaystyle \mathbb{L}({\mathcal{Y}}_T|\mathit{\theta })}$ of Equation 20 by the expression of Equation 9.

The Bayesian framework (as demonstrated in this article) has various properties which makes it superior to ML estimation (MLE) (McElreath, 2018). Those properties are in particular useful for the analysis of biophysical data since very often the dynamics of interest are hidden or latent in the data. Models with a hidden structure are called singular. For regular (non-singular) statistical models, maxima ${\mathit{\theta }}_{\mathrm{ML}}$ of the posterior or likelihood converge in distribution

to the true value ${\mathit{\theta }}_{\mathrm{true}}$,where ${\mathbf{F}}^{-1}({\mathit{\theta }}_{\mathrm{true}})$ is the inverse Fisher information matrix. Under those conditions it is justified to derive from the curvature of the likelihood at ${\mathit{\theta }}_{\mathrm{ML}}$ via the Cramer-Rao-bound theorem

a confidence volume for the inferred parameters. In contrast, consider for example the type of data investigated in this study which probes the protein dynamics by current and light. Singularity means that the Fisher information matrix of a model is not invertible leading to the breakdown of the Cramer-Rao Bound theorem. Due to the breakdown, it cannot be guaranteed that even in the asymptotic limit the log-likelihood function can be approximated by a quadratic form Watanabe, 2007. Thus, usually the MLE does not obey Equation 22. Consequently, the posterior distribution is usually not a normal distribution either (Watanabe, 2007). Using the full posterior distribution without further approximations detects the resulting problems such as deviation from normality or non-identifiability of parameters, related to the singularity. In conclusion, the posterior is still a valid representation of parameter plausibility while ML fails.

In the following, we write the time $t$ as function argument rather than a subscript. Following standard approaches, we attribute to each state of the Markov model an element of a vector space with dimension $M$. At a time, a channel can only be in a single state. This implies that the set of possible states is S:=$\{(1,0,0,\mathrm{\dots }),(0,1,0,\mathrm{\dots }),\mathrm{\dots },(\mathrm{\dots },0,1)\}\subset {\{0,\mathrm{\hspace{0.17em}1}\}}^M$. In the following, Greek subscripts refer to different states while Latin subscripts refer to different channels. By $\mathbf{s}(t)={\mathbf{e}}_{\alpha }$ we specify that the channel is in state α at time $t$. Mathematically, ${\mathbf{e}}_{\alpha }$ stands for the α-th canonical unit Cartesian vector (Table 1).

Assuming that the state transitions can be modeled by a first order Markov process, the path probability can be decomposed as the product of conditional probabilities as follows:

Markov models (MMs) and rate models are widely used for modeling molecular kinetics (Appendix 2). They provide an interpretation of the data in terms of a set of conformational states and the transition rates between these states. For exactness it remains indispensable to model the dynamics with a HMMs (Noé et al., 2013). The core of a hidden Markov model is a conventional Markov model, which is supplemented with a an additional observation model. We will therefore first focus on a conventional Markov model. State-to-state transitions can be equivalently described with a transition matrix $\mathbf{T}$ in discrete time or with a rate matrix $\mathbf{K}$ in continuous time, as follows:

where $\exp$ is the matrix exponential. We aim to infer the elements of the rate matrix $\mathbf{K}$, constituting a kinetic model or reaction network of the channel. Realizations of sequences of states can be produced by the Doob-Gillespie algorithm Gillespie Daniel T. (1977). To derive succinct equations for the stochastic dynamics of a system, it is beneficial to consider the time propagation of an ensemble of virtual system copies. This allows to ascribe a probability vector $\mathbf{p}(t)$ to the system, in which each element $p_{\alpha }(t)$ is the probability to find the system at $t$ in state α. One can interpret the probability vector $\mathbf{p}$ as the instantaneous expectation value of the state vector $\mathbf{s}$.

The probability vector obeys the discrete-time Master equation

We model the experimentally observed system as a collection of non-interacting channels. A single channel can be modeled with a first-order MM. The same applies to the ensemble of non-interacting channels. We focus on modeling the time course of extensive macroscopic observables such as the mean current and fluorescence signals as well as their fluctuations. A central quantity is the vector $\mathbf{n}(t)$ which is the occupancy of the channel states at time $t$:

This quantity, like $\mathbf{s}(t)$, is a random variate. Unlike $\mathbf{s}(t)$, its domain is not confined to canonical unit vectors but to $\mathbf{n}\in {\mathbb{N}}^M$. From the linearity of Equation 28 in the channel dimension and from the single-channel CME Equation 27b one can immediately derive the equation for the time evolution of the mean occupancy $\overline{\mathbf{n}}(t)=E[\mathbf{n}(t)]$:

with the transition matrix $\mathbf{T}$. The full distribution $\mathbb{P}(\mathbf{n}(t+1)|\mathbf{n}(t))$ is a generalized multinomial distribution. To understand the generalized multinomial distribution and how it can be constructed from the (conventional) multinomial distribution, consider the simplified case where all channels are assumed to be in the same state α. Already after one time step, the channels will have spread out over the state space. The channel distribution after one time step is parametrized by the transition probabilities in row number α of the single-channel transition matrix $\mathbf{T}$. According to the theory of Markov models, the final distribution of channels originating from state α is the multinomial distribution

In general, the initial ensemble will not have only one but multiple occupied channel states. Because of the independence of the channels, one can imagine each initial sub-population spreading out over the state space independently. Each sub-population with initial state α gives rise to its own final multinomial distribution that contributes $n_{\beta }^{(\alpha )}$ transitions into state β to the total final distribution. The total number of channels at $t+1$ in each state can then be simply found by adding the number of channels transitioning out of the different states α.

Evidently, the total number of channels is conserved during propagation. The distribution of $\mathbf{n}(t+1)$, defined by Equations 30; 31, is called the generalized multinomial distribution:

While no simple expression exists for the generalized multinomial distribution, closed form expressions for its moments can be readily derived. For large $N_{\mathrm{ch}}$ each $\mathbb{P}({\mathbf{n}}^{(\alpha )}(t+1)\mid n_{\alpha }{\mathbf{e}}_{\alpha })$ can be approximated by a multivariate-normal distribution such that also $\mathrm{general}-\mathrm{multinomial}(\mathbf{n}(t),\mathbf{T})$ has a multivariate-normal approximation. In the next section, we combine the kinetics of channel ensembles with the KF by a moment expansion of the governing equations for the ensemble probability evolution.

The multinomial distribution (Fredkin and Rice, 1992) has the following mean and covariance matrix

where ${\mathbf{T}}_{:,\alpha }$ denotes the column number α of the transition matrix and $\mathrm{diag}({\mathbf{T}}_{:,\alpha })$ describes the diagonal matrix with ${\mathbf{T}}_{;,\alpha }$ on its diagonal. Combining Equation 31 with Equations 33; 34 we deduce the mean and variance of the generalized multinomial distribution:

Note that Equations 35; 36 are conditional expectations that depend on the random state $\mathbf{n}$ at the previous time ${\displaystyle t}$ and not only on the previous mean ${\displaystyle \overline{\mathbf{n}}}$. To find the absolute mean, the law of total expectation is applied to Equation 35, giving

in agreement with the simple derivation of Equation 29. We introduce a shorthand ${\displaystyle \mathbf{P}(t):=\mathrm{c}\mathrm{o}\mathrm{v}(\mathbf{n}(t),\mathbf{n}(t))}$ for the absolute covariance matrix of ${\displaystyle \mathbf{n}(t+1)}$. Similarly, ${\displaystyle \mathbf{P}(t)}$ can be found by applying the law of total variance decomposition (Weiss, 2005 to Equations 35; 36), giving

Equations 37, 38d dare compact analytical expressions for the mean and the covariance matrix of the occupancy vector ${\displaystyle \mathbf{n}}$ at ${\displaystyle t+1}$ that depend on the mean ${\displaystyle \overline{\mathbf{n}}}$ and covariance matrix ${\displaystyle \mathbf{P}}$ at the previous time step $t$. Chaining these equations for different time steps ${\displaystyle t=0,\dots ,T}$ allows to model the whole evolution of a channel ensemble. Moreover, these two equations together with the output statistics of ${\displaystyle \mathbb{O}(\mathbf{y}|\mathbf{n}(t))}$ are sufficient to formulate correction equations Equation 59 of the KF (Moffatt, 2007; Anderson and Moore, 2012). These equations will be used in a Bayesian context to sample the posterior distribution of the model parameters. The sampling entails repeated numerical evaluation of the model likelihood. Therefore, analytical equations for the ensemble evolution that can be quickly evaluated on a computer millions of times are indispensable. This was achieved by deriving Equation 37, Equation 38d. Comparing Equation 38d with the KF prediction equation (Anderson and Moore, 2012) for ${\displaystyle \mathbf{P}(t)}$, we obtain the state-dependent covariance matrix of Equation 3 as

In the following section on properties of measured data and the KF, we no longer need to refer to the random variate ${\displaystyle \mathbf{n}(t)}$. All subsequent equations can be formulated by only using the mean hidden state ${\displaystyle \overline{\mathbf{n}}(t)}$ and the variance-covariance matrix of the hidden state $\mathbf{P}(t)$. We therefore drop the overbar in $\overline{\mathbf{n}}(t)$ so that the symbol $\mathbf{n}(t)$ refers from now on to the mean hidden state.

In the following, we develop a model for the conditional observation distribution ${\displaystyle \mathbb{O}(\mathbf{y}|\mathbf{n}(t))}$ (Appendix 5 for experimental details). Together with the hidden ensemble dynamics this will enable us to derive the output statistics of the KF (see, below). Let $\mathbf{y}(t)$ be the vector of all observations at $t$. Components of the vector are the ion current and fluorescence intensity.

As outlined in the introduction part, in Equation 4 we model the observation by using a conditional probability distribution ${\displaystyle \mathbb{O}(\mathbf{y}(t)|\mathbf{n}(t))}$ that only depends on the mean hidden state $\mathbf{n}(t)$, as well as on fixed channel and other measurement parameters. ${\displaystyle \mathbb{O}(\mathbf{y}(t)|\mathbf{n}(t))}$ is modeled as a multivariate normal distribution with mean $\mathbf{Hn}(t)$ and variance-covariance matrix $\mathbf{\Sigma }(t)$, that can in general depend on the mean state vector $\mathbf{n}(t)$ (much like the covariance matrix of the kinetics in (Equation 38d) ). The observation matrix $\mathbf{H}\in {ℝ}^{N_{\mathrm{obs}}\times M}$ projects the hidden state vector $\mathbf{n}(t)$ onto $\mathbf{Hn}(t)\in {ℝ}^{N_{\mathrm{obs}}}$, the observation space. The observation distribution is

This measurement model is very flexible and allows to include different types of signals and error sources arising from both the molecules and the instruments. A summary of the signals and sources of measurement error and their contributions to the parameters of ${\displaystyle \mathbb{O}(\mathbf{y}(t)|\mathbf{n}(t))}$ is provided by Table 2. Below we address the two types of signals and four noise sources one by one. For this, we decompose the observation matrix and the observation noise covariance matrix into the individual terms:

In the following, we report the individual matrices for the exemplary CCCO model with one open state $\alpha =4$ and three closed states $\alpha =1,2,3$. Matrices can be constructed analogously for the other models. For the definition of ${\mathbf{\Sigma }}_{\mathrm{back}}$ refer to (Appendix 5).