Studying dyadic structure-function relationships: a review of current modeling approaches and new insights into Ca2+ (mis)handling.

Excitation-contraction coupling in cardiac myocytes requires calcium influx through L-type calcium channels in the sarcolemma, which gates calcium release through sarcoplasmic reticulum ryanodine receptors in a process known as calcium-induced calcium release, producing a myoplasmic calcium transient and enabling cardiomyocyte contraction. The spatio-temporal dynamics of calcium release, buffering, and reuptake into the sarcoplasmic reticulum play a central role in excitation-contraction coupling in both normal and diseased cardiac myocytes. However, further quantitative understanding of these cells' calcium machinery and the study of mechanisms that underlie both normal cardiac function and calcium-dependent etiologies in heart disease requires accurate knowledge of cardiac ultrastructure, protein distribution and subcellular function. As current imaging techniques are limited in spatial resolution, limiting insight into changes in calcium handling, computational models of excitation-contraction coupling have been increasingly employed to probe these structure-function relationships. This review will focus on the development of structural models of cardiac calcium dynamics at the subcellular level, orienting the reader broadly towards the development of models of subcellular calcium handling in cardiomyocytes. Specific focus will be given to progress in recent years in terms of multi-scale modeling employing resolved spatial models of subcellular calcium machinery. A review of the state-of-the-art will be followed by a review of emergent insights into calcium-dependent etiologies in heart disease and, finally, we will offer a perspective on future directions for related computational modeling and simulation efforts.

An Evaluation of the Accuracy of Classical Models for Computing the Membrane Potential and Extracellular Potential for Neurons.

Two mathematical models are part of the foundation of Computational neurophysiology; (a) the Cable equation is used to compute the membrane potential of neurons, and, (b) volume-conductor theory describes the extracellular potential around neurons. In the standard procedure for computing extracellular potentials, the transmembrane currents are computed by means of (a) and the extracellular potentials are computed using an explicit sum over analytical point-current source solutions as prescribed by volume conductor theory. Both models are extremely useful as they allow huge simplifications of the computational efforts involved in computing extracellular potentials. However, there are more accurate, though computationally very expensive, models available where the potentials inside and outside the neurons are computed simultaneously in a self-consistent scheme. In the present work we explore the accuracy of the classical models (a) and (b) by comparing them to these more accurate schemes. The main assumption of (a) is that the ephaptic current can be ignored in the derivation of the Cable equation. We find, however, for our examples with stylized neurons, that the ephaptic current is comparable in magnitude to other currents involved in the computations, suggesting that it may be significant-at least in parts of the simulation. The magnitude of the error introduced in the membrane potential is several millivolts, and this error also translates into errors in the predicted extracellular potentials. While the error becomes negligible if we assume the extracellular conductivity to be very large, this assumption is, unfortunately, not easy to justify a priori for all situations of interest.

Pleiotropic effects of schizophrenia-associated genetic variants in neuron firing and cardiac pacemaking revealed by computational modeling.

Schizophrenia patients have an increased risk of cardiac dysfunction. A possible factor underlying this comorbidity are the common variants in the large set of genes that have recently been discovered in genome-wide association studies (GWASs) as risk genes of schizophrenia. Many of these genes control the cell electrogenesis and calcium homeostasis. We applied biophysically detailed models of layer V pyramidal cells and sinoatrial node cells to study the contribution of schizophrenia-associated genes on cellular excitability. By including data from functional genomics literature to simulate the effects of common variants of these genes, we showed that variants of voltage-gated Na+ channel or hyperpolarization-activated cation channel-encoding genes cause qualitatively similar effects on layer V pyramidal cell and sinoatrial node cell excitability. By contrast, variants of Ca2+ channel or transporter-encoding genes mostly have opposite effects on cellular excitability in the two cell types. We also show that the variants may crucially affect the propagation of the cardiac action potential in the sinus node. These results may help explain some of the cardiac comorbidity in schizophrenia, and may facilitate generation of effective antipsychotic medications without cardiac side-effects such as arrhythmia.

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

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

