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

I have read the journal’s policy and the authors of this manuscript have the following competing interests: KJ, SW and AT have financial relationships with Organos Inc, and the company may benefit from commercialization of the results of this research.

Article Type: research-article Article History
Abstract

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.

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.

Jæger,Wall,Tveito,and Beard: Computational prediction of drug response in short QT syndrome type 1 based on measurements of compound effect in stem cell-derived cardiomyocytes

1 Introduction

Short QT (SQT) syndrome is a cardiac channelopathy characterized by an abnormally short duration of the QT interval of the patient’s electrocardiogram (ECG) [13]. The syndrome is associated with increased risk of atrial fibrillation, ventricular arrhythmias and sudden cardiac death [4, 5] and was first described by Gussak et al. in 2000 [1]. The first identified SQT subtype, termed SQT1, results from an increase in the transmembrane potassium current, IKr, caused by a mutation in the IKr encoding gene KCNH2 [6]. Later, several additional subtypes of SQT syndrome have been identified, originating from mutations in other genes, including gain of function alterations in genes encoding the potassium channels responsible for the IKs and IK1 currents and loss of function alterations in genes encoding calcium channels [2, 3].

Because of the lethal risks associated with the syndrome, there is an urgent need for effective pharmacological therapies. One modern approach is to identify compounds that can make the electrophysiological properties of cardiomyocytes (CMs) populated with SQT mutated channels more similar to the electrophysiological properties of CMs populated with wild type (WT), unmutated, channels (see, e.g., [711]). For SQT1, which is characterized by an increased IKr current, drugs inhibiting the IKr current have been proposed as possible candidates (see, e.g., [11, 12]). However, the effect of a drug on WT IKr channels can be very different from the effect on mutated channels. For example, the IKr blockers sotalol and ibutilide have been shown to be ineffective for SQT1 patients [12]. This may be explained by the fact that these drugs primarily affect the inactivated state of the IKr channels, and the SQT1 mutation impairs the inactivation of the channels [6, 13, 14]. On the other hand, quinidine, which affects both the open and inactivated states of IKr channels has proven to be more effective for SQT1 patients [15]. These examples demonstrate that investigations into drug effects on specific SQT mutations are needed in order to find appropriate pharmacological treatments.

A promising prospect in the quest for suitable drugs for SQT syndrome is the development of human induced pluripotent stem cells (hiPSCs) (see, e.g., [1619]). These cells can be generated from individual patients and differentiated into a large number of different cell types, including CMs. Thus, cells can be generated from SQT patients, allowing for investigations of drug effects for CMs populated by ion channels affected by the patients’ specific mutations. Indeed, in [9, 11], several drugs have been tested on hiPSC-derived CMs (hiPSC-CMs) from an SQT1 patient, revealing possible promising drugs.

However, a difficulty with using hiPSC-CMs in drug testing applications is that the electrophysiological properties of hiPSC-CMs differ significantly from the electrophysiological properties of adult CMs (see, e.g., [20, 21]). In general, hiPSC-CMs are recognized as electrophysiologically immature, with properties more similar to those of fetal CMs. These differences imply that the drug response observed for hiPSC-CMs might not be the same as the corresponding drug response for adult native CMs.

Mathematical modeling has been proposed as a possible tool to help translate the drug response of hiPSC-CMs to the drug response of adult CMs (see, e.g., [2225]). The development of mathematical models of the dynamics underlying human cardiac action potentials is an active field of research, and a large number of models have been developed, including models for adult undiseased ventricular CMs [26, 27], atrial CMs [28, 29], hiPSC-CMs [30, 31] and CMs affected by mutations [32, 33]. In these models, changes in the membrane potential are represented by individual transmembrane ionic currents (see, e.g., (1) below), and the models can therefore be useful for investigating how drug effects on individual ion channels affect the composite dynamics of the full action potential. Furthermore, the action potential models can be combined with spatial models of electric conduction (see, e.g., [3438]) to provide insight into drug effects on cardiac conduction properties and mechanisms of arrhythmia.

In previous studies, we have applied a procedure based on mathematical action potential modeling to estimate drug effects on individual ion channels from measurements of hiPSC-CMs and predict adult CM drug responses [22, 23]. This procedure is based on the assumption that the function of individual proteins, e.g. ion channels, is the same for hiPSC-CMs and adult CMs, and that only the density of the proteins differs between hiPSC-CMs and adult CMs (see Fig 2 below). From this assumption, it follows that the effect of a drug on an individual ion channel is the same for hiPSC-CMs and adult CMs. Assuming that we have correctly identified the drug effect on individual channels in the hiPSC-CM case, we can therefore directly translate this effect to the adult case by inserting the inferred mechanisms into a model for adult cells.

The aim of the present study is to use this computational procedure to predict drug effects for adult SQT1 CMs based on measurements of drug effects on the action potential of SQT1 hiPSC-CMs from [9, 11] and to extend these results into prediction of patient QT changes. Our overall computational pipeline is depicted in Fig 1. We consider the four drugs quinidine, ajmaline, mexilietine and ivabradine, shown to potentially be useful for SQT1 patients in [9, 11] and show predicted drug responses for the drugs on the ventricular action potential and QT interval of adult patients. We validate our pipeline using data for quinidine, where measurements of the drug effect on the QT interval have been conducted for adult patients [15]. The predicted drug response turns out to be in good agreement with the measured drug effect, indicating that the computational procedure could be useful for predicting adult drug responses from measurements of hiPSC-CMs.

Illustration of the computational pipeline.
Fig 1

Illustration of the computational pipeline.

1) Biomarkers from the cardiac AP are taken from hiPSC-CMs under drug testing. 2) These biomarkers from dose escalation studies are inverted into an SQT1 model of the AP of hiPSC-CMs. Inversion into a matched model provides determination of drug effects on specific channels [23]. 3) The drug effects determined in 2 are inserted into a model of adult CMs with the same SQT1 mutation to give a prediction of drug effect in mature CMs [22]. 4) The adult CM model is converted into pseudo-ECG waveforms for prediction of QT segment changes in SQT1 patients under the estimated effect of the drug.

2 Methods

In this section, we describe the methods applied in this study. We start by describing the basic modeling assumptions underlying the approaches used for computational identification of drug effects and mapping of drug effects from hiPSC-CMs to adult CMs. Then, we describe details of the applied action potential model and inversion procedure. Finally, we describe the approach used to estimate the QT interval in the adult case. Note that the majority of these methods are to a large extent based on the methods described in [23].

2.1 Action potential models

In this section, we describe the framework for modeling the action potentials of hiPSC-CMs and adult CMs, with and without the SQT1 mutation, and the relationship between these models. Three of the main modeling assumptions are illustrated in Fig 2, building on the approaches of [22, 23]. In the next subsections we will explain these assumptions and demonstrate how they affect the action potential modeling.

Illustration of the assumptions underlying the computational maturation approach.
Fig 2

Illustration of the assumptions underlying the computational maturation approach.

1) The density of different types of ion channels (and other membrane or intracellular proteins) may differ between hiPSC-CMs and adult CMs, but the function of the individual channels is the same. In the model, the density difference is represented by the parameter λ. 2) The SQT1 mutation affects the individual IKr channels in exactly the same manner for hiPSC-CMs and adult CMs. In the model, the mutation is represented by an adjusted model for the open probability, oKr. 3) The effect of a drug on a single protein is the same for hiPSC-CMs and adult CMs. In the model, the drug effect is represented by the parameter ε.

2.1.1 Action potential modeling framework

In the action potential model, the membrane potential, v (in mV), is governed by an equation of the form

where Ij are different membrane current densities. These current densities can be expressed on the form
where Nj is the number of proteins of type j on the membrane, A is the area of the membrane, Cm is the specific membrane capacitance, and ij is the average current through a single protein of type j. For currents through voltage-gated ion channels, this ij is typically given on the form
where g0,j is the conductance through a single open channel, oj is the open probability of the channel, and Ej is the equilibrium potential of the channel. In this case, it is common to introduce combined parameters of the form
and write (2) on the form

2.1.2 Assumption 1: Functional invariance of wild type (WT) ion channels during maturation

A number of electrophysiological properties have been shown to differ between hiPSC-CMs and adult CMs (see, e.g., [20, 21]). In addition, the electrophysiological properties of different samples of WT hiPSC-CMs often vary significantly (see, e.g., [39, 40]). We assume that these differences in electrophysiological properties (both between samples of hiPSC-CMs and between adult CMs and hiPSC-CMs) are due to:

    Differences in the geometry of the cells,

    Differences in the number of membrane proteins, like ion channels, pumps and exchangers,

    Differences in the number of intracellular proteins, like calcium buffers, ryanodine receptors and SERCA pumps.

However, we assume that the function of the individual membrane and intracellular proteins is the same for the different WT cases.

Considering the model (1) and (2), these assumptions imply that the density ρj=NjA may differ between cells, but that ij remains the same. Parameterizing the model to a specific set of measurements can then be accomplished by adjustments of the densities represented by adjustment factors λj, such that

where ρj* is the default density of proteins of type j. Incorporating these adjustment factors into the model (1), specific adult and hiPSC-CM versions of the model are given by
respectively, where λjA and λjhiPSC specify the adjustment factors for the protein densities in the two cases. For currents through ion channels, these adjustment factors can be incorporated by adjusting the conductances, gj (see (4) and (5)). Similar adjustment factors can be set up for the density of intracellular proteins and the cell geometry, see [23]. A list of the maturation-dependent parameter values of the base model are given in Table IV of the S1 Appendix.

2.1.3 Assumption 2: Functional invariance of mutated ion channels during maturation

