An in silico hiPSC-Derived Cardiomyocyte Model Built With Genetic Algorithm

The formulation of in silico biophysical models generally requires optimization strategies for reproducing experimentally observed phenomena. In electrophysiological modeling, robust nonlinear regressive methods are often crucial for guaranteeing high fidelity models. Human induced pluripotent stem cell-derived cardiomyocytes (hiPSC-CMs), though nascent, have proven to be useful in cardiac safety pharmacology, regenerative medicine, and in the implementation of patient-specific test benches for investigating inherited cardiac disorders. This study demonstrates the potency of heuristic techniques at formulating biophysical models, with emphasis on a hiPSC-CM model using a novel genetic algorithm (GA) recipe we proposed. The proposed GA protocol was used to develop a hiPSC-CM biophysical computer model by fitting mathematical formulations to experimental data for five ionic currents recorded in hiPSC-CMs. The maximum conductances of the remaining ionic channels were scaled based on recommendations from literature to accurately reproduce the experimentally observed hiPSC-CM action potential (AP) metrics. Near-optimal parameter fitting was achieved for the GA-fitted ionic currents. The resulting model recapitulated experimental AP parameters such as AP durations (APD₅₀, APD₇₅, and APD₉₀), maximum diastolic potential, and frequency of automaticity. The outcome of this work has implications for validating the biophysics of hiPSC-CMs in their use as viable substitutes for human cardiomyocytes, particularly in cardiac safety pharmacology and in the study of inherited cardiac disorders. This study presents a novel GA protocol useful for formulating robust numerical biophysical models. The proposed protocol is used to develop a hiPSC-CM model with implications for cardiac safety pharmacology.

Introduction

High fidelity numerical biophysical models have the potential to provide the missing mechanistic link between the experimental observations and their clinical implications. Formulating such models are primarily optimization problems with the experimental data as targets. Metaheuristic algorithms, such as genetic algorithms (GAs; Smith et al., 1995), are well suited for this task due to their inherent stochastic and judicious exploration of the solution space. In electrophysiological studies, single- and multi-cell mathematical models have aided efforts at elucidating various biological processes, such as perception, cognition, and cardiac function, which is the focus of this study.

Human induced pluripotent stem cells (hiPSCs), since their discovery by Takahashi et al. (2007), have been instrumental at the development of ethically sound patient-specific human cell and tissue models. These models have, in turn, allowed scientists to investigate the underpinnings of congenital and drug-induced disorders. A prominent derivative of this cell type is the human induced pluripotent stem cell-derived cardiomyocyte (hiPSC-CM). hiPSC-CMs have recently proven to be both useful and promising in cardiac safety pharmacology, where these cells are adopted as test benches for studying drug effects on cardiac function (Gintant et al., 2019). hiPSC-CMs have also been used to better understand drug-induced and inherited cardiac disorders such as long QT syndrome (Moretti et al., 2010; Egashira et al., 2012) and catecholaminergic polymorphic ventricular tachycardia (Fatima et al., 2011; Jung et al., 2012). Prior to the advent of hiPSC-CMs, researchers often relied on animal models to make extrapolations of disease effects in human cells. This route possesses a significant caveat of dissimilar genotypic representation by these models. While freshly excised human cardiac cells are the ideal candidates for inductive human cardiac studies, obtaining these uncommon cells is highly invasive relative to hiPSC-CMs, which can even be derived from superficial somatic cells such as skin cells. Although hiPSC-CM is apparently a convenient choice for extensive long- and short-term cardiac studies, there are some deficiencies in the electrophysiological properties. There are ongoing attempts by scientists to better understand these cells and validate them as viable substitutes for human cardiomyocytes. These attempts are mainly in vitro and in silico. For instance, the Comprehensive in vitro Proarrhythmia Assay initiative, first presented in 2013, stipulates research directions around the use of hiPSC-CM experimental data and mathematical modeling in proarrhythmic risk assessment (Gintant et al., 2017).

Since the pioneering work by Hodgkin and Huxley (1952) involving the formulation of a biophysical mathematical model of the squid giant axon, there have been numerous attempts at creating biophysical models for different cell types (Di Francesco and Noble, 1985; Luo and Rudy, 1994; Courtemanche et al., 1998; Kurata et al., 2002). These models have been helpful at elucidating the electrical dynamics of native cardiomyocytes. Recently, there have been few attempts at formulating in silico hiPSC-CM models (Paci et al., 2013, 2018; Koivumäki et al., 2018; Kernik et al., 2019). However, the inherent variability in these cells pose challenges to deriving formulations that could reproduce the wide range of behaviors of hiPSC-CMs. One recent approach is to combine experimental data from multiple labs to derive model formulations on averaged data (Kernik et al., 2019). However, care must be taken while combining data from disparate sources because the inconsistencies in cell origins (Hwang et al., 2015), culturing environments, experimental protocols and conditions, as well as cell maturation levels (Narsinh et al., 2011) may introduce unwanted deviations in the model. This study joins in the efforts at formulating a robust in silico hiPSC-CM model based mostly on the experimental data from a single lab for maintaining phenotypical consistency. In achieving this, we present a novel customizable GA protocol employed for a near-optimal fitting of model ionic current formulations to the experimental data.

Genetic algorithm is a heuristic optimization method inspired by Darwinian evolution and natural selection (Smith et al., 1995). The process, similar to the biological counterpart, involves population (i.e., sets of possible solutions termed chromosomes) initialization to crossover, which involves a combinatorial shuffling of parameters (termed genes) about a specified number of pivot points (termed crossover points) between two parent solutions to generate offspring. Offspring solutions may then undergo mutation, typically done in an additive or multiplicative fashion. A fitness criterion must be defined to facilitate the propagation of near-optimal solutions over the specified number of iterations termed generations. Metaheuristic optimization methods offer resilience against local optima and saddle point convergence relative to gradient-based nonlinear optimization methods, such as the Newton–Raphson and Levenberg–Marquardt methods (Rhinehart, 2016). GA-based parametrization has been used to fine-tune mathematical models of murine myocytes (Bot et al., 2012; Groenendaal et al., 2015) and canine atrial cells (Syed et al., 2005) to incorporate cell-specific experimental data.

Our previous work (Akwaboah et al., 2020) demonstrated the feasibility of using GA to fit cardiac cell biophysical model formulations. In this study, we demonstrate how the GA-based ionic current fitting outcomes can be further integrated into the biophysical numerical model development process. An improved GA-based parametrization protocol was developed to build a complete in silico biophysical model of a hiPSC-CM.

Materials and Methods

Preparations of hiPSC-CMs

Human induced pluripotent stem cells usage was approved by the Institutional Stem Cell Research Oversight committee at Masonic Medical Research Institute. hiPSCs (from WiCell) were maintained on growth factor-reduced Matrigel coated plates in E8 medium with E8 supplement. Cardiac differentiation was induced by 6 μM CHIR99021, 10 ng/ml Activin A, and 5 nM rapamycin in RPMI1640 medium containing B27 (minus insulin) and 50 μg/ml L-ascorbic acid (cell culture tested powder) as the basal medium. After 24 h, media was changed in the same basal medium. The following day, cells were kept in RPMI1640/B27 (-insulin) with the addition of 5 μM XAV939 and 10 μM KY0211 for days 2–6, changing the media every other day. Afterward, cells were kept in RPMI1640/B27 (+insulin) with 50 μg/mL L-ascorbic acid for days 8–10, followed by a purification medium, RPMI 1640 without glucose, and supplemented with 4 mM sodium L-lactate. Cells were kept in the purification medium for days 12–14, then back to the basal medium of RPMI/B27(+insulin) supplemented with 20 ng/ml triiodothyronine (T3) and 1 μM dexamethasone until day 30. For single myocytes, the monolayers were then dissociated around day 25 with 0.05% trypsin, 1 mg/ml collagenase II, and plated onto Matrigel coated dishes. All voltage clamp recordings were made 5 days after recovery. Experiments were typically performed on hiPSC-CMs at least 20 days post-differentiation.

Voltage Clamp Recordings

Whole-cell patch clamp recordings were obtained for five ionic currents, namely: fast sodium current, I_Na; transient outward potassium current, I_to; L-type calcium current, I_CaL; rapid delayed rectifier potassium current, I_Kr; and hyperpolarization-activated current, I_f. Figure 1 presents the experimental current-voltage (IV) plots for all five currents.

FIGURE 1