NS5806 partially restores action potential duration but fails to ameliorate calcium transient dysfunction in a computational model of canine heart failure.

AIMS: The study investigates how increased Ito, as mediated by the activator NS5806, affects excitation-contraction coupling in chronic heart failure (HF). We hypothesized that restoring spike-and-dome morphology of the action potential (AP) to a healthy phenotype would be insufficient to restore the intracellular Ca(2) (+) transient (CaT), due to HF-induced remodelling of Ca(2+) handling. METHODS AND RESULTS: An existing mathematical model of the canine ventricular myocyte was modified to incorporate recent experimental data from healthy and failing myocytes, resulting in models of both healthy and HF epicardial, midmyocardial, and endocardial cell variants. Affects of NS5806 were also included in HF models through its direct interaction with Kv4.3 and Kv1.4. Single-cell simulations performed in all models (control, HF, and HF + drug) and variants (epi, mid, and endo) assessed AP morphology and underlying ionic processes with a focus on calcium transients (CaT), how these were altered in HF across the ventricular wall, and the subsequent effects of varying compound concentration in HF. Heart failure model variants recapitulated a characteristic increase in AP duration (APD) in the disease. The qualitative effects of application of half-maximal effective concentration (EC50) of NS5806 on APs and CaT are heterogeneous and non-linear. Deepening in the AP notch with drug is a direct effect of the activation of Ito; both Ito and consequent alteration of IK1 kinetics cause decrease in AP plateau potential. Decreased APD50 and APD90 are both due to altered IK1. Analysis revealed that drug effects depend on transmurality. Ca(2+) transient morphology changes-increased amplitude and shorter time to peak-are due to direct increase in ICa,L and indirect larger SR Ca(2+) release subsequent to Ito activation. CONCLUSIONS: Downstream effects of a compound acting exclusively on sarcolemmal ion channels are difficult to predict. Remediation of APD to pre-failing state does not ameliorate dysfunction in CaT; however, restoration of notch depth appears to impart modest benefit and a likelihood of therapeutic value in modulating early repolarization.

Computing the stability of steady-state solutions of mathematical models of the electrical activity in the heart.

Instabilities in the electro-chemical resting state of the heart can generate ectopic waves that in turn can initiate arrhythmias. We derive methods for computing the resting state for mathematical models of the electro-chemical process underpinning a heartbeat, and we estimate the stability of the resting state by invoking the largest real part of the eigenvalues of a linearized model. The implementation of the methods is described and a number of numerical experiments illustrate the feasibility of the methods. In particular, we test the methods for problems where we can compare the solutions with analytical results, and problems where we have solutions computed by independent software. The software is also tested for a fairly realistic 3D model.

Defining candidate drug characteristics for Long-QT (LQT3) syndrome.

Mutations of the SCN5A gene can significantly alter the function of cardiac myocyte sodium channels leading to increased risk of ventricular arrhythmia. Over the past decade, detailed Markov models of the action potential of cardiac cells have been developed. In such models, the effects of a drug can be treated as alterations in on- and off rates between open and inactivated states on one hand, and blocked states on the other hand. Our aim is to compute the rates specifying a drug in order to: (a) restore the steady-state open probability of the mutant channel to that of normal wild type channels; and (b) minimize the difference between whole cell currents in drugged mutant and wild type cells. The difference in the electrochemical state vector of the cell can be measured in a norm taking all components and their dynamical properties into account. Measured with this norm, the difference between the state of the mutant and wild-type cell was reduced by a factor of 36 after the drug was introduced and by factors of 4 over mexitiline and 25 over lidocaine. The results suggest the potential to synthesize more effective drugs based on mechanisms of action of existing compounds.

Control of Ca2+ release by action potential configuration in normal and failing murine cardiomyocytes.

