Bayesian Statistical Inference in Ion-Channel Models with Exact Missed Event Correction

Identifying the generator matrix of a stationary Markov chain using partially observable data

Given that most states in real-world systems are inaccessible, it is critical to study the inverse problem of an irreversibly stationary Markov chain regarding how a generator matrix can be identified using minimal observations. The hitting-time distribution of an irreversibly stationary Markov chain is first generalized from a reversible case. The hitting-time distribution is then decoded via the taboo rate, and the results show remarkably that under mild conditions, the generator matrix of a reversible Markov chain or a specific case of irreversibly stationary ones can be identified by utilizing observations from all leaves and two adjacent states in each cycle. Several algorithms are proposed for calculating the generator matrix accurately, and numerical examples are presented to confirm their validity and efficiency. An application to neurophysiology is provided to demonstrate the applicability of such statistics to real-world data. This means that partially observable data can be used to identify the generator matrix of a stationary Markov chain.

Permissive and nonpermissive channel closings in CFTR revealed by a factor graph inference algorithm

The closing of the gated ion channel in the cystic fibrosis transmembrane conductance regulator can be categorized as nonpermissive to reopening, which involves the unbinding of ADP or ATP, or permissive, which does not. Identifying the type of closing is of interest as interactions with nucleotides can be affected in mutants or by introducing agonists. However, all closings are electrically silent and difficult to differentiate. For single-channel patch-clamp traces, we show that the type of the closing can be accurately determined by an inference algorithm implemented on a factor graph, which we demonstrate using both simulated and lab-obtained patch-clamp traces.

Bayesian inference of kinetic schemes for ion channels by Kalman filtering

Inferring adequate kinetic schemes for ion channel gating from ensemble currents is a daunting task due to limited information in the data. We address this problem by using a parallelized Bayesian filter to specify hidden Markov models for current and fluorescence data. We demonstrate the flexibility of this algorithm by including different noise distributions. Our generalized Kalman filter outperforms both a classical Kalman filter and a rate equation approach when applied to patch-clamp data exhibiting realistic open-channel noise. The derived generalization also enables inclusion of orthogonal fluorescence data, making unidentifiable parameters identifiable and increasing the accuracy of the parameter estimates by an order of magnitude. By using Bayesian highest credibility volumes, we found that our approach, in contrast to the rate equation approach, yields a realistic uncertainty quantification. Furthermore, the Bayesian filter delivers negligibly biased estimates for a wider range of data quality. For some data sets, it identifies more parameters than the rate equation approach. These results also demonstrate the power of assessing the validity of algorithms by Bayesian credibility volumes in general. Finally, we show that our Bayesian filter is more robust against errors induced by either analog filtering before analog-to-digital conversion or by limited time resolution of fluorescence data than a rate equation approach.

Parameter Optimization for Ion Channel Models: Integrating New Data with Known Channel Properties

Ion channels play a central role in membrane physiology, but to fully understand how they operate, one must have accurate kinetic mechanisms. Estimating kinetics is not trivial when the mechanism is complex, and a large number of parameters must be extracted from data. Furthermore, the information contained in the data is often limited, and the model may not be fully determined. The solution is to reduce the number of parameters and to estimate them in such a way that they not only describe well the new data but also agree with the existing knowledge. In a previous study, we presented a comprehensive formalism for estimating kinetic parameters subject to a variety of explicit and implicit constraints that define quantitative relationships between parameters and describe specific mechanism properties. Here, we introduce the reader to the QuB software, which implements this constraining formalism. QuB features a powerful visual interface and a high-level scripting language that can be used to formulate kinetic models and constraints of arbitrary complexity, and to efficiently estimate the parameters from a variety of experimental data.

Toward a More General Understanding of Bohr's Complementarity: Insights from Modeling of Ion Channels

Some contemporary theorists such as Mazzocchi, Theise and Kafatos are convinced that the reformed complementarity may redefine how we might exploit the complexity theory in 21st-century life sciences research. However, the motives behind the profound re-invention of "biological complementarity" need to be substantiated with concrete shreds of evidence about this principle's applicability in real-life science experimentation, which we found missing in the literature. This paper discusses such pieces of evidence by confronting Bohr's complementarity and ion channel modeling practice. We examine whether and to what extent this principle might assist in developing ion channel models incorporating both deterministic and stochastic solutions. According to the "mutual exclusiveness of experimental setups" version of Bohr's complementarity, this principle is needed when two mutually exclusive perspectives or approaches are right, necessary in a particular context, and are not contradictory as they arise in mutually exclusive conditions (mutually exclusive experimental or modeling setups). A detailed examination of the modeling practice reveals that both solutions are often used simultaneously in a single ion channel model, suggesting that the opposite conceptual frameworks can coexist in the same modeling setup. We concluded that Bohr's complementarity might find applications in these complex modeling setups but only through its realistic phenomenological interpretation that allows applying different modes of description regardless of the nature of the underlying ion channel opening process. Also, we propose the combined use of complementarity and Complex thinking in building the multifaceted ion channel models. Overall, this paper's results support the efforts to establish a more general form of complementarity to meet today's complexity theory-inspired life sciences modeling demands.