Figure 1. Experimental data IV plots. (A) fast sodium current (I_Na); (B) transient outward potassium current (I_to); (C) delayed rectifier potassium current (I_Kr); (D) L-type calcium current (I_CaL); and (E) Hyperpolarization-activated current (I_f). The errors bars represent the experimental standard deviations from multiple cell recordings.

I_Na was measured as described previously (Goodrow et al., 2018). Briefly, the bath solution consisted of 2 mM CaCl₂, 10 mM glucose, 1 mM MgCl₂, 105 mM N-Methyl D Glucamine, 40 mM NaCl, and 10 mM HEPES free acid. The pH was adjusted to 7.4 with HCl. Also, a 300 μM CdCl₂ is added to block the calcium currents which may interfere with I_Na recording. The pipette solution was composed of 1 mM MgCl₂, 15 mM NaCl, 5 mM KCl, 120 mM CsF, 10 mM HEPES, and 10 mM EGTA. Before the experiments, 5 mM Na₂ATP was added, and the pH was adjusted to 7.2 by the addition of CsOH. I_Na was recorded by applying command voltages (25 ms-long pulses) in steps of 5 mV over the range of –80 – +35 mV from a holding potential of –120 mV. All I_Na measurements (n = 8–15) were taken at 20°C and a lowered extracellular (bath) sodium concentration of 40 mM to ensure an adequate voltage control. Temperature and concentration extrapolation facilitated by the Q₁₀ (temperature adjustment factor) correction and Goldman–Katz constant field equation was predicted by Goodrow et al. (2018) to be a factor of 7 and was subsequently adopted in the I_Na curve fitting. More details about the I_Na experimental protocol can be found in Goodrow et al. (2018). I_Kr was recorded as described previously in Treat et al. (2019). In recording the potassium currents, the conventional K⁺ pipette solution of 90 mM K⁺-aspartate, 45 mM KCl, 10 mM NaCl, 1 mM MgCl₂, 10 mM HEPES, 5 mM EGTA, and 5 mM MgATP, with a pH of 7.2 (maintained by adding KOH), was used. I_Kr was recorded by applying 300 ms-long test pulses between –40 and +60 mV in steps of 20 mV from a –80 mV holding potential as described in Doss et al. (2012). I_CaL was recorded using 300 ms-long test pulses applied between –40 and +60 mV in steps of 10 mV from a holding potential of –40 mV. I_to was measured as described in Cordeiro et al. (2013). Briefly, the voltage clamp protocol consisted of a holding potential of –80 mV, followed by a brief –50 mV potential to ensure that all sodium channels were inactivated. This is important as I_to closely follows the depolarization facilitated by Na⁺ influx. Test pulses were then applied in steps of 10 mV from –40 to +50 mV for each clamp voltage. The half-inactivation voltage of I_to, V_1/2 = −41.1 ± 0.2mV, and a slope, k = 6.68± 0.19, were used in our model based on previous studies (Cordeiro et al., 2013). I_f was recorded using a holding potential of –40 mV, followed by pulses from –110 to –40 mV in steps of 10 mV.

All voltage clamp recordings were made using a MultiClamp 700A (Molecular Devices, Sunnyvale, CA, United States). Whole cell patch pipettes were made from glass capillary tubes (1.5 mm O.D., Fisher Scientific, Pittsburg, PA, United States) and pulled on a gravity puller (Narishige Corporation, East Meadow, NY, United States). The resistance ranged from 1.0 to 3.0 MΩ and electronic compensation of series resistance was applied (∼60–70%). Capacitance of the hiPSC-CMs was measured by applying –5 mV voltage steps. Signals were acquired at sampling rate of 10–25 kHz, filtered at 4–6 kHz, and digitized with a Digidata 1322 converter (Molecular Devices, Sunnyvale, CA, United States). Drugs were rapidly applied to the cell using a four-barrel quartz micromanifold (ALA Scientific Instruments Inc., Westbury, NY, United States) placed 500 μm from the cell. Acquisition and analysis were performed using the pClamp9 software. All hiPSC-CM experiments were performed at 36°C, except peak I_Na which was recorded at room temperature (20°C).

Pooled data are presented as Mean ± SEM. Statistical analysis was performed using an ANOVA followed by a Student–Newman–Keuls test using the SigmaStat software.

Mathematical Model of hiPSC-CM

Formulations from six existing cardiac cell models were adopted based on being representative of human cardiomyocyte electrophysiology and/or reproducing spontaneous activity. The formulations of five key currents (I_Na, I_to, I_CaL, I_Kr, and I_f) were optimized by GA based on the experimental data acquired in our lab. The rest of the current formulations, namely, ultrarapid, and slow delayed potassium rectifier currents (I_Kur and I_Ks), sodium-calcium exchanger current (I_NCX or I_NaCa), sodium-potassium exchange pump current (I_NaK), calcium pump current (I_pCa), background sodium and calcium currents (I_bNa and I_bCa), acetylcholine-activated inward-rectifying potassium current (I_KAch), and inward rectifier current (I_K1), were adjusted through proportional scaling based on published literature. The patch clamp experiments in our and other laboratories have revealed a very low to negligible I_K1 in hiPSC-CMs (Doss et al., 2012; Bett et al., 2013; Vaidyanathan et al., 2016). Similarly, most of the hiPSC-CMs in our lab do not express I_Ks. The experimental study by Ma et al. (2011) also reported I_Ks in only five out of 16 cells, and its role in AP repolarization was insignificant. Therefore, we did not include I_K1 and I_Ks in the GA-based optimization. The I_Na formulation was adopted from the Luo–Rudy model (Luo and Rudy, 1994), which is a widely adopted mammalian ventricular myocyte model formulated based on the Hodgkin–Huxley formalism. The I_to formulation was adopted from the Grandi–Pandit human atrial cell model (Grandi et al., 2011) and modified based on the experimental data acquired in hiPSC-CMs by Cordeiro et al. (2013). The I_f formulation was adopted from the human cardiac Purkinje cell model by Stewart et al. (2009). I_CaL and I_Kr were formulated based on the mammalian sinoatrial nodal cell model by Kurata et al. (2002). The choice of this model is due to the inherent automaticity exhibited by hiPSC-CMs similar to that of the nodal cells. The intracellular calcium dynamics in hiPSC-CMs is similar to the nodal cells due to the lack of a mature T-tubules structure (Li et al., 2013). Therefore, the intracellular calcium handing [involving the Ca²⁺ release and uptake fluxes between the sarcoplasmic reticulum (SR) and the cytosol as well as the intra-SR Ca²⁺ transfer flux occurring between the junctional SR and network SR] was adopted from the Kurata model as well. Figure 2 presents a schematic of our hiPSC-CM model depicting the constituent ionic currents and fluxes. Table 1 lists the references to the adopted ionic current formulations. The transmembrane voltage, V_m, is given by the following first order differential equation:

TABLE 1

Table 1. Summary of sources (References) to adopted model formulations.

FIGURE 2

Figure 2. Schematic of the hiPSC-CM model depicting the major ionic currents and fluxes.

$C_m\frac{d{V}_m}{dt}=-\sum _i{I}_{{\text{ion}}_i}(1)$

where, C_m is the membrane capacitance and I_{ion_i} is an instance of the 14 ionic current formulations constituting the model in this study.

Numerical Implementations: Ionic Current Formulations and Whole Cell Integration

