Markov models for ion channels: versatility versus identifiability and speed

Mathematical models of C-type and N-type inactivating heteromeric voltage gated potassium channels

Voltage gated potassium channels can be composed of either four identical, or different, pore-forming protein subunits. While the voltage gated channels with identical subunits have been extensively studied both physiologically and mathematically, those with multiple subunit types, termed heteromeric channels, have not been. Here we construct, and explore the predictive outputs of, mechanistic models for heteromeric voltage gated potassium channels that possess either N-type or C-type inactivation kinetics. For both types of inactivation, we first build Markov models of four identical pore-forming inactivating subunits. Combining this with previous results regarding non-inactivating heteromeric channels, we are able to define models for heteromeric channels containing both non-inactivating and inactivating subunits of any ratio. We simulate each model through three unique voltage clamp protocols to identify steady state properties. In doing so, we generate predictions about the impact of adding additional inactivating subunits on a total channel's kinetics. We show that while N-type inactivating subunits appear to have a non-linear impact on the level of inactivation the channel experiences, the effect of C-type inactivating subunits is almost linear. Finally, to combat the computational issues of working with a large number of state variables we define model reductions for both types of heteromeric channels. For the N-type heteromers we derive a quasi-steady-state approximation and indicate where the approximation is appropriate. With the C-type heteromers we are able to write an explicit model reduction bringing models of greater than 10 dimensions down to 2.

Complex biophysical changes and reduced neuronal firing in an SCN8A variant associated with developmental delay and epilepsy

Mutations in the SCN8A gene, encoding the voltage-gated sodium channel NaV1.6, are associated with a range of neurodevelopmental syndromes. The p.(Gly1625Arg) (G1625R) mutation was identified in a patient diagnosed with developmental epileptic encephalopathy (DEE). While most of the characterized DEE-associated SCN8A mutations were shown to cause a gain-of-channel function, we show that the G1625R variant, positioned within the S4 segment of domain IV, results in complex effects. Voltage-clamp analyses of NaV1.6G1625R demonstrated a mixture of gain- and loss-of-function properties, including reduced current amplitudes, increased time constant of fast voltage-dependent inactivation, a depolarizing shift in the voltage dependence of activation and inactivation, and increased channel availability with high-frequency repeated depolarization. Current-clamp analyses in transfected cultured neurons revealed that these biophysical properties caused a marked reduction in the number of action potentials when firing was driven by the transfected mutant NaV1.6. Accordingly, computational modeling of mature cortical neurons demonstrated a mild decrease in neuronal firing when mimicking the patients' heterozygous SCN8A expression. Structural modeling of NaV1.6G1625R suggested the formation of a cation-π interaction between R1625 and F1588 within domain IV. Double-mutant cycle analysis revealed that this interaction affects the voltage dependence of inactivation in NaV1.6G1625R. Together, our studies demonstrate that the G1625R variant leads to a complex combination of gain and loss of function biophysical changes that result in an overall mild reduction in neuronal firing, related to the perturbed interaction network within the voltage sensor domain, necessitating personalized multi-tiered analysis for SCN8A mutations for optimal treatment selection.

Quantifying the efficacy of voltage protocols in characterising ion channel kinetics: A novel information-theoretic approach

Voltage-clamp experiments are commonly utilised to characterise cellular ion channel kinetics. In these experiments, cells are stimulated using a known time-varying voltage, referred to as the voltage protocol, and the resulting cellular response, typically in the form of current, is measured. Parameters of models that describe ion channel kinetics are then estimated by solving an inverse problem which aims to minimise the discrepancy between the predicted response of the model and the actual measured cell response. In this paper, a novel framework to evaluate the information content of voltage-clamp protocols in relation to ion channel model parameters is presented. Additional quantitative information metrics that allow for comparisons among various voltage protocols are proposed. These metrics offer a foundation for future optimal design frameworks to devise novel, information-rich protocols. The efficacy of the proposed framework is evidenced through the analysis of seven voltage protocols from the literature. By comparing known numerical results for inverse problems using these protocols with the information-theoretic metrics, the proposed approach is validated. The essential steps of the framework are: (i) generate random samples of the parameters from chosen prior distributions; (ii) run the model to generate model output (current) for all samples; (iii) construct reduced-dimensional representations of the time-varying current output using proper orthogonal decomposition (POD); (iv) estimate information-theoretic metrics such as mutual information, entropy equivalent variance, and conditional mutual information using non-parametric methods; (v) interpret the metrics; for example, a higher mutual information between a parameter and the current output suggests the protocol yields greater information about that parameter, resulting in improved identifiability; and (vi) integrate the information-theoretic metrics into a single quantitative criterion, encapsulating the protocol's efficacy in estimating model parameters.

Empirical Quantification of Predictive Uncertainty Due to Model Discrepancy by Training with an Ensemble of Experimental Designs: An Application to Ion Channel Kinetics

When using mathematical models to make quantitative predictions for clinical or industrial use, it is important that predictions come with a reliable estimate of their accuracy (uncertainty quantification). Because models of complex biological systems are always large simplifications, model discrepancy arises-models fail to perfectly recapitulate the true data generating process. This presents a particular challenge for making accurate predictions, and especially for accurately quantifying uncertainty in these predictions. Experimentalists and modellers must choose which experimental procedures (protocols) are used to produce data used to train models. We propose to characterise uncertainty owing to model discrepancy with an ensemble of parameter sets, each of which results from training to data from a different protocol. The variability in predictions from this ensemble provides an empirical estimate of predictive uncertainty owing to model discrepancy, even for unseen protocols. We use the example of electrophysiology experiments that investigate the properties of hERG potassium channels. Here, 'information-rich' protocols allow mathematical models to be trained using numerous short experiments performed on the same cell. In this case, we simulate data with one model and fit it with a different (discrepant) one. For any individual experimental protocol, parameter estimates vary little under repeated samples from the assumed additive independent Gaussian noise model. Yet parameter sets arising from the same model applied to different experiments conflict-highlighting model discrepancy. Our methods will help select more suitable ion channel models for future studies, and will be widely applicable to a range of biological modelling problems.

The physics of heart rhythm disorders

The global burden caused by cardiovascular disease is substantial, with heart disease representing the most common cause of death around the world. There remains a need to develop better mechanistic models of cardiac function in order to combat this health concern. Heart rhythm disorders, or arrhythmias, are one particular type of disease which has been amenable to quantitative investigation. Here we review the application of quantitative methodologies to explore dynamical questions pertaining to arrhythmias. We begin by describing single-cell models of cardiac myocytes, from which two and three dimensional models can be constructed. Special focus is placed on results relating to pattern formation across these spatially-distributed systems, especially the formation of spiral waves of activation. Next, we discuss mechanisms which can lead to the initiation of arrhythmias, focusing on the dynamical state of spatially discordant alternans, and outline proposed mechanisms perpetuating arrhythmias such as fibrillation. We then review experimental and clinical results related to the spatio-temporal mapping of heart rhythm disorders. Finally, we describe treatment options for heart rhythm disorders and demonstrate how statistical physics tools can provide insights into the dynamics of heart rhythm disorders.

Ion channel model reduction using manifold boundaries