We assume that a mutation affects an individual ion channel in exactly the same manner for hiPSC-CMs and adult CMs. SQT1 is known to affect the IKr current [3, 41], and we assume that the effect on the current can be represented by adjusting the model for the open probability, oKr (see (3)). More specifically, we adjust the voltage dependence of the steady state inactivation gate, xKr2, by shifting the inactivation towards more positive potentials, as described for the SQT1 mutation N588K in [3]. Because we assume that the mutation affects an individual IKr channel in the same manner for hiPSC-CMs and adult CMs, we apply the same adjustment of the model for oKr in the hiPSC and adult SQT1 cases. In addition, we assume that the SQT1 mutation can be represented solely by this adjustment of oKr, and that the density of proteins, including the density of IKr channels, is the same in the WT and SQT1 cases.

2.1.4 Assumption 3: Identical drug effects for identical proteins

The third assumption illustrated in Fig 2 is that since the function of a single ion channel is the same for hiPSC-CMs and adult CMs, the effect of a drug on an individual ion channel will also be identical for hiPSC-CMs and adult CMs. This means that if we are able to determine the effect of a drug on an individual channel in the hiPSC-CM case, we can find the effect of the same drug in the adult case by incorporating the same single channel drug effect in the adult version of the model.

Following [23], we use a simple IC50-based modeling approach to represent the effect of a drug. That is, we assume that the average single channel current through a channel j in the presence of the drug dose D is given by

Here, ij(0) is the average single channel current when no drug is present and εj (in μM−1) represents the effect of the drug, defined as

where IC50j (in μM) is the drug concentration that blocks the current j by 50%. Incorporating drug effects into the hiPSC-CM and adult CM models (7) and (8), we obtain
where we note that εj is the same in the two versions of the model.

In this study, we will use this modeling framework to estimate the drug effect (in the form of ε values) of drugs based on action potential measurements of hiPSC-CMs with the SQT1 mutation N588K. Next, we will insert the estimated ε values into a model for adult ventricular cells with the same mutation to estimate the effect of the drug for an adult patient.

2.1.5 Adult and hiPSC-CM base model variability

In our computations, we use action potential measurements of hiPSC-CMs to predict drug effects for adult patients. In this procedure, we need to fit the hiPSC-CM base model (8) to match data of specific samples of hiPSC-CMs. The properties of these hiPSC-CMs tend to vary between experiments, also in the control case (see, e.g., the biomarkers for the zero dose cases in Fig 7). Therefore, we let the values of the adjustment factors, λjhiPSC, vary for each experiment (that is, for each considered drug). The adjustment factors are, however, assumed to be the same in the control case and for all doses of a specific drug. Following [23], in the adult case, we consider only one default case defined by the parameters λjA. In other words, the adult base model is assumed to be the same in all cases, regardless of the variability in the parameters λjhiPSC.

2.1.6 Base model formulation

In order to represent the action potentials of adult ventricular CMs and hiPSC-CMs (both with and without the SQT1 mutation N588K), we use a slightly modified version of the base model formulation from [23], following the form (1) and (2). Three main adjustments have been incorporated into the current version of the model. First, we have modified the steady-state values of the gating variables of the IKr current to better fit measured IKr currents in the WT and SQT1 cases from [14] (see Fig 5). Second, we have extended the temperature-dependence of the model by including Q10 values for a number of the currents, following [42, 43]. This was done because we consider two different temperatures in our computations. In the adult case, we assume body temperature (310 K), while for the hiPSC-CM case, we assume room temperature (296 K) since the considered hiPSC-CM data (from [9, 11]) are obtained at room temperature. Third, we have incorporated a dynamic model for the intracellular Na+ concentration. This was done in order to make the frequency dependent QT interval changes more physiologically realistic (see Fig 11). The full base model formulation is found in the S1 Appendix. The system of ordinary differential equations (ODEs) defined by the base model is solved using the ode15s solver in MATLAB. The MATLAB code for the base model is found in S1 Code.

2.1.7 Stimulation protocol

In the adult case, we use 1 Hz pacing unless otherwise specified. Furthermore, we run each simulation for at least 500 pacing cycles after each parameter change in order to obtain new stable solutions before recording the action potentials (and pseudo-ECGs, see Section 2.4 below).

In the hiPSC-CMs case, we use 0.2 Hz pacing, as specified in [11]. As a compromise between the need to obtain new stable solutions after each parameter change and the need for reducing the computing time of the inversion procedure, we run the simulation for 30 pacing cycles before measuring the action potential for each iteration in the continuation method used for optimization (see Section 2.3.2). However, the parameter changes between each iteration of the continuation method are expected to be quite small and we update the initial conditions between each of the 20 iterations. Therefore, the final solutions of the inversion are expected to be quite close to the steady state solutions for the applied parameters. This has also been confirmed by numerical experiments. For example, the APD90 value for the inversion of data for 10 μM of quinidine was 323.7 ms using 30 pacing cycles from the states saved in the second-to-last continuation iteration. When the simulation was allowed to continue for 500 pacing cycles, the APD90 value changed by only 0.5% to 325.5 ms.

2.1.8 Parameterization of default adult and hiPSC-CM base models

In order to parameterize the default adult and hiPSC-CM versions of the base model for WT and SQT1, the model parameters are adjusted by hand. For the hiPSC-CM case, adjustments are made to the conductance of all membrane currents. The purpose of the adjustments is to obtain a model in rough agreement with some measured current densities and AP characteristics for WT and SQT1 hiPSC-CMs from [9] to use as a suitable starting point for the inversions of hiPSC-CM data. A comparison of the AP characteristics and current densities from [9] to the corresponding values for the default WT and SQT1 versions of the hiPSC-CM model are given in S1 Fig of the Supporting Information.

In the adult case, the base model parameters are fitted to information about the heart rate dependent value of the QTp interval (see Section 2.4) measured for adults with and without the SQT1 mutation from [15]. Because the intracellular Na+ and Ca+ dynamics are important for the frequency response of action potential models [27], the conductance of the currents responsible for the balance of intracellular Ca2+ and Na+ concentrations are adjusted (i.e., the conductance of ICaL, IbCa, IpCa, INaCa, INa, and INaK) in order to fit the model to the heart rate dependent QTp interval. In addition, to achieve the measured difference between WT and SQT1 QTp intervals, the conductance of the different repolarizing currents IKr, IKs, and IbCl are adjusted. We also attempt to make the current densities during an action potential relatively close to the current densities of the well-established adult ventricular AP models of Grandi et al. and O’Hara et al. [26, 27]. A comparison of some of the main current densities between the adult WT base model and the Grandi et al. and the O’Hara et al. models are given in S2 Fig of the Supporting Information.

2.2 Inversion of IKr measurements

In order to represent the WT and SQT1 IKr currents, parameters of the model for the IKr open probability are fitted to data of WT and SQT1 IKr currents from [14]. In this data, it is revealed that the steady state inactivation of the current is shifted towards higher potentials in the SQT1 (N588K) case compared to the WT case (see the right panel of Fig 5 below). Therefore, we assume that the mutation can be represented in the model by an adjusted model for the steady state value of the inactivation gate, xKr2. To find a model matching the data as well as possible, we introduce six free parameters c1c6 in the steady state activation and inactivation gates of the IKr current on the form:

See S1 Appendix for a description of the full IKr model formulation. To find optimal values of the six parameters c1c6, we run simulations of the voltage clamp protocol applied in [14]. That is, we first fix the membrane potential at -80 mV and run a simulation to steady state. Then, we increase the membrane potential to a value between -50 mV and 100 mV and compare the current after 2 seconds of simulation to the corresponding current reported in [14]. In addition, we compare the steady state values of the inactivation gate xKr2 directly to the steady state inactivation measurements reported in [14].

In the simulations trying to adjust the IKr model to measurements from [14], we adjust the temperature and intracellular and extracellular potassium concentrations in the model to match those reported for the experiments, that is T = 37°C, [K+]i = 130 mM, [K+]e = 4 mM.

We find the optimal parameters c = (c1, …, c6) by minimizing a cost function of the form

using the Nelder-Mead algorithm [44]. Here, vi are each of the values of the membrane potential considered in the data set, Id(vi) are each of the measured data points for relative IKr currents, and Ib(c, vi) are the corresponding currents computed for the model specified by the parameters c.

2.3 Inversion of action potential measurements

The inversion of data of SQT1 hiPSC-CMs exposed to various drugs is performed using the inversion procedure described in [23]. In the inversion, the adjustment factors λKr, λCaL, λNa, λK1, λf, λNaCa, λNaK, and λbCl for the currents IKr, ICaL, INa, IK1, If, INaCa, INaK, and IbCl, and the combined adjustment factor λB, adjusting all the intracellular calcium buffers as in [23] are treated as free parameters. In addition, we assume that the drugs affect the IKr, ICaL or INa currents for quinidine, ajmaline and mexiletine and the IKr, If or INa currents for ivabradine. Therefore, we introduce the drug parameters εKr, εf and εNa for ivabradine and εKr, εCaL, and εNa for the remaining drugs. For a given drug, these twelve parameters are fitted to data simultaneously using information about both the control case and each drug dose included in the data set as explained below.

2.3.1 Cost function definition

In the inversion procedure, we wish to minimize a cost function of the form

where Dd represent each of the drug doses included in the data set (including the control case, D0 = 0), Hi represent the different cost function terms included in the cost function, wd,i are weights for each of the cost function terms and doses and Hreg is a regularization term (see (19) below). The cost function consists of terms of the form
where Ri(λ, ε, Dd) is a biomarker computed from the solution of the model specified by λ and ε for the drug dose Dd, and Ri*(Dd) is the corresponding measured biomarker of the data we are trying to invert.

