Search

The cable equation was first derived to model transport of electrical signals in telegraphic cables. But it later gained enormous popularity as a model of transport of electrical signals along a neuronal axon. In Chapter 9, we will discuss how this equation is derived and how the different terms in the equation come about. But here, we will just take a simple version of the equations for granted and then try to solve them. We will observe that the few techniques we learned above are actually sufficient to solve the non-linear reaction-diffusion equations we consider here.

In mathematics, a common approach to solving a new problem is to break down the problem into sub-problems that you know how to solve and then try to glue the pieces together to yield a solution of the new problem. This approach is also very useful in software development; well-tested pieces of software can be glued together in order to obtain solutions to a wider class of problems. Operator splitting is a technique that illustrates this principle very well. Rather complex equations can be broken down into more familiar problems and solved individually. It’s a miracle that it works, but it does. And it’s no miracle because there are proofs of convergence. Anyway, we will illustrate operator splitting with two examples and then come back to this technique when the equations become more challenging.

In the examples we have considered so far, the error introduced by the numerical scheme has been O(∆t). This can easily be improved as we will show through two simple examples.

In the previous chapter, we saw that the simple explicit numerical scheme resulted in an instability problem. We also saw that the problem could be resolved by using sufficiently short time steps. But in many situations, short time steps become exceedingly short, as can be seen, e.g., in the stability criterion (3.56). This means that we have to perform computations for a very large number of time steps to reach the final time and it is therefore tempting to look for alternatives. The most common alternative is to use an implicit scheme, which generally allows for much longer time steps.

In Chapter 8, we introduced mathematical models of the action potential across a cell membrane. Understanding the properties of the cell membrane is absolutely essential in order to understand the electrophysiology of excitable cells. However, some essential properties can only be studied in spatially resolved models; i.e., in models representing spatial variation across a single cell or a collection of cells.

In these notes, we have considered models of electrophysiology across several scales. The first was the membrane model. It assumes that the action potential is similar across the whole cell membrane, and the model represents the action potential as a function of time alone. No spatial variable is involved in the pure membrane models, so a length scale of these models does not make sense.

The diffusion equation appears in many applications in science and engineering, and computational physiology is no exception. In its most basic form, the diffusion equation is also useful as an example of how to deal with a PDE using numerical methods. We will start by considering it as a stand-alone model, but in the next chapterswe will study it in combination with non-linear ODEs. This chapter therefore serves as a warm-up for the more complex models. We will also follow the path we started above. In the very simplest case of an ODE, we found a formula for the solution of both the differential equation and the numerical scheme approximating the equation. One nice consequence of this is that, as we saw above, we can explicitly study the error introduced by the numerical approximation. In this chapter, we will use the same approach to study the error of the numerical approximation of the diffusion equation.

We will introduce numerical methods for systems of ODEs by considering the celebrated FitzHugh-Nagumo model published by FitzHugh [1] in 1961 and, independently, by Nagumo et. al. [2] in 1962. The model is a system of ordinary differential equations with two unknowns, and is commonly used as a simple model for the action potentials of excitable pacemaker cells.

In Chapter 6, we studied a simple version of the cable equation, where a diffusion term was added to the FitzHugh-Nagumo equations. In this chapter, we will revisit the cable equation and go through a simple derivation of the model. In addition, we will consider the numerical solution of the cable equation for a neuronal axon with membrane dynamics modeled by the Hodgkin-Huxley model.

As mentioned earlier, the bidomain system is currently the standard mathematical model of cardiac electrophysiology. This system is nowroutinely solved and provides valuable insights into the conduction of electrical signals in cardiac tissue. However, the model has one glaring limitation: The cardiomyocyte is nowhere to be found in the model, since the extracellular space, the intracellular space and the cell membrane are all assumed to be everywhere in the computational domain. The cell was lost in homogenization! There is a tremendous advantage to this because the model becomes much simpler and thus solvable for the whole human heart. And it works! But the downside is of course that the cell is the essential building block of the tissue and leaving it out of the model has consequences. For instance, it becomes impossible to investigate the detailed dynamics of the electrochemical processes in the vicinity of a small collection of cells.

