Markov chain Monte Carlo fitting of single-channel data from inositol trisphosphate receptors

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.

Computing Optimal Properties of Drugs Using Mathematical Models of Single Channel Dynamics

Single channel dynamics can be modeled using stochastic differential equations, and the dynamics of the state of the channel (e.g. open, closed, inactivated) can be represented using Markov models. Such models can also be used to represent the effect of mutations as well as the effect of drugs used to alleviate deleterious effects of mutations. Based on the Markov model and the stochastic models of the single channel, it is possible to derive deterministic partial differential equations (PDEs) giving the probability density functions (PDFs) of the states of the Markov model. In this study, we have analyzed PDEs modeling wild type (WT) channels, mutant channels (MT) and mutant channels for which a drug has been applied (MTD). Our aim is to show that it is possible to optimize the parameters of a given drug such that the solution of theMTD model is very close to that of the WT: the mutation’s effect is, theoretically, reduced significantly.We will present the mathematical framework underpinning this methodology and apply it to several examples. In particular, we will show that it is possible to use the method to, theoretically, improve the properties of some well-known existing drugs.

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

The stochastic behavior of single ion channels is most often described as an aggregated continuous-time Markov process with discrete states. For ligand-gated channels each state can represent a different conformation of the channel protein or a different number of bound ligands. Single-channel recordings show only whether the channel is open or shut: states of equal conductance are aggregated, so transitions between them have to be inferred indirectly. The requirement to filter noise from the raw signal further complicates the modeling process, as it limits the time resolution of the data. The consequence of the reduced bandwidth is that openings or shuttings that are shorter than the resolution cannot be observed; these are known as missed events. Postulated models fitted using filtered data must therefore explicitly account for missed events to avoid bias in the estimation of rate parameters and therefore assess parameter identifiability accurately. In this article, we present the first, to our knowledge, Bayesian modeling of ion-channels with exact missed events correction. Bayesian analysis represents uncertain knowledge of the true value of model parameters by considering these parameters as random variables. This allows us to gain a full appreciation of parameter identifiability and uncertainty when estimating values for model parameters. However, Bayesian inference is particularly challenging in this context as the correction for missed events increases the computational complexity of the model likelihood. Nonetheless, we successfully implemented a two-step Markov chain Monte Carlo method that we called "BICME", which performs Bayesian inference in models of realistic complexity. The method is demonstrated on synthetic and real single-channel data from muscle nicotinic acetylcholine channels. We show that parameter uncertainty can be characterized more accurately than with maximum-likelihood methods. Our code for performing inference in these ion channel models is publicly available.

Decoding Single Molecule Time Traces with Dynamic Disorder

Single molecule time trajectories of biomolecules provide glimpses into complex folding landscapes that are difficult to visualize using conventional ensemble measurements. Recent experiments and theoretical analyses have highlighted dynamic disorder in certain classes of biomolecules, whose dynamic pattern of conformational transitions is affected by slower transition dynamics of internal state hidden in a low dimensional projection. A systematic means to analyze such data is, however, currently not well developed. Here we report a new algorithm—Variational Bayes-double chain Markov model (VB-DCMM)—to analyze single molecule time trajectories that display dynamic disorder. The proposed analysis employing VB-DCMM allows us to detect the presence of dynamic disorder, if any, in each trajectory, identify the number of internal states, and estimate transition rates between the internal states as well as the rates of conformational transition within each internal state. Applying VB-DCMM algorithm to single molecule FRET data of H-DNA in 100 mM-Na⁺ solution, followed by data clustering, we show that at least 6 kinetic paths linking 4 distinct internal states are required to correctly interpret the duplex-triplex transitions of H-DNA.

We have developed a new algorithm to better decode single molecule data with dynamic disorder. Our new algorithm, which represents a substantial improvement over other methodologies, can detect the presence of dynamic disorder in each trajectory and quantify the kinetic characteristics of underlying energy landscape. As a model system, we applied our algorithm to the single molecule FRET time traces of H-DNA. While duplex-triplex transitions of H-DNA are conventionally interpreted in terms of two-state kinetics, slowly varying dynamic patterns corresponding to hidden internal states can also be identified from the individual time traces. Our algorithm reveals that at least 4 distinct internal states are required to correctly interpret the data.

Computing rates of Markov models of voltage-gated ion channels by inverting partial differential equations governing the probability density functions of the conducting and non-conducting states