Considered biomarkers. We consider each of the biomarkers, Ri, included in the data sets from [9, 11], i.e., APD50, APD90, APA, dvdt and RMP. These biomarkers are illustrated in the left panel of Fig 3. Here, the maximal upstroke velocity, dvdt (in mV/ms), is defined as the maximum value of the derivative of the membrane potential with respect to time, the resting membrane potential, RMP (in mV), is defined as the minimum value of the membrane potential, the action potential amplitude, APA (in mV), is defined as the difference between maximum and minimum values of the membrane potential, and the APD50 and APD90 values are defined as the time (in ms) from the time point of the maximum upstroke velocity to the time point where the membrane potential first reaches a value below 50% or 90%, respectively, of the action potential amplitude. All biomarkers are computed directly from the model solution, found using an adaptive time step in MATLAB’s ODE solver ode15s.

Illustration of considered biomarkers.
Fig 3

Illustration of considered biomarkers.

Left: Illustration of the five biomarkers included in the cost function of the inversion procedure: The action potential durations (in ms) at 90% and 50% repolarization (APD90 and APD50, respectively), the maximal upstroke velocity (dvdt, in mV/ms), the action potential amplitude (APA, in mV) and the resting membrane potential (RMP, in mV). Right: Illustration of the QTp interval (in ms) computed from a pseudo-ECG waveform. The QTp interval is defined as the time from the onset of the QRS complex to the peak of the T-wave.

Cost function weights. We use the weight 3 for HAPA, the weight 6 for HAPD90 and the weight 1 for the remaining terms. In addition, the weights for the control case, w0,j, and the weight for HAPD90 for the largest dose are multiplied by the total number of doses included in the data set. The term HRMP is only included for the drug quinidine because it was not specified in the data set for the remaining drugs [9, 11]. The specific weights used for the cost function terms are chosen to prevent large errors for single biomarkers and to give an overall acceptable fit.

Regularization term. The regularization term, Hreg(λ, ε), is defined as

Here, the first term is defined to make the inversion procedure avoid solutions with unrealistically large perturbations of some of the currents, and the second term is defined to make the inversion select small drug perturbations if small and large perturbations result in almost indistinguishable solutions (D˜ is the median of the considered drug doses of the data set). In our computations, we set the weights for the terms to wλ = 1 and wε = 0.001. We let the set Sε consist of all the currents we assume may be affected by the drug and the set Sλ consist of the IKr and ICaL currents. Note that we here assume that λj = 0 defines the starting point of the inversion procedure (i.e., the default hiPSC-CM base model parameterization).

2.3.2 Optimization method

We apply the continuation-based optimization algorithm from [23] to minimize the cost function (17). A simple example application of the continuation algorithm with two free parameters fitted to a single AP waveform is illustrated in Fig 4. In short, the algorithm consists of introducing a parameter, θ, which is gradually increased from 0 to 1 in M iterations (or θ-steps). For each θ-step, we seek the optimal parameters fitting a temporary objective that is, as θ is increased, gradually adjusted from the default model solution to the actual data we are trying to invert. More specifically, for a given θ ∈ [0, 1], we try to minimize the cost function (17) and (18) with the biomarkers given by

where Ri(0) are the biomarkers of the default model used as a starting point for the optimization, Ri* are the biomarkers of the data and Ri*(θ) define a temporary objective. The purpose of gradually stepping from the default model solution to the data in this manner is that we know the optimal solution for θ = 0, and for each θ-step, we assume that the optimal parameters are quite close to the optimal parameters of the last step. Therefore, we can keep the space in which to search for optimal parameters quite small in each iteration.

Illustration of the continuation algorithm used for optimization in the inversion method.
Fig 4

Illustration of the continuation algorithm used for optimization in the inversion method.

A) The problem is defined by a default model and some data we are trying to invert by finding an optimal model parameterization fitting the data. B) In the continuation algorithm, we seek temporary optimal parameters in M iterations (θ-steps). The objective for each θ-step is gradually changed from the default model to the data we are trying to invert. C) In each θ-step, we look for optimal parameters for the temporary objective by drawing NG random guesses in the vicinity of the optimal parameters from the previous θ-step. For each random guess, we run NNM Nelder-Mead iterations, and from the result, we select the best fit as the new optimal parameters. D) The final parameterization is given by the optimal parameters found in the last θ-step.

In our computations, we use M = 20 continuation iterations (θ-steps) with NG = 100 or 200 randomly chosen initial guesses for the first fifteen and the last five iterations, respectively. The initial guesses for λ are chosen within 20% above or below the optimal values from the previous iteration, and the initial guesses for εm are chosen within [εm−1/5, 5εm−1], where εm−1 is the optimal ε from the previous iteration. From these initial guesses we run NNM = 30 or 60 iterations of the Nelder-Mead algorithm [44] for the first fifteen and the last five iterations, respectively.

We have chosen to use the continuation-based optimization algorithm because we have previous experience from [23] using this algorithm for this type of optimization problems. Furthermore, in S4 Fig of the Supporting Information, we show a comparison of the continuation method with other optimization methods. The other methods are provided in MATLAB’s Global Optimization Toolbox [45]. We compare the methods for cases using simulated data and thus the exact minimum is known. This allows for comparison of the methods and it turns out that the continuation method compares well with the other methods, especially when the number of free parameters increases.

2.4 Computation of pseudo-ECG and QT interval

In order to estimate drug effects on adult QT intervals, we apply a simple pseudo-ECG calculation applied in a number of earlier studies (e.g., [4650]). The calculation follows a two-step procedure.

Step 1

The membrane potential and membrane currents along a strand of cylindrical cells are computed using the cable equation as explained i detail in [35]. In this step, cells are connected to each other by gap junctions, each cell is assumed to be isopotential, and the extracellular potential is assumed to be constant in space. The membrane potential vk in each cell k is modeled by

where r is the cell radius, L is the cell length, RCG is the ratio between the capacitive and the geometrical cell areas (see, e.g., [35]), Cm is the specific cell capacitance and σik-1/2 is the averaged intracellular conductivity between cells k − 1 and k. This averaged intracellular conductivity is given by
where Rmyo is the myoplasmic resistance and Rg is the gap junction resistance. Furthermore, Iionk is the sum of the ionic current densities of the cell, modeled by the base model (i.e., ∑j Ij in (1)).

The ODE system defined by (21) is solved using the ode15s solver in MATLAB, and the solution is used to compute the membrane currents

originating from each cell i at each time point n.

Step 2

The extracellular potential originating from the cell strand is computed for each time step n using the so-called point-source approximation (see, e.g., [46, 51]),

where σe is the extracellular conductivity and |rrk| is the Euclidean distance between the point at which we are measuring the extracellular potential, r, and the center of cell number k, rk.

Parameter values

The parameters of the pseudo-ECG simulations are based on the parameters of earlier computations of pseudo-ECGs (e.g., [35, 46, 49]). We consider a cell stand of 100 cells of length L = 150μm and radius r = 10μm. We let the cell strand exhibit transmural heterogeneity of ion channel density, with an endocardial region consisting of the first 25 cells, a midmyocardial region for the next 35 cells and an epicardial region in the last 40 cells. The default adult base model defines the ion channel density in the epicardial region. In the endocardial region, the Ito and IKs channel densities are reduced to 1% and 31%, respectively, compared to the epicardial region. Similarly, in the midmyocardial region, the Ito and IKs channel densities are reduced to 85% and 11%, respectively, compared to the epicardial region.

Furthermore, we set Cm = 1μF/cm2, RCG = 2, Rmyo = 0.15 kΩcm, and Rg = 0.0015 kΩcm2, resulting in a conduction velocity along the cell strand of approximately 52 cm/s, close to experimentally measured conduction velocities from literature (∼50 cm/s [52]). Stimulation is applied for the first two cells of the endocardial region, and the extracellular potential is measured 2 cm from the end of the epicardial region. Because we in this study only are interested in the time course of the ECG and not the amplitude in mV, we do not assign a specific value for σe and report the computed pseudo-ECG without numbers on the ue-axis.

Definition of the QT interval

After computing a pseudo-ECG using the described approach, the QT interval is computed from the ECG waveforms. The QT interval is often defined as the time from the start of the QRS interval to the end of the T-wave (see Fig 3). However, in the study [15], which we will use for comparison to the computational results, an alternative definition of the so-called QTp interval is applied, defined as the time from the start of the QRS complex to the peak of the T-wave, as illustrated in the right panel of Fig 3. We will therefore use this definition in this study. Note that for the SQT1 case, the computed T-wave is inverted to have the opposite sign of the QRS complex (see Fig 6). In that case, we define the peak of the T-wave as the time point when the minimum, i.e. the maximum absolute deviation of the T-wave, is reached.

3 Results

In this section, we describe the results of our applications of the above mentioned methods. First, we set up the default hiPSC-CM and adult base models for WT and SQT1 based on measurements of WT and SQT1 IKr currents from [14]. Next, we apply the inversion procedure to identify the effect of four drugs on individual ion channels based on action potential measurements of SQT1 hiPSC-CMs from [9, 11]. Finally, we estimate the corresponding drug effect for adult patients with short QT syndrome using these identified drug effects.

3.1 Computational representation of the SQT1 mutation

In order to represent the SQT1 mutation N588K in the model, we first adjust the model of the open probability of the IKr current to data for WT and SQT1 IKr currents from [14] as described in Section 2.2. The optimal values for the model given on the form (13)(15) returned by the inversion procedure are c1 = −2.7, c2 = −15.3, c3 = 70, c4 = 20.9, c5 = −62, c6 = 1.85. In Fig 5, the fitted model solutions of IKr in the WT and SQT1 cases are compared to the data from [14]. In the right panel, we observe that the inactivation is shifted towards higher values of the membrane potential in the SQT1 case.

Representation of the SQT1 mutation N588K in the IKr model.
Fig 5

Representation of the SQT1 mutation N588K in the IKr model.

Left panel: Steady state IKr currents obtained at different fixed values of the membrane potential divided by the currents obtained at v = 0 mV. The results obtained for the fitted WT and SQT1 IKr models are compared to corresponding data from [14]. Right panel: Comparison of the steady state inactivation gate in the WT and SQT1 models of IKr and steady state inactivation data from [14]. In the SQT1 case, the inactivation is shifted towards higher values of the membrane potential. Data used in this figure can be found in S1 Data.

