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

A Reinterpretation of the Relationship Between Persistent and Resurgent Sodium Currents

The resurgent sodium current (INaR) activates on membrane repolarization, such as during the downstroke of neuronal action potentials. Due to its unique activation properties, INaR is thought to drive high rates of repetitive neuronal firing. However, INaR is often studied in combination with the persistent or non-inactivating portion of sodium currents (INaP). We used dynamic clamp to test how INaR and INaP individually affect repetitive firing in adult cerebellar Purkinje neurons from male and female mice. We learned INaR does not scale repetitive firing rates due to its rapid decay at subthreshold voltages, and that subthreshold INaP is critical in regulating neuronal firing rate. Adjustments to the Nav conductance model used in these studies revealed INaP and INaR can be inversely scaled by adjusting occupancy in the slow inactivated kinetic state. Together with additional dynamic clamp experiments, these data suggest the regulation of sodium channel slow inactivation can fine-tune INaP and Purkinje neuron repetitive firing rates.

Evaluating sequential and allosteric activation models in IKs channels with mutated voltage sensors

The ion-conducting IKs channel complex, important in cardiac repolarization and arrhythmias, comprises tetramers of KCNQ1 α-subunits along with 1-4 KCNE1 accessory subunits and calmodulin regulatory molecules. The E160R mutation in individual KCNQ1 subunits was used to prevent activation of voltage sensors and allow direct determination of transition rate data from complexes opening with a fixed number of 1, 2, or 4 activatable voltage sensors. Markov models were used to test the suitability of sequential versus allosteric models of IKs activation by comparing simulations with experimental steady-state and transient activation kinetics, voltage-sensor fluorescence from channels with two or four activatable domains, and limiting slope currents at negative potentials. Sequential Hodgkin-Huxley-type models approximately describe IKs currents but cannot explain an activation delay in channels with only one activatable subunit or the hyperpolarizing shift in the conductance-voltage relationship with more activatable voltage sensors. Incorporating two voltage sensor activation steps in sequential models and a concerted step in opening via rates derived from fluorescence measurements improves models but does not resolve fundamental differences with experimental data. Limiting slope current data that show the opening of channels at negative potentials and very low open probability are better simulated using allosteric models of activation with one transition per voltage sensor, which implies that movement of all four sensors is not required for IKs conductance. Tiered allosteric models with two activating transitions per voltage sensor can fully account for IKs current and fluorescence activation kinetics in constructs with different numbers of activatable voltage sensors.

Computational Modeling of Cardiac Electrophysiology

Mathematical modeling and simulation are well-established and powerful tools to integrate experimental data of individual components of cardiac electrophysiology, excitation-contraction coupling, and regulatory signaling pathways, to gain quantitative and mechanistic insight into pathophysiological processes and guide therapeutic strategies. Here, we briefly describe the processes governing cardiac myocyte electrophysiology and Ca2+ handling and their regulation, as well as action potential propagation in tissue. We discuss the models and methods used to describe these phenomena, including procedures for model parameterization and validation, in addition to protocols for model interrogation and analysis and techniques that account for phenotypic variability and parameter uncertainty. Our objective is to provide a summary of basic concepts and approaches as a resource for scientists training in this discipline and for all researchers aiming to gain an understanding of cardiac modeling studies.

Emerging methods to model cardiac ion channel and myocyte electrophysiology

In the field of cardiac electrophysiology, modeling has played a central role for many decades. However, even though the effort is well-established, it has recently seen a rapid and sustained evolution in the complexity and predictive power of the models being created. In particular, new approaches to modeling have allowed the tracking of parallel and interconnected processes that span from the nanometers and femtoseconds that determine ion channel gating to the centimeters and minutes needed to describe an arrhythmia. The connection between scales has brought unprecedented insight into cardiac arrhythmia mechanisms and drug therapies. This review focuses on the generation of these models from first principles, generation of detailed models to describe ion channel kinetics, algorithms to create and numerically solve kinetic models, and new approaches toward data gathering that parameterize these models. While we focus on application of these models for cardiac arrhythmia, these concepts are widely applicable to model the physiology and pathophysiology of any excitable cell.