In the in silico implementation of current formulations, computational equivalents of the experimental voltage clamping protocols for the five ionic currents mentioned earlier were implemented. An initial forward Euler implementation was done for all five currents; various time steps were tried, where a sufficiently small time-step for computational stability meant a prolonged runtime, as in the case of the rapid I_Na and I_CaL currents, we deferred to an adaptive solver. The I_Na model was simulated using the implicit backward difference (BDF) integration (Iserles, 2008) by implementing voltage test pulses from –80 to +35 mV in the steps of 5 mV with a duration of 25 ms (similar to the experimental protocol). The peak inward current values were then identified and used to produce IV plots. The BDF employed an adaptive computational time-step; where signal gradients were high (as in the case of the rapid depolarization), a suitably small time-step is automatically adopted and vice versa. This allowed for an accelerated computational runtime with a minimal approximation error in single run and multi-generational GA fitting processes. The I_Kr voltage clamp simulations relied on a forward Euler integration with a computational time-step of 0.2 ms. The chosen time values ensured a trade-off between high resolution time discretization and extended time spans. Voltage pulse intensities applied from a –80-mV holding potential ranged from –40 to +60 mV in steps of 20 mV at a 300 ms pulse duration. Peak I_Kr tail currents were measured at –50 mV. Voltage clamping for I_CaL was simulated using the implicit, adaptive time step, BDF integration. The characteristic equations were modified based on ten Tusscher et al. (2004) to allow the dependence on a dynamic Ca²⁺ reversal potential, rather than the fixed-Ca²⁺-potential formalism adopted in the original Kurata model (Kurata et al., 2002). Application of 500 ms voltage pulses ranging from –40 to +60 mV were applied from a holding potential of –30 mV, similar to the experimental protocol. I_to voltage clamp simulations were performed by the explicit forward Euler integration with a fixed time step of 0.5 ms. Voltage pulse intensities applied from a holding potential of –80 mV followed a brief –50 mV potential (as described in the experimental protocol) ranged from –40 to +50 mV in steps of 10 mV. Pulse duration of 500 ms was adopted. I_f implementation was executed by the forward Euler integration with a fixed time step of 10 ms. Voltage pulses that are 500 ms long were applied from a holding potential of –40 mV over a range of –110 to –40 mV in 10 mV increments.

We performed whole cell simulations using the implicit Radau adaptive integration method provided by the solve_ivp module in the SciPy python package. Alternatively, we implemented a faster C/C++ version, where a forward Euler integration was employed.

To enable the pacing of the model by an external stimulus of varying frequencies, the spontaneous firing of APs was disabled by decreasing the maximum conductance of I_f by 50%. The effects of adrenergic stimulation using isoproterenol were simulated by altering the maximum conductance of five currents as described previously (Shah et al., 2019). Briefly, the conductance of L-type Ca²⁺ channel (I_CaL), Na⁺/K⁺ pump (I_NaK), Na⁺/Ca²⁺ exchanger (I_NaCa), and SR Ca²⁺-ATPase (SERCA) were up-regulated by 100, 30, 30, and 20%, respectively, while I_K1 was down-regulated by 20%. The effects of cholinergic stimulation using acetylcholine (ACh or CCh) were simulated by scaling the maximum conductance of I_KACh by 200%. The effects of 4-aminopyridine (4-AP) were simulated by varying levels of I_Kur block.

Model Parametrization

Characteristic equations for the five currents were parametrized, over which, the optimization was done. For each ionic current optimization, parametrization was executed by converting scaling coefficients (a_xi) and half-activation/inactivation voltages (b_xi) and slopes (c_xi) of the sigmoidal gating equations into free parameters as presented in Eq. (2) below.

${\mathrm{\alpha }}_{xi}\left(V_m\right)|{\mathrm{\beta }}_{xi}\left(V_m\right)=\frac{a_{\mathit{xi}}\left(V_m-b_{\mathit{xi}}\right)}{\text{exp}\left(\frac{V_m-b_{\mathit{xi}}}{c_{\mathit{xi}}}\right)}(2)$

Parameter sets for each current were bundled as chromosomes and optimized heuristically by the GA (described next). Accordingly, 20, 18, 20, 6, and 7 free parameters were created for the equations of I_Na, I_to, I_Kr, I_CaL, and I_f, respectively. Detailed parametric equations for each current are given in the Supplementary Material.

The Genetic Algorithm Protocol

Myriads of GA protocols can be formulated by adjusting the population size, fitness criterion, number of generations, crossover, and mutation schemes. A description of the GA protocol adopted and its implementation are discussed in this section.

A starting population was initialized by generating individuals composed of genes selected randomly from a uniformly distributed interval. Constraints (interval limits) were chosen such that the existing model parameter values were contained in the range. The rationale behind this was to initialize the search from a physiologically feasible solution space, as there exists the caveat of multiple individual solutions producing the same/similar model effects. The initial population is, thus, governed mathematically by:

${\mathrm{\rho }}^{\{0\}}\in [\left(1-{\mathrm{\lambda }}_a\right)\mathrm{\theta },\left(1+{\mathrm{\lambda }}_a\right)\mathrm{\theta }](3)$

where, ρ^{0} is the stochastic-drawn initial population parameter set, λ_a is the initial population gene range width determination constant, and θ is the original model parameter. The condition was applied to genes in all chromosomes in the generation. Superscripts in format {x} throughout this paper indicates the generation number.

Population size of 10 × Chromosome size was chosen based on the recommendation by Storn (1996). However, an exception of reduced population size was made in the case of an extended computational runtime, which, in turn, was limited by the number of processors available [28 cores per node for the Ohio Supercomputer HPC clusters (Ohio Supercomputer Center (OSC), 1987) and 24 cores per node for the Extreme Science and Engineering Discovery Environment HPC (Towns et al., 2014)]. A detailed algorithmic description of the GA protocol used in fitting the ionic currents is presented in Figure 3.

FIGURE 3

Figure 3. Algorithmic description of the genetic algorithm (GA) protocol used in fitting the ionic currents. θ: original model chromosome; θ[i]: original model gene; N_x: chromosomal length/number of genes; m: population size; s: offspring proportion (parent-offspring ratio =  1:s); J(): loss function (SSE or RMSE); λ_a: initial population constraint; λ_b: mutation constraint; and χ: an iterable comprising lists of crossover point(s). Length(χ[i]) ∈ (0,s]. All intervals [i.e., [a, b], (a, b] or (a, b)] are uniformly distributed. Parameters belonging to these intervals were randomly drawn.

The SSE loss or RMSE loss for the IV plot points between the model output and experimental data was adopted as the objective function to be minimized to ensure a robust curve fitting. For adaptive time-step, ordinary differential equation solvers available in the SciPy solve_ivp python module—implicit BDF and Radau (Guglielmi and Hairer, 2001) methods as well as explicit Runge–Kutta methods (Dormand and Prince, 1980; Petzold, 1983)—were used. The ionic current models were implemented as importable modules to allow parallel computation using the python multiprocessing module.

Mathematically, the objective function for individual fitness computation, J(θ_n,X) can be expressed as:

$J_{\text{SSE}}({\mathrm{\theta }}_n,X)=\sum _{i=0}^T{\left[X-I\left({\mathrm{\theta }}_n\right)\right]}^2(4)$

$J_{\text{RMSE}}({\mathrm{\theta }}_n,X)=\sqrt{\frac{1}{T}\cdot \sum _{i=0}^T{\left[X-I\left({\mathrm{\theta }}_n\right)\right]}^2}(5)$

where, θ_n is symbolic of a parameter set, X is a set of experimental data points, I(θ_n) represent the ionic current modeled as a function of the parameter set, and T is the number of data points (IV plots).

Single- and multi-point crossover scheme were adopted in fitting the five currents. For each offspring production, crossover points were stochastically drawn from a uniformly distributed interval of positive integers not exceeding the chromosomal length. The convention was to use a larger number of crossover points to compensate for a small population size relative to the number of parameters. Crossover involved the selection of multiple mating pairs from the pool of the best performing individuals. Genes in these mating pairs were then shuffled about the selected crossover-points. The convention chosen here was to generate the same number of offspring as the number of crossover points. This convention, though arbitrary, was sufficient at producing a reasonable fitting. It, however, does not indicate the maximum allowable number of unique offspring. The mutation protocol adopted in this work is an additive scheme involving addition of positive or diminutive proportions of the present offspring parameters to the offspring. The value of these proportions, like in the case of the initial population generation, are randomly selected from a uniformly distributed range. This can be expressed in a matrix form as:

$M=O+\left({\mathrm{\Lambda }}^{°}O\right)(6)$

where, M is the resulting mutant matrix, O is the offspring matrix, and Λ is the proportion matrix. The constraints (mutation coefficient range width determination constant, λ_b) for the elements in the proportion matrix is given mathematically as:

${\mathrm{\Lambda }}_{uv}^{\left(j\right)}\in [-{\mathrm{\lambda }}_b,{\mathrm{\lambda }}_b](7)$

The proportion-offspring multiplication in Eq. (6) is not a traditional matrix multiplication, but rather the Hadamard element-wise multiplication (u and v are row and column indices, respectively, whereas j is generation/iteration number). This way, multi-gene mutations involving all genes per chromosome are executed with random proportions simultaneously. After the two operations, element-wise multiplication and matrix addition, the mutant matrix picks on the offspring size, η = sr; where r is the parent population size. The parent-offspring-mutant ratio is therefore 1:s:s and, consequently, the population size, m = r(1 + 2s). A summary of the various GA protocols adopted in curve fitting for the five experimental-data-complemented ionic currents are presented in Table 2. In all fittings, the GA was run for 100 generations.