Mathematical models of voltage-gated ion channels are used in basic research, industrial and clinical settings. These models range in complexity, but typically contain numerous variables representing the proportion of channels in a given state, and parameters describing the voltage-dependent rates of transition between states. An open problem is selecting the appropriate degree of complexity and structure for an ion channel model given data availability. Here, we simplify a model of the cardiac human Ether-à-go-go related gene (hERG) potassium ion channel, which carries cardiac IKr, using the manifold boundary approximation method (MBAM). The MBAM approximates high-dimensional model-output manifolds by reduced models describing their boundaries, resulting in models with fewer parameters (and often variables). We produced a series of models of reducing complexity starting from an established five-state hERG model with 15 parameters. Models with up to three fewer states and eight fewer parameters were shown to retain much of the predictive capability of the full model and were validated using experimental hERG1a data collected in HEK293 cells at 37°C. The method provides a way to simplify complex models of ion channels that improves parameter identifiability and will aid in future model development.

A Quantitative Systems Pharmacology Perspective on the Importance of Parameter Identifiability

There is an inherent tension in Quantitative Systems Pharmacology (QSP) between the need to incorporate mathematical descriptions of complex physiology and drug targets with the necessity of developing robust, predictive and well-constrained models. In addition to this, there is no “gold standard” for model development and assessment in QSP. Moreover, there can be confusion over terminology such as model and parameter identifiability; complex and simple models; virtual populations; and other concepts, which leads to potential miscommunication and misapplication of methodologies within modeling communities, both the QSP community and related disciplines. This perspective article highlights the pros and cons of using simple (often identifiable) vs. complex (more physiologically detailed but often non-identifiable) models, as well as aspects of parameter identifiability, sensitivity and inference methodologies for model development and analysis. The paper distills the central themes of the issue of identifiability and optimal model size and discusses open challenges.

Creating Ion Channel Kinetic Models Using Cloud Computing

Computational modeling of ion channels provides key insight into experimental electrophysiology results and can be used to connect channel dynamics to emergent phenomena observed at the tissue and organ levels. However, creation of these models requires substantial mathematical and computational background. This tutorial seeks to lower the barrier to creating these models by providing an automated pipeline for creating Markov models of an ion channel kinetics dataset. We start by detailing how to encode sample voltage-clamp protocols and experimental data into the program and its implementation in a cloud computing environment. We guide the reader on how to build a containerized instance, push the machine image, and finally run the routine on cluster nodes. While providing open-source code has become more standard in computational studies, this tutorial provides unprecedented detail on the use of the program and the creation of channel models, starting from inputting the raw experimental data. © 2022 Wiley Periodicals LLC. Basic Protocol: Creation of ion channel kinetic models with a cloud computing environment Alternate Protocol: Instructions for use in a standard high-performance compute cluster.

Ion Channel Modeling beyond State of the Art: A Comparison with a System Theory-Based Model of the Shaker-Related Voltage-Gated Potassium Channel Kv1.1

The mathematical modeling of ion channel kinetics is an important tool for studying the electrophysiological mechanisms of the nerves, heart, or cancer, from a single cell to an organ. Common approaches use either a Hodgkin-Huxley (HH) or a hidden Markov model (HMM) description, depending on the level of detail of the functionality and structural changes of the underlying channel gating, and taking into account the computational effort for model simulations. Here, we introduce for the first time a novel system theory-based approach for ion channel modeling based on the concept of transfer function characterization, without a priori knowledge of the biological system, using patch clamp measurements. Using the shaker-related voltage-gated potassium channel Kv1.1 (KCNA1) as an example, we compare the established approaches, HH and HMM, with the system theory-based concept in terms of model accuracy, computational effort, the degree of electrophysiological interpretability, and methodological limitations. This highly data-driven modeling concept offers a new opportunity for the phenomenological kinetic modeling of ion channels, exhibiting exceptional accuracy and computational efficiency compared to the conventional methods. The method has a high potential to further improve the quality and computational performance of complex cell and organ model simulations, and could provide a valuable new tool in the field of next-generation in silico electrophysiology.

Engaging biological oscillators through second messenger pathways permits emergence of a robust gastric slow-wave during peristalsis

Peristalsis, the coordinated contraction—relaxation of the muscles of the stomach is important for normal gastric motility and is impaired in motility disorders. Coordinated electrical depolarizations that originate and propagate within a network of interconnected layers of interstitial cells of Cajal (ICC) and smooth muscle (SM) cells of the stomach wall as a slow-wave, underly peristalsis. Normally, the gastric slow-wave oscillates with a single period and uniform rostrocaudal lag, exhibiting network entrainment. Understanding of the integrative role of neurotransmission and intercellular coupling in the propagation of an entrained gastric slow-wave, important for understanding motility disorders, however, remains incomplete. Using a computational framework constituted of a novel gastric motility network (GMN) model we address the hypothesis that engaging biological oscillators (i.e., ICCs) by constitutive gap junction coupling mechanisms and enteric neural innervation activated signals can confer a robust entrained gastric slow-wave. We demonstrate that while a decreasing enteric neural innervation gradient that modulates the intracellular IP₃ concentration in the ICCs can guide the aboral slow-wave propagation essential for peristalsis, engaging ICCs by recruiting the exchange of second messengers (inositol trisphosphate (IP₃) and Ca²⁺) ensures a robust entrained longitudinal slow-wave, even in the presence of biological variability in electrical coupling strengths. Our GMN with the distinct intercellular coupling in conjunction with the intracellular feedback pathways and a rostrocaudal enteric neural innervation gradient allows gastric slow waves to oscillate with a moderate range of frequencies and to propagate with a broad range of velocities, thus preventing decoupling observed in motility disorders. Overall, the findings provide a mechanistic explanation for the emergence of decoupled slow waves associated with motility impairments of the stomach, offer directions for future experiments and theoretical work, and can potentially aid in the design of new interventional pharmacological and neuromodulation device treatments for addressing gastric motility disorders.

The coordinated contraction and relaxation of the muscles of the stomach, known as peristalsis is important for normal gastric motility and primarily governed by electrical depolarizations that originate and propagate within a network of interconnected layers of interstitial cells of Cajal (ICCs) and smooth muscle cells of the stomach wall as a slow-wave. Under normal conditions, a gastric slow-wave oscillates with a single period and uniform rostrocaudal lag, exhibiting network entrainment. However, the understanding of intrinsic and extrinsic mechanisms that ensure propagation of a robust entrained slow-wave remains incomplete. Here, using a computational framework, we show that in conjunction with an enteric neural innervation gradient along the rostrocaudal ICC chain, and intercellular electrical coupling, the intercellular exchange of inositol trisphosphate between ICCs prevents decoupling by extending the longitudinal entrainment range along the stomach wall, even when variability in intercellular coupling exists. The findings from our study indicate ways that ensure the rostrocaudal spread of a robust gastric slow-wave and provide a mechanistic explanation for the emergence of decoupled slow waves associated with motility impairments of the stomach.

Identification of structures for ion channel kinetic models

Markov models of ion channel dynamics have evolved as experimental advances have improved our understanding of channel function. Past studies have examined limited sets of various topologies for Markov models of channel dynamics. We present a systematic method for identification of all possible Markov model topologies using experimental data for two types of native voltage-gated ion channel currents: mouse atrial sodium currents and human left ventricular fast transient outward potassium currents. Successful models identified with this approach have certain characteristics in common, suggesting that aspects of the model topology are determined by the experimental data. Incorporating these channel models into cell and tissue simulations to assess model performance within protocols that were not used for training provided validation and further narrowing of the number of acceptable models. The success of this approach suggests a channel model creation pipeline may be feasible where the structure of the model is not specified a priori.