Modeling the kinetics of heteromeric potassium channels

Mechanistic mathematical modeling has long been used as a tool for answering questions in cellular physiology. To mathematically describe cellular processes such as cell excitability, volume regulation, neurotransmitter release, and hormone secretion requires accurate descriptions of ion channel kinetics. One class of ion channels currently lacking a physiological model framework is the class of channels built with multiple different potassium protein subunits called heteromeric voltage gated potassium channels. Here we present a novel mathematical model for heteromeric potassium channels that captures both the number and type of protein subunits present in each channel. Key model assumptions are validated by showing our model is the reduction of a Markov model and through observations about voltage clamp data. We then show our model's success in replicating kinetic properties of concatemeric channels with different numbers of Kv1.1 and Kv1.2 subunits. Finally, through comparisons with multiple expression experiments across multiple voltage gated potassium families, we use the model to make predictions about the importance and effect of genetic mutations in heteromeric channel formation.

Automatic modeling of dynamic drug-hERG channel interactions using three voltage protocols and machine learning techniques: A simulation study

BACKGROUND

Assessment of drug cardiac safety is critical in the development of new compounds and is commonly addressed by evaluating the half-maximal blocking concentration of the potassium human ether-à-go-go related gene (hERG) channels. However, recent works have evidenced that the modelling of drug-binding dynamics to hERG can help to improve early cardiac safety assessment. Our goal is to develop a methodology to automatically generate Markovian models of the drug-hERG channel interactions.

METHODS

The training and the test sets consisted of 20800 and 5200 virtual drugs, respectively, distributed into 104 groups with different affinities and kinetics to the conformational states of the channel. In our system, drugs may bind to any state (individually or simultaneously), with different degrees of preference for a conformational state and the change of the conformational state of the drug bound channels may be restricted or allowed. To model such a wide range of possibilities, 12 Markovian chains are considered. Our approach uses the response of the drugs to our three previously developed voltage clamp protocols, which enhance the differences in the probabilities of occupying a certain conformational state of the channel (open, closed and inactivated). The computing tool is comprised of a classifier and a parameter optimizer and uses linear interpolation, support vector machines and a simplex method for function minimization.

RESULTS

We propose a novel methodology that automatically generates dynamic drug models using Markov model formulations and that elucidates the states where the drug binds and unbinds and the preferential binding state using data obtained from simple voltage clamp protocols that captures the preferential state-dependent binding properties, the relative affinities, trapping and non-trapping dynamics and the onset of I_Kr block. Overall, the tool correctly predicted the class of 92.04% of the drugs and the model provided by the tool accurately fitted the response of the target compound, the mean accuracy being 97.53%. Moreover, generation of the dynamic model of an I_Kr blocker from its response to our voltage clamp protocols usually takes less than an hour on a common desktop computer.

CONCLUSION

Our methodology could be very useful to model and simulate dynamic drug-hERG channel interactions. It would contribute to the improvement of the preclinical assessment of the proarrhythmic risk of drugs that inhibit I_Kr and the efficacy of antiarrhythmic I_Kr blockers.

Intrinsic mechanisms in the gating of resurgent Na+ currents