TABLE 2

Table 2. Summary of the genetic algorithm protocols for ionic currents.

To quantitatively ascertain the goodness of fit at the end of the GA, the coefficient of determination, R², was used. This statistical measure gives the degree of variability in the data accounted for by the fitted model. The mathematical formula for R² is given in Eq. (8). This value typically ranges between 0 and 1 (it may also be negative in the case of nonlinear regression), with 1 interpreted as an absolute fit accounting for 100% of the variability in the data captured by the fitted model:

$R^2=1-\frac{{\sum }^N{\left(y_i-\widehat{y}\right)}^2}{{\sum }^N{\left(y_i-\overline{y}\right)}^2}=1-\frac{\text{SSE}}{\text{SST}}(8)$

where, y_i is experimental (actual) current values, $\widehat{y}$ is the model fitted current value, and $\overline{y}$ is the mean of the experimental values. SSE and SST are the model residual sum of squared error and total sum of squared errors for the experimental data, respectively.

Once the five abovementioned currents were optimized by GA, the rest of the currents were manually scaled based on published literature to obtain the experimentally recorded AP morphology. The scaling used for the ionic currents and cell related constants are listed in Supplementary Tables 7, 8. Code implementation of the GA parameter optimization and the hiPSC-CM biophysical model are available at https://github.com/Adakwaboah/hiPSC-CM_Computational_Model.

Results

GA-Based Optimizations

Using the protocols described in the Methods section, parameters of the five ionic current formulations (I_Na, I_to, I_CaL, I_Kr, and I_f) were optimized to reproduce the experimental data. The GA process for each ionic current being fitted were performed multiple times to ensure consistency and reproducibility of this meta-heuristic. For each of the fittings, the initial and final model IV plots (representative) are shown in Figure 4 (left panels). In addition, representative time course plots for each fitted current model obtained using simulated patch clamp protocols are shown in the middle panels. The fitted formulations were able to reproduce the experimental mean current recordings (see Supplementary Figure 7 for representative recordings). The right panels in Figure 4 show the near-optimal parameter convergence in the form of losses over 100 generations for the various GA trails per current. The apparent non-uniform decreasing trend in loss plots is evident of the stochastic yet perpetual search for the global minimum. Table 3 summarizes the initial and final fitted loss metrics (mean and standard deviations for multiple trials, n = 5) for all currents. Improvement in fitting can be seen in the increasing R² values toward unity. The original and fitted parameter values for all currents along with the corresponding detailed formulations are given in the Supplementary Tables 1–6. Supplementary Figures 1–4, 6 compare the current activation/inactivation kinetics and corresponding time constants in our optimized model with two recent hiPSC-CM numerical models, namely, models from Paci et al. (2018) and Kernik et al. (2019). Computational runtimes for a single GA trial over the 100-generation period with parallel fitness computation over 28 cores were 11, 13, 7, 12, and 12 min for I_Na, I_CaL, I_f, I_to, and I_Kr, respectively.

TABLE 3

Table 3. Loss and fitness values of the GA-based optimization of currents.

FIGURE 4

Figure 4. GA-based parameter fitting of ionic currents. (A) fast sodium current (I_Na); (B) L-type calcium current (I_CaL); (C) Hyperpolarization-activated current (I_f); (D) transient outward potassium current (I_to); and (E) delayed rectifier potassium current (I_Kr). The left plots show initial and final fitted I-V plots for each current compared with the corresponding experimental values. Middle plots show the simulated time course of current activations based on the corresponding voltage clamp protocols (insets). Right plots show loss (RMSE or SSE) over 100 generations showing convergence during the GA fitting process.

Simulated Action Potentials

The model was simulated for 10 min to achieve the steady states for all currents, fluxes, and ionic concentrations. The steady state (initial conditions) values of the model variables have been reported in Supplementary Table 9. Figure 5 shows a comparison of the spontaneous APs simulated by the hiPSC-CM model as compared to the experimentally recorded ones. Our model was able to accurately reproduce the experimentally observed AP morphology as well as the automaticity of hiPSC-CMs. Time course metrics, such as the AP duration at 50, 75, and 90% repolarization, that is, APD₅₀, APD₇₅, and APD₉₀, respectively; cycle length (CL), that is, AP peak-to-peak duration which includes the diastolic resting phase duration; maximum diastolic potential (MDP); and beats per minute (BPM), were used to assess the model AP morphology compared to the experimental counterpart. Table 4 presents the comparison of the AP metrics produced by the model to those recorded in the experiments, whereas Table 5 lists the calcium transient parameters, such as rise time from 10 to 50% and 10 to 90% of peak value (RT₁₀₅₀ and RT₁₀₉₀, respectively), decay time from 90 to 10% of peak value (DT₉₀₁₀), peak time (T_peak), and frequency of spontaneous activation. Supplementary Table 10 compares the AP metrics and calcium transients of our model to those of Paci et al. (2018) and Kernik et al. (2019). It should be noted that the experimental AP traces used in our model differ considerably than those used in the other models. Wide variability in AP morphologies and calcium transients in hiPSC-CMs has been reported in literature indicating a large population heterogeneity in these cells. For example, Doss et al. (2012) documented CL variations from 327 to 7,063 ms; APD₉₀ variations from 70 to 789 ms; AP amplitudes ranging from 58 to 121 mV; and V_max ranging from 5 to 86 V/s among the hiPSC-CM populations. Figure 6 shows the spontaneous AP generation and corresponding transient plots of the constituent ionic currents. Figure 7 shows the intracellular ionic concentrations ([Na⁺]_i, [K⁺]_i, and [Ca²⁺]_i), whereas Figure 8 shows the Ca²⁺ concentrations in NSR ([Ca²⁺]_up), JSR ([Ca²⁺]_rel), and subspace ([Ca²⁺]_sub) during spontaneous APs. Supplementary Figure 8 shows various transient ionic currents during multiple spontaneous AP firing. Supplementary Figure 9 presents a comparison of spontaneous APs vs. stimulus-elicited APs and the corresponding calcium transients in our model. The stimulus-elicited APs have a steeper Phase 0 depolarization and higher AP magnitude due to a higher I_Na amplitude. The diastolic depolarization seen in spontaneous APs is absent in paced APs due to a partial block of I_f. However, the amplitude of calcium transients is higher in spontaneous APs than the paced APs.

TABLE 4

Table 4. Comparison of action potential parameters.

TABLE 5

Table 5. Intracellular calcium transient parameters in our model during spontaneous AP.

FIGURE 5

Figure 5. Comparison of simulated action potential morphology (blue) to experimentally recorded (orange).

FIGURE 6

Figure 6. Model AP (Panels A and B) and corresponding ionic current transients (Panels C–J) during spontaneous AP generation.

FIGURE 7

Figure 7. Intracellular ionic concentrations. (A) [Na⁺]_i, (B) [K⁺]_i, and (C) [Ca²⁺]_i during spontaneous APs.

FIGURE 8

Figure 8. Ca²⁺ concentrations in (A) NSR ([Ca²⁺]_up), (B) JSR ([Ca²⁺]_rel), and (C) subspace ([Ca²⁺]_sub) during spontaneous APs.

Sensitivity Analysis

To analyze the sensitivity of the baseline hiPSC-CM model, the maximum conductances of the various ionic currents were varied from 0 (complete block) to 200% (2 × enhancement), and their effects on the AP parameters such as APD (APD₅₀, APD₇₅, and APD₉₀), AP magnitude, dV/dt_max, MDP, CL, and spontaneous firing rate (BPM) were studied. Scaling factors from 0 to 200% of the baseline channel conductance values were used in computing the correlation coefficients between various AP parameters and the corresponding ionic current alterations. Figure 9 shows the outcome of the systematic sensitivity analysis of the model. The figure shows a strong positive correlation coefficient between MDP and the ionic currents I_f and I_to, which indicates that an increase in either of these currents causes an increase in the MDP. I_KACh, on the other hand, shows a strong negative correlation with MDP, while showing a positive correlation with CL. I_Kr expectedly shows a strong negative correlation with APD (APD₉₀, APD₇₅, and APD₅₀). This also has implications on CL, which is revealed by a strong negative correlation. Varying the I_Na reveals its strong effect on the AP amplitude corroborated by a strong positive correlation coefficient. It, however, shows an almost neutral correlation with APD. I_f shows a strong positive correlation with MDP, implying a direct relationship between the two. A similar relationship exists between I_f and BPM. I_f also has a strong negative correlation with CL. The I_KACh shows a strong negative correlation with MDP and APD while showing a positive correlation with CL.