Markov models of ion channel dynamics have evolved as experimental advances have improved our understanding of channel function. Past studies have examined limited sets of various structures for Markov models of channel dynamics. Here, we present a computational routine designed to thoroughly search for Markov model topologies for simulating whole-cell currents. We tested this method on two distinct types of voltage-gated cardiac ion channels and found the number of states and connectivity required to recapitulate experimentally observed kinetics. Successful models identified with this approach have certain characteristics in common, suggesting that model structures are determined by the experimental data. Incorporation of these models into higher scale action potential and cable (an approximation of one-dimensional action potential propagation) simulations, identified key channel phenomena that were required for proper function. These methods provide a route to create functional channel models that can be used for action potential simulation without pre-defining their structure ahead of time.

Independent Markov decomposition: Toward modeling kinetics of biomolecular complexes

To advance the mission of in silico cell biology, modeling the interactions of large and complex biological systems becomes increasingly relevant. The combination of molecular dynamics (MD) simulations and Markov state models (MSMs) has enabled the construction of simplified models of molecular kinetics on long timescales. Despite its success, this approach is inherently limited by the size of the molecular system. With increasing size of macromolecular complexes, the number of independent or weakly coupled subsystems increases, and the number of global system states increases exponentially, making the sampling of all distinct global states unfeasible. In this work, we present a technique called independent Markov decomposition (IMD) that leverages weak coupling between subsystems to compute a global kinetic model without requiring the sampling of all combinatorial states of subsystems. We give a theoretical basis for IMD and propose an approach for finding and validating such a decomposition. Using empirical few-state MSMs of ion channel models that are well established in electrophysiology, we demonstrate that IMD models can reproduce experimental conductance measurements with a major reduction in sampling compared with a standard MSM approach. We further show how to find the optimal partition of all-atom protein simulations into weakly coupled subunits.

Computational prediction of drug response in short QT syndrome type 1 based on measurements of compound effect in stem cell-derived cardiomyocytes

Short QT (SQT) syndrome is a genetic cardiac disorder characterized by an abbreviated QT interval of the patient's electrocardiogram. The syndrome is associated with increased risk of arrhythmia and sudden cardiac death and can arise from a number of ion channel mutations. Cardiomyocytes derived from induced pluripotent stem cells generated from SQT patients (SQT hiPSC-CMs) provide promising platforms for testing pharmacological treatments directly in human cardiac cells exhibiting mutations specific for the syndrome. However, a difficulty is posed by the relative immaturity of hiPSC-CMs, with the possibility that drug effects observed in SQT hiPSC-CMs could be very different from the corresponding drug effect in vivo. In this paper, we apply a multistep computational procedure for translating measured drug effects from these cells to human QT response. This process first detects drug effects on individual ion channels based on measurements of SQT hiPSC-CMs and then uses these results to estimate the drug effects on ventricular action potentials and QT intervals of adult SQT patients. We find that the procedure is able to identify IC50 values in line with measured values for the four drugs quinidine, ivabradine, ajmaline and mexiletine. In addition, the predicted effect of quinidine on the adult QT interval is in good agreement with measured effects of quinidine for adult patients. Consequently, the computational procedure appears to be a useful tool for helping predicting adult drug responses from pure in vitro measurements of patient derived cell lines.

Modeling the Interactions Between Sodium Channels Provides Insight Into the Negative Dominance of Certain Channel Mutations

BACKGROUND

Nav1.5 cardiac Na+ channel mutations can cause arrhythmogenic syndromes. Some of these mutations exert a dominant negative effect on wild-type channels. Recent studies showed that Na+ channels can dimerize, allowing coupled gating. This leads to the hypothesis that allosteric interactions between Na+ channels modulate their function and that these interactions may contribute to the negative dominance of certain mutations.

METHODS

To investigate how allosteric interactions affect microscopic and macroscopic channel function, we developed a modeling paradigm in which Markovian models of two channels are combined. Allosteric interactions are incorporated by modifying the free energies of the composite states and/or barriers between states.

RESULTS

Simulations using two generic 2-state models (C-O, closed-open) revealed that increasing the free energy of the composite states CO/OC leads to coupled gating. Simulations using two 3-state models (closed-open-inactivated) revealed that coupled closings must also involve interactions between further composite states. Using two 6-state cardiac Na+ channel models, we replicated previous experimental results mainly by increasing the energies of the CO/OC states and lowering the energy barriers between the CO/OC and the CO/OO states. The channel model was then modified to simulate a negative dominant mutation (Nav1.5 p.L325R). Simulations of homodimers and heterodimers in the presence and absence of interactions showed that the interactions with the variant channel impair the opening of the wild-type channel and thus contribute to negative dominance.

CONCLUSION

Our new modeling framework recapitulates qualitatively previous experimental observations and helps identifying possible interaction mechanisms between ion channels.

An Understandable, Extensible, and Reusable Implementation of the Hodgkin-Huxley Equations Using Modelica

The Hodgkin-Huxley model of the squid giant axon has been used for decades as the basis of many action potential models. These models are usually communicated using just a list of equations or a circuit diagram, which makes them unnecessarily complicated both for novices and for experts. We present a modular version of the Hodgkin-Huxley model that is more understandable than the usual monolithic implementations and that can be easily reused and extended. Our model is written in Modelica using software engineering concepts, such as object orientation and inheritance. It retains the electrical analogy, but names and explains individual components in biological terms. We use cognitive load theory to measure understandability as the amount of items that have to be kept in working memory simultaneously. The model is broken down into small self-contained components in human-readable code with extensive documentation. Additionally, it features a hybrid diagram that uses biological symbols in an electrical circuit and that is directly tied to the model code. The new model design avoids many redundancies and reduces the cognitive load associated with understanding the model by a factor of 6. Extensions can be easily applied due to an unifying interface and inheritance from shared base classes. The model can be used in an educational context as a more approachable introduction to mathematical modeling in electrophysiology. Additionally the modeling approach and the base components can be used to make complex Hodgkin-Huxley-type models more understandable and reusable.

Deterministic modeling of single-channel and whole-cell currents

As a complement to the experimental work, mathematical models are extensively used to study the functional properties of ionic channels. Even though it is generally assumed that the gating of ionic channels is a Markovian phenomenon, reports based on non-traditional analyses of experimental recordings suggest that non-Markovian processes might be also present. While the stochastic Markov models are by far the most adopted approach for the modeling of ionic channels, a model based on the idea of a deterministic process underlying the gating of ionic channels was proposed by Liebovitch and Toth (Liebovitch, L.S. and Toth, T.I., 1991. Journal of Theoretical Biology, 148(2), pp.243-267.) Here, by using a voltage-dependent K⁺ channel as a first approximation, we propose a modified version of the deterministic model of Liebovitch and Toth that, in addition to reproducing the single-channel currents simulated by a two-states Markov model, it is capable of reproducing the whole-cell currents produced by a population of K⁺ channels.

A Mathematical Model of the Mouse Atrial Myocyte With Inter-Atrial Electrophysiological Heterogeneity