The resurgent component of the voltage-gated sodium current (I_NaR) is a depolarizing conductance, revealed on membrane hyperpolarizations following brief depolarizing voltage steps, which has been shown to contribute to regulating the firing properties of numerous neuronal cell types throughout the central and peripheral nervous systems. Although mediated by the same voltage-gated sodium (Nav) channels that underlie the transient and persistent Nav current components, the gating mechanisms that contribute to the generation of I_NaR remain unclear. Here, we characterized Nav currents in mouse cerebellar Purkinje neurons, and used tailored voltage-clamp protocols to define how the voltage and the duration of the initial membrane depolarization affect the amplitudes and kinetics of I_NaR. Using the acquired voltage-clamp data, we developed a novel Markov kinetic state model with parallel (fast and slow) inactivation pathways and, we show that this model reproduces the properties of the resurgent, as well as the transient and persistent, Nav currents recorded in (mouse) cerebellar Purkinje neurons. Based on the acquired experimental data and the simulations, we propose that resurgent Na⁺ influx occurs as a result of fast inactivating Nav channels transitioning into an open/conducting state on membrane hyperpolarization, and that the decay of I_NaR reflects the slow accumulation of recovered/opened Nav channels into a second, alternative and more slowly populated, inactivated state. Additional simulations reveal that extrinsic factors that affect the kinetics of fast or slow Nav channel inactivation and/or impact the relative distribution of Nav channels in the fast- and slow-inactivated states, such as the accessory Navβ4 channel subunit, can modulate the amplitude of I_NaR.

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.

Rate-dependent effects of lidocaine on cardiac dynamics: Development and analysis of a low-dimensional drug-channel interaction model

State-dependent sodium channel blockers are often prescribed to treat cardiac arrhythmias, but many sodium channel blockers are known to have pro-arrhythmic side effects. While the anti and proarrhythmic potential of a sodium channel blocker is thought to depend on the characteristics of its rate-dependent block, the mechanisms linking these two attributes are unclear. Furthermore, how specific properties of rate-dependent block arise from the binding kinetics of a particular drug is poorly understood. Here, we examine the rate-dependent effects of the sodium channel blocker lidocaine by constructing and analyzing a novel drug-channel interaction model. First, we identify the predominant mode of lidocaine binding in a 24 variable Markov model for lidocaine-sodium channel interaction by Moreno et al. Specifically, we find that (1) the vast majority of lidocaine bound to sodium channels is in the neutral form, i.e., the binding of charged lidocaine to sodium channels is negligible, and (2) neutral lidocaine binds almost exclusively to inactivated channels and, upon binding, immobilizes channels in the inactivated state. We then develop a novel 3-variable lidocaine-sodium channel interaction model that incorporates only the predominant mode of drug binding. Our low-dimensional model replicates an extensive amount of the voltage-clamp data used to parameterize the Moreno et al. model. Furthermore, the effects of lidocaine on action potential upstroke velocity and conduction velocity in our model are similar to those predicted by the Moreno et al. model. By exploiting the low-dimensionality of our model, we derive an algebraic expression for level of rate-dependent block as a function of pacing frequency, restitution properties, diastolic and plateau potentials, and drug binding rate constants. Our model predicts that the level of rate-dependent block is sensitive to alterations in restitution properties and increases in diastolic potential, but it is insensitive to variations in the shape of the action potential waveform and lidocaine binding rates.

Cardiac arrhythmias are often treated with drugs that block and alter the kinetics of membrane sodium channels. However, different drugs interact with sodium channels in different ways, and the complexity of the drug-channel interactions makes it difficult to predict whether a particular sodium channel blocker will reduce or increase the probability of cardiac arrhythmias. Here, we characterize the binding kinetics and effects on electrical signal propagation of the antiarrhythmic drug lidocaine, which is an archetypical example of a safe sodium channel blocker. Through analysis of a high-dimensional biophysically-detailed model of lidocaine-sodium channel interaction, we identify the predominant lidocaine binding pathway. We then incorporate only the key features of the predominant binding pathway into a novel low-dimensional model of lidocaine-sodium channel interaction. Our analysis of the low-dimensional model characterizes how the key binding properties of lidocaine affect electrical signal generation and propagation in the heart, and therefore our results are a step towards understanding the features that differentiate pro- and antiarrhythmic sodium channel blockers.

Ranolazine as an Alternative Therapy to Flecainide for SCN5A V411M Long QT Syndrome Type 3 Patients