Markov models are ubiquitously used to represent the function of single ion channels. However, solving the inverse problem to construct a Markov model of single channel dynamics from bilayer or patch-clamp recordings remains challenging, particularly for channels involving complex gating processes. Methods for solving the inverse problem are generally based on data from voltage clamp measurements. Here, we describe an alternative approach to this problem based on measurements of voltage traces. The voltage traces define probability density functions of the functional states of an ion channel. These probability density functions can also be computed by solving a deterministic system of partial differential equations. The inversion is based on tuning the rates of the Markov models used in the deterministic system of partial differential equations such that the solution mimics the properties of the probability density function gathered from (pseudo) experimental data as well as possible. The optimization is done by defining a cost function to measure the difference between the deterministic solution and the solution based on experimental data. By evoking the properties of this function, it is possible to infer whether the rates of the Markov model are identifiable by our method. We present applications to Markov model well-known from the literature.

Statistical analysis of modal gating in ion channels

Ion channels regulate the concentrations of ions within cells. By stochastically opening and closing its pore, they enable or prevent ions from crossing the cell membrane. However, rather than opening with a constant probability, many ion channels switch between several different levels of activity even if the experimental conditions are unchanged. This phenomenon is known as modal gating: instead of directly adapting its activity, the channel seems to mix sojourns in active and inactive modes in order to exhibit intermediate open probabilities. Evidence is accumulating that modal gating rather than modulation of opening and closing at a faster time scale is the primary regulatory mechanism of ion channels. However, currently, no method is available for reliably calculating sojourns in different modes. In order to address this challenge, we develop a statistical framework for segmenting single-channel datasets into segments that are characteristic for particular modes. The algorithm finds the number of mode changes, detects their locations and infers the open probabilities of the modes. We apply our approach to data from the inositol-trisphosphate receptor. Based upon these results, we propose that mode changes originate from alternative conformational states of the channel protein that determine a certain level of channel activity.

Quantitative properties and receptor reserve of the DAG and PKC branch of G(q)-coupled receptor signaling