Biophysically detailed mathematical models of cardiac electrophysiology provide an alternative to experimental approaches for investigating possible ionic mechanisms underlying the genesis of electrical action potentials and their propagation through the heart. The aim of this study was to develop a biophysically detailed mathematical model of the action potentials of mouse atrial myocytes, a popular experimental model for elucidating molecular and cellular mechanisms of arrhythmogenesis. Based on experimental data from isolated mouse atrial cardiomyocytes, a set of mathematical equations for describing the biophysical properties of membrane ion channel currents, intracellular Ca²⁺ handling, and Ca²⁺-calmodulin activated protein kinase II and β-adrenergic signaling pathways were developed. Wherever possible, membrane ion channel currents were modeled using Markov chain formalisms, allowing detailed representation of channel kinetics. The model also considered heterogeneous electrophysiological properties between the left and the right atrial cardiomyocytes. The developed model was validated by its ability to reproduce the characteristics of action potentials and Ca²⁺ transients, matching quantitatively to experimental data. Using the model, the functional roles of four K⁺ channel currents in atrial action potential were evaluated by channel block simulations, results of which were quantitatively in agreement with existent experimental data. To conclude, this newly developed model of mouse atrial cardiomyocytes provides a powerful tool for investigating possible ion channel mechanisms of atrial electrical activity at the cellular level and can be further used to investigate mechanisms underlying atrial arrhythmogenesis.

Combination of "generalized Trotter operator splitting" and "quadratic adaptive algorithm" method for tradeoff among speedup, stability, and accuracy in the Markov chain model of sodium ion channels in the ventricular cell model

The fast hybrid operator splitting (HOS) and stable uniformization (UNI) methods have been proposed to save computation cost and enhance stability for Markov chain model in cardiac cell simulations. Moreover, Chen-Chen-Luo's quadratic adaptive algorithm (CCL) combined with HOS or UNI was used to improve the tradeoff between speedup and stability, but without considering accuracy. To compromise among stability, acceleration, and accuracy, we propose a generalized Trotter operator splitting (GTOS) method combined with CCL independent of the asymptotic property of a particular ion-channel model. Due to the accuracy underestimation of the mixed root mean square error (MRMSE) method, threshold root mean square error (TRMSE) is proposed to evaluate computation accuracy. With the fixed time-step RK4 as a reference, the second-order GTOS combined with CCL (30.8-fold speedup) for the wild-type Markov chain model with nine states (WT-9 model) or (7.4-fold) for the wild-type Markov chain model with eight states (WT-8 model) is faster than UNI combined with CCL (15.6-fold) for WT-9 model or (1.2-fold) for WT-8 model, separately. Besides, the second-order GTOS combined with CCL has 3.81% TRMSE for WT-9 model or 4.32% TRMSE for WT-8 model more accurate than 72.43% TRMSE for WT-9 model or 136.17% TRMSE for WT-8 model of HOS combined with CCL. To compromise speedup and accuracy, low-order GTOS combined with CCL is suggested to have the advantages of high precision and low computation cost. For high-accuracy requirements, high-order GTOS combined with CCL is recommended. Graphical abstract.

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.

Creation and application of virtual patient cohorts of heart models

Patient-specific cardiac models are now being used to guide therapies. The increased use of patient-specific cardiac simulations in clinical care will give rise to the development of virtual cohorts of cardiac models. These cohorts will allow cardiac simulations to capture and quantify inter-patient variability. However, the development of virtual cohorts of cardiac models will require the transformation of cardiac modelling from small numbers of bespoke models to robust and rapid workflows that can create large numbers of models. In this review, we describe the state of the art in virtual cohorts of cardiac models, the process of creating virtual cohorts of cardiac models, and how to generate the individual cohort member models, followed by a discussion of the potential and future applications of virtual cohorts of cardiac models. This article is part of the theme issue 'Uncertainty quantification in cardiac and cardiovascular modelling and simulation'.

A distinct molecular mechanism by which phenytoin rescues a novel long QT 3 variant

BACKGROUND

Genetic variants in SCN5A can result in channelopathies such as the long QT syndrome type 3 (LQT3), but the therapeutic response to Na⁺ channel blockers can vary. We previously reported a case of an infant with malignant LQT3 and a missense Q1475P SCN5A variant, who was effectively treated with phenytoin, but only partially with mexiletine. Here, we functionally characterized this variant and investigated possible mechanisms for the differential drug actions.

METHODS

Wild-type or mutant Na_v1.5 cDNAs were examined in transfected HEK293 cells with patch clamping and biochemical assays. We used computational modeling to provide insights into altered channel kinetics and to predict effects on the action potential.

RESULTS

The Q1475P variant in Na_v1.5 reduced the current density and channel surface expression, characteristic of a trafficking defect. The variant also led to positive shifts in the voltage dependence of steady-state activation and inactivation, faster inactivation and recovery from inactivation, and increased the "late" Na⁺ current. Simulations of Na_v1.5 gating with a 9-state Markov model suggested that transitions from inactivated to closed states were accelerated in Q1475P channels, leading to accumulation of channels in non-inactivated closed states. Simulations with a human ventricular myocyte model predicted action potential prolongation with Q1475P, compared with wild type, channels. Patch clamp data showed that mexiletine and phenytoin similarly rescued some of the gating defects. Chronic incubation with mexiletine, but not phenytoin, rescued the Na_v1.5-Q1475P trafficking defect, thus increasing mutant channel expression.

CONCLUSIONS

The gain-of-function effects of Na_v1.5-Q1475P predominate to cause a malignant long QT phenotype. Phenytoin partially corrects the gating defect without restoring surface expression of the mutant channel, whereas mexiletine restores surface expression of the mutant channel, which may explain the lack of efficacy of mexiletine when compared to phenytoin. Our data makes a case for experimental studies before embarking on a one-for-all therapy of arrhythmias.

Slowing of the Time Course of Acidification Decreases the Acid-Sensing Ion Channel 1a Current Amplitude and Modulates Action Potential Firing in Neurons

Acid-sensing ion channels (ASICs) are H⁺-activated neuronal Na⁺ channels. They are involved in fear behavior, learning, neurodegeneration after ischemic stroke and in pain sensation. ASIC activation has so far been studied only with fast pH changes, although the pH changes associated with many roles of ASICs are slow. It is currently not known whether slow pH changes can open ASICs at all. Here, we investigated to which extent slow pH changes can activate ASIC1a channels and induce action potential signaling. To this end, ASIC1a current amplitudes and charge transport in transfected Chinese hamster ovary cells, and ASIC-mediated action potential signaling in cultured cortical neurons were measured in response to defined pH ramps of 1-40 s duration from pH 7.4 to pH 6.6 or 6.0. A kinetic model of the ASIC1a current was developed and integrated into the Hodgkin-Huxley action potential model. Interestingly, whereas the ASIC1a current amplitude decreased with slower pH ramps, action potential firing was higher upon intermediate than fast acidification in cortical neurons. Indeed, fast pH changes (<4 s) induced short action potential bursts, while pH changes of intermediate speed (4-10 s) induced longer bursts. Slower pH changes (>10 s) did in many experiments not generate action potentials. Computer simulations corroborated these observations. We provide here the first description of ASIC function in response to defined slow pH changes. Our study shows that ASIC1a currents, and neuronal activity induced by ASIC1a currents, strongly depend on the speed of pH changes. Importantly, with pH changes that take >10 s to complete, ASIC1a activation is inefficient. Therefore, it is likely that currently unknown modulatory mechanisms allow ASIC activity in situations such as ischemia and inflammation.

Multi-scale approaches for the simulation of cardiac electrophysiology: II - Tissue-level structure and function