The prolongation of the QT interval represents the main feature of the long QT syndrome (LQTS), a life-threatening genetic disease. The heterozygous SCN5A V411M mutation of the human sodium channel leads to a LQTS type 3 with severe proarrhythmic effects due to an increase in the late component of the sodium current (INaL). The two sodium blockers flecainide and ranolazine are equally recommended by the current 2015 ESC guidelines to treat patients with LQTS type 3 and persistently prolonged QT intervals. However, awareness of pro-arrhythmic effects of flecainide in LQTS type 3 patients arose upon the study of the SCN5A E1784K mutation. Regarding SCN5A V411M individuals, flecainide showed good results albeit in a reduced number of patients and no evidence supporting the use of ranolazine has ever been released. Therefore, we ought to compare the effect of ranolazine and flecainide in a SCN5A V411M model using an in-silico modeling and simulation approach. We collected clinical data of four patients. Then, we fitted four Markovian models of the human sodium current (INa) to experimental and clinical data. Two of them correspond to the wild type and the heterozygous SCN5A V411M scenarios, and the other two mimic the effects of flecainide and ranolazine on INa. Next, we inserted them into three isolated cell action potential (AP) models for endocardial, midmyocardial and epicardial cells and in a one-dimensional tissue model. The SCN5A V411M mutation produced a 15.9% APD90 prolongation in the isolated endocardial cell model, which corresponded to a 14.3% of the QT interval prolongation in a one-dimensional strand model, in keeping with clinical observations. Although with different underlying mechanisms, flecainide and ranolazine partially countered this prolongation at the isolated endocardial model by reducing the APD90 by 8.7 and 4.3%, and the QT interval by 7.2 and 3.2%, respectively. While flecainide specifically targeted the mutation-induced increase in peak INaL, ranolazine reduced it during the entire AP. Our simulations also suggest that ranolazine could prevent early afterdepolarizations triggered by the SCN5A V411M mutation during bradycardia, as flecainide. We conclude that ranolazine could be used to treat SCN5A V411M patients, specifically when flecainide is contraindicated.

Electrophysiological and anatomical factors determine arrhythmic risk in acute myocardial ischaemia and its modulation by sodium current availability

Acute myocardial ischaemia caused by coronary artery disease is one of the main causes of sudden cardiac death. Even though sodium current blockers are used as anti-arrhythmic drugs, decreased sodium current availability, also caused by mutations, has been shown to increase arrhythmic risk in ischaemic patients. The mechanisms are still unclear. Our goal is to exploit perfect control and data transparency of over 300 high-performance computing simulations to investigate arrhythmia mechanisms in acute myocardial ischaemia with variable sodium current availability. The human anatomically based torso-biventricular electrophysiological model used includes representation of realistic ventricular anatomy and fibre architecture, as well as ionic to electrocardiographic properties. Simulations show that reduced sodium current availability increased arrhythmic risk in acute regional ischaemia due to both electrophysiological (increased dispersion of refractoriness across the ischaemic border zone) and anatomical factors (conduction block from the thin right ventricle to thick left ventricle). The asymmetric ventricular anatomy caused high arrhythmic risk specifically for ectopic stimuli originating from the right ventricle and ventricular base. Increased sodium current availability was ineffective in reducing arrhythmic risk for septo-basal ectopic excitation. Human-based multiscale modelling and simulations reveal key electrophysiological and anatomical factors determining arrhythmic risk in acute ischaemia with variable sodium current availability.

A Review on Applications of Computational Methods in Drug Screening and Design

Drug development is one of the most significant processes in the pharmaceutical industry. Various computational methods have dramatically reduced the time and cost of drug discovery. In this review, we firstly discussed roles of multiscale biomolecular simulations in identifying drug binding sites on the target macromolecule and elucidating drug action mechanisms. Then, virtual screening methods (e.g., molecular docking, pharmacophore modeling, and QSAR) as well as structure- and ligand-based classical/de novo drug design were introduced and discussed. Last, we explored the development of machine learning methods and their applications in aforementioned computational methods to speed up the drug discovery process. Also, several application examples of combining various methods was discussed. A combination of different methods to jointly solve the tough problem at different scales and dimensions will be an inevitable trend in drug screening and design.