FIGURE 9

Figure 9. Sensitivity analysis of model showing color-coded correlation coefficients corresponding to the ionic current variations and their effects on the AP parameters.

Figure 10 shows specific cases of ionic current blockades and their effects on the AP morphology (Figures 10A–F), AP durations (Figures 10G–I), and CL (Figures 10J–L). I_Kr is the principal repolarization current in the rapid repolarization phase in hiPSC-CMs and its blockade prolongs the AP in the Phase 3 as shown in Panel A. With a varying extent of the I_Kr block, APD successively increases (Figure 10G) and cycle length decreases (Figure 10J). Blocking I_CaL significantly shortens the AP (Figures 10B,H), thereby reducing the CL (Figure 10K). Figure 10C shows the reduction in the AP magnitude as a result of the I_Na block. Excessive blockage of I_Na (beyond 40%), however, prevents the initiation of a spontaneous AP. The pacemaker current, I_f, plays a major role in maintaining spontaneous activity. A slight block in this current reduces the frequency of automaticity, and lowers the MDP and magnitude of AP (Figure 10D). More severe blocks (beyond 50%) of I_f inhibit the spontaneous AP generation. Blocking I_KACh prolongs the AP in Phases 2 and 3 as shown in Figures 10E,I. Subsequently, there is an elevation of MDP (not shown) with a higher extent of I_KACh block. Interestingly, although AP was prolonged, the CL was seen decreasing monotonically with the extent of I_KACh block (Figure 10L). It was observed that in the diastolic phase (Phase 4), a reduced repolarization effect of I_KACh favors the diastolic depolarization offered by I_f.

FIGURE 10

Figure 10. Effect of blocking various ionic currents on AP characteristics in the hiPSC-CM model such as AP morphology (A–F), AP durations (G–I), and cycle lengths (J–L).

We further investigated the contribution of two atria-specific currents, namely, I_Kur and I_KACh, which have been recorded in atrial- and nodal-like hiPSC-CMs (Lodrini et al., 2020). Blocking I_Kur as a result of simulating the effects of 4-AP (50 μM) reduced Phase 1 repolarization resulting in AP prolongation and increased AP magnitude as shown in Figure 11A. The extent of AP prolongation was proportional to the extent of I_Kur block. The spontaneous firing rate (BPM), however, remained unchanged. The AP prolongation in all phases of repolarization in our model is in agreement with the recent experimental findings (Devalla et al., 2015). Figure 11B shows the effect of vagal stimulation by carbachol (CCh; 10 μM). The 200% enhancement of I_KACh slowed down the spontaneous activity by ∼10% without any significant effect on the APD.

FIGURE 11

Figure 11. (A) Effect of varying extent of I_Kur block on model AP and CL. (B) Effect of I_KACh enhancement on model AP and CL.

Figure 12A shows the model behavior in hyperkalemic conditions. Our model showed an increased spontaneous firing rate and diastolic depolarization when extracellular K⁺ was increased from 5.4 to 8 mM. This behavior is in agreement with several experimental studies in nodal cells (DiFrancesco et al., 1986; Frace et al., 1992; Hoppe and Beuckelmann, 1998). Figure 12B shows the effects of completely blocking I_NaCa. Inhibition of I_NaCa reduces the spontaneous firing rate and hyperpolarizes the membrane potential in our model. However, it does not abolish the spontaneous activity as reported in (Blinova et al., 2018).

FIGURE 12

Figure 12. (A) Effect of hyperkalemia ([K⁺]_o = 8 mM) conditions on model AP and frequency of automaticity. (B) Effect of I_NaCa blockade on model AP and frequency of automaticity.

Effects of Adrenergic and Cholinergic Stimulation

The model was challenged with stressors such as adrenergic stimulation using isoproterenol and cholinergic stimulation using ACh. For isoproterenol stimulation, the spontaneous APs were suppressed by blocking I_f by 50%. The model was burst paced at 5 Hz for 5 s to overload the SR with Ca²⁺ in the presence of isoproterenol effects. The model exhibited several spontaneous DAD-triggered APs post burst-pacing as shown in Figure 13 (black arrows). Figure 14A shows experimental AP traces when the spontaneously beating hiPSC-CMs were treated with 1 μM ACh. The ACh was administered for 20–40 s which depolarized the membrane and suppressed the spontaneous AP firing. This was unexpected because one would expect the spontaneous activity to slow down similar to that in atrial cells (or as seen in Figure 11B). After the washout of ACh (post 40 s), the spontaneous AP was resumed, albeit at a slightly higher frequency. Our model was able to reproduce the experimental behavior (shown in Figure 14A) when I_KACh was enhanced to 300%, as shown in Figure 14B. Additionally, our model was able to provide useful insights into the unusual behavior observed. The ACh exposure causes a persistently increased I_KACh during diastolic interval, which balances the inward I_f, thereby suppressing the AP generation. After washout, the I_KACh returned to normal levels and the I_f-facilitated spontaneous APs resumed. The elevated intracellular K⁺ levels cause the frequency of spontaneous AP to be higher after washout which gradually returned to a normal firing rate.

FIGURE 13

Figure 13. Model response to 5 Hz burst pacing in presence of isoproterenol effects. Black arrows indicate spontaneous DAD-triggered APs post burst-pacing.

FIGURE 14

Figure 14. (A) A representative experimental AP trace in hiPSC-CM treated with 1 μM ACh. (B) Model behavior with 300% enhancement of I_KACh simulating the effects of ACh. The ACh was administered from 20–40 s followed by washout.

Discussion

We present a robust approach to fit the experimental electrophysiological data to theoretical formulations in order to develop high-fidelity numerical biophysical models. We demonstrate the effectiveness of a GA-based optimization method to develop an in silico model of hiPSC-CMs that accurately reproduces the experimental measurements. The model behavior was extensively validated based on AP morphology, sensitivity analysis, and various ion channel blocking mechanisms.

The experimental electrophysiological data in hiPSC-CMs show a wide range of variability, presumably due to the different techniques used to direct the hiPSCs to the cardiac lineage (Biendarra-Tiegs et al., 2020). This imposes additional challenges in fitting the experimental data to model equations. We utilized the GA method to attain model optimization, which is an evolutionary metaheuristic method inspired by Darwinian evolution and natural selection. Optimization, by this approach, is executed in a stochastic combinatorial manner which makes it less susceptible to getting stuck at local minima, unlike the gradient-based methods, and tend to converge at the global optimized solution (Smith et al., 1995). Previous attempts at implementing GA-based fitting of cardiac models mostly focused on the simple fitting of maximum channel conductance to recapitulate the desired AP morphology (Syed et al., 2005; Bot et al., 2012; Groenendaal et al., 2015; Tomek et al., 2019; Smirnov et al., 2020). We, on the other hand, performed GA optimization at the underlying ionic current level, which is a more realistic and robust approach. Our approach ensures that the model ionic formulations adhere correctly to the experimental recordings and is especially more suitable for modeling hiPSC cells which exhibit a wide range of phenomenological variation. Smirnov et al. (2020) modified the GA protocol by Bot et al. (2012) by adopting vector mutation terms from the Cauchy distribution that promote drastic variance between the mutants and their uncorrelated parent proportions. Tomek et al. (2019), on the other hand, used a multi-objective GA approach for parameter fitting. These works, however, do not mention constraining the range of each gene to a defined neighborhood of the respective original model parameters. Since the uniqueness of a solution is often not guaranteed in such non-convex optimization problems, physiological relevance must be upheld. We address this by imposing the customizable constraint, λ_a, which symmetrically bounds the range from which the initial parameter values (from published physiologically validated models) are drawn. The significance of this constraint, in conjunction with the population size, is the ability to define the initial search resolution. A wider unconstrained parameter range is likely to be initially underrepresented if the population size is not large enough, which, in turn, may favor the speed of convergence (offered by extreme chromosomal variance) over physiological relevance. Another novelty in our GA protocol, to the best of our knowledge, is the correlation between the parent proportion and the number of crossover points. To maintain sufficient diversity, more offspring and mutants are produced as the number of crossover points increase. Furthermore, an increase in the crossover points offer an extensive combinatorial gene shuffling during crossover per generation; we therefore propose a proportionate increase in the population size to accommodate the resulting diverse offspring and mutants.