Computational models of the heart, from cell-level models, through one-, two- and three-dimensional tissue-level simplifications, to biophysically-detailed three-dimensional models of the ventricles, atria or whole heart, allow the simulation of excitation and propagation of this excitation, and have provided remarkable insight into the normal and pathological functioning of the heart. In this article we present equations for modelling cellular excitation (i.e. the cell action potential) from both a phenomenological and a biophysical perspective. Hodgkin-Huxley formalism is discussed, along with the current generation of biophysically-detailed cardiac cell models. Alternative Markovian formulations for modelling ionic currents are also presented. Equations describing propagation of this cellular excitation, through one-, two- and three-dimensional idealised or realistic tissues, are then presented. For all types of model, from cell to tissue, methods for discretisation and integration of the underlying equations are discussed. The article finishes with a discussion of two tissue-level experimental imaging techniques - diffusion tensor magnetic resonance imaging and optical imaging - that can be used to provide data for parameterisation and validation of cell- and tissue-level cardiac models.

Inference-based assessment of parameter identifiability in nonlinear biological models

As systems approaches to the development of biological models become more mature, attention is increasingly focusing on the problem of inferring parameter values within those models from experimental data. However, particularly for nonlinear models, it is not obvious, either from inspection of the model or from the experimental data, that the inverse problem of parameter fitting will have a unique solution, or even a non-unique solution that constrains the parameters to lie within a plausible physiological range. Where parameters cannot be constrained they are termed 'unidentifiable'. We focus on gaining insight into the causes of unidentifiability using inference-based methods, and compare a recently developed measure-theoretic approach to inverse sensitivity analysis to the popular Markov chain Monte Carlo and approximate Bayesian computation techniques for Bayesian inference. All three approaches map the uncertainty in quantities of interest in the output space to the probability of sets of parameters in the input space. The geometry of these sets demonstrates how unidentifiability can be caused by parameter compensation and provides an intuitive approach to inference-based experimental design.

Four Ways to Fit an Ion Channel Model

Mathematical models of ionic currents are used to study the electrophysiology of the heart, brain, gut, and several other organs. Increasingly, these models are being used predictively in the clinic, for example, to predict the risks and results of genetic mutations, pharmacological treatments, or surgical procedures. These safety-critical applications depend on accurate characterization of the underlying ionic currents. Four different methods can be found in the literature to fit voltage-sensitive ion channel models to whole-cell current measurements: method 1, fitting model equations directly to time-constant, steady-state, and I-V summary curves; method 2, fitting by comparing simulated versions of these summary curves to their experimental counterparts; method 3, fitting to the current traces themselves from a range of protocols; and method 4, fitting to a single current trace from a short and rapidly fluctuating voltage-clamp protocol. We compare these methods using a set of experiments in which hERG1a current was measured in nine Chinese hamster ovary cells. In each cell, the same sequence of fitting protocols was applied, as well as an independent validation protocol. We show that methods 3 and 4 provide the best predictions on the independent validation set and that short, rapidly fluctuating protocols like that used in method 4 can replace much longer conventional protocols without loss of predictive ability. Although data for method 2 are most readily available from the literature, we find it performs poorly compared to methods 3 and 4 both in accuracy of predictions and computational efficiency. Our results demonstrate how novel experimental and computational approaches can improve the quality of model predictions in safety-critical applications.

Rapid Characterization of hERG Channel Kinetics II: Temperature Dependence

Ion channel behavior can depend strongly on temperature, with faster kinetics at physiological temperatures leading to considerable changes in currents relative to room temperature. These temperature-dependent changes in voltage-dependent ion channel kinetics (rates of opening, closing, inactivating, and recovery) are commonly represented with Q₁₀ coefficients or an Eyring relationship. In this article, we assess the validity of these representations by characterizing channel kinetics at multiple temperatures. We focus on the human Ether-à-go-go-Related Gene (hERG) channel, which is important in drug safety assessment and commonly screened at room temperature so that results require extrapolation to physiological temperature. In Part I of this study, we established a reliable method for high-throughput characterization of hERG1a (Kv11.1) kinetics, using a 15-second information-rich optimized protocol. In this Part II, we use this protocol to study the temperature dependence of hERG kinetics using Chinese hamster ovary cells overexpressing hERG1a on the Nanion SyncroPatch 384PE, a 384-well automated patch-clamp platform, with temperature control. We characterize the temperature dependence of hERG gating by fitting the parameters of a mathematical model of hERG kinetics to data obtained at five distinct temperatures between 25 and 37°C and validate the models using different protocols. Our models reveal that activation is far more temperature sensitive than inactivation, and we observe that the temperature dependency of the kinetic parameters is not represented well by Q₁₀ coefficients; it broadly follows a generalized, but not the standardly-used, Eyring relationship. We also demonstrate that experimental estimations of Q₁₀ coefficients are protocol dependent. Our results show that a direct fit using our 15-s protocol best represents hERG kinetics at any given temperature and suggests that using the Generalized Eyring theory is preferable if no experimental data are available to derive model parameters at a given temperature.

Rapid Characterization of hERG Channel Kinetics I: Using an Automated High-Throughput System

Predicting how pharmaceuticals may affect heart rhythm is a crucial step in drug development and requires a deep understanding of a compound’s action on ion channels. In vitro hERG channel current recordings are an important step in evaluating the proarrhythmic potential of small molecules and are now routinely performed using automated high-throughput patch-clamp platforms. These machines can execute traditional voltage-clamp protocols aimed at specific gating processes, but the array of protocols needed to fully characterize a current is typically too long to be applied in a single cell. Shorter high-information protocols have recently been introduced that have this capability, but they are not typically compatible with high-throughput platforms. We present a new 15 second protocol to characterize hERG (Kv11.1) kinetics, suitable for both manual and high-throughput systems. We demonstrate its use on the Nanion SyncroPatch 384PE, a 384-well automated patch-clamp platform, by applying it to Chinese hamster ovary cells stably expressing hERG1a. From these recordings, we construct 124 cell-specific variants/parameterizations of a hERG model at 25°C. A further eight independent protocols are run in each cell and are used to validate the model predictions. We then combine the experimental recordings using a hierarchical Bayesian model, which we use to quantify the uncertainty in the model parameters, and their variability from cell-to-cell; we use this model to suggest reasons for the variability. This study demonstrates a robust method to measure and quantify uncertainty and shows that it is possible and practical to use high-throughput systems to capture full hERG channel kinetics quantitatively and rapidly.

Detecting undetectables: Can conductances of action potential models be changed without appreciable change in the transmembrane potential?

Mathematical models describing the dynamics of the cardiac action potential are of great value for understanding how changes to the system can disrupt the normal electrical activity of cells and tissue in the heart. However, to represent specific data, these models must be parameterized, and adjustment of the maximum conductances of the individual contributing ionic currents is a commonly used method. Here, we present a method for investigating the uniqueness of such resulting parameterizations. Our key question is: Can the maximum conductances of a model be changed without giving any appreciable changes in the action potential? If so, the model parameters are not unique and this poses a major problem in using the models to identify changes in parameters from data, for instance, to evaluate potential drug effects. We propose a method for evaluating this uniqueness, founded on the singular value decomposition of a matrix consisting of the individual ionic currents. Small singular values of this matrix signify lack of parameter uniqueness and we show that the conclusion from linear analysis of the matrix carries over to provide insight into the uniqueness of the parameters in the nonlinear case. Using numerical experiments, we quantify the identifiability of the maximum conductances of well-known models of the cardiac action potential. Furthermore, we show how the identifiability depends on the time step used in the observation of the currents, how the application of drugs may change identifiability, and, finally, how the stimulation protocol can be used to improve the identifiability of a model.

A Mathematical Model of the Human Cardiac Na+ Channel