Delay differential equation-based models of cardiac tissue: Efficient implementation and effects on spiral-wave dynamics

Delay differential equations (DDEs) recently have been used in models of cardiac electrophysiology, particularly in studies focusing on electrical alternans, instabilities, and chaos. A number of processes within cardiac cells involve delays, and DDEs can potentially represent mechanisms that result in complex dynamics both at the cellular level and at the tissue level, including cardiac arrhythmias. However, DDE-based formulations introduce new computational challenges due to the need for storing and retrieving past values of variables at each spatial location. Cardiac tissue simulations that use DDEs may require over 28 GB of memory if the history of variables is not managed carefully. This paper addresses both computational and dynamical issues. First, we present new methods for the numerical solution of DDEs in tissue to mitigate the memory requirements associated with the history of variables. The new methods exploit the different time scales of an action potential to dynamically optimize history size. We find that the proposed methods decrease memory usage by up to 95% in cardiac tissue simulations compared to straightforward history-management algorithms. Second, we use the optimized methods to analyze for the first time the dynamics of wave propagation in two-dimensional cardiac tissue for models that include DDEs. In particular, we study the effects of DDEs on spiral-wave dynamics, including wave breakup and chaos, using a canine myocyte model. We find that by introducing delays to the gating variables governing the calcium current, DDEs can induce spiral-wave breakup in 2D cardiac tissue domains.

A Molecularly Detailed NaV1.5 Model Reveals a New Class I Antiarrhythmic Target

•

Antiarrhythmic therapies remain suboptimal due to our inability to predict how drug interactions with ion channels will affect the ability of the tissue to initiate and sustain an arrhythmia.

•

We built a computational framework that allows for in silico design of precision-targeted therapeutic agents that simultaneously assesses antiarrhythmic markers of success and failure at multiple spatial and time scales. Using this framework, a novel in silico mexiletine “booster” was designed that may dramatically improve the efficacy of mexiletine in suppression of arrhythmia triggers.

•

These results provide a roadmap for the design of novel molecular-based therapy to treat myriad arrhythmia syndromes, including ventricular tachycardia, heart failure arrhythmias, and inherited arrhythmia syndromes.

•

In summary, computational modeling approaches to drug discovery represent a novel tool to design and test precision-targeted therapeutic agents. By exploiting nontraditional ion channel drug targets, an entirely new dimension can be added to the wide parameter space of traditional antiarrhythmic drugs to develop more precision-targeted and potent Class I therapeutic agents.

Antiarrhythmic treatment strategies remain suboptimal due to our inability to predict how drug interactions with ion channels will affect the ability of the tissues to initiate and sustain an arrhythmia. We built a multiscale molecular model of the Na⁺ channel domain III (domain III voltage-sensing domain) to highlight the molecular underpinnings responsible for mexiletine drug efficacy. This model predicts that a hyperpolarizing shift in the domain III voltage-sensing domain is critical for drug efficacy and may be leveraged to design more potent Class I molecules. The model was therefore used to design, in silico, a theoretical mexiletine booster that can dramatically rescue a mutant resistant to the potent antiarrhythmic effects of mexiletine. Our framework provides a strategy for in silico design of precision-targeted therapeutic agents that simultaneously assesses antiarrhythmic markers of success and failure at multiple spatial and time scales. This approach provides a roadmap for the design of novel molecular-based therapy to treat myriad arrhythmia syndromes, including ventricular tachycardia, heart failure arrhythmias, and inherited arrhythmia syndromes.

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.

A computational model of induced pluripotent stem-cell derived cardiomyocytes incorporating experimental variability from multiple data sources