There have been very few attempts at implementing hiPSC-CM biophysical models due to inconsistent experimental data. An earlier model by Paci et al. (2013) was formulated based on the limited data at that time. More recent models by Paci et al. (2018, 2020), and Koivumäki et al. (2018) incorporated a more realistic calcium handling. The study by Kernik et al. (2019) adopted experimental data from multiple sources in an attempt to cover the range of variability seen in these cells. However, the unresolved inconsistencies in the recording and clamping protocols used by the various sources have the potential to introduce unwarranted deviations in the hiPSC-CM electrophysiological parameter range as well as the generalizability of the baseline model. Our model recapitulates the experimentally recorded hiPSC-CM AP morphology with high fidelity. Moreover, it is able to qualitatively reproduce the experimentally reported (Ma et al., 2011) effects of: (i) APD prolongation caused by I_Kr block, (ii) reduction in AP magnitude and rate of change of upstroke voltage as a result of Na⁺ channel block, (iii) loss of notch-shaped AP (Phase 1 repolarization) due to I_to blockade, (iv) triangulation of AP due to I_Kr or I_to blockades, (v) significant AP shortening due to I_CaL block, (vi) loss of automaticity as a result of I_f blockade, and (vii) alterations in spontaneous firing rate as a result of I_NCX block and hyperkalemic conditions. Our model was also able to produce the experimentally observed effects of channel blocks on the frequency of spontaneous APs (cycle length) and MDPs in hiPSC cells. In the presence of adrenergic stimulation and hypercalcemia, our model was able to generate a DAD-induced triggered activity as a result of SR Ca²⁺ overload.

One of the advantages of our model over the existing ones is the inclusion of atria-specific ionic channels, I_Kur and I_KACh, which have recently been found to be present and functional in hiPSC-CMs (Zhao et al., 2018). The inclusion of I_KACh allows for the investigation of the variability in the spontaneous beating frequency influenced by parasympathetic influences and/or the presence of ACh which have been found to reduce the heart rate. Indeed, our model was able to reproduce and provide explanation on the experimentally observed suppression of automaticity in hiPSC-CMs upon treatment with ACh. Moreover, our model can qualitatively reproduce the effects of atria-specific drugs such as carbachol and 4-AP. I_KACh is actively involved in the maintenance of atrial fibrillation, including chronic atrial fibrillation (Dobrev et al., 2005). Recent advances in cell differentiation techniques use I_KACh and I_Kur as markers to identify atrial-like hiPSC-CMs, which were preferentially produced by retinoic acid treatment (Argenziano et al., 2018). The atrial-like CMs are being increasingly used for disease modeling and pre-clinical screening of antiarrhythmic drugs (Devalla et al., 2015). As such, our model is valuable in disease modeling and simulations of atrial phenotypes.

We adopted simpler Hodgkin-Huxley type current formulations to limit the number of optimization parameters and to keep the simulations computationally tractable. It also avoided overfitting of the experimental data. The main aim of this study was to demonstrate the feasibility of GA-based parametrization of model equations which was done by optimizing five ionic current formulations based on the experimental data. The remaining components of the model, including the intracellular calcium handling, were adopted from a nodal cell model, which has a very close morphological resemblance with the hiPSC-CMs. Nonetheless, the GA-based optimization can be seamlessly extended to the whole cell level as more and more experimental data becomes available. The calcium handling in our model is based on nodal cell formulations and may not represent true calcium transients in hiPSC-CMs. Our model is able to reproduce DAD-triggered APs, but does not produce other arrhythmogenic events such as EADs and alternans. Our model also does not reproduce the hyperkalemia-induced slowing of spontaneous activity as shown by Blinova et al. (2018). The hiPSC-CM population used in Blinova et al. (2018) showed a minor role of I_f in sustaining the spontaneous electrical activity, as revealed by their tests with ivabradine [see the original Figure 2A(i)], which could explain the different response of hiPSC-CMs to hyperkalemia reported in their study. Notwithstanding these limitations, our model behavior is consistent with the findings of various in vitro drug tests, in which such arrhythmic markers were observed only in a portion of the pluripotent cells used for testing (Blinova et al., 2018).

Conclusion

Human induced pluripotent stem cell-derived cardiomyocytes have received significant attention lately; with applications in regenerative medicine, cardiac safety pharmacology, and the implementation of patient specific models for studying drug-induced and inherited cardiac diseases. We present a GA-based model parametrization methodology to incorporate experimental data into numerical models. The proposed method was utilized to formulate a biophysical computer model of hiPSC-CMs based on the experimental data and available literature, as a potential tool for studying and simulating the dynamics of hiPSC-CM electrophysiology. Ionic current formulations of five key currents, namely, fast sodium current (I_Na), transient outward potassium current (I_to), L-type calcium current (I_CaL), rapid delayed rectifier current (I_Kr), and hyperpolarization-activated current (I_f), were optimized by the GA protocol to fit their experimental data. These were then combined with adjusted formulations for nine other currents imported from existing models to faithfully reproduce experimentally obtained hiPSC-CM AP morphology and spontaneous activity. The model was able to accurately reproduce the experimentally recorded AP characteristics and channel blocking drug effects. The model was able to provide insights into the causes of the experimentally observed suppression of automaticity in hiPSC-CMs during cholinergic stimulation. The outcome of this work has implications on validating the biophysics of hiPSC-CMs in their use as viable substitutes for human cardiomyocytes. Specifically, in the study of inherited cardiac disorders and in cardiac safety pharmacology, where drug-induced cardiac disorders are investigated.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

AA performed the model formulations, GA-based parameterization, and model simulations. BT performed the single-cell simulations and sensitivity analysis. PY performed the model assessment and data analysis. JT and JC performed the electrophysiology experiments and experimental data analysis. JC also provided experimental and clinical guidance wherever needed. MB-H reviewed and modified the mathematical formulations. MD conceived the study, assessed the model, guided the students, and confirmed the model outcomes. All authors contributed to the article and approved the submitted version.

Funding

The study was supported by the National Institutes of Health (NIH) award no. 1R15HL145530-01A1 (to MD) and by the Free and Accepted Masons of New York, Florida, Massachusetts, Connecticut, Maryland, Wisconsin, Washington, and Rhode Island (to JC).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2021.675867/full#supplementary-material

References

Akwaboah, A. D., Yamlome, P., Treat, J. A., Cordeiro, J. M., and Deo, M. (2020). “Genetic algorithm for fitting cardiac cell biophysical model formulations,” in Proceedings of the 2020 42nd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC), Montreal, QC, 2463–2466. doi: 10.1109/EMBC44109.2020.9175707

PubMed Abstract | CrossRef Full Text | Google Scholar

Argenziano, M., Lambers, E., Hong, L., Sridhar, A., Zhang, M., Chalazan, B., et al. (2018). Electrophysiologic characterization of calcium handling in human induced pluripotent stem cell-derived atrial cardiomyocytes. Stem Cell Rep. 10, 1867–1878. doi: 10.1016/j.stemcr.2018.04.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Bett, G. C. L., Kaplan, A. D., Lis, A., Cimato, T. R., Tzanakakis, E. S., Zhou, Q., et al. (2013). Electronic “expression” of the inward rectifier in cardiocytes derived from human-induced pluripotent stem cells. Heart Rhythm 10, 1903–1910. doi: 10.1016/j.hrthm.2013.09.061

PubMed Abstract | CrossRef Full Text | Google Scholar

Biendarra-Tiegs, S. M., Secreto, F. J., and Nelson, T. J. (2020). “Addressing variability and heterogeneity of induced pluripotent stem cell-derived cardiomyocytes,” in Advances in Experimental Medicine and Biology, ed. K. Turksen (Cham: Springer), 1–29. doi: 10.1007/5584_2019_350

CrossRef Full Text | Google Scholar

Blinova, K., Dang, Q., Millard, D., Smith, G., Pierson, J., Guo, L., et al. (2018). International multisite study of human-induced pluripotent stem cell-derived cardiomyocytes for drug proarrhythmic potential assessment. Cell Rep. 24, 3582–3592. doi: 10.1016/j.celrep.2018.08.079

PubMed Abstract | CrossRef Full Text | Google Scholar