Sodium ion channel is a membrane protein that plays an important role in excitable cells, as it is responsible for the initiation of action potentials. Understanding the electrical characteristics of sodium channels is essential in predicting their behavior under different physiological conditions. We investigated several Markov models for the human cardiac sodium channel Na_V1.5 to derive a minimal mathematical model that describes the reported experimental data obtained using major voltage clamp protocols. We obtained simulation results for peak current-voltage relationships, the voltage dependence of normalized ion channel conductance, steady-state inactivation, activation and deactivation kinetics, fast and slow inactivation kinetics, and recovery from inactivation kinetics. Good agreement with the experimental data provides us with the mechanisms of the fast and slow inactivation of the human sodium channel and the coupling of its inactivation states to the closed and open states in the activation pathway.

Complexity reduction in human atrial modeling using extended Kalman filter

Human atrial tissue electrophysiology is modeled upon biophysical details obtained from cellular level measurements. Data collected for this purpose typically represent a unique state of the tissue. As reproducing dynamic cases such as subject-varied and/or disorder-varied electrophysiological properties is in question, such complex models are typically hard to use. Hence, there is a need for simpler yet biophysically accurate and mathematically tractable models to be used for case-specific reproductions and simulations. In this study, a scheme for parameter estimation of a phenomenological cardiac model to match a targeted behavior generated from a complex model is used. Specifically, an algorithm incorporating extended Kalman filter (EKF) into the scheme is proposed. Its performance is then compared to that of particle swarm optimization (PSO) and sequential quadratic programming (SQP), algorithms that have been widely used for parameter optimization. Both robustness and adaptability performance of the algorithms are tested through various designs. For this, reproducing action potential (AP) waveforms of varying remodeling states of atrial fibrillation (AF) at different stimulus protocols was targeted. Also, randomly generated initial parameter sets are included in the tests. In addition, AP duration (APD) restitution curve (RC) is used for a multiscale evaluation of fitting performance. Finally, wavefront propagation on 2D of a selected AF remodeling state using parameter solutions from each of the algorithms is simulated for a qualitative evaluation. In general, PSO yielded superior performances than EKF and SQP with respect to fitting AP waveforms. Considering both AP and APD RC, however, EKF yielded the best accuracies. Also, more accurate spiral wave reentry is obtained with EKF. Overall, EKF algorithm yielded the best performance in robustness and adaptability. Graphical Abstract.

Bond graph modelling of the cardiac action potential: implications for drift and non-unique steady states

Mathematical models of cardiac action potentials have become increasingly important in the study of heart disease and pharmacology, but concerns linger over their robustness during long periods of simulation, in particular due to issues such as model drift and non-unique steady states. Previous studies have linked these to violation of conservation laws, but only explored those issues with respect to charge conservation in specific models. Here, we propose a general and systematic method of identifying conservation laws hidden in models of cardiac electrophysiology by using bond graphs, and develop a bond graph model of the cardiac action potential to study long-term behaviour. Bond graphs provide an explicit energy-based framework for modelling physical systems, which makes them well suited for examining conservation within electrophysiological models. We find that the charge conservation laws derived in previous studies are examples of the more general concept of a 'conserved moiety'. Conserved moieties explain model drift and non-unique steady states, generalizing the results from previous studies. The bond graph approach provides a rigorous method to check for drift and non-unique steady states in a wide range of cardiac action potential models, and can be extended to examine behaviours of other excitable systems.

Reproducible model development in the cardiac electrophysiology Web Lab

The modelling of the electrophysiology of cardiac cells is one of the most mature areas of systems biology. This extended concentration of research effort brings with it new challenges, foremost among which is that of choosing which of these models is most suitable for addressing a particular scientific question. In a previous paper, we presented our initial work in developing an online resource for the characterisation and comparison of electrophysiological cell models in a wide range of experimental scenarios. In that work, we described how we had developed a novel protocol language that allowed us to separate the details of the mathematical model (the majority of cardiac cell models take the form of ordinary differential equations) from the experimental protocol being simulated. We developed a fully-open online repository (which we termed the Cardiac Electrophysiology Web Lab) which allows users to store and compare the results of applying the same experimental protocol to competing models. In the current paper we describe the most recent and planned extensions of this work, focused on supporting the process of model building from experimental data. We outline the necessary work to develop a machine-readable language to describe the process of inferring parameters from wet lab datasets, and illustrate our approach through a detailed example of fitting a model of the hERG channel using experimental data. We conclude by discussing the future challenges in making further progress in this domain towards our goal of facilitating a fully reproducible approach to the development of cardiac cell models.

Comparing two sequential Monte Carlo samplers for exact and approximate Bayesian inference on biological models

Bayesian methods are advantageous for biological modelling studies due to their ability to quantify and characterize posterior variability in model parameters. When Bayesian methods cannot be applied, due either to non-determinism in the model or limitations on system observability, approximate Bayesian computation (ABC) methods can be used to similar effect, despite producing inflated estimates of the true posterior variance. Owing to generally differing application domains, there are few studies comparing Bayesian and ABC methods, and thus there is little understanding of the properties and magnitude of this uncertainty inflation. To address this problem, we present two popular strategies for ABC sampling that we have adapted to perform exact Bayesian inference, and compare them on several model problems. We find that one sampler was impractical for exact inference due to its sensitivity to a key normalizing constant, and additionally highlight sensitivities of both samplers to various algorithmic parameters and model conditions. We conclude with a study of the O'Hara-Rudy cardiac action potential model to quantify the uncertainty amplification resulting from employing ABC using a set of clinically relevant biomarkers. We hope that this work serves to guide the implementation and comparative assessment of Bayesian and ABC sampling techniques in biological models.

Sinusoidal voltage protocols for rapid characterisation of ion channel kinetics

Key points

Ion current kinetics are commonly represented by current–voltage relationships, time constant–voltage relationships and subsequently mathematical models fitted to these. These experiments take substantial time, which means they are rarely performed in the same cell.

Rather than traditional square‐wave voltage clamps, we fitted a model to the current evoked by a novel sum‐of‐sinusoids voltage clamp that was only 8 s long.

Short protocols that can be performed multiple times within a single cell will offer many new opportunities to measure how ion current kinetics are affected by changing conditions.

The new model predicts the current under traditional square‐wave protocols well, with better predictions of underlying currents than literature models. The current under a novel physiologically relevant series of action potential clamps is predicted extremely well.

The short sinusoidal protocols allow a model to be fully fitted to individual cells, allowing us to examine cell–cell variability in current kinetics for the first time.

Abstract

Understanding the roles of ion currents is crucial to predict the action of pharmaceuticals and mutations in different scenarios, and thereby to guide clinical interventions in the heart, brain and other electrophysiological systems. Our ability to predict how ion currents contribute to cellular electrophysiology is in turn critically dependent on our characterisation of ion channel kinetics – the voltage‐dependent rates of transition between open, closed and inactivated channel states. We present a new method for rapidly exploring and characterising ion channel kinetics, applying it to the hERG potassium channel as an example, with the aim of generating a quantitatively predictive representation of the ion current. We fitted a mathematical model to currents evoked by a novel 8 second sinusoidal voltage clamp in CHO cells overexpressing hERG1a. The model was then used to predict over 5 minutes of recordings in the same cell in response to further protocols: a series of traditional square step voltage clamps, and also a novel voltage clamp comprising a collection of physiologically relevant action potentials. We demonstrate that we can make predictive cell‐specific models that outperform the use of averaged data from a number of different cells, and thereby examine which changes in gating are responsible for cell–cell variability in current kinetics. Our technique allows rapid collection of consistent and high quality data, from single cells, and produces more predictive mathematical ion channel models than traditional approaches.