Key points

Induced pluripotent stem cell‐derived cardiomyocytes (iPSC‐CMs) capture patient‐specific genotype–phenotype relationships, as well as cell‐to‐cell variability of cardiac electrical activity

Computational modelling and simulation provide a high throughput approach to reconcile multiple datasets describing physiological variability, and also identify vulnerable parameter regimes

We have developed a whole‐cell model of iPSC‐CMs, composed of single exponential voltage‐dependent gating variable rate constants, parameterized to fit experimental iPSC‐CM outputs

We have utilized experimental data across multiple laboratories to model experimental variability and investigate subcellular phenotypic mechanisms in iPSC‐CMs

This framework links molecular mechanisms to cellular‐level outputs by revealing unique subsets of model parameters linked to known iPSC‐CM phenotypes

Abstract

There is a profound need to develop a strategy for predicting patient‐to‐patient vulnerability in the emergence of cardiac arrhythmia. A promising in vitro method to address patient‐specific proclivity to cardiac disease utilizes induced pluripotent stem cell‐derived cardiomyocytes (iPSC‐CMs). A major strength of this approach is that iPSC‐CMs contain donor genetic information and therefore capture patient‐specific genotype–phenotype relationships. A cited detriment of iPSC‐CMs is the cell‐to‐cell variability observed in electrical activity. We postulated, however, that cell‐to‐cell variability may constitute a strength when appropriately utilized in a computational framework to build cell populations that can be employed to identify phenotypic mechanisms and pinpoint key sensitive parameters. Thus, we have exploited variation in experimental data across multiple laboratories to develop a computational framework for investigating subcellular phenotypic mechanisms. We have developed a whole‐cell model of iPSC‐CMs composed of simple model components comprising ion channel models with single exponential voltage‐dependent gating variable rate constants, parameterized to fit experimental iPSC‐CM data for all major ionic currents. By optimizing ionic current model parameters to multiple experimental datasets, we incorporate experimentally‐observed variability in the ionic currents. The resulting population of cellular models predicts robust inter‐subject variability in iPSC‐CMs. This approach links molecular mechanisms to known cellular‐level iPSC‐CM phenotypes, as shown by comparing immature and mature subpopulations of models to analyse the contributing factors underlying each phenotype. In the future, the presented models can be readily expanded to include genetic mutations and pharmacological interventions for studying the mechanisms of rare events, such as arrhythmia triggers.

Induced pluripotent stem cell‐derived cardiomyocytes (iPSC‐CMs) capture patient‐specific genotype–phenotype relationships, as well as cell‐to‐cell variability of cardiac electrical activity

Computational modelling and simulation provide a high throughput approach to reconcile multiple datasets describing physiological variability, and also identify vulnerable parameter regimes

We have developed a whole‐cell model of iPSC‐CMs, composed of single exponential voltage‐dependent gating variable rate constants, parameterized to fit experimental iPSC‐CM outputs

We have utilized experimental data across multiple laboratories to model experimental variability and investigate subcellular phenotypic mechanisms in iPSC‐CMs

This framework links molecular mechanisms to cellular‐level outputs by revealing unique subsets of model parameters linked to known iPSC‐CM phenotypes

In silico Assessment of Pharmacotherapy for Human Atrial Patho-Electrophysiology Associated With hERG-Linked Short QT Syndrome