Bot, C. T., Kherlopian, A. R., Ortega, F. A., Christini, D. J., and Krogh-Madsen, T. (2012). Rapid genetic algorithm optimization of a mouse computational model: Benefits for anthropomorphization of neonatal mouse cardiomyocytes. Front. Physiol. 3:421. doi: 10.3389/fphys.2012.00421

PubMed Abstract | CrossRef Full Text | Google Scholar

Cordeiro, J. M., Nesterenko, V. V., Sicouri, S., Goodrow, R. J., Treat, J. A., Desai, M., et al. (2013). Identification and characterization of a transient outward K+ current in human induced pluripotent stem cell-derived cardiomyocytes. J. Mol. Cell. Cardiol. 60, 36–46. doi: 10.1016/j.yjmcc.2013.03.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Courtemanche, M., Ramirez, R. J., and Nattel, S. (1998). Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model. Am. J. Physiol. Heart Circ. Physiol. 275, H301–H321. doi: 10.1152/ajpheart.1998.275.1.h301

PubMed Abstract | CrossRef Full Text | Google Scholar

Devalla, H. D., Schwach, V., Ford, J. W., Milnes, J. T., El-Haou, S., Jackson, C., et al. (2015). Atrial-like cardiomyocytes from human pluripotent stem cells are a robust preclinical model for assessing atrial-selective pharmacology. EMBO Mol. Med. 7, 394–410. doi: 10.15252/emmm.201404757

PubMed Abstract | CrossRef Full Text | Google Scholar

Di Francesco, D., and Noble, D. (1985). A model of cardiac electrical activity incorporating ionic pumps and concentration changes. Philos. Trans. R. Soc. Lond. B Biol. Sci. 307, 353–398.

Google Scholar

DiFrancesco, D., Ferroni, A., Mazzanti, M., and Tromba, C. (1986). Properties of the hyperpolarizing-activated current (if) in cells isolated from the rabbit sino-atrial node. J. Physiol. 377, 61–88. doi: 10.1113/jphysiol.1986.sp016177

PubMed Abstract | CrossRef Full Text | Google Scholar

Dobrev, D., Friedrich, A., Voigt, N., Jost, N., Wettwer, E., Christ, T., et al. (2005). The G protein-gated potassium current IK, ACh is constitutively active in patients with chronic atrial fibrillation. Circulation 112, 3697–3706. doi: 10.1161/CIRCULATIONAHA.105.575332

PubMed Abstract | CrossRef Full Text | Google Scholar

Dormand, J. R., and Prince, P. J. (1980). A family of embedded Runge-Kutta formulae. J. Comput. Appl. Math. 6, 19–26. doi: 10.1016/0771-050X(80)90013-3

CrossRef Full Text | Google Scholar

Doss, M. X., Di Diego, J. M., Goodrow, R. J., Wu, Y., Cordeiro, J. M., Nesterenko, V. V., et al. (2012). Maximum diastolic potential of human induced pluripotent stem cell-derived cardiomyocytes depends critically on IKr. PLoS One 7:e40288. doi: 10.1371/journal.pone.0040288

PubMed Abstract | CrossRef Full Text | Google Scholar

Egashira, T., Yuasa, S., Suzuki, T., Aizawa, Y., Yamakawa, H., Matsuhashi, T., et al. (2012). Disease characterization using LQTS-specific induced pluripotent stem cells. Cardiovasc. Res. 95, 419–429. doi: 10.1093/cvr/cvs206

PubMed Abstract | CrossRef Full Text | Google Scholar

Fatima, A., Xu, G., Shao, K., Papadopoulos, S., Lehmann, M., Arnáiz-Cot, J. J., et al. (2011). In vitro modeling of ryanodine receptor 2 dysfunction using human induced pluripotent stem cells. Cell. Physiol. Biochem. 28, 579–592. doi: 10.1159/000335753

PubMed Abstract | CrossRef Full Text | Google Scholar

Frace, A. M., Maruoka, F., and Noma, A. (1992). External K+ increases Na+ conductance of the hyperpolarization-activated current in rabbit cardiac pacemaker cells. Pflügers Arch. Eur. J. Physiol. 421, 94–96. doi: 10.1007/BF00374739

CrossRef Full Text | Google Scholar

Gintant, G., Burridge, P., Gepstein, L., Harding, S., Herron, T., Hong, C., et al. (2019). Use of human induced pluripotent stem cell-derived cardiomyocytes in preclinical cancer drug cardiotoxicity testing: a scientific statement from the American heart association. Circ. Res. 125, e75–e92. doi: 10.1161/RES.0000000000000291

PubMed Abstract | CrossRef Full Text | Google Scholar

Gintant, G., Fermini, B., Stockbridge, N., and Strauss, D. (2017). The evolving roles of human iPSC-derived cardiomyocytes in drug safety and discovery. Cell Stem Cell 21, 14–17. doi: 10.1016/j.stem.2017.06.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Goodrow, R. J., Desai, S., Treat, J. A., Panama, B. K., Desai, M., Nesterenko, V. V., et al. (2018). Biophysical comparison of sodium currents in native cardiac myocytes and human induced pluripotent stem cell-derived cardiomyocytes. J. Pharmacol. Toxicol. Methods 90, 19–30. doi: 10.1016/j.vascn.2017.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Grandi, E., Pandit, S. V., Voigt, N., Workman, A. J., Dobrev, D., Jalife, J., et al. (2011). Human atrial action potential and Ca 2+ model: Sinus rhythm and chronic atrial fibrillation. Circ. Res. 109, 1055–1066. doi: 10.1161/CIRCRESAHA.111.253955

PubMed Abstract | CrossRef Full Text | Google Scholar

Grandi, E., Pasqualini, F. S., and Bers, D. M. (2010). A novel computational model of the human ventricular action potential and Ca transient. J. Mol. Cell. Cardiol. 48, 112–121. doi: 10.1016/j.yjmcc.2009.09.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Groenendaal, W., Ortega, F. A., Kherlopian, A. R., Zygmunt, A. C., Krogh-Madsen, T., and Christini, D. J. (2015). Cell-specific cardiac electrophysiology models. PLoS Comput. Biol. 11:e1004242. doi: 10.1371/journal.pcbi.1004242

PubMed Abstract | CrossRef Full Text | Google Scholar

Guglielmi, N., and Hairer, E. (2001). Implementing Radau IIA methods for stiff delay differential equations. Computing (Vienna/New York) 67, 1–12. doi: 10.1007/s006070170013

CrossRef Full Text | Google Scholar

Hodgkin, A. L., and Huxley, A. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117, 500–544.

Google Scholar

Hoppe, C. U., and Beuckelmann, D. J. (1998). Characterization of the hyperpolarization-activated inward current in isolated human atrial myocytes. Cardiovasc. Res. 38, 788–801. doi: 10.1016/S0008-6363(98)00047-9

CrossRef Full Text | Google Scholar

Hwang, H. S., Kryshtal, D. O., Feaster, T. K., Sánchez-Freire, V., Zhang, J., Kamp, T. J., et al. (2015). Comparable calcium handling of human iPSC-derived cardiomyocytes generated by multiple laboratories. J. Mol. Cell. Cardiol. 85, 79–88. doi: 10.1016/j.yjmcc.2015.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Iserles, A. (2008). A First Course in the Numerical Analysis of Differential Equations: Cambridge Texts in Applied Mathematics, 2nd Edn. Cambridge: Cambridge University Press. doi: 10.1017/CBO9780511995569

CrossRef Full Text | Google Scholar

Jung, C. B., Moretti, A., Mederos y Schnitzler, M., Iop, L., Storch, U., Bellin, M., et al. (2012). Dantrolene rescues arrhythmogenic RYR2 defect in a patient-specific stem cell model of catecholaminergic polymorphic ventricular tachycardia. EMBO Mol. Med. 4, 180–191. doi: 10.1002/emmm.201100194

PubMed Abstract | CrossRef Full Text | Google Scholar

Kernik, D. C., Morotti, S., Wu, H., Garg, P., Duff, H. J., Kurokawa, J., et al. (2019). A computational model of induced pluripotent stem-cell derived cardiomyocytes incorporating experimental variability from multiple data sources. J. Physiol. 597, 4533–4564. doi: 10.1113/JP277724

PubMed Abstract | CrossRef Full Text | Google Scholar