The WT and SQT1 versions of the IKr model fitted to measurements from [14] are inserted into the base model formulation. The parameters, λhiPSC, for the hiPSC-CM versions of the model are then fitted to action potential measurements of WT and SQT1 hiPSC-CMs from [9]. In addition, the parameters of adult versions of the base model are fitted to information about the QTp interval for adults with and without the SQT1 mutation from [15]. These adjustments are made by hand-tuning the parameters of the default hiPSC-CM and adult versions of the base model from [23], as described in Section 2.1.8.

The action potentials of the resulting models are plotted in the upper panel of Fig 6, along with the computed pseudo-ECG for the adult case. The pseudo-ECG is computed using the approach described in Section 2.4. Note that the only difference between the WT and SQT1 versions of the models is a difference in the steady state inactivation gate of the IKr current (see (15)), and that the remaining model parameters (including the density of IKr channels) are the same in the WT and SQT1 cases. Note also that the adjustment of the IKr steady state inactivation gate is exactly the same in the hiPSC-CM and adult versions of the model (following assumption 2 in Fig 2), but that the geometry of the cell and the density of different types of membrane and intracellular proteins are different between the models for hiPSC-CMs and adult CMs (following assumption 1 in Fig 2).

Properties of the base models for hiPSC-CMs and adult ventricular CMs in the WT and SQT1 cases.
Fig 6

Properties of the base models for hiPSC-CMs and adult ventricular CMs in the WT and SQT1 cases.

Upper panel: Comparison of the action potentials computed for the WT and SQT1 versions of the hiPSC-CM (left) and adult (center) models, in addition to the pseudo-ECG for the adult model (right). The only difference between the formulations of the WT and SQT1 models is a shift in the inactivation gate of IKr as illustrated in the right panel of Fig 5. Lower panel: APD90 values and QTp intervals computed using the WT and SQT1 versions of the models. The computed APD90 values for hiPSC-CMs are compared to data from [9], and the computed QTp intervals for adults are compared to data from [15]. Data used in this figure can be found in S1 Data.

In the lower panel of Fig 6, we report the APD90 values of the WT and SQT1 models in the hiPSC-CM and adult cases. In addition, we report the QTp interval computed from the pseudo-ECG (see Fig 3). In the hiPSC-CM case, we compare the APD90 values of the model to the corresponding values reported in [9]. The APD90 value is reduced from 324 ms in the WT case to 188 ms in the SQT1 case, in good agreement with the values reported in [9]. In the adult case, the APD90 value is similarly reduced from 272 ms to 189 ms. Furthermore, the QTp interval is reduced from 292 ms in the WT case to 200 ms in the SQT1 case, close to reported QTp values from [15]. Moreover, in Figure 1 of [9], an ECG from the SQT1 patient whose skin fibroblasts were used to generate the SQT1 hiPSC-CMs of the study is reported. Considering the signal recorded in lead V3 (as in [15]), the QTp interval can be estimated to be about 170 ms. Correcting for heart rate using Bazett’s formula and assuming a heart rate of 60 beats per minute (as is used in the simulation) the QTp interval is about 190 ms, similar to the value reported in [15] and close to the value given by the model (200 ms). We also note that in the computed pseudo-ECG, the T-wave is positive in the WT case and negative in the SQT1 case. A note on this property is found in the S2 Appendix.

3.2 Computational identification of drug response from membrane potential measurements of SQT1 hiPSC-CMs

In [9, 11], action potential biomarkers are reported from measurements of hiPSC-CMs with the SQT1 mutation N588K exposed to drugs attempting to make the properties of the SQT1 cells more similar to those of WT cells. In this paper, we wish to use this data to estimate the corresponding drug responses for adult SQT1 patients. In order to make these predictions, we first need to identify the effect of the drugs on individual ion channels in the hiPSC-CM case, based on the data provided in [9, 11]. These drug effects are predicted using the inversion procedure described in Section 2.3 using the SQT1 hiPSC-CM model from Fig 6 as a starting point for the inversion.

Fig 7 shows how well the biomarkers of the fitted SQT1 hiPSC-CM model match the corresponding biomarkers reported in the data. We observe that the biomarkers of the fitted models are quite similar to the values of the data. In particular, the model seems to capture the APD90 increase resulting from the drugs quite well. In the lower panel of Fig 8, the action potentials of the SQT1 hiPSC-CM models fitted to the data are plotted for the control case and for each of the drug doses, along with the action potential of the default base model for WT hiPSC-CM. We observe that each of the drugs increases the action potential duration, but not enough to fully recapture the WT hiPSC-CM action potential length.

Comparison of measured and computed action potential biomarkers.
Fig 7

Comparison of measured and computed action potential biomarkers.

We consider the biomarkers reported in the data from [9, 11] (green) and computed from the fitted SQT1 hiPSC-CM models returned by the inversion procedure described in Section 2.3 (purple) for the drugs quinidine, ivabradine, ajmaline and mexiletine. Note that the definition of each of the biomarkers are illustrated in Fig 3 and that RMP data are only included for the quinidine case. Data are shown as the mean ± SEM (standard error of the mean). Data used in this figure can be found in S1 Data.

Result of the inversion procedure applied to data for the drugs quinidine, ivabradine, ajmaline and mexiletine from [9, 11].
Fig 8

Result of the inversion procedure applied to data for the drugs quinidine, ivabradine, ajmaline and mexiletine from [9, 11].

Upper panel: Predicted drug effect on individual ion channels in the form of ε-values returned by the inversion procedure. Lower panel: Action potentials of the fitted SQT1 hiPSC-CM models in the control case and for doses of each drug. For comparison, we also show the action potential of the WT hiPSC-CM model. Action potential characteristics of the fitted models in the lower panel are compared to data from [9, 11] in Fig 7. Data used in this figure can be found in S1 Data.

The upper panel of Fig 8 reports the drug effects identified by the inversion procedure in the form of ε values. Recall that εj is defined as 1/IC50j, where IC50j is the drug concentration that blocks the current j by 50%. In other words, a high ε value corresponds to a significant drug effect. We observe that the inversion procedure predicts that all four drugs have a significant effect on the IKr current. In addition, ajmaline is estimated to considerably block the ICaL current, and mexiletine is estimated to block INa.

In [11], measured drug effects on SQT1 IKr currents are reported for the maximum dose of each of the considered drugs. In order to get an impression of the accuracy of the IC50 values estimated by the inversion procedure, we compare the block percentages reported in [11] to the corresponding block percentages estimated by the computational procedure. The results are displayed in Fig 9, and we observe that the estimated block percentages seem to be in good agreement with the measured values.

Comparison of reported and estimated block percentages for IKr.
Fig 9

Comparison of reported and estimated block percentages for IKr.

We consider the reported block percentages for 10 μM quinidine, 10 μM ivabradine, 30 μM ajmaline and 100 μM mexiletine from [11] and the block percentages estimated based on the IC50 values identified in the inversion procedure. Data are shown as the mean ± SEM. Data used in this figure can be found in S1 Data.

3.3 Estimation of drug response for adult SQT1 patients

Following assumption 3 in Fig 2, we assume that the effect of a drug on a single ion channel is the same for hiPSC-CMs and adult CMs. Using this assumption, we can estimate the effect of the drugs quinidine, ivabradine, ajmaline and mexiletine for adult SQT1 CMs by inserting the ε values identified in Fig 8 into the adult SQT1 base model. In Fig 10, we report the resulting estimated adult drug effects. In the upper panel, we show the ventricular action potentials for the SQT1 case with no drug and with each of the considered drug doses present. In addition, we show the action potential in the WT case. The lower panel of Fig 10 displays the adult QTp intervals computed from pseudo-ECG simulations as explained in Section 2.4. We observe that the drugs increase the QTp interval for the SQT1 case to be more similar to the WT QTp interval. For 10 μM of quinidine, the QTp interval is estimated to increase from 200 ms to 257 ms. Similarly, the QTp interval is predicted to increase to 238 ms for 10 μM of ivabradine, and 221 ms for 100 μM of mexiletine. For comparison, the WT QTp interval is 292 ms. For 30 μM of ajmaline, the QTp interval is predicted to increase to only 209 ms. This is probably caused by the fact that ajmaline is estimated to have considerable effect on both the depolarizing current ICaL and the repolarizing current IKr, and that these two drug effects cancel each other out in the adult model.

Estimated drug effects for adult patients with the SQT1 mutation.
Fig 10

Estimated drug effects for adult patients with the SQT1 mutation.

The figure displays the estimated drug effect of quinidine, ivabradine, ajmaline and mexiletine for adult patients with the SQT1 mutation (and a heart rate of 60 beats/min) based on the result of the inversion of hiPSC-CM data reported in Fig 8. Upper panel: Estimated drug effect on the adult ventricular action potential. Lower panel: Estimated drug effect on the QTp interval. See Fig 3 for an illustration of the definition of the QTp interval. Note that the blue lines represent the adult WT case. Data used in this figure can be found in S1 Data.

The effect of the drug quinidine on the QTp interval of adult SQT1 patients is reported for a number of different heart rates in [15]. In Fig 11, we compare the predicted adult drug response of quinidine from our computational procedure to the measured responses in [15]. In the study [15], mean serum concentration was reported as 2.1 mg/L. Applying the molecular weight of quinidine of 324.4 g/mol [53], we consider the drug concentration 6.5 μM in the computational model. Moreover, we use the ε values identified based on measurements of SQT1 hiPSC-CMs from [9] (see Fig 8). The QTp intervals are computed for different pacing frequencies using the pseudo-ECG approach described in Section 2.4. In Fig 11, we observe that the QTp intervals predicted by the model for different heart rates fit well with the measured data in both the WT case, in the SQT1 case and in the SQT1 case with quinidine.