In the previous chapters, we introduced techniques that can be used to find numerical solutions of differential equations. The examples we considered were simple, theoretical differential equations without units. From this point forwards, we will look at how to apply the methods introduced in the previous chapters to differential equations that are set up to model aspects of electrophysiology.

In this introductory chapter we will introduce two essential concepts: differential equations and numerical methods. If you want to spend as little energy as possible (don’t be ashamed of that - energy preservation is both fashionable and a fundamental property of many biological mechanisms - it’s fine) you can get a good overview from this chapter alone.

Evaluation of arrhythmogenic drugs is required by regulatory agencies before any new compound can obtain market approval. Despite rigorous review, cardiac disorders remain the second most common cause for safety-related market withdrawal. On the other hand, false-positive preclinical findings prohibit potentially beneficial candidates from moving forward in the development pipeline. Complex

Atrial fibrillation (AF) is a common health problem with substantial individual and societal costs. The origin of AF has been debated for more than a century, and the precise, biophysical mechanisms that are responsible for the initiation and maintenance of the chaotic electrochemical waves that define AF, remains unclear. It is well accepted that the outlet of the pulmonary veins is the primary anatomical site of AF initiation, and that electrical isolation of these regions remains the most effective treatment for AF. Furthermore, it is well known that certain ion channel or transporter mutations can significantly increase the likelihood of AF. Here, we present a computational model capable of characterizing functionally important features of the microanatomical and electrophysiological substrate that represents the transition from the pulmonary veins (PV) to the left atrium (LA) of the human heart. This model is based on a finite element representation of every myocyte in a segment of this (PV/LA) region. Thus, it allows for investigation a mix of typical PV and LA myocytes. We use the model to investigate the likelihood of ectopic beats and re-entrant waves in a cylindrical geometry representing the transition from PV to LA. In particular, we investigate and illustrate how six different AF- associated mutations can alter the probability of ectopic beats and re-entry in this region.