Short QT syndrome variant 1 (SQT1) arises due to gain-of-function mutations to the human Ether-à-go-go-Related Gene (hERG), which encodes the α subunit of channels carrying rapid delayed rectifier potassium current, I_Kr. In addition to QT interval shortening and ventricular arrhythmias, SQT1 is associated with increased risk of atrial fibrillation (AF), which is often the only clinical presentation. However, the underlying basis of AF and its pharmacological treatment remain incompletely understood in the context of SQT1. In this study, computational modeling was used to investigate mechanisms of human atrial arrhythmogenesis consequent to a SQT1 mutation, as well as pharmacotherapeutic effects of selected class I drugs–disopyramide, quinidine, and propafenone. A Markov chain formulation describing wild type (WT) and N588K-hERG mutant I_Kr was incorporated into a contemporary human atrial action potential (AP) model, which was integrated into one-dimensional (1D) tissue strands, idealized 2D sheets, and a 3D heterogeneous, anatomical human atria model. Multi-channel pharmacological effects of disopyramide, quinidine, and propafenone, including binding kinetics for I_Kr/hERG and sodium current, I_Na, were considered. Heterozygous and homozygous formulations of the N588K-hERG mutation shortened the AP duration (APD) by 53 and 86 ms, respectively, which abbreviated the effective refractory period (ERP) and excitation wavelength in tissue, increasing the lifespan and dominant frequency (DF) of scroll waves in the 3D anatomical human atria. At the concentrations tested in this study, quinidine most effectively prolonged the APD and ERP in the setting of SQT1, followed by disopyramide and propafenone. In 2D simulations, disopyramide and quinidine promoted re-entry termination by increasing the re-entry wavelength, whereas propafenone induced secondary waves which destabilized the re-entrant circuit. In 3D simulations, the DF of re-entry was reduced in a dose-dependent manner for disopyramide and quinidine, and propafenone to a lesser extent. All of the anti-arrhythmic agents promoted pharmacological conversion, most frequently terminating re-entry in the order quinidine > propafenone = disopyramide. Our findings provide further insight into mechanisms of SQT1-related AF and a rational basis for the pursuit of combined I_Kr and I_Na block based pharmacological strategies in the treatment of SQT1-linked AF.

Human Atrial Arrhythmogenesis and Sinus Bradycardia in KCNQ1-Linked Short QT Syndrome: Insights From Computational Modelling

Atrial fibrillation (AF) and sinus bradycardia have been reported in patients with short QT syndrome variant 2 (SQT2), which is underlain by gain-of-function mutations in KCNQ1 encoding the α subunit of channels carrying slow delayed rectifier potassium current, I_Ks. However, the mechanism(s) underlying the increased atrial arrhythmogenesis and impaired cardiac pacemaking activity arising from increased I_Ks remain unclear. Possible pharmacological interventions of AF in the SQT2 condition also remain to be elucidated. Using computational modelling, we assessed the functional impact of SQT2 mutations on human sinoatrial node (SAN) pacemaking, atrial repolarisation and arrhythmogenesis, and efficacy of the anti-arrhythmic drug quinidine. Markov chain formulations of I_Ks describing two KCNQ1 mutations – V141M and V307L – were developed from voltage-clamp experimental data and then incorporated into contemporary action potential (AP) models of human atrial and SAN cells, the former of which were integrated into idealised and anatomically detailed tissue models. Both mutations shortened atrial AP duration (APD) through distinct I_Ks ‘gain-of-function’ mechanisms, whereas SAN pacemaking rate was slowed markedly only by the V141M mutation. Differences in APD restitution steepness influenced re-entry dynamics in tissue – the V141M mutation promoted stationary and stable spiral waves whereas the V307L mutation promoted non-stationary and unstable re-entrant waves. Both mutations shortened tissue excitation wavelength through reduced effective refractory period but not conduction velocity, which served to increase the lifespan of re-entrant excitation in a 3D anatomical human atria model, as well as the dominant frequency (DF), which was higher for the V141M mutation. Quinidine was effective at terminating arrhythmic excitation waves associated with the V307L but not V141M mutation, and reduced the DF in a dose-dependent manner under both mutation conditions. This study provides mechanistic insights into different AF/bradycardia phenotypes in SQT2 and the efficacy of quinidine pharmacotherapy.

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.

Computational motility models of neurogastroenterology and neuromodulation