Comparison of measured and estimated adult QTp intervals.
Fig 11

Comparison of measured and estimated adult QTp intervals.

We compare the QTp intervals (see Fig 3) estimated by the computational procedure and the measured QTp intervals from [15] for different heart rates. We consider the WT case, the SQT1 case and the SQT1 case with cells exposed to 6.5 μM quinidine. The drug effect estimated by the procedure is based on measured drug effects for SQT1 hiPSC-CMs from [9]. Data used in this figure can be found in S1 Data.

4 Discussion

In this paper, we have applied the computational procedure introduced in [22, 23] to data of hiPSC-CMs derived from an SQT1 patient, found in [9, 11]. In this procedure, we first identified the drug response on individual ion channels based on measurements of drug effects for SQT1 hiPSC-CMs and then estimated the drug effect for an adult SQT1 patient by inserting these single channel drug effects into a model for adult cells.

4.1 Computational representation of the SQT mutation

Short QT syndrome is characterized by an abnormally short duration of the QT interval, and several different genetic mutations have proven to give rise to different subtypes of the syndrome [13]. SQT1–SQT3 are caused by gain of function mutations in the potassium channels responsible for the IKr, IKs and IK1 currents, respectively, while SQT4–SQT6 are associated with loss of function mutations in genes encoding calcium channels. Furthermore, SQT7 and SQT8 have been characterized by a loss of function mutation in a gene encoding sodium channels and a mutation in a gene encoding the cardiac Cl/HCO3 exchangers AE3, respectively [54]. Because of the different effects of different mutations, it is likely that pharmacological treatment of SQT syndrome must be tailored to the individual subtypes. Furthermore, in order to study the effects of mutations computationally, different mathematical representations are needed for the mechanistic causes of the shortening.

In this study, we have considered the SQT1 subtype, caused by a mutation (N588K) leading to increased IKr currents. Based on measurements of WT and mutated IKr currents [14], we represented the mutation by shifting the inactivation of the channels towards more positive potentials (see (15) and Fig 5). This computational representation of the mutation follows the same lines as [49, 55, 56], which used similar adjustments of the model for the IKr, IKs and IK1 open probabilities to represent the SQT1, SQT2 and SQT3 subtypes, respectively.

In the base model used in this study (see S1 Appendix), as well as in other commonly used action potential models (e.g., [26, 27, 30, 57]), the open probability of voltage-gated ion channels are governed by gating variables following the Hodgkin-Huxley formalism. However, more detailed and versatile models for the open probability have been introduced in the form of more complicated Markov models (see, e.g., [58, 59] for a comparison of the two formalisms). Such Markov models have been suggested to give a more realistic representation of both the effect of mutations and the effect of drugs (see, e.g., [6062]), and was applied in computational studies of SQT1 in [48, 63, 64] and SQT3 in [65]. On the other hand, complex Markov models are associated with a significant increase in the number of free parameters to determine in the models, and we have therefore chosen to use a more simplified Hodgkin-Huxley formalism in this study.

The WT and SQT1 versions of the IKr current were inserted into the base models for hiPSC-CMs and adult CMs to define WT and SQT1 versions of these action potential models. The resulting computed action potentials displayed a clear action potential shortening for the SQT1 case for both hiPSC-CMs and adult CMs (see Fig 6). The computed pseudo-ECGs for the adult case also displayed a clear reduction in the duration of the QT interval for SQT1. Furthermore, the APD90 values of the hiPSC-CM model and the QT intervals of the adult model appeared to fit well with measured data for WT and SQT1 from [9, 15].

4.2 Computational identification of drug effects from membrane potential measurements of SQT1 hiPSC-CMs

Because of the possibly lethal risks associated with SQT syndrome, there is an urgent need for efficient therapies, including pharmacological treatment. In [9, 11], a number of drug candidates were tested on hiPSC-CMs generated from an SQT1 patient, and the drugs quinidine, ivabradine, ajmailine and mexiletine seemed to be the most promising for increasing the action potential duration of the SQT1 hiPSC-CMs. Other drugs shown to block WT IKr (sotalol, ranolazine, flecainide, amiodarone) did not seem to work as well for hiPSC-CMs affected by the SQT1 mutation. This highlights the need for mutation-specific investigations of drug effects. Here we have tried to extend such investigations computationally, and have used the measurements from [9, 11] to estimate the effect of quinidine, ivabradine, ajmaline and mexiletine for adult SQT1 patients using the computational procedure introduced in [22, 23].

4.2.1 Identifiability of model parameters

In [22, 23], drug effects were identified in an inversion procedure based on optical measurements of the membrane potential and cytosolic calcium concentration of hiPSC-CMs in a microphysiological system [66]. In the present study, we have used a different type of data in the inversion procedure, i.e., membrane potential measurements from patch-clamp recordings. This alternative data type could have both advantages and disadvantages compared to the optical measurements used in [22, 23]. A disadvantage is that we only have measurements of the membrane potential and not any information about the calcium transient. Information about the calcium transient has been shown to potentially improve the identification of parameters in mathematical action potential models [22, 67, 68]. However, an advantage of the patch-clamp recordings is that we can obtain information about the actual value of the membrane potential instead of just relative pixel intensities from optical measurements of voltage sensitive dyes. Therefore, we can get information about the resting membrane potential, the action potential amplitude and the upstroke velocity, which could all be useful for the identification of parameters. For example, the upstroke velocity is difficult to obtain from optical measurements of the action potential, which makes it very hard to identify drug effects on the fast sodium current, INa, unless additional data, e.g. extracellular measurements are included (see, e.g., [69]). Note also that identification of drug effects on INa from measurements of hiPSC-CMs could be hampered by the fact that the diastolic membrane potential of hiPSC-CMs often is depolarized compared to that of adult CMs (see, e.g., [21, 39]). This can lead to inactivation of the Na+ channels, making the upstroke of the action potential to a lager extent driven by ICaL as opposed to INa, like for pacemaker cells [21, 70]. In the current study, however, we consider measurements of hiPSC-CMs with a resting membrane potential of about −80 mV, indicating that INa is important for determining the upstroke velocity of these cells.

With information about the upstroke velocity, it should therefore be possible to identify INa. However, full identification of all model parameters based on membrane potential measurements is very challenging, and in some cases several combinations of parameter values result in virtually identical action potentials (see, e.g., [68, 71, 72]). In order to investigate which currents were identifiable in the SQT1 hiPSC-CM base model, we applied the singular value decomposition approach from [73]. The result of this analysis is given in S3 Fig of the Supporting Information. In short, the analysis reveals that the INa, ICaL, and IKr currents are the most identifiable currents of the model. This means that we expect the identifiability of the three main currents investigated for drug effects to be relatively high. In S4 Fig of the Supporting Information, we also show examples of the inversion method applied to simulated data. This enables us to assess the accuracy of our approach, and we observe that we are able to find the correct solution using the continuation method in examples where adjustment factors, λ, for these three currents and up to two additional currents are treated as free parameters. Similar experiments have shown similar results (see, e.g., Figure 6 in [23]) and although this does not prove that we are able to get the correct parameterization of the model in all cases, it indicates that the method is useful. A general proof of uniqueness of a minimum is most likely not possible unless strong regularization terms are added.

For the If-blocking drug ivabradine, we also estimate drug effects on If, but we note that the identifiability of this current is not very high according to the singular value decomposition analysis and that the estimated drug effect for If consequently is associated with a large degree of uncertainty. The λ-values found for the remaining parameters in the inversion procedure are also associated with a degree of uncertainty. However, these parameters are not used directly in the prediction of drug effects in the adult case (see Section 2.1).

It should also be noted that the hand-tuned parameterizations chosen as a starting point for the inversions of hiPSC-CM data and as the default adult base models are associated with uncertainty. In the hiPSC-CM case, some of the current densities are adjusted to measurements of hiPSC-CMs from [9] (see S1 Fig), but for the remaining parameter values, it is possible that other combinations of parameters would give equally good or better fits to the data. Similarly, for the adult case, we have fitted the model to information about the heart rate dependent QTp interval from [15]. In addition, we have attempted to make the size of some of the main currents similar to the size of the currents in the models [26, 27] (see S2 Fig). However, it is possible that the chosen parameterization is not unique or that better parameterizations exist.

4.2.2 Identified drug effects

Using the computational inversion procedure described in Section 2.3, we estimated drug effects on major ion channels based on action potential measurements of SQT1 hiPSC-CMs. For quinidine, ajmaline and mexiletine, we considered drug effects on the currents IKr, ICaL and INa, and for ivabradine, we considered drug effects on the currents IKr, If and INa, because these currents have been shown to be affected by ivabradine [7476].

Each of the considered drugs had a considerable effect on the IKr current, with predicted IC50 values of 8.1 μM for quinidine, 12.6 μM for ivabradine, 69 μM for ajmaline and 280 μM for mexiletine (see Fig 8). The effect of a drug can be very different for WT IKr channels compared to that of mutated IKr channels (see, e.g., [77, 78]). For instance, in [77], the IC50 value of quinidine was found to be 0.6 μM for WT IKr channels and 2.2 μM for the SQT1 mutation N588K. Consequently, it is not reasonable to compare the IC50 values identified for the SQT1 case to IC50 values found in WT IKr experiments. However, we observe that the identified IC50 value for quinidine is quite close to the IC50 value reported for the SQT1 case in [77] (8.1 μM vs 2.2 μM). Furthermore, the predicted block percentage for 10 μM of quinidine is very close to the measured block percentage reported in [11] (see Fig 9). The block percentages for the remaining drugs are also quite close to the percentages reported in [11], indicating reasonable identified IC50 values for SQT1 IKr. For ivabradine, the identified IC50 value for IKr (12.6 μM) is also in good agreement with the IC50 value of 10.3 μM reported for IKr affected by the SQT1 (N588K) mutation in [74].