Ion current kinetics are commonly represented by current–voltage relationships, time constant–voltage relationships and subsequently mathematical models fitted to these. These experiments take substantial time, which means they are rarely performed in the same cell.

Rather than traditional square‐wave voltage clamps, we fitted a model to the current evoked by a novel sum‐of‐sinusoids voltage clamp that was only 8 s long.

Short protocols that can be performed multiple times within a single cell will offer many new opportunities to measure how ion current kinetics are affected by changing conditions.

The new model predicts the current under traditional square‐wave protocols well, with better predictions of underlying currents than literature models. The current under a novel physiologically relevant series of action potential clamps is predicted extremely well.

The short sinusoidal protocols allow a model to be fully fitted to individual cells, allowing us to examine cell–cell variability in current kinetics for the first time.

Multiscale Kinetic Modeling Reveals an Ensemble of Cl-/H+ Exchange Pathways in ClC-ec1 Antiporter

Despite several years of research, the ion exchange mechanisms in chloride/proton antiporters and many other coupled transporters are not yet understood at the molecular level. Here, we present a novel approach to kinetic modeling and apply it to ion exchange in ClC-ec1. Our multiscale kinetic model is developed by (1) calculating the state-to-state rate coefficients with reactive and polarizable molecular dynamics simulations, (2) optimizing these rates in a global kinetic network, and (3) predicting new electrophysiological results. The model shows that the robust Cl:H exchange ratio (2.2:1) can indeed arise from kinetic coupling without large protein conformational changes, indicating a possible facile evolutionary connection to chloride channels. The E148 amino acid residue is shown to couple chloride and proton transport through protonation-dependent blockage of the central anion binding site and an anion-dependent pK_a value, which influences proton transport. The results demonstrate how an ensemble of different exchange pathways, as opposed to a single series of transitions, culminates in the macroscopic observables of the antiporter, such as transport rates, chloride/proton stoichiometry, and pH dependence.

Tailoring Mathematical Models to Stem-Cell Derived Cardiomyocyte Lines Can Improve Predictions of Drug-Induced Changes to Their Electrophysiology

Human induced pluripotent stem cell derived cardiomyocytes (iPSC-CMs) have applications in disease modeling, cell therapy, drug screening and personalized medicine. Computational models can be used to interpret experimental findings in iPSC-CMs, provide mechanistic insights, and translate these findings to adult cardiomyocyte (CM) electrophysiology. However, different cell lines display different expression of ion channels, pumps and receptors, and show differences in electrophysiology. In this exploratory study, we use a mathematical model based on iPSC-CMs from Cellular Dynamic International (CDI, iCell), and compare its predictions to novel experimental recordings made with the Axiogenesis Cor.4U line. We show that tailoring this model to the specific cell line, even using limited data and a relatively simple approach, leads to improved predictions of baseline behavior and response to drugs. This demonstrates the need and the feasibility to tailor models to individual cell lines, although a more refined approach will be needed to characterize individual currents, address differences in ion current kinetics, and further improve these results.

Observability of Complex Systems: Finding the Gap

For a reconstruction of state and parameter values in a dynamic system model, first the question whether these values can be uniquely determined from the data must be answered. This structural model property is known as observability or, in case of parameter calibration only, identifiability. Testing a given model for observability is a well studied problem in the systems and control sciences. However, it is increasingly difficult, if not impossible, to address this property for large size models that, nowadays, are frequently used. We demonstrate the application of a recently developed algorithm that overcomes this problem and is remarkably efficient. As an illustration we show how an observability analysis for a Chinese Hamster Ovary Cell model (34 states, 117 parameters), a JAKSTAT signalling model (31 states, 51 parameters), and a MAP Kinase model (100 states, 88 parameters) can be established in a very short time.

Uncertainty Quantification Reveals the Importance of Data Variability and Experimental Design Considerations for in Silico Proarrhythmia Risk Assessment

The Comprehensive in vitro Proarrhythmia Assay (CiPA) is a global initiative intended to improve drug proarrhythmia risk assessment using a new paradigm of mechanistic assays. Under the CiPA paradigm, the relative risk of drug-induced Torsade de Pointes (TdP) is assessed using an in silico model of the human ventricular action potential (AP) that integrates in vitro pharmacology data from multiple ion channels. Thus, modeling predictions of cardiac risk liability will depend critically on the variability in pharmacology data, and uncertainty quantification (UQ) must comprise an essential component of the in silico assay. This study explores UQ methods that may be incorporated into the CiPA framework. Recently, we proposed a promising in silico TdP risk metric (qNet), which is derived from AP simulations and allows separation of a set of CiPA training compounds into Low, Intermediate, and High TdP risk categories. The purpose of this study was to use UQ to evaluate the robustness of TdP risk separation by qNet. Uncertainty in the model parameters used to describe drug binding and ionic current block was estimated using the non-parametric bootstrap method and a Bayesian inference approach. Uncertainty was then propagated through AP simulations to quantify uncertainty in qNet for each drug. UQ revealed lower uncertainty and more accurate TdP risk stratification by qNet when simulations were run at concentrations below 5× the maximum therapeutic exposure (C_max). However, when drug effects were extrapolated above 10× C_max, UQ showed that qNet could no longer clearly separate drugs by TdP risk. This was because for most of the pharmacology data, the amount of current block measured was <60%, preventing reliable estimation of IC₅₀-values. The results of this study demonstrate that the accuracy of TdP risk prediction depends both on the intrinsic variability in ion channel pharmacology data as well as on experimental design considerations that preclude an accurate determination of drug IC₅₀-values in vitro. Thus, we demonstrate that UQ provides valuable information about in silico modeling predictions that can inform future proarrhythmic risk evaluation of drugs under the CiPA paradigm.

Nonequilibrium response of a voltage gated sodium ion channel and biophysical characterization of dynamic hysteresis

Here we have studied the dynamic as well as the non-equilibrium thermodynamic response properties of voltage-gated Na-ion channel. Using sinusoidally oscillating external voltage protocol we have both kinetically and energetically studied the non-equilibrium steady state properties of dynamic hysteresis in details. We have introduced a method of estimating the work done associated with the dynamic memory due to a cycle of oscillating voltage. We have quantitatively characterised the loop area of ionic current which gives information about the work done to sustain the dynamic memory only for ion conduction, while the loop area of total entropy production rate gives the estimate of work done for overall gating dynamics. The maximum dynamic memory of Na-channel not only depends on the frequency and amplitude but it also depends sensitively on the mean of the oscillating voltage and here we have shown how the system optimize the dynamic memory itself in the biophysical range of field parameters. The relation between the average ionic current with increasing frequency corresponds to the nature of the average dissipative work done at steady state. It is also important to understand that the utilization of the energy from the external field can not be directly obtained only from the measurement of ionic current but also requires nonequilibrium thermodynamic study.

A Parsimonious Model of the Rabbit Action Potential Elucidates the Minimal Physiological Requirements for Alternans and Spiral Wave Breakup