Analysis of patchclamp recordings: model-free multiscale methods and software

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.

Calibration of ionic and cellular cardiac electrophysiology models

Cardiac electrophysiology models are among the most mature and well-studied mathematical models of biological systems. This maturity is bringing new challenges as models are being used increasingly to make quantitative rather than qualitative predictions. As such, calibrating the parameters within ion current and action potential (AP) models to experimental data sets is a crucial step in constructing a predictive model. This review highlights some of the fundamental concepts in cardiac model calibration and is intended to be readily understood by computational and mathematical modelers working in other fields of biology. We discuss the classic and latest approaches to calibration in the electrophysiology field, at both the ion channel and cellular AP scales. We end with a discussion of the many challenges that work to date has raised and the need for reproducible descriptions of the calibration process to enable models to be recalibrated to new data sets and built upon for new studies. This article is categorized under: Analytical and Computational Methods > Computational Methods Physiology > Mammalian Physiology in Health and Disease Models of Systems Properties and Processes > Cellular Models.

Exposing Hidden High-Affinity RNA Conformational States

RNA recognition frequently results in conformational changes that optimize intermolecular binding. As a consequence, the overall binding affinity of RNA to its binding partners depends not only on the intermolecular interactions formed in the bound state but also on the energy cost associated with changing the RNA conformational distribution. Measuring these "conformational penalties" is, however, challenging because bound RNA conformations tend to have equilibrium populations in the absence of the binding partner that fall outside detection by conventional biophysical methods. In this study we employ as a model system HIV-1 TAR RNA and its interaction with the ligand argininamide (ARG), a mimic of TAR's cognate protein binding partner, the transactivator Tat. We use NMR chemical shift perturbations and relaxation dispersion in combination with Bayesian inference to develop a detailed thermodynamic model of coupled conformational change and ligand binding. Starting from a comprehensive 12-state model of the equilibrium, we estimate the energies of six distinct detectable thermodynamic states that are not accessible by currently available methods. Our approach identifies a minimum of four RNA intermediates that differ in terms of the TAR conformation and ARG occupancy. The dominant bound TAR conformation features two bound ARG ligands and has an equilibrium population in the absence of ARG that is below detection limit. Consequently, even though ARG binds to TAR with an apparent overall weak affinity (Kdapp ≈ 0.2 mM), it binds the prefolded conformation with a Kd in the nM range. Our results show that conformational penalties can be major determinants of RNA-ligand binding affinity as well as a source of binding cooperativity, with important implications for a predictive understanding of how RNA is recognized and for RNA-targeted drug discovery.

A HIDDEN MARKOV MODEL APPROACH TO CHARACTERIZING THE PHOTO-SWITCHING BEHAVIOR OF FLUOROPHORES

Fluorescing molecules (fluorophores) that stochastically switch between photon-emitting and dark states underpin some of the most celebrated advancements in super-resolution microscopy. While this stochastic behavior has been heavily exploited, full characterization of the underlying models can potentially drive forward further imaging methodologies. Under the assumption that fluorophores move between fluorescing and dark states as continuous time Markov processes, the goal is to use a sequence of images to select a model and estimate the transition rates. We use a hidden Markov model to relate the observed discrete time signal to the hidden continuous time process. With imaging involving several repeat exposures of the fluorophore, we show the observed signal depends on both the current and past states of the hidden process, producing emission probabilities that depend on the transition rate parameters to be estimated. To tackle this unusual coupling of the transition and emission probabilities, we conceive transmission (transition-emission) matrices that capture all dependencies of the model. We provide a scheme of computing these matrices and adapt the forward-backward algorithm to compute a likelihood which is readily optimized to provide rate estimates. When confronted with several model proposals, combining this procedure with the Bayesian Information Criterion provides accurate model selection.

Quantitative Kinetic Models from Intravital Microscopy: A Case Study Using Hepatic Transport