In addition, for ajmaline, the inversion procedure identified an IC50 value of 47 μM for ICaL. This is comparable to what was measured in [79] (70.8 μM). The inversion procedure also identified a significant block of INa for the drug mexiletine (see Fig 8). This is reasonable, since mexiletine is known as an INa blocker [80]. However, the identified IC50 value of 200 μM is larger than IC50 values found in literature (7.6 μM–43 μM) [8183]. Similarly, ivabradine is an If blocker with an IC50 value found to be about 2 μM in [75], but the computational procedure underestimates the drug effect to an IC50 value of about 42 μM. This discrepancy is not unexpected, since the identifiability of If is predicted to be low (see S3 Fig). The IC50 value for INa block of ivabradine, was estimated to 86 μM, closer to, but still a bit higher than the measured IC50 value of 30 μM reported in [76]. Possible effects of ajmaline on INa, of quinidine on ICaL and INa and of mexiletine on ICaL are also underestimated by the inversion procedure compared to the IC50 values reported in [83]. This could indicate that the biomarkers considered in the inversion procedure (see Figs 3 and 7) are not enough to properly characterize these less pronounced drug effects or that the applied base model formulation or simplified drug modeling are not sufficient to represent these effects. It also could be that the regularization term used for the ε values can cause underestimation due to its pressure on finding the minimum effect in the chosen fit. It should also be noted that IC50 values identified in different experiments tend to vary considerably.

4.3 Computational maturation of drug response

Using the single channel drug effects identified from membrane potential measurements of hiPSC-CM, we estimated the drug effects in the adult case based on the three assumptions explained in Fig 2. First, we assumed that the density of different types of ion channels (and other membrane and intracellular proteins) differ between hiPSC-CMs and adult CMs, but that the function of a single channel is the same for adult CMs and hiPSC-CMs. This means that a set of parameters representing the density of different types of proteins, e.g., ion channel conductances (see (4)), are different between the hiPSC-CM and adult CM models, but that the models for, e.g., the open probability of the channels are the same in both versions of the model. Second, we assumed that a mutation affects the function of an ion channel in exactly the same manner for hiPSC-CMs and adult CMs. Therefore, we used the same adjustment of the open probability of the IKr current (see Fig 5) to represent the SQT1 mutation N588K for both hiPSC-CMs and adult CMs (see Fig 6). Third, we assumed that a drug affects a single ion channel in exactly the same manner for hiPSC-CMs and adult CMs, and we could therefore estimate the drug effect in the adult case based on the drug effect identified for hiPSC-CMs.

This modeling framework was used in [22, 23] to predict side effects of drugs for adult WT CMs based on measurements of WT hiPSC-CMs. However, the procedure can also take the effect of mutations into account (as specified in the second assumption of Fig 2), and in the present study, we have used the procedure to estimate drug effects on the ventricular action potential and the QT interval of adult SQT1 patients based on measurements of hiPSC-CMs generated from an SQT1 patient.

4.3.1 Estimation of the QT interval

Because the effect of one of the drugs considered in this study has been investigated for adult SQT1 patients in the form of measured QT intervals of patients’ ECGs, we wished to be able to compare the drug responses predicted by our computational procedure to these measured QT intervals. However, the mathematical base model from [23] only represents the adult ventricular action potential and cannot be directly used to compute the adult QT interval. On the other hand, it is generally recognized that there is a clear link between the duration of the ventricular action potential and the QT interval of a patient’s ECG. For example, in [84], a simple linear relationship between the APD90 value and the QT interval was observed. Conversely, in the computational study [85], it was shown that there was not a direct linear realtionship between drug effects on the action potential duration and drug effects on the QT interval. For example, drugs blocking the INa current gave rise to prolongations of the computed QT intervals, even though no prolongation was observed for the action potential duration. In this study, we have chosen to predict the QT interval by computing a pseudo-ECG using a very simple mathematical representation, as applied in several previous studies of drug effects on the QT interval (e.g., [46, 49, 50]). The applied method appeared to give rise to reasonable pseudo-ECG waveforms and QT intervals (see, e.g., Fig 6), but it is a very simplified representation of the dynamics underlying the ECG, and more realistic spatial modeling, possibly extended to whole-heart simulations including a surrounding torso like in [8588], could potentially improve the reliability of the computed QT intervals.

4.3.2 Predicted adult drug responses for SQT1 patients

The estimated effects of the four drugs quinidine, ivabradine, ajmaline and mexiletine on the ventricular action potential and QT interval of adult SQT1 patients revealed that all four drugs are expected to lead to an increase in the action potential duration and the QT interval of SQT1 patients (see Fig 10). However, the extent of the prolongation of the QT interval was estimated to be quite different for the different drugs. The most significant prolongation was observed for 10 μM of quinidine, which was predicted to increase the QTp interval from 200 ms to 257 ms, closer to the WT value of 292 ms. For the drug ivabradine, the procedure also estimated a significant increase in the adult QT interval, and for the drug mexilitine, a smaller, but clearly visible QT interval increase was predicted. For ajmaline, however, the estimated drug effect in the adult case was quite small, even though a significant APD prolongation was observed in the hiPSC-CM case (compare Figs 8 and 10). This is probably explained by the fact that the inversion procedure estimates a significant block of both IKr and ICaL for ajmaline which due to the difference in protein densities between the hiPSC-CM and adult cases cancel each other out in the adult model, but not in the hiPSC-CM model.

It is generally hard to validate the drug responses predicted by the computational procedure because data of drug responses for adult SQT1 patients would be required. However, for the drug quinidine, such experiments are reported in [15] in the form of measured QT intervals for a number of different heart rates. These measurements have been conducted for SQT1 patients without pharmacological treatment, SQT1 patients using the drug quinidine and for healthy adults without the SQT1 mutation [15]. We found that the effect of quinidine for adult patients estimated by the computational procedure fitted well with these measured drug responses (see Fig 11), indicating that the applied computational procedure shows promise for reliably predicting adult drug responses based on measurements of hiPSC-CMs.

4.4 Conclusions

In this paper, we have used a computational procedure to predict the effect of the drugs quinidine, ivabradine, ajmaline and mexilitine for clinical biomarkers of adult SQT1 patients, based on measured drug responses for hiPSC-CMs generated from an SQT1 patient. The procedure consists of identifying drug effects on individual ion channels based on membrane potential measurements of drug effects for hiPSC-CMs and then mapping this effect up to the adult case based on the assumption of functional invariance of ion channels during maturation. All four drugs were predicted to lead to increases in the ventricular action potential duration and QT interval, and the drugs quinidine and ivabradine were estimated to be most effective. Furthermore, the effect of quinidine predicted by the computational procedure appeared to be in good agreement with measured effects of the drug. Consequently, we conclude that the computational procedure applied in this study could be a useful tool for estimating adult drug responses based on measurements of hiPSC-CMs, also as new measurements become available, including measurements of alternative drugs or mutations.

References

IGussak, PBrugada, JBrugada, RSWright, SLKopecky, BRChaitman, et al. Idiopathic short QT interval: a new clinical syndrome? Cardiology. 2000;94(2):99102. 10.1159/000047299

OCampuzano, GSarquella-Brugada, SCesar, EArbelo, JBrugada, RBrugada. Recent advances in short QT syndrome. Frontiers in Cardiovascular Medicine. 2018;5:149. 10.3389/fcvm.2018.00149

CPatel, GXYan, CAntzelevitch. Short QT syndrome: from bench to bedside. Circulation: Arrhythmia and Electrophysiology. 2010;3(4):401408. 10.1161/CIRCEP.109.921056

FGaita, CGiustetto, FBianchi, CWolpert, RSchimpf, RRiccardi, et al. Short QT syndrome: a familial cause of sudden death. Circulation. 2003;108(8):965970. 10.1161/01.CIR.0000085071.28695.C4

RSchimpf, CWolpert, FGaita, CGiustetto, MBorggrefe. Short QT syndrome. Cardiovascular Research. 2005;67(3):357366. 10.1016/j.cardiores.2005.03.026

RBrugada, KHong, RDumaine, JCordeiro, FGaita, MBorggrefe, et al. Sudden death associated with short-QT syndrome linked to mutations in HERG. Circulation. 2004;109(1):3035. 10.1161/01.CIR.0000109482.92774.3A

HAbriel, JSRougier. β-Blockers in Congenital Short-QT Syndrome as Ion Channel Blockers. Journal of Cardiovascular Electrophysiology. 2013;24(10):11721174. 10.1111/jce.12204

AMazzanti, RMaragna, GVacanti, AKostopoulou, MMarino, NMonteforte, et al. Hydroquinidine prevents life-threatening arrhythmic events in patients with short QT syndrome. Journal of the American College of Cardiology. 2017;70(24):30103015. 10.1016/j.jacc.2017.10.025

IEl-Battrawy, HLan, LCyganek, ZZhao, XLi, FBuljubasic, et al. Modeling Short QT Syndrome Using Human-Induced Pluripotent Stem Cell–Derived Cardiomyocytes. Journal of the American Heart Association. 2018;7(7):e007394. 10.1161/JAHA.117.007394

10 

IEl-Battrawy, JBesler, VLiebe, RSchimpf, ETülümen, BRudic, et al. Long-Term Follow-Up of Patients With Short QT Syndrome: Clinical Profile and Outcome. Journal of the American Heart Association. 2018;7(23):e010073. 10.1161/JAHA.118.010073

11 

ZZhao, XLi, IEl-Battrawy, HLan, RZhong, QXu, et al. Drug Testing in Human-Induced Pluripotent Stem Cell–Derived Cardiomyocytes From a Patient With Short QT Syndrome Type 1. Clinical Pharmacology & Therapeutics. 2019;106(3):642651. 10.1002/cpt.1449

12 