Mutations are known to cause perturbations in essential functional features of integral membrane proteins, including ion channels. Even restricted or point mutations can result in substantially changed properties of ion currents. The additive effect of these alterations for a specific ion channel can result in significantly changed properties of the action potential (AP). Both AP shortening and AP prolongation can result from known mutations, and the consequences can be life-threatening. Here, we present a computational method for identifying new drugs utilizing combinations of existing drugs. Based on the knowledge of theoretical effects of existing drugs on individual ion currents, our aim is to compute optimal combinations that can ‘repair’ the mutant AP waveforms so that the baseline AP-properties are restored. More specifically, we compute optimal, combined, drug concentrations such that the waveforms of the transmembrane potential and the cytosolic calcium concentration of the mutant cardiomyocytes (CMs) becomes as similar as possible to their wild type counterparts after the drug has been applied. In order to demonstrate the utility of this method, we address the question of computing an optimal drug for the short QT syndrome type 1 (SQT1). For the SQT1 mutation N588K, there are available data sets that describe the effect of various drugs on the mutated K+ channel. These published findings are the basis for our computational analysis which can identify optimal compounds in the sense that the AP of the mutant CMs resembles essential biomarkers of the wild type CMs. Using recently developed insights regarding electrophysiological properties among myocytes from different species, we compute optimal drug combinations for hiPSC-CMs, rabbit ventricular CMs and adult human ventricular CMs with the SQT1 mutation. Since the ‘composition’ of ion channels that form the AP is different for the three types of myocytes under consideration, so is the composition of the optimal drug., Author summary Poly-pharmacology (using multiple drugs to treat disease) has been proposed for improving cardiac anti-arrhythmic therapy for at least two decades. However, the specific arrhythmia contexts in which polytherapy is likely to be both safe and effective have remained elusive. Type 1 short QT syndrome (SQT1) is a rare form of cardiac arrhythmia that results from mutations to the human Ether-á-go-go Related Gene (hERG) potassium channel. Functionally, these mutations are remarkably consistent in that they permit the channel to open earlier during each heart beat. While hundreds of compounds are known to inhibit hERG channels, the specific effect of SQT1 mutations that allows for early channel opening also limits the ability of most of those compounds to correct SQT1 dysfunction. Here, we have applied a suite of ventricular cardiomyocyte computational models to ask whether polytherapy may offer a more effective therapeutic strategy in SQT1, and if so, what the likely characteristics of that strategy are. Our analyses suggest that simultaneous induction of late sodium current and partial hERG blockade offers a promising strategy. While no activators of late sodium current have been clinically approved, several experimental compounds are available and may provide a basis for interrogating this strategy. The method presented here can be used to compute optimal drug combinations provided that the effect of each drug on every relevant ion channel is known.

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., Author summary A number of cardiac disorders originate from genetic mutations affecting the function of ion channels populating the membrane of cardiomyocytes. One example is short QT syndrome, associated with increased risk of arrhythmias and sudden cardiac death. Cardiomyocytes derived from human induced pluripotent stem cells (hiPSC-CMs) provide a promising platform for testing potential pharmacological treatments for such disorders, as human cardiomyocytes exhibiting specific mutations can be generated and exposed to drugs in vitro. However, the electrophysiological properties of hiPSC-CMs differ significantly from those of adult native cardiomyocytes. Therefore, drug effects observed for hiPSC-CMs could possibly be very different from corresponding drug effects for adult cells in vivo. In this study, we apply a computational framework for translating drug effects observed for hiPSC-CMs derived from a short QT patient to drug effects for adult short QT cardiomyocytes. For one of the considered drugs, the effect on adult QT intervals has been measured and these measurements turn out to be in good agreement with the response estimated by the computational procedure. Thus, the computational framework shows promise for being a useful tool for predicting adult drug responses from measurements of hiPSC-CMs, allowing earlier identification of compounds to accurately treat cardiac diseases.

We want to be able to perform accurate simulations of a large number of cardiac cells based on mathematical models where each individual cell is represented in the model. This implies that the computational mesh has to have a typical resolution of a fewµm leading to huge computational challenges. In this paper we use a certain operator splitting of the coupled equations and showthat this leads to systems that can be solved in parallel. This opens up for the possibility of simulating large numbers of coupled cardiac cells.

Excitable cells are of vital importance in biology, and mathematical models have contributed significantly to understand their basic mechanisms. However, classical models of excitable cells are based on severe assumptions that may limit the accuracy of the simulation results. Here, we derive a more detailed approach to modeling that has recently been applied to study the electrical properties of both neurons and cardiomyocytes. The model is derived from first principles and opens up possibilities for studying detailed properties of excitable cells.We refer to the model as the EMI model because both the extracellular space (E), the cell membrane (M) and the intracellular space (I) are explicitly represented in the model, in contrast to classical spatial models of excitable cells. Later chapters of the present text will focus on numerical methods and software for solving the model. Also, in the next chapter, the model will be extended to account for ionic concentrations in the intracellular and extracellular spaces.

Using animal cells and tissues as precise measuring devices for developing new drugs presents a long-standing challenge for the pharmaceutical industry. Despite the very significant resources that continue to be dedicated to animal testing of new compounds, only qualitative results can be obtained. This often results in both false positives and false negatives. Here, we show how the effect of drugs applied to animal ventricular myocytes can be translated, quantitatively, to estimate a number of different effects of the same drug on human cardiomyocytes. We illustrate and validate our methodology by translating, from animal to human, the effect of dofetilide applied to dog cardiomyocytes, the effect of E-4031 applied to zebrafish cardiomyocytes, and, finally, the effect of sotalol applied to rabbit cardiomyocytes. In all cases, the accuracy of our quantitative estimates are demonstrated. Our computations reveal that, in principle, electrophysiological data from testing using animal ventricular myocytes, can give precise, quantitative estimates of the effect of new compounds on human cardiomyocytes.