The success of neuromodulation therapies, particularly in the brain, spinal cord, and peripheral nerves, has been greatly aided by computational, biophysical models. However, treating gastrointestinal disorders with electrical stimulation has been much less explored, partly because the mode of action of such treatments is unclear, and selection of stimulation parameters is often empirical. Progress in gut neuromodulation is limited by the comparative lack of biophysical models capable of simulating neuromodulation of gastrointestinal function. Here, we review the recently developed biophysical models of electrically-active cells in the gastrointestinal system that contribute to motility. Biophysical models are replacing phenomenologically-defined models due to advancements in electrophysiological characterization of key players in the gut: enteric neurons, smooth muscle fibers, and interstitial cells of Cajal. In this review, we explore existing biophysically-defined cellular and network models that contribute to gastrointestinal motility. We focus on recent models that are laying the groundwork for modeling electrical stimulation of the gastrointestinal system. Developing models of gut neuromodulation will improve our mechanistic understanding of these treatments, leading to better parameterization, selectivity, and efficacy of neuromodulation to treat gastrointestinal disorders. Such models may have direct clinical translation to current neuromodulation therapies, such as sacral nerve stimulation.

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.

How to Connect Cardiac Excitation to the Atomic Interactions of Ion Channels

Many have worked to create cardiac action potential models that explicitly represent atomic-level details of ion channel structure. Such models have the potential to define new therapeutic directions and to show how nanoscale perturbations to channel function predispose patients to deadly cardiac arrhythmia. However, there have been significant experimental and theoretical barriers that have limited model usefulness. Recently, many of these barriers have come down, suggesting that considerable progress toward creating these long-sought models may be possible in the near term.

Computational Analysis of the Mode of Action of Disopyramide and Quinidine on hERG-Linked Short QT Syndrome in Human Ventricles

The short QT syndrome (SQTS) is a rare cardiac disorder associated with arrhythmias and sudden death. Gain-of-function mutations to potassium channels mediating the rapid delayed rectifier current, I_Kr, underlie SQTS variant 1 (SQT1), in which treatment with Na⁺ and K⁺ channel blocking class Ia anti-arrhythmic agents has demonstrated some efficacy. This study used computational modeling to gain mechanistic insights into the actions of two such drugs, disopyramide and quinidine, in the setting of SQT1. The O'Hara-Rudy (ORd) human ventricle model was modified to incorporate a Markov chain formulation of I_Kr describing wild type (WT) and SQT1 mutant conditions. Effects of multi-channel block by disopyramide and quinidine, including binding kinetics and altered potency of I_Kr/hERG channel block in SQT1 and state-dependent block of sodium channels, were simulated on action potential and multicellular tissue models. A one-dimensional (1D) transmural ventricular strand model was used to assess prolongation of the QT interval, effective refractory period (ERP), and re-entry wavelength (WL) by both drugs. Dynamics of re-entrant excitation waves were investigated using a 3D human left ventricular wedge model. In the setting of SQT1, disopyramide, and quinidine both produced a dose-dependent prolongation in (i) the QT interval, which was primarily due to I_Kr block, and (ii) the ERP, which was mediated by a synergistic combination of I_Kr and I_Na block. Over the same range of concentrations quinidine was more effective in restoring the QT interval, due to more potent block of I_Kr. Both drugs demonstrated an anti-arrhythmic increase in the WL of re-entrant circuits. In the 3D wedge, disopyramide and quinidine at clinically-relevant concentrations decreased the dominant frequency of re-entrant excitations and exhibited anti-fibrillatory effects; preventing formation of multiple, chaotic wavelets which developed in SQT1, and could terminate arrhythmias. This computational modeling study provides novel insights into the clinical efficacy of disopyramide and quinidine in the setting of SQT1; it also dissects ionic mechanisms underlying QT and ERP prolongation. Our findings show that both drugs demonstrate efficacy in reversing the SQT1 phenotype, and indicate that disopyramide warrants further investigation as an alternative to quinidine in the treatment of SQT1.