FGaita, CGiustetto, FBianchi, RSchimpf, MHaissaguerre, LCalò, et al. Short QT syndrome: pharmacological treatment. Journal of the American College of Cardiology. 2004;43(8):14941499. 10.1016/j.jacc.2004.02.034

13 

JCordeiro, RBrugada, YWu, KHong, RDumaine. Modulation of IKr inactivation by mutation N588K in KCNH2: a link to arrhythmogenesis in short QT syndrome. Cardiovascular Research. 2005;67(3):498509. 10.1016/j.cardiores.2005.02.018

14 

MJMcPate, RSDuncan, JTMilnes, HJWitchel, JCHancox. The N588K-HERG K+ channel mutation in the’short QT syndrome’: mechanism of gain-in-function determined at 37°C. Biochemical and Biophysical Research Communications. 2005;334(2):441449. 10.1016/j.bbrc.2005.06.112

15 

CWolpert, RSchimpf, CGiustetto, CAntzelevitch, JCordeiro, RDumaine, et al. Further insights into the effect of quinidine in short QT syndrome caused by a mutation in HERG. Journal of Cardiovascular Electrophysiology. 2005;16(1):5458. 10.1046/j.1540-8167.2005.04470.x

16 

KTakahashi, SYamanaka. Induction of pluripotent stem cells from mouse embryonic and adult fibroblast cultures by defined factors. Cell. 2006;126(4):663676. 10.1016/j.cell.2006.07.024

17 

YYoshida, SYamanaka. Induced Pluripotent Stem Cells 10 Years Later. Circulation Research. 2017;120(12):19581968. 10.1161/CIRCRESAHA.117.311080

18 

ADi Baldassarre, ECimetta, SBollini, GGaggi, BGhinassi. Human-Induced Pluripotent Stem Cell Technology and Cardiomyocyte Generation: Progress and Clinical Applications. Cells. 2018;7(6):48. 10.3390/cells7060048

19 

LYe, XNi, ZAZhao, WLei, SHu. The Application of Induced Pluripotent Stem Cells in Cardiac Disease Modeling and Drug Testing. Journal of Cardiovascular Translational Research. 2018;11(5):366374. 10.1007/s12265-018-9811-3

20 

MEHartman, DFDai, MALaflamme. Human pluripotent stem cells: Prospects and challenges as a source of cardiomyocytes for in vitro modeling and cell-based cardiac repair. Advanced Drug Delivery Reviews. 2016;96:317. 10.1016/j.addr.2015.05.004

21 

PGarg, VGarg, RShrestha, MCSanguinetti, TJKamp, JCWu. Human Induced Pluripotent Stem Cell–Derived Cardiomyocytes as Models for Cardiac Channelopathies. Circulation Research. 2018;123(2):224243. 10.1161/CIRCRESAHA.118.311209

22 

ATveito, KHJæger, NHuebsch, BCharrez, AGEdwards, SWall, et al. Inversion and computational maturation of drug response using human stem cell derived cardiomyocytes in microphysiological systems. Scientific Reports. 2018;8(1):17626. 10.1038/s41598-018-35858-7

23 

KHJæger, VCharwat, BCharrez, HFinsberg, MMMaleckar, SWall, et al. Improved computational identification of drug response using optical measurements of human stem cell derived cardiomyocytes in microphysiological systems. Frontiers in Pharmacology. 2020;10:1648. 10.3389/fphar.2019.01648

24 

JQGong, EASobie. Population-based mechanistic modeling allows for quantitative predictions of drug responses across cell types. NPJ Systems Biology and Applications. 2018;4(1):11. 10.1038/s41540-018-0047-2

25 

MPaci, JHyttinen, BRodriguez, SSeveri. Human induced pluripotent stem cell-derived versus adult cardiomyocytes: an in silico electrophysiological study on effects of ionic current block. British Journal of Pharmacology. 2015;172(21):51475160. 10.1111/bph.13282

26 

TO’Hara, LVirág, AVarró, YRudy. Simulation of the undiseased human cardiac ventricular action potential: Model formulation and experimental validation. PLoS Computational Biology. 2011;7(5):e1002061. 10.1371/journal.pcbi.1002061

27 

EGrandi, FSPasqualini, DMBers. A novel computational model of the human ventricular action potential and Ca transient. Journal of Molecular and Cellular Cardiology. 2010;48(1):112121. 10.1016/j.yjmcc.2009.09.019

28 

MMMaleckar, JLGreenstein, NATrayanova, WRGiles. Mathematical simulations of ligand-gated and cell-type specific effects on the action potential of human atrium. Progress in Biophysics and Molecular Biology. 2008;98(2-3):161170. 10.1016/j.pbiomolbio.2009.01.010

29 

MCourtemanche, RJRamirez, SNattel. Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model. American Journal of Physiology-Heart and Circulatory Physiology. 1998;275(1):H301H321. 10.1152/ajpheart.1998.275.1.H301

30 

MPaci, JHyttinen, KAalto-Setälä, SSeveri. Computational models of ventricular-and atrial-like human induced pluripotent stem cell derived cardiomyocytes. Annals of Biomedical Engineering. 2013;41(11):23342348. 10.1007/s10439-013-0833-3

31 

DCKernik, SMorotti, HWu, PGarg, HJDuff, JKurokawa, et al. A computational model of induced pluripotent stem-cell derived cardiomyocytes incorporating experimental variability from multiple data sources. The Journal of Physiology. 2019;. 10.1113/JP277724

32 

SVecchietti, EGrandi, SSeveri, IRivolta, CNapolitano, SGPriori, et al. In silico assessment of Y1795C and Y1795H SCN5A mutations: implication for inherited arrhythmogenic syndromes. American Journal of Physiology-Heart and Circulatory Physiology. 2007;292(1):H56H65. 10.1152/ajpheart.00270.2006

33 

MPaci, EPassini, SSeveri, JHyttinen, BRodriguez. Phenotypic variability in LQT3 human induced pluripotent stem cell-derived cardiomyocytes and their response to anti-arrhythmic pharmacological therapy: an in silico approach. Heart Rhythm. 2017;. 10.1016/j.hrthm.2017.07.026

34 

LTung. A bi-domain model for describing ischemic myocardial dc potentials. Massachusetts Institute of Technology; 1978.

35 

RMShaw, YRudy. Ionic mechanisms of propagation in cardiac tissue: roles of the sodium and L-type calcium currents during reduced excitability and decreased gap junction coupling. Circulation Research. 1997;81(5):727741. 10.1161/01.RES.81.5.727

36 

NATrayanova. Whole-heart modeling: Applications to cardiac electrophysiology and electromechanics. Circulation Research. 2011;108(1):11328. 10.1161/CIRCRESAHA.110.223610

37 

ATveito, KHJæger, MKuchta, KAMardal, MERognes. A cell-based framework for numerical modeling of electrical conduction in cardiac tissue. Frontiers in Physics. 2017;5:48. 10.3389/fphy.2017.00048

38 

KHJæger, AGEdwards, AMcCulloch, ATveito. Properties of cardiac conduction in a cell-based computational model. PLoS Computational Biology. 2019;15(5):e1007042. 10.1371/journal.pcbi.1007042

39 

MHoekstra, CLMummery, AWilde, CRBezzina, AOVerkerk. Induced pluripotent stem cell derived cardiomyocytes as models for cardiac arrhythmias. Frontiers in Physiology. 2012;3:346. 10.3389/fphys.2012.00346

40 

JKGibson, YYue, JBronson, CPalmer, RNumann. Human stem cell-derived cardiomyocytes detect drug-mediated changes in action potentials and ion currents. Journal of Pharmacological and Toxicological Methods. 2014;70(3):255267. 10.1016/j.vascn.2014.09.005

41 

JCHancox, DGWhittaker, CDu, AGStuart, HZhang. Emerging therapeutic targets in the short QT syndrome. Expert Opinion on Therapeutic Targets. 2018;22(5):439451. 10.1080/14728222.2018.1470621

42 

Paci M, Passini E, Klimas A, Severi S, Hyttinen J, Rodriguez B, et al. In silico populations optimized on optogenetic recordings predict drug effects in human induced pluripotent stem cell-derived cardiomyocytes. In: 2018 Computing in Cardiology Conference (CinC). vol. 45. IEEE; 2018. p. 1–4.

43 

TRShannon, FWang, JPuglisi, CWeber, DMBers. A mathematical treatment of integrated Ca dynamics within the ventricular myocyte. Biophysical Journal. 2004;87(5):33513371. 10.1529/biophysj.104.047449

44 

JANelder, RMead. A simplex method for function minimization. The Computer Journal. 1965;7(4):308313. 10.1093/comjnl/7.4.308

45 

MATLAB Global Optimization Toolbox; Version 4.3.

46 

KGima, YRudy. Ionic current basis of electrocardiographic waveforms: a model study. Circulation Research. 2002;90(8):889896. 10.1161/01.RES.0000016960.61087.86

47 

EPueyo, ZHusti, THornyik, IBaczkó, PLaguna, AVarró, et al. Mechanisms of ventricular rate adaptation as a predictor of arrhythmic risk. American Journal of Physiology-Heart and Circulatory Physiology. 2010;298(5):H1577H1587. 10.1152/ajpheart.00936.2009

48 

IAdeniran, MJMcPate, HJWitchel, JCHancox, HZhang. Increased vulnerability of human ventricle to re-entrant excitation in hERG-linked variant 1 short QT syndrome. PLoS Computational Biology. 2011;7(12). 10.1371/journal.pcbi.1002313

49 

CLuo, KWang, HZhang. Effects of amiodarone on short QT syndrome variant 3 in human ventricles: a simulation study. Biomedical Engineering Online. 2017;16(1):69. 10.1186/s12938-017-0369-0

50 

BTrenor, JGomis-Tena, KCardona, LRomero, SRajamani, LBelardinelli, et al. In silico assessment of drug safety in human heart applied to late sodium current blockers. Channels. 2013;7(4):249262. 10.4161/chan.24905