Cardiomyocytes from failing hearts exhibit spatially nonuniform or dyssynchronous sarcoplasmic reticulum (SR) Ca(2+) release. We investigated the contribution of action potential (AP) prolongation in mice with congestive heart failure (CHF) after myocardial infarction. AP recordings from CHF and control myocytes were included in a computational model of the dyad, which predicted more dyssynchronous ryanodine receptor opening during stimulation with the CHF AP. This prediction was confirmed in cardiomyocyte experiments, when cells were alternately stimulated by control and CHF AP voltage-clamp waveforms. However, when a train of like APs was used as the voltage stimulus, the control and CHF AP produced a similar Ca(2+) release pattern. In this steady-state condition, greater integrated Ca(2+) entry during the CHF AP lead to increased SR Ca(2+) content. A resulting increase in ryanodine receptor sensitivity synchronized SR Ca(2+) release in the mathematical model, thus offsetting the desynchronizing effects of reduced driving force for Ca(2+) entry. A modest nondyssynchronous prolongation of Ca(2+) release was nevertheless observed during the steady-state CHF AP, which contributed to increased time-to-peak measurements for Ca(2+) transients in failing cells. Thus, dyssynchronous Ca(2+) release in failing mouse myocytes does not result from electrical remodeling, but rather other alterations such as T-tubule reorganization.

A method for analyzing the stability of the resting state for a model of pacemaker cells surrounded by stable cells.

The purpose of this paper is to derive and analyze methods for examining the stability of solutions of partial differential equations modeling collections of excitable cells. In particular, we derive methods for estimating the principal eigenvalue of a linearized version of the Luo-Rudy I model close to an equilibrium solution. It has been suggested that the stability of a collection of unstable cells surrounded by a large collection of stable cells can be studied by considering only a collection of unstable cells equipped with a Dirichlet type boundary condition. This method has earlier been applied to analytically assess the stability of a reduced version the Luo-Rudy I model. In this paper we analyze the accuracy of this technique and apply it to the full Luo-Rudy I model. Furthermore, we extend the method to provide analytical results for the FitzHugh-Nagumo model in the case where a collection of unstable cells is surrounded by a collection of stable cells. All our analytical findings are complemented by numerical computations computing the principal eigenvalue of a discrete version of linearized models.

Synchronizing computer simulations with measurement data for a case of atrial flutter.

Atrial flutter is a common supraventricular tachycardia that can be treated using radiofrequency catheter ablation, a procedure that is guided by electroanatomical mapping systems. In this paper, we propose an algorithm for incorporating mapping data into computer simulations of atrial electrical activity with the purpose of creating a more accurate map of electrical activation. The algorithm takes as input the extracellular potential values recorded at a number of sites throughout the atria and estimates the activation time for the entire atrial domain. We test the algorithm using synthetic mapping data and an anatomically detailed atrial geometry with an activation pattern typical of atrial flutter. The results show that the algorithm performs well with synthetic mapping data with information from relatively few mapping sites and in the presence of modeling and measurement error.

A note on a method for determining advantageous properties of an anti-arrhythmic drug based on a mathematical model of cardiac cells.

Regional hyperkalemia during acute ischemia may provoke cardiac arrhythmias such as ventricular fibrillation. Despite intense research efforts over the last decades, the problem of finding an efficient anti-arrhythmic drug without dangerous side effects is still open. One approach to analyze the effect of anti-arrhythmic drugs is to do simulations based on mathematical models of collections of cardiomyocytes. Such simulations have recently illuminated the pro-arrhythmic capability of well-established anti-arrhythmic drugs. The purpose of the present note is to introduce a method intended for computing advantageous properties of an anti-arrhythmic drug. For a given model of a normal and an ischemic cell, we introduce a drug as a vector of non-negative real numbers whose components are multiplied by individual terms representing specific ionic currents. The drug vector is computed such that the action potentials of the resulting drugged cells are as close as possible to the action potential of a normal (not drugged) cell. Numerical simulations based on the Luo-Rudy I model and the Hund-Rudy model show that the classical shortened action potential obtained due to hyperkalemia is prolonged by using the drug computed by this method. Furthermore, for both models a 2D collection of spatially coupled ischemic cells give arrhythmogenic solutions before the drug is applied, and stable solutions after the drug is applied. It is emphasized that we do not address the possibility of realizing a drug with the properties computed in this note.

A condition for setting off ectopic waves in computational models of excitable cells.