The liver performs critical physiological functions, including metabolizing and removing substances, such as toxins and drugs, from the bloodstream. Hepatotoxicity itself is intimately linked to abnormal hepatic transport, and hepatotoxicity remains the primary reason drugs in development fail and approved drugs are withdrawn from the market. For this reason, we propose to analyze, across liver compartments, the transport kinetics of fluorescein-a fluorescent marker used as a proxy for drug molecules-using intravital microscopy data. To resolve the transport kinetics quantitatively from fluorescence data, we account for the effect that different liver compartments (with different chemical properties) have on fluorescein's emission rate. To do so, we develop ordinary differential equation transport models from the data where the kinetics is related to the observable fluorescence levels by "measurement parameters" that vary across different liver compartments. On account of the steep non-linearities in the kinetics and stochasticity inherent to the model, we infer kinetic and measurement parameters by generalizing the method of parameter cascades. For this application, the method of parameter cascades ensures fast and precise parameter estimates from noisy time traces.

Moving average filtering with deconvolution (MAD) for hidden Markov model with filtering and correlated noise

Ion channel data recorded using the patch clamp technique are low-pass filtered to remove high-frequency noise. Almanjahie et al. (Eur Biophys J 44:545-556, 2015) based statistical analysis of such data on a hidden Markov model (HMM) with a moving average adjustment for the filter but without correlated noise, and used the EM algorithm for parameter estimation. In this paper, we extend their model to include correlated noise, using signal processing methods and deconvolution to pre-whiten the noise. The resulting data can be modelled as a standard HMM and parameter estimates are again obtained using the EM algorithm. We evaluate this approach using simulated data and also apply it to real data obtained from the mechanosensitive channel of large conductance (MscL) in Escherichia coli. Estimates of mean conductances are comparable to literature values. The key advantages of this method are that it is much simpler and computationally considerably more efficient than currently used HMM methods that include filtering and correlated noise.

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

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

Estimating kinetic mechanisms with prior knowledge I: Linear parameter constraints

New mathematical tools are needed to incorporate existing knowledge into kinetic models of ion channels and other proteins. Salari et al. describe an algebraic transformation that can enforce linearly interdependent parameters into kinetic models in order to test new hypotheses.

To understand how ion channels and other proteins function at the molecular and cellular levels, one must decrypt their kinetic mechanisms. Sophisticated algorithms have been developed that can be used to extract kinetic parameters from a variety of experimental data types. However, formulating models that not only explain new data, but are also consistent with existing knowledge, remains a challenge. Here, we present a two-part study describing a mathematical and computational formalism that can be used to enforce prior knowledge into the model using constraints. In this first part, we focus on constraints that enforce explicit linear relationships involving rate constants or other model parameters. We develop a simple, linear algebra–based transformation that can be applied to enforce many types of model properties and assumptions, such as microscopic reversibility, allosteric gating, and equality and inequality parameter relationships. This transformation converts the set of linearly interdependent model parameters into a reduced set of independent parameters, which can be passed to an automated search engine for model optimization. In the companion article, we introduce a complementary method that can be used to enforce arbitrary parameter relationships and any constraints that quantify the behavior of the model under certain conditions. The procedures described in this study can, in principle, be coupled to any of the existing methods for solving molecular kinetics for ion channels or other proteins. These concepts can be used not only to enforce existing knowledge but also to formulate and test new hypotheses.

The reproducibility of research and the misinterpretation of p-values

We wish to answer this question: If you observe a 'significant' p-value after doing a single unbiased experiment, what is the probability that your result is a false positive? The weak evidence provided by p-values between 0.01 and 0.05 is explored by exact calculations of false positive risks. When you observe p = 0.05, the odds in favour of there being a real effect (given by the likelihood ratio) are about 3 : 1. This is far weaker evidence than the odds of 19 to 1 that might, wrongly, be inferred from the p-value. And if you want to limit the false positive risk to 5%, you would have to assume that you were 87% sure that there was a real effect before the experiment was done. If you observe p = 0.001 in a well-powered experiment, it gives a likelihood ratio of almost 100 : 1 odds on there being a real effect. That would usually be regarded as conclusive. But the false positive risk would still be 8% if the prior probability of a real effect were only 0.1. And, in this case, if you wanted to achieve a false positive risk of 5% you would need to observe p = 0.00045. It is recommended that the terms 'significant' and 'non-significant' should never be used. Rather, p-values should be supplemented by specifying the prior probability that would be needed to produce a specified (e.g. 5%) false positive risk. It may also be helpful to specify the minimum false positive risk associated with the observed p-value. Despite decades of warnings, many areas of science still insist on labelling a result of p < 0.05 as 'statistically significant'. This practice must contribute to the lack of reproducibility in some areas of science. This is before you get to the many other well-known problems, like multiple comparisons, lack of randomization and p-hacking. Precise inductive inference is impossible and replication is the only way to be sure. Science is endangered by statistical misunderstanding, and by senior people who impose perverse incentives on scientists.