G_q protein–coupled receptors (G_qPCRs) of the plasma membrane activate the phospholipase C (PLC) signaling cascade. PLC cleaves the membrane lipid phosphatidylinositol 4,5-bisphosphate (PIP₂) into the second messengers diacylgycerol (DAG) and inositol 1,4,5-trisphosphate (IP₃), leading to calcium release, protein kinase C (PKC) activation, and in some cases, PIP₂ depletion. We determine the kinetics of each of these downstream endpoints and also ask which is responsible for the inhibition of KCNQ2/3 (K_V7.2/7.3) potassium channels in single living tsA-201 cells. We measure DAG production and PKC activity by Förster resonance energy transfer–based sensors, and PIP₂ by KCNQ2/3 channels. Fully activating endogenous purinergic receptors by uridine 5′triphosphate (UTP) leads to calcium release, DAG production, and PKC activation, but no net PIP₂ depletion. Fully activating high-density transfected muscarinic receptors (M₁Rs) by oxotremorine-M (Oxo-M) leads to similar calcium, DAG, and PKC signals, but PIP₂ is depleted. KCNQ2/3 channels are inhibited by the Oxo-M treatment (85%) and not by UTP (<1%), indicating that depletion of PIP₂ is required to inhibit KCNQ2/3 in response to receptor activation. Overexpression of A kinase–anchoring protein (AKAP)79 or calmodulin (CaM) does not increase KCNQ2/3 inhibition by UTP. From these results and measurements of IP₃ and calcium presented in our companion paper (Dickson et al. 2013. J. Gen. Physiol.http://dx.doi.org/10.1085/jgp.201210886), we extend our kinetic model for signaling from M₁Rs to DAG/PKC and IP₃/calcium signaling. We conclude that calcium/CaM and PKC-mediated phosphorylation do not underlie dynamic KCNQ2/3 channel inhibition during G_qPCR activation in tsA-201 cells. Finally, our experimental data provide indirect evidence for cleavage of PI(4)P by PLC in living cells, and our modeling revisits/explains the concept of receptor reserve with measurements from all steps of G_qPCR signaling.

Dynamic metabolic control of an ion channel

G-protein-coupled receptors mediate responses to external stimuli in various cell types. We are interested in the modulation of KCNQ2/3 potassium channels by the Gq-coupled M1 muscarinic (acetylcholine) receptor (M1R). Here, we describe development of a mathematical model that incorporates all known steps along the M1R signaling cascade and accurately reproduces the macroscopic behavior we observe when KCNQ2/3 currents are inhibited following M1R activation. Gq protein-coupled receptors of the plasma membrane activate phospholipase C (PLC) which cleaves the minor plasma membrane lipid phosphatidylinositol 4,5-bisphosphate (PI(4,5)P2) into the second messengers diacylgycerol and inositol 1,4,5-trisphosphate, leading to calcium release, protein kinase C (PKC) activation, and PI(4,5)P2 depletion. Combining optical and electrical techniques with knowledge of relative abundance of each signaling component has allowed us to develop a kinetic model and determine that (i) M1R activation and M1R/Gβ interaction are fast; (ii) Gαq/Gβ separation and Gαq/PLC interaction have intermediate time constants; (iii) the amount of activated PLC limits the rate of KCNQ2/3 suppression; (iv) weak PLC activation can elicit robust calcium signals without net PI(4,5)P2 depletion or KCNQ2/3 channel inhibition; and (v) depletion of PI(4,5)P2, and not calcium/CaM or PKC-mediated phosphorylation, closes KCNQ2/3 potassium channels, thereby increasing neuronal excitability.

Bayesian approaches for mechanistic ion channel modeling

We consider the Bayesian analysis of mechanistic models describing the dynamic behavior of ligand-gated ion channels. The opening of the transmembrane pore in an ion channel is brought about by conformational changes in the protein, which results in a flow of ions through the pore. Remarkably, given the diameter of the pore, the flow of ions from a small number of channels or indeed from a single ion channel molecule can be recorded experimentally. This produces a large time-series of high-resolution experimental data, which can be used to investigate the gating process of these channels. We give a brief overview of the achievements and limitations of alternative maximum-likelihood approaches to this type of modeling, before investigating the statistical issues associated with analyzing stochastic model reaction mechanisms from a Bayesian perspective. Finally, we compare a number of Markov chain Monte Carlo algorithms that may be used to tackle this challenging inference problem.

MCMC can detect nonidentifiable models

Continuous-time Markov models have been considered the best representation for the stochastic dynamics of ion channels for more than thirty years. For most single-channel data sets, several open and closed states are required for accurately representing the dynamics. However, each data point only shows if the channel is open or closed but not in which state it is. Consequently, some model structures are inherently overparameterized and therefore, in principle, unsuitable for representing any data--those models are called "nonidentifiable". As of this writing, it seems to be poorly understood which continuous-time Markov models are identifiable and which are not, therefore the unconscious use of a nonidentifiable model is a considerable concern. To address this problem, an improved variant of a recently published Markov-chain Monte Carlo method is presented. The algorithm is tested using test data as well as experimental data. We demonstrate that, opposed to a widely used maximum-likelihood estimator, it gives clear warning signs when a nonidentifiable model is used for fitting. Furthermore, for test data that was generated from a nonidentifiable model, the Markov-chain Monte Carlo results recover much more information from the data than maximum-likelihood estimation.

Comparison of models for IP3 receptor kinetics using stochastic simulations

Inositol 1,4,5-trisphosphate receptor (IP₃R) is a ubiquitous intracellular calcium (Ca²⁺) channel which has a major role in controlling Ca²⁺ levels in neurons. A variety of computational models have been developed to describe the kinetic function of IP₃R under different conditions. In the field of computational neuroscience, it is of great interest to apply the existing models of IP₃R when modeling local Ca²⁺ transients in dendrites or overall Ca²⁺ dynamics in large neuronal models. The goal of this study was to evaluate existing IP₃R models, based on electrophysiological data. This was done in order to be able to suggest suitable models for neuronal modeling. Altogether four models (Othmer and Tang, 1993; Dawson etal., 2003; Fraiman and Dawson, 2004; Doi etal., 2005) were selected for a more detailed comparison. The selection was based on the computational efficiency of the models and the type of experimental data that was used in developing the model. The kinetics of all four models were simulated by stochastic means, using the simulation software STEPS, which implements the Gillespie stochastic simulation algorithm. The results show major differences in the statistical properties of model functionality. Of the four compared models, the one by Fraiman and Dawson (2004) proved most satisfactory in producing the specific features of experimental findings reported in literature. To our knowledge, the present study is the first detailed evaluation of IP₃R models using stochastic simulation methods, thus providing an important setting for constructing a new, realistic model of IP₃R channel kinetics for compartmental modeling of neuronal functions. We conclude that the kinetics of IP₃R with different concentrations of Ca²⁺ and IP₃ should be more carefully addressed when new models for IP₃R are developed.

A multi-scale approach to airway hyperresponsiveness: from molecule to organ

Airway hyperresponsiveness (AHR), a characteristic of asthma that involves an excessive reduction in airway caliber, is a complex mechanism reflecting multiple processes that manifest over a large range of length and time scales. At one extreme, molecular interactions determine the force generated by airway smooth muscle (ASM). At the other, the spatially distributed constriction of the branching airways leads to breathing difficulties. Similarly, asthma therapies act at the molecular scale while clinical outcomes are determined by lung function. These extremes are linked by events operating over intermediate scales of length and time. Thus, AHR is an emergent phenomenon that limits our understanding of asthma and confounds the interpretation of studies that address physiological mechanisms over a limited range of scales. A solution is a modular computational model that integrates experimental and mathematical data from multiple scales. This includes, at the molecular scale, kinetics, and force production of actin-myosin contractile proteins during cross-bridge and latch-state cycling; at the cellular scale, Ca²⁺ signaling mechanisms that regulate ASM force production; at the tissue scale, forces acting between contracting ASM and opposing viscoelastic tissue that determine airway narrowing; at the organ scale, the topographic distribution of ASM contraction dynamics that determine mechanical impedance of the lung. At each scale, models are constructed with iterations between theory and experimentation to identify the parameters that link adjacent scales. This modular model establishes algorithms for modeling over a wide range of scales and provides a framework for the inclusion of other responses such as inflammation or therapeutic regimes. The goal is to develop this lung model so that it can make predictions about bronchoconstriction and identify the pathophysiologic mechanisms having the greatest impact on AHR and its therapy.

Potassium-selective block of barium permeation through single KcsA channels

Ba²⁺, a doubly charged analogue of K⁺, specifically blocks K⁺ channels by virtue of electrostatic stabilization in the permeation pathway. Ba²⁺ block is used here as a tool to determine the equilibrium binding affinity for various monovalent cations at specific sites in the selectivity filter of a noninactivating mutant of KcsA. At high concentrations of external K⁺, the block-time distribution is double exponential, marking at least two Ba²⁺ sites in the selectivity filter, in accord with a Ba²⁺-containing crystal structure of KcsA. By analyzing block as a function of extracellular K⁺, we determined the equilibrium dissociation constant of K⁺ and of other monovalent cations at an extracellular site, presumably S1, to arrive at a selectivity sequence for binding at this site: Rb⁺ (3 µM) > Cs⁺ (23 µM) > K⁺ (29 µM) > NH₄⁺ (440 µM) >> Na⁺ and Li⁺ (>1 M). This represents an unusually high selectivity for K⁺ over Na⁺, with |ΔΔG⁰| of at least 7 kcal mol⁻¹. These results fit well with other kinetic measurements of selectivity as well as with the many crystal structures of KcsA in various ionic conditions.

MCMC estimation of Markov models for ion channels

Ion channels are characterized by inherently stochastic behavior which can be represented by continuous-time Markov models (CTMM). Although methods for collecting data from single ion channels are available, translating a time series of open and closed channels to a CTMM remains a challenge. Bayesian statistics combined with Markov chain Monte Carlo (MCMC) sampling provide means for estimating the rate constants of a CTMM directly from single channel data. In this article, different approaches for the MCMC sampling of Markov models are combined. This method, new to our knowledge, detects overparameterizations and gives more accurate results than existing MCMC methods. It shows similar performance as QuB-MIL, which indicates that it also compares well with maximum likelihood estimators. Data collected from an inositol trisphosphate receptor is used to demonstrate how the best model for a given data set can be found in practice.

Cardiac cell modelling: observations from the heart of the cardiac physiome project

In this manuscript we review the state of cardiac cell modelling in the context of international initiatives such as the IUPS Physiome and Virtual Physiological Human Projects, which aim to integrate computational models across scales and physics. In particular we focus on the relationship between experimental data and model parameterisation across a range of model types and cellular physiological systems. Finally, in the context of parameter identification and model reuse within the Cardiac Physiome, we suggest some future priority areas for this field.

Inositol trisphosphate receptor and ion channel models based on single-channel data

The inositol trisphosphate receptor (IPR) plays an important role in controlling the dynamics of intracellular Ca(2+). Single-channel patch-clamp recordings are a typical way to study these receptors as well as other ion channels. Methods for analyzing and using this type of data have been developed to fit Markov models of the receptor. The usual method of parameter fitting is based on maximum-likelihood techniques. However, Bayesian inference and Markov chain Monte Carlo techniques are becoming more popular. We describe the application of the Bayesian methods to real experimental single-channel data in three ion channels: the ryanodine receptor, the K(+) channel, and the IPR. One of the main aims of all three studies was that of model selection with different approaches taken. We also discuss the modeling implications for single-channel data that display different levels of channel activity within one recording.

A kinetic model of the inositol trisphosphate receptor based on single-channel data

In many cell types, the inositol trisphosphate receptor is one of the important components controlling intracellular calcium dynamics, and an understanding of this receptor is necessary for an understanding of calcium oscillations and waves. Based on single-channel data from the type-I inositol trisphosphate receptor, and using a Markov chain Monte Carlo approach, we show that the most complex time-dependent model that can be unambiguously determined from steady-state data is one with three closed states and one open state, and we determine how the rate constants depend on calcium. Because the transitions between these states are complex functions of calcium concentration, each model state must correspond to a group of physical states. We fit two different topologies and find that both models predict that the main effect of [Ca(2+)] is to modulate the probability that the receptor is in a state that is able to open, rather than to modulate the transition rate to the open state.