The immature physiology of cardiomyocytes derived from human induced pluripotent stem cells (hiPSCs) limits their utility for drug screening and disease modelling. Here we show that suitable combinations of mechanical stimuli and metabolic cues can enhance the maturation of hiPSC-derived cardiomyocytes, and that the maturation-inducing cues have phenotype-dependent effects on the cells' action-potential morphology and calcium handling. By using microfluidic chips that enhanced the alignment and extracellular-matrix production of cardiac microtissues derived from genetically distinct sources of hiPSC-derived cardiomyocytes, we identified fatty-acid-enriched maturation media that improved the cells' mitochondrial structure and calcium handling, and observed divergent cell-source-dependent effects on action-potential duration (APD). Specifically, in the presence of maturation media, tissues with abnormally prolonged APDs exhibited shorter APDs, and tissues with aberrantly short APDs displayed prolonged APDs. Regardless of cell source, tissue maturation reduced variabilities in spontaneous beat rate and in APD, and led to converging cell phenotypes (with APDs within the 300-450 ms range characteristic of human left ventricular cardiomyocytes) that improved the modelling of the effects of pro-arrhythmic drugs on cardiac tissue.

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.

ObjectiveMechanistic modeling of neurons is an essential component of computational neuroscience that enables scientists to simulate, explain, and explore neural activity. The conventional approach to simulation of extracellular neural recordings first computes transmembrane currents using the cable equation and then sums their contribution to model the extracellular potential. This two-step approach relies on the assumption that the extracellular space is an infinite and homogeneous conductive medium, while measurements are performed using neural probes. The main purpose of this paper is to assess to what extent the presence of the neural probes of varying shape and size impacts the extracellular field and how to correct for them.ApproachWe apply a detailed modeling framework allowing explicit representation of the neuron and the probe to study the effect of the probes and thereby estimate the effect of ignoring it. We use meshes with simplified neurons and different types of probe and compare the extracellular action potentials with and without the probe in the extracellular space. We then compare various solutions to account for the probes’ presence and introduce an efficient probe correction method to include theprobe effectin modeling of extracellular potentials.Main resultsOur computations show that microwires hardly influence the extracellular electric field and their effect can therefore be ignored. In contrast, Multi-Electrode Arrays (MEAs) significantly affect the extracellular field by magnifying the recorded potential. While MEAs behave similarly to infinite insulated planes, we find that their effect strongly depends on the neuron-probe alignment and probe orientation.SignificanceIgnoring theprobe effectmight be deleterious in some applications, such as neural localization and parameterization of neural models from extracellular recordings. Moreover, the presence of the probe can improve the interpretation of extracellular recordings, by providing a more accurate estimation of the extracellular potential generated by neuronal models.