Calcium binding and voltage gating in Cx46 hemichannels

The opening of connexin (Cx) hemichannels in the membrane is tightly regulated by calcium (Ca2+) and membrane voltage. Electrophysiological and atomic force microscopy experiments indicate that Ca2+ stabilizes the hemichannel closed state. However, structural data show that Ca2+ binding induces an electrostatic seal preventing ion transport without significant structural rearrangements. In agreement with the closed-state stabilization hypothesis, we found that the apparent Ca2+ sensitivity is increased as the voltage is made more negative. Moreover, the voltage and Ca2+ dependence of the channel kinetics indicate that the voltage sensor movement and Ca2+ binding are allosterically coupled. An allosteric kinetic model in which the Ca2+ decreases the energy necessary to deactivate the voltage sensor reproduces the effects of Ca2+ and voltage in Cx46 hemichannels. In agreement with the model and suggesting a conformational change that narrows the pore, Ca2+ inhibits the water flux through Cx hemichannels. We conclude that Ca2+ and voltage act allosterically to stabilize the closed conformation of Cx46 hemichannels.

A single Markov-type kinetic model accounting for the macroscopic currents of all human voltage-gated sodium channel isoforms

Modelling ionic channels represents a fundamental step towards developing biologically detailed neuron models. Until recently, the voltage-gated ion channels have been mainly modelled according to the formalism introduced by the seminal works of Hodgkin and Huxley (HH). However, following the continuing achievements in the biophysical and molecular comprehension of these pore-forming transmembrane proteins, the HH formalism turned out to carry limitations and inconsistencies in reproducing the ion-channels electrophysiological behaviour. At the same time, Markov-type kinetic models have been increasingly proven to successfully replicate both the electrophysiological and biophysical features of different ion channels. However, in order to model even the finest non-conducting molecular conformational change, they are often equipped with a considerable number of states and related transitions, which make them computationally heavy and less suitable for implementation in conductance-based neurons and large networks of those. In this purely modelling study we develop a Markov-type kinetic model for all human voltage-gated sodium channels (VGSCs). The model framework is detailed, unifying (i.e., it accounts for all ion-channel isoforms) and computationally efficient (i.e. with a minimal set of states and transitions). The electrophysiological data to be modelled are gathered from previously published studies on whole-cell patch-clamp experiments in mammalian cell lines heterologously expressing the human VGSC subtypes (from NaV1.1 to NaV1.9). By adopting a minimum sequence of states, and using the same state diagram for all the distinct isoforms, the model ensures the lightest computational load when used in neuron models and neural networks of increasing complexity. The transitions between the states are described by original ordinary differential equations, which represent the rate of the state transitions as a function of voltage (i.e., membrane potential). The kinetic model, developed in the NEURON simulation environment, appears to be the simplest and most parsimonious way for a detailed phenomenological description of the human VGSCs electrophysiological behaviour.

Structural identifiability of equilibrium ligand-binding parameters

Precise mathematical descriptions of ligand–protein interactions are hindered by the inability to experimentally measure affinity and cooperativity, although these parameters can be estimated from agonist binding models. Middendorf and Aldrich present a method to determine the accuracy of parameters estimated in this way.

Understanding the interactions of proteins with their ligands requires knowledge of molecular properties, such as binding site affinities and the effects that binding at one site exerts on binding at other sites (cooperativity). These properties cannot be measured directly and are usually estimated by fitting binding data with models that contain these quantities as parameters. In this study, we present a general method for answering the critical question of whether these parameters are identifiable (i.e., whether their estimates are accurate and unique). In cases in which parameter estimates are not unique, our analysis provides insight into the fundamental causes of nonidentifiability. This approach can thus serve as a guide for the proper design and analysis of protein–ligand binding experiments. We show that the equilibrium total binding relation can be reduced to a conserved mathematical form for all models composed solely of bimolecular association reactions and to a related, conserved form for all models composed of arbitrary combinations of binding and conformational equilibria. This canonical mathematical structure implies a universal parameterization of the binding relation that is consistent with virtually any physically reasonable binding model, for proteins with any number of binding sites. Matrix algebraic methods are used to prove that these universal parameter sets are structurally identifiable (SI; i.e., identifiable under conditions of noiseless data). A general approach for assessing and understanding the factors governing practical identifiability (i.e., the identifiability under conditions of real, noisy data) of these SI parameter sets is presented in the companion paper by Middendorf and Aldrich (2017. J. Gen. Physiol.https://doi.org/10.1085/jgp.201611703).