The purpose of this paper is to study the stability of steady state solutions of the Monodomain model equipped with Luo-Rudy I kinetics. It is well established that re-entrant arrhythmias can be created in computational models of excitable cells. Such arrhythmias can be initiated by applying an external stimulus that interacts with a partially refractory region, and spawn breaking waves that can eventually generate extremely complex wave patterns commonly referred to as fibrillation. An ectopic wave is one possible stimulus that may initiate fibrillation. Physiologically, it is well known that ectopic waves exist, but the mechanism for initiating ectopic waves in a large collection of cells is poorly understood. In the present paper we consider computational models of collections of excitable cells in one and two spatial dimensions. The cells are modeled by Luo-Rudy I kinetics, and we assume that the spatial dynamics is governed by the Monodomain model. The mathematical analysis is carried out for a reduced model that is known to provide good approximations of the initial phase of solutions of the Luo-Rudy I model. A further simplification is also introduced to motivate and explain the results for the more complicated models. In the analysis the cells are divided into two regions; one region (N) consists of normal cells as model by the standard Luo-Rudy I model, and another region (A) where the cells are automatic in the sense that they would act as pacemaker cells if they where isolated from their surroundings. We let delta denote the spatial diffusion and a denote a characteristic length of the automatic region. It has previously been shown that reducing diffusion or increasing the automatic region enhances ectopic activity. Here we derive a condition for the transition from stable resting state to ectopic wave spread. Under suitable assumptions on the model we provide mathematical and computational arguments indicating that there is a constant eta such that a steady state solution of this system is stable whenever delta approximately > etaa(2), and unstable whenever delta approximately < etaa(2).

Stochastic binding of Ca2+ ions in the dyadic cleft; continuous versus random walk description of diffusion.

Ca(2+) signaling in the dyadic cleft in ventricular myocytes is fundamentally discrete and stochastic. We study the stochastic binding of single Ca(2+) ions to receptors in the cleft using two different models of diffusion: a stochastic and discrete Random Walk (RW) model, and a deterministic continuous model. We investigate whether the latter model, together with a stochastic receptor model, can reproduce binding events registered in fully stochastic RW simulations. By evaluating the continuous model goodness-of-fit for a large range of parameters, we present evidence that it can. Further, we show that the large fluctuations in binding rate observed at the level of single time-steps are integrated and smoothed at the larger timescale of binding events, which explains the continuous model goodness-of-fit. With these results we demonstrate that the stochasticity and discreteness of the Ca(2+) signaling in the dyadic cleft, determined by single binding events, can be described using a deterministic model of Ca(2+) diffusion together with a stochastic model of the binding events, for a specific range of physiological relevant parameters. Time-consuming RW simulations can thus be avoided. We also present a new analytical model of bimolecular binding probabilities, which we use in the RW simulations and the statistical analysis.

Contribution of the Na+/Ca2+ exchanger to rapid Ca2+ release in cardiomyocytes.

Trigger Ca(2+) is considered to be the Ca(2+) current through the L-type Ca(2+) channel (LTCC) that causes release of Ca(2+) from the sarcoplasmic reticulum. However, cell contraction also occurs in the absence of the LTCC current (I(Ca)). In this article, we investigate the contribution of the Na(+)/Ca(2+) exchanger (NCX) to the trigger Ca(2+). Experimental data from rat cardiomyocytes using confocal microscopy indicating that inhibition of reverse mode Na(+)/Ca(2+) exchange delays the Ca(2+) transient by 3-4 ms served as a basis for the mathematical model. A detailed computational model of the dyadic cleft (fuzzy space) is presented where the diffusion of both Na(+) and Ca(2+) is taken into account. Ionic channels are included at discrete locations, making it possible to study the effect of channel position and colocalization. The simulations indicate that if a Na(+) channel is present in the fuzzy space, the NCX is able to bring enough Ca(2+) into the cell to affect the timing of release. However, this critically depends on channel placement and local diffusion properties. With fuzzy space diffusion in the order of four orders of magnitude lower than in water, triggering through LTCC alone was up to 5 ms slower than with the presence of a Na(+) channel and NCX.