Human induced pluripotent stem cell derived cardiomyocytes (hiPSC-CM) are a promising in vitro tool for drug development and disease modeling, but their immature electrophysiology limits diagnostic utility. Tissue engineering approaches involving aligned 3D cultures enhance hiPSC-CM structural maturation but are insufficient to induce mature electrophysiology. We hypothesized that mimicking post-natal switching of the heart’s primary ATP source from glycolysis to fatty acid oxidation could enhance electrophysiological maturation of hiPSC-CM. We combined hiPSC-CM with microfabricated culture chambers to form 3D cardiac microphysiological systems (MPS) that enhanced immediate microtissue alignment and tissue specific extracellular matrix (ECM) production. Using Robust Experimental design, we identified a maturation media that improved calcium handling in MPS derived from two genetically distinct hiPSC sources. Although calcium handling and metabolic maturation were improved in both genotypes, there was a divergent effect on action potential duration (APD): MPS that started with abnormally prolonged APD exhibited shorter APD in response to maturation media, whereas the same media prolonged the APD in MPS that started with aberrantly short APD. Importantly, the APD of both genotypes was brought near the range of 270-300ms observed in human left ventricular cardiomyocytes. Mathematical modeling explained these divergent phenotypes, and further predicted the response of matured MPS to drugs with known pro-arrhythmic effects. These results suggest that systematic combination of biophysical stimuli and metabolic cues can enhance the electrophysiological maturation of hiPSC-derived cardiomyocytes. However, they also reveal that maturation-inducing cues can have differential effects on electrophysiology depending on the baseline phenotype of hiPSC-CM. In silico models provide a valuable tool for predicting how changes in cellular maturation will manifest in drug responsiveness.

While cardiomyocytes differentiated from human induced pluripotent stems cells (hiPSCs) hold great promise for drug screening, the electrophys-iological properties of these cells can be variable and immature, producing results that are significantly different from their human adult counterparts. Here, we describe a computational framework to address this limitation, and show how in silico methods, applied to measurements on immature cardiomyocytes, can be used to both identify drug action and to predict its effect in mature cells. Our synthetic and experimental results indicate that optically obtained waveforms of voltage and calcium from microphysiological systems can be inverted into information on drug ion channel blockage, and then, through assuming functional invariance of proteins during maturation, this data can be used to predict drug induced changes in mature ventricular cells. Together, this pipeline of measurements and computational analysis could significantly improve the ability of hiPSC derived cardiomycocytes to predict dangerous drug side effects.

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.

The conduction of electrical signals through cardiac tissue is essential for maintaining the function of the heart, and conduction abnormalities are known to potentially lead to life-threatening arrhythmias. The properties of cardiac conduction have therefore been the topic of intense study for decades, but a number of questions related to the mechanisms of conduction still remain unresolved. In this paper, we demonstrate how the so-called EMI model may be used to study some of these open questions. In the EMI model, the extracellular space, the cell membrane, the intracellular space and the cell connections are all represented as separate parts of the computational domain, and the model therefore allows for study of local properties that are hard to represent in the classical homogenized bidomain or monodomain models commonly used to study cardiac conduction. We conclude that a non-uniform sodium channel distribution increases the conduction velocity and decreases the time delays over gap junctions of reduced coupling in the EMI model simulations. We also present a theoretical optimal cell length with respect to conduction velocity and consider the possibility of ephaptic coupling (i.e. cell-to-cell coupling through the extracellular potential) acting as an alternative or supporting mechanism to gap junction coupling. We conclude that for a non-uniform distribution of sodium channels and a sufficiently small intercellular distance, ephaptic coupling can influence the dynamics of the sodium channels and potentially provide cell-to-cell coupling when the gap junction connection is absent., Author summary The electrochemical wave traversing the heart during every beat is essential for cardiac pumping function and supply of blood to the body. Understanding the stability of this wave is crucial to understanding how lethal arrhythmias are generated. Despite this importance, our knowledge of the physical determinants of wave propagation are still evolving. One particular challenge has been the lack of accurate mathematical models of conduction at the cellular level. Because cardiac muscle is an electrical syncytium, in which direct charge transfer between cells drives wave propagation, classical bidomain and monodomain tissue models employ a homogenized approximation of this process. This approximation is not valid at the length scale of single cells, and prevents any analysis of how cellular structures impact cardiac conduction. Instead, so-called microdomain models must be used for these questions. Here we utilize a recently developed modelling framework that is well suited to represent small collections of cells. By applying this framework, we show that concentration of sodium channels at the longitudinal borders of myocytes accelerates cardiac conduction. We also demonstrate that when juxtaposed cells are sufficiently close, this non-uniform distribution induces large ephaptic currents, which contribute to intercellular coupling.