Elucidating the underlying mechanisms of fatal cardiac arrhythmias requires a tight integration of electrophysiological experiments, models, and theory. Existing models of transmembrane action potential (AP) are complex (resulting in over parameterization) and varied (leading to dissimilar predictions). Thus, simpler models are needed to elucidate the "minimal physiological requirements" to reproduce significant observable phenomena using as few parameters as possible. Moreover, models have been derived from experimental studies from a variety of species under a range of environmental conditions (for example, all existing rabbit AP models incorporate a formulation of the rapid sodium current, INa, based on 30 year old data from chick embryo cell aggregates). Here we develop a simple "parsimonious" rabbit AP model that is mathematically identifiable (i.e., not over parameterized) by combining a novel Hodgkin-Huxley formulation of INa with a phenomenological model of repolarization similar to the voltage dependent, time-independent rectifying outward potassium current (IK). The model was calibrated using the following experimental data sets measured from the same species (rabbit) under physiological conditions: dynamic current-voltage (I-V) relationships during the AP upstroke; rapid recovery of AP excitability during the relative refractory period; and steady-state INa inactivation via voltage clamp. Simulations reproduced several important "emergent" phenomena including cellular alternans at rates > 250 bpm as observed in rabbit myocytes, reentrant spiral waves as observed on the surface of the rabbit heart, and spiral wave breakup. Model variants were studied which elucidated the minimal requirements for alternans and spiral wave break up, namely the kinetics of INa inactivation and the non-linear rectification of IK.The simplicity of the model, and the fact that its parameters have physiological meaning, make it ideal for engendering generalizable mechanistic insight and should provide a solid "building-block" to generate more detailed ionic models to represent complex rabbit electrophysiology.

Computational cardiology and risk stratification for sudden cardiac death: one of the grand challenges for cardiology in the 21st century

Risk stratification in the context of sudden cardiac death has been acknowledged as one of the major challenges facing cardiology for the past four decades. In recent years, the advent of high performance computing has facilitated organ-level simulation of the heart, meaning we can now examine the causes, mechanisms and impact of cardiac dysfunction in silico. As a result, computational cardiology, largely driven by the Physiome project, now stands at the threshold of clinical utility in regards to risk stratification and treatment of patients at risk of sudden cardiac death. In this white paper, we outline a roadmap of what needs to be done to make this translational step, using the relatively well-developed case of acquired or drug-induced long QT syndrome as an exemplar case.

Uncertainty and variability in computational and mathematical models of cardiac physiology

KEY POINTS

Mathematical and computational models of cardiac physiology have been an integral component of cardiac electrophysiology since its inception, and are collectively known as the Cardiac Physiome. We identify and classify the numerous sources of variability and uncertainty in model formulation, parameters and other inputs that arise from both natural variation in experimental data and lack of knowledge. The impact of uncertainty on the outputs of Cardiac Physiome models is not well understood, and this limits their utility as clinical tools. We argue that incorporating variability and uncertainty should be a high priority for the future of the Cardiac Physiome. We suggest investigating the adoption of approaches developed in other areas of science and engineering while recognising unique challenges for the Cardiac Physiome; it is likely that novel methods will be necessary that require engagement with the mathematics and statistics community.

ABSTRACT

The Cardiac Physiome effort is one of the most mature and successful applications of mathematical and computational modelling for describing and advancing the understanding of physiology. After five decades of development, physiological cardiac models are poised to realise the promise of translational research via clinical applications such as drug development and patient-specific approaches as well as ablation, cardiac resynchronisation and contractility modulation therapies. For models to be included as a vital component of the decision process in safety-critical applications, rigorous assessment of model credibility will be required. This White Paper describes one aspect of this process by identifying and classifying sources of variability and uncertainty in models as well as their implications for the application and development of cardiac models. We stress the need to understand and quantify the sources of variability and uncertainty in model inputs, and the impact of model structure and complexity and their consequences for predictive model outputs. We propose that the future of the Cardiac Physiome should include a probabilistic approach to quantify the relationship of variability and uncertainty of model inputs and outputs.

Reconstruction of Cell Surface Densities of Ion Pumps, Exchangers, and Channels from mRNA Expression, Conductance Kinetics, Whole-Cell Calcium, and Current-Clamp Voltage Recordings, with an Application to Human Uterine Smooth Muscle Cells

Uterine smooth muscle cells remain quiescent throughout most of gestation, only generating spontaneous action potentials immediately prior to, and during, labor. This study presents a method that combines transcriptomics with biophysical recordings to characterise the conductance repertoire of these cells, the ‘conductance repertoire’ being the total complement of ion channels and transporters expressed by an electrically active cell. Transcriptomic analysis provides a set of potential electrogenic entities, of which the conductance repertoire is a subset. Each entity within the conductance repertoire was modeled independently and its gating parameter values were fixed using the available biophysical data. The only remaining free parameters were the surface densities for each entity. We characterise the space of combinations of surface densities (density vectors) consistent with experimentally observed membrane potential and calcium waveforms. This yields insights on the functional redundancy of the system as well as its behavioral versatility. Our approach couples high-throughput transcriptomic data with physiological behaviors in health and disease, and provides a formal method to link genotype to phenotype in excitable systems. We accurately predict current densities and chart functional redundancy. For example, we find that to evoke the observed voltage waveform, the BK channel is functionally redundant whereas hERG is essential. Furthermore, our analysis suggests that activation of calcium-activated chloride conductances by intracellular calcium release is the key factor underlying spontaneous depolarisations.

A well-known problem in electrophysiologal modeling is that the parameters of the gating kinetics of the ion channels cannot be uniquely determined from observed behavior at the cellular level. One solution is to employ simplified “macroscopic” currents that mimic the behavior of aggregates of distinct entities at the protein level. The gating parameters of each channel or pump can be determined by studying it in isolation, leaving the general problem of finding the densities at which the channels occur in the plasma membrane. We propose an approach, which we apply to uterine smooth muscle cells, whereby we constrain the list of possible entities by means of transcriptomics and chart the indeterminacy of the problem in terms of the kernel of the corresponding linear transformation. A graphical representation of this kernel visualises the functional redundancy of the system. We show that the role of certain conductances can be fulfilled, or compensated for, by suitable combinations of other conductances; this is not always the case, and such “non-substitutable” conductances can be regarded as functionally non-redundant. Electrogenic entities belonging to the latter category are suitable putative clinical targets.

Parameterization for In-Silico Modeling of Ion Channel Interactions with Drugs

Since the first Hodgkin and Huxley ion channel model was described in the 1950s, there has been an explosion in mathematical models to describe ion channel function. As experimental data has become richer, models have concomitantly been improved to better represent ion channel kinetic processes, although these improvements have generally resulted in more model complexity and an increase in the number of parameters necessary to populate the models. Models have also been developed to explicitly model drug interactions with ion channels. Recent models of drug-channel interactions account for the discrete kinetics of drug interaction with distinct ion channel state conformations, as it has become clear that such interactions underlie complex emergent kinetics such as use-dependent block. Here, we describe an approach for developing a model for ion channel drug interactions. The method describes the process of extracting rate constants from experimental electrophysiological function data to use as initial conditions for the model parameters. We then describe implementation of a parameter optimization method to refine the model rate constants describing ion channel drug kinetics. The algorithm takes advantage of readily available parallel computing tools to speed up the optimization. Finally, we describe some potential applications of the platform including the potential for gaining fundamental mechanistic insights into ion channel function and applications to in silico drug screening and development.

ModFossa: A library for modeling ion channels using Python

The creation and simulation of ion channel models using continuous-time Markov processes is a powerful and well-used tool in the field of electrophysiology and ion channel research. While several software packages exist for the purpose of ion channel modeling, most are GUI based, and none are available as a Python library. In an attempt to provide an easy-to-use, yet powerful Markov model-based ion channel simulator, we have developed ModFossa, a Python library supporting easy model creation and stimulus definition, complete with a fast numerical solver, and attractive vector graphics plotting.