Koivumäki, J. T., Naumenko, N., Tuomainen, T., Takalo, J., Oksanen, M., Puttonen, K. A., et al. (2018). Structural immaturity of human iPSC-derived cardiomyocytes: In silico investigation of effects on function and disease modeling. Front. Physiol. 9:80. doi: 10.3389/fphys.2018.00080

PubMed Abstract | CrossRef Full Text | Google Scholar

Kurata, Y., Hisatome, I., Imanishi, S., and Shibamoto, T. (2002). Dynamical description of sinoatrial node pacemaking: Improved mathematical model for primary pacemaker cell. Am J. Physiol. Heart Circ. Physiol. 283, 2074–2101. doi: 10.1152/ajpheart.00900.2001

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, S., Chen, G., and Li, R. A. (2013). Calcium signalling of human pluripotent stem cell-derived cardiomyocytes. J. Physiol. 591, 5279–5290. doi: 10.1113/jphysiol.2013.256495

PubMed Abstract | CrossRef Full Text | Google Scholar

Lodrini, A. M., Barile, L., Rocchetti, M., and Altomare, C. (2020). Human induced pluripotent stem cells derived from a cardiac somatic source: insights for an in-vitro cardiomyocyte platform. Int. J. Mol. Sci. 21:507. doi: 10.3390/ijms21020507

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, C. H., and Rudy, Y. (1994). A dynamic model of the cardiac ventricular action potential: I. Simulations of ionic currents and concentration changes. Circ. Res. 74, 1071–1096. doi: 10.1161/01.RES.74.6.1071

CrossRef Full Text | Google Scholar

Ma, J., Guo, L., Fiene, S. J., Anson, B. D., Thomson, J. A., Kamp, T. J., et al. (2011). High purity human-induced pluripotent stem cell-derived cardiomyocytes: electrophysiological properties of action potentials and ionic currents. Am. J. Physiol. Heart Circ. Physiol. 301, 2006–2017. doi: 10.1152/ajpheart.00694.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Moretti, A., Bellin, M., Welling, A., Jung, C. B., Lam, J. T., Bott-Flügel, L., et al. (2010). Patient-specific induced pluripotent stem-cell models for long-QT syndrome. N. Engl. J. Med. 363, 1397–1409.

Google Scholar

Narsinh, K. H., Sun, N., Sanchez-Freire, V., Lee, A. S., Almeida, P., Hu, S., et al. (2011). Single cell transcriptional profiling reveals heterogeneity of human induced pluripotent stem cells. J. Clin. Investig. 121, 1217–1221. doi: 10.1172/JCI44635

PubMed Abstract | CrossRef Full Text | Google Scholar

Ohio Supercomputer Center (OSC) (1987). Ohio Supercomputer Center. Columbus, OH: OSC.

Google Scholar

Paci, M., Hyttinen, J., Aalto-Setälä, K., and Severi, S. (2013). Computational models of ventricular-and atrial-like human induced pluripotent stem cell derived cardiomyocytes. Ann. Biomed. Eng. 41, 2334–2348.

Google Scholar

Paci, M., Passini, E., Klimas, A., Severi, S., Hyttinen, J., Rodriguez, B., et al. (2020). All-optical electrophysiology refines populations of in silico human iPSC-CMs for drug evaluation. Biophys. J. 118, 2596–2611. doi: 10.1016/j.bpj.2020.03.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Paci, M., Pölönen, R. P., Cori, D., Penttinen, K., Aalto-Setälä, K., Severi, S., et al. (2018). Automatic optimization of an in silico model of human iPSC derived cardiomyocytes recapitulating calcium handling abnormalities. Front. Physiol. 9:709. doi: 10.3389/fphys.2018.00709

PubMed Abstract | CrossRef Full Text | Google Scholar

Petzold, L. (1983). Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations. SIAM J. Sci. Stat. Comput. 4, 136–148. doi: 10.1137/0904010

CrossRef Full Text | Google Scholar

Rhinehart, R. R. (2016). “Finite difference schemes,” in Nonlinear Regression Modeling for Engineering Applications: Modeling, Model Validation, and Enabling Design of Experiments, (New York, NY: Wiley), 157–175. doi: 10.1002/9781118597972

CrossRef Full Text | Google Scholar

Shah, C., Jiwani, S., Limbu, B., Weinberg, S., and Deo, M. (2019). Delayed afterdepolarization-induced triggered activity in cardiac purkinje cells mediated through cytosolic calcium diffusion waves. Physiol. Rep. 7:e14296. doi: 10.14814/phy2.14296

PubMed Abstract | CrossRef Full Text | Google Scholar

Smirnov, D., Pikunov, A., Syunyaev, R., Deviatiiarov, R., Gusev, O., Aras, K., et al. (2020). Genetic algorithm-based personalized models of human cardiac action potential. PLoS One 15:e0231695. doi: 10.1371/journal.pone.0231695

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, A. E., Gulsen, M., and Tate, D. M. (1995). A genetic algorithm approach to curve fitting. Int. J. Product. Res. 33, 1911–1923. doi: 10.1080/00207549508904789

CrossRef Full Text | Google Scholar

Stewart, P., Aslanidi, O. V., Noble, D., Noble, P. J., Boyett, M. R., and Zhang, H. (2009). Mathematical models of the electrical action potential of Purkinje fibre cells. Philos. Trans. Ser. A Math. Phys. Eng. Sci. 367, 2225–2255. doi: 10.1098/rsta.2008.0283

PubMed Abstract | CrossRef Full Text | Google Scholar

Storn, R. (1996). “On the usage of differential evolution for function optimization,” in Proceedings of the Biennial Conference of the North American Fuzzy Information Processing Society - NAFIPS, (Berkeley, CA: IEEE), 519–523. doi: 10.1109/nafips.1996.534789

CrossRef Full Text | Google Scholar

Syed, Z., Vigmond, E., Nattel, S., and Leon, L. J. (2005). Atrial cell action potential parameter fitting using genetic algorithms. Med. Biol. Eng. Comput. 43, 561–571. doi: 10.1007/BF02351029

PubMed Abstract | CrossRef Full Text | Google Scholar

Takahashi, K., Tanabe, K., Ohnuki, M., Narita, M., Ichisaka, T., Tomoda, K., et al. (2007). Induction of pluripotent stem cells from adult human fibroblasts by defined factors. Cell 131, 861–872. doi: 10.1016/j.cell.2007.11.019

PubMed Abstract | CrossRef Full Text | Google Scholar

ten Tusscher, K. H. W. J., Noble, D., Noble, P. J., and Panfilov, A. V. (2004). A model for human ventricular tissue. Am. J. Physiol. Heart Circ. Physiol. 286, H1573–H1589. doi: 10.1152/ajpheart.00794.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Tomek, J., Bueno-Orovio, A., Passini, E., Zhou, X., Minchole, A., Britton, O., et al. (2019). Development, calibration, and validation of a novel human ventricular myocyte model in health, disease, and drug block. eLife 8:e48890. doi: 10.7554/eLife.48890

PubMed Abstract | CrossRef Full Text | Google Scholar

Towns, J., Cockerill, T., Dahan, M., Foster, I., Gaither, K., Grimshaw, A., et al. (2014). XSEDE: accelerating scientific discovery. Comput. Sci. Eng. 16, 62–74. doi: 10.1109/MCSE.2014.80

CrossRef Full Text | Google Scholar

Treat, J. A., Goodrow, R. J., Bot, C. T., Haedo, R. J., and Cordeiro, J. M. (2019). Pharmacological enhancement of repolarization reserve in human induced pluripotent stem cells derived cardiomyocytes. Biochem. Pharmacol. 169:113608. doi: 10.1016/j.bcp.2019.08.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Vaidyanathan, R., Markandeya, Y. S., Kamp, T. J., Makielski, J. C., January, C. T., and Eckhardt, L. L. (2016). IK1-enhanced human-induced pluripotent stem cell-derived cardiomyocytes: An improved cardiomyocyte model to investigate inherited arrhythmia syndromes. Am. J. Physiol. Heart Circ. Physiol. 310, H1611–H1621. doi: 10.1152/ajpheart.00481.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Z., Lan, H., El-Battrawy, I., Li, X., Buljubasic, F., Sattler, K., et al. (2018). Ion channel expression and characterization in human induced pluripotent stem cell-derived cardiomyocytes. Stem Cells Int. 2018:6067096. doi: 10.1155/2018/6067096

PubMed Abstract | CrossRef Full Text | Google Scholar