51 

ATveito, KHJæger, GTLines, ŁPaszkowski, JSundnes, AGEdwards, et al. An evaluation of the accuracy of classical models for computing the membrane potential and extracellular potential for neurons. Frontiers in Computational Neuroscience. 2017;11:27. 10.3389/fncom.2017.00027

52 

RWeingart. The actions of ouabain on intercellular coupling and conduction velocity in mammalian ventricular muscle. The Journal of Physiology. 1977;264(2):341365. 10.1113/jphysiol.1977.sp011672

53 

National Center for Biotechnology Information. PubChem Database. Quinidine, CID = 441074;. https://pubchem.ncbi.nlm.nih.gov/compound/Quinidine.

54 

IEl-Battrawy, JBesler, XLi, HLan, ZZhao, VLiebe, et al. Impact of antiarrhythmic drugs on the outcome of Short QT Syndrome. Frontiers in Pharmacology. 2019;10:771. 10.3389/fphar.2019.00771

55 

DLWeiss, GSeemann, FBSachse, ODössel. Modelling of short QT syndrome in a heterogeneous model of the human ventricular wall. EP Europace. 2005;7(s2):S105S117. 10.1016/j.eupc.2005.04.008

56 

Bartolucci C, Moreno C, de la Cruz A, Lambiase P, Severi S, Valenzuela C. Linking a novel mutation to its short QT phenotype through multiscale computational modelling. In: Computing in Cardiology 2014. IEEE; 2014. p. 1017–1020.

57 

KHten Tusscher, AVPanfilov. Cell model for efficient simulation of wave propagation in human ventricular tissue under normal and pathological conditions. Physics in Medicine and Biology. 2006;51(23):6141. 10.1088/0031-9155/51/23/014

58 

YRudy, JRSilva. Computational biology in the study of cardiac ion channels and cell electrophysiology. Quarterly Reviews of Biophysics. 2006;39(1):57116. 10.1017/S0033583506004227

59 

MFink, DNoble. Markov models for ion channels: versatility versus identifiability and speed. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences. 2009;367(1896):21612179. 10.1098/rsta.2008.0301

60 

CEClancy, ZIZhu, YRudy. Pharmacogenetics and anti-arrhythmic drug therapy: A theoretical investigation. AJP: Heart and Circulatory Physiology. 2007;292(1):H66H75. 10.1152/ajpheart.00312.2006

61 

ATveito, GTLines, PLi, AMcCulloch. Defining candidate drug characteristics for Long-QT (LQT3) syndrome. Mathematical Biosciences and Engineering. 2011;8(3):86173. 10.3934/mbe.2011.8.861

62 

Tveito A, Lines GT. Computing Characterizations of Drugs for Ion Channels and Receptors Using Markov Models. Springer-Verlag, Lecture Notes, vol. 111; 2016.

63 

DGWhittaker, HNi, APBenson, JCHancox, HZhang. Computational analysis of the mode of action of disopyramide and quinidine on hERG-linked short QT syndrome in human ventricles. Frontiers in Physiology. 2017;8:759. 10.3389/fphys.2017.00759

64 

Luo C, Wang K, Liu T, Zhang H. Computational Analysis of the Action of Chloroquine on Short QT Syndrome Variant 1 and Variant 3 in Human Ventricles. In: 2018 40th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC). IEEE; 2018:5462–5465.

65 

IAdeniran, DGWhittaker, AEl Harchi, JCHancox, HZhang. In silico investigation of a KCNQ1 mutation associated with short QT syndrome. Scientific Reports. 2017;7(1):114. 10.1038/s41598-017-08367-2

66 

AMathur, PLoskill, KShao, NHuebsch, SHong, SGMarcus, et al. Human iPSC-based cardiac microphysiological system for drug screening applications. Scientific Reports. 2015;5:8883. 10.1038/srep08883

67 

CMRees, JHYang, MSantolini, AJLusis, JNWeiss, AKarma. The Ca2+ transient as a feedback sensor controlling cardiomyocyte ionic conductances in mouse populations. eLife. 2018;7:e36717. 10.7554/eLife.36717

68 

AXSarkar, EASobie. Regression analysis for constraining free parameters in electrophysiological models of cardiac cells. PLoS Computational Biology. 2010;6(9). 10.1371/journal.pcbi.1000914

69 

KHJæger, VCharwat, SWall, KHealy, ATveito. Identifying drug response by combining measurements of the membrane potential, the cytosolic calcium concentration, and the extracellular potential in microphysiological systems. Frontiers in Pharmacology. 2021; 11:569489. 10.3389/fphar.2020.569489

70 

QWen-ning, LShao-ru. Effect of tetrodotoxin on the transmembrane potential of the pacemaker cell in the sinus venosus of toad. Journal of Tongji Medical University. 1988;8(2):101105. 10.1007/BF02887804

71 

AXSarkar, DJChristini, EASobie. Exploiting mathematical models to illuminate electrophysiological variability between individuals. The Journal of Physiology. 2012;590(11):25552567. 10.1113/jphysiol.2011.223313

72 

SDokos, NHLovell. Parameter estimation in cardiac ionic models. Progress in Biophysics and Molecular Biology. 2004;85(2-3):407431. 10.1016/j.pbiomolbio.2004.02.002

73 

KHJæger, SWall, ATveito. Detecting undetectables: Can conductances of action potential models be changed without appreciable change in the transmembrane potential? Chaos: An Interdisciplinary Journal of Nonlinear Science. 2019;29(7):073102. 10.1063/1.5087629

74 

DMelgari, KEBrack, CZhang, YZhang, AEl Harchi, JSMitcheson, et al. hERG potassium channel blockade by the HCN channel inhibitor bradycardic agent ivabradine. Journal of the American Heart Association. 2015;4(4):e001813. 10.1161/JAHA.115.001813

75 

ABucchi, ATognati, RMilanesi, MBaruscotti, DDiFrancesco. Properties of ivabradine-induced block of HCN1 and HCN4 pacemaker channels. The Journal of Physiology. 2006;572(2):335346. 10.1113/jphysiol.2005.100776

76 

NHaechl, JEbner, KHilber, HTodt, XKoenig. Pharmacological Profile of the Bradycardic Agent Ivabradine on Human Cardiac Ion Channels. Cellular Physiology and Biochemistry. 2019;53:3648. 10.33594/000000119

77 

MMcPate, RDuncan, JHancox, HWitchel. Pharmacology of the short QT syndrome N588K-hERG K+ channel mutation: differential impact on selected class I and class III antiarrhythmic drugs. British Journal of Pharmacology. 2008;155(6):957966. 10.1038/bjp.2008.325

78 

DMelgari, AEl Harchi, CEDempsey, JCHancox. Sensitivity of Flecainide Inhibition of hERG Channels to Channel Inactivation. Biophysical Journal. 2014;106(2):138a. 10.1016/j.bpj.2013.11.807

79 

MBébarová, PMatejovic, MPásek, MSimurdova, JSimurda. Effect of ajmaline on action potential and ionic currents in rat ventricular myocytes. General Physiology and Biophysics. 2005;24(3):311.

80 

KRCourney. Comparatibe actions of mexiletine on sodium channels in nerve, skeletal and cardiac muscle. European Journal of Pharmacology. 1981;74(1):918. 10.1016/0014-2999(81)90317-4

81 

YQu, HMVargas. Proarrhythmia risk assessment in human induced pluripotent stem cell-derived cardiomyocytes using the maestro MEA platform. Toxicological Sciences. 2015;147(1):286295. 10.1093/toxsci/kfv128

82 

XDHu, JQQian. DDPH inhibited L-type calcium current and sodium current in single ventricular myocyte of guinea pig. Acta pharmacologica Sinica. 2001;22(5):415419.

83 

GRMirams, YCui, ASher, MFink, JCooper, BMHeath, et al. Simulation of multiple ion channel block provides improved early prediction of compounds’ clinical torsadogenic risk. Cardiovascular Research. 2011;91(1):5361. 10.1093/cvr/cvr044

84 

TMing, ZGuang-lan, ZXiao-rong, ZJun-yuan, ZYi, LShao-ru. Relationship between action potential duration of ventricular cells and heart rate in dog under natural breathing. Journal of Tongji Medical University. 1987;7(3):148152. 10.1007/BF02888208

85 

NZemzemi, MOBernabeu, JSaiz, JCooper, PPathmanathan, GRMirams, et al. Computational assessment of drug-induced effects on the electrocardiogram: from ion channel to body surface potentials. British Journal of Pharmacology. 2013;168(3):718733. 10.1111/j.1476-5381.2012.02200.x

86 

GTLines, PGrottum, ATveito. Modeling the electrical activity of the heart: a bidomain model of the ventricles embedded in a torso. Computing and Visualization in Science. 2003;5(4):195213. 10.1007/s00791-003-0100-5

87 

MPotse, BDubé, AVinet. Cardiac anisotropy in boundary-element models for the electrocardiogram. Medical & Biological Engineering & Computing. 2009;47(7):719729. 10.1007/s11517-009-0472-x

88 

MBoulakia, SCazeau, MAFernández, JFGerbeau, NZemzemi. Mathematical modeling of electrocardiograms: a numerical study. Annals of Biomedical Engineering. 2010;38(3):10711097. 10.1007/s10439-009-9873-0

89 

SZhang. Isolation and characterization of IKr in cardiac myocytes by Cs+ permeation. American Journal of Physiology-Heart and Circulatory Physiology. 2006;290(3):H1038H1049. 10.1152/ajpheart.00679.2005

90 

TO’Hara, LVirág, AVarró, YRudy. Simulation of the undiseased human cardiac ventricular action potential: Model formulation and experimental validation. PLoS Computational Biology. 2011;7(5):e1002061. 10.1371/journal.pcbi.1002061