Impact Factor 4.566 | CiteScore 5.6
More on impact ›

ORIGINAL RESEARCH article

Front. Physiol., 26 June 2018 | https://doi.org/10.3389/fphys.2018.00709

Automatic Optimization of an in Silico Model of Human iPSC Derived Cardiomyocytes Recapitulating Calcium Handling Abnormalities

  • 1Faculty of Biomedical Sciences and Engineering, BioMediTech Institute, Tampere University of Technology, Tampere, Finland
  • 2Faculty of Medicine and Life Sciences, BioMediTech Institute, University of Tampere, Tampere, Finland
  • 3Department of Electrical, Electronic and Information Engineering “Guglielmo Marconi”, University of Bologna, Cesena, Italy
  • 4Heart Hospital, Tampere University Hospital, Tampere, Finland

The growing importance of human induced pluripotent stem cell-derived cardiomyoyctes (hiPSC-CMs), as patient-specific and disease-specific models for studying cellular cardiac electrophysiology or for preliminary cardiotoxicity tests, generated better understanding of hiPSC-CM biophysical mechanisms and great amount of action potential and calcium transient data. In this paper, we propose a new hiPSC-CM in silico model, with particular attention to Ca2+ handling. We used (i) the hiPSC-CM Paci2013 model as starting point, (ii) a new dataset of Ca2+ transient measurements to tune the parameters of the inward and outward Ca2+ fluxes of sarcoplasmic reticulum, and (iii) an automatic parameter optimization to fit action potentials and Ca2+ transients. The Paci2018 model simulates, together with the typical hiPSC-CM spontaneous action potentials, more refined Ca2+ transients and delayed afterdepolarizations-like abnormalities, which the old Paci2013 was not able to predict due to its mathematical formulation. The Paci2018 model was validated against (i) the same current blocking experiments used to validate the Paci2013 model, and (ii) recently published data about effects of different extracellular ionic concentrations. In conclusion, we present a new and more versatile in silico model, which will provide a platform for modeling the effects of drugs or mutations that affect Ca2+ handling in hiPSC-CMs.

Introduction

Human induced pluripotent stem cell-derived cardiomyocytes (hiPSC-CMs) are cardiac cells derived from stem cells, which have been produced by donor's differentiated cells by means of reprogramming (Takahashi et al., 2007). The role of hiPSC-CMs, has become more and more pervasive in basic electrophysiological studies as well as in applied research, such as pharmacological tests, since their discover in 2007. As an in vitro human model, hiPSC-CMs strongly impacted the study of biophysical mechanisms underlying cardiac electrophysiology at cellular level, both in control and diseased conditions. Especially hiPSC-CMs' patient- and disease-specificity is fundamental to assess the effects of genetic mutations, such as Long QT (LQT) (Moretti et al., 2010; Lahti et al., 2012; Ma et al., 2013), catecholaminergic polymorphic ventricular tachycardia (CPVT) (Kujala et al., 2012) and hypertrophic cardiomyopathy (HCM) (Ojala et al., 2016), on the functionality of cardiomyocytes. Ever since the beginning of the Comprehensive In vitro Proarrhythmic Assay (CIPA) (http://cipaproject.org/) in 2013, hiPSC-CMs have had a dramatic impact on pharmacology, serving as a powerful in vitro model to test the in silico model predictions regarding cardiac safety or drug toxicity at cellular level.

During the last 10 years, many progresses were done in terms of efficiency of hiPSC-CM production and availability of commercial cell lines. This enabled also the hiPSC-CM electrophysiological and pharmacological evaluation by means of medium-throughput (Rajamohan et al., 2016) or even high-throughput systems (Entcheva and Bub, 2016; Klimas et al., 2016), where the use of voltage- and calcium-sensitive dyes has been combined with hiPSC-CM optogenetic stimulation. Together with the availability of these experimental data, new methods have been developed to process them (Björk et al., 2017; Ahola et al., 2018).

The importance of Ca2+ cycling to basic cardiac functionality, together with the growing availability of hiPSC-CM Ca2+ cycling data, makes Ca2+ transients, and biomarkers computed onto them, as interesting as action potential (AP) measurements. Ca2+ is fundamental in the heart excitation-contraction (EC) coupling, i.e., how the electrical and the mechanical properties of the heart are linked together and how the AP leads to the cardiomyocyte contraction. The elements involved in this phenomenon are the L-type Ca2+ channels, the sarcoplasmic reticulum (SR) and the sarcomeres, i.e., the contractile unit of the cardiomyocyte. During the AP upstroke, the L-type Ca2+ channels open and Ca2+ flows into the cardiomyocyte. This Ca2+ influx is sufficient to trigger the Ca2+ release from SR through the ryanodine-sensitive Ca2+ channels, which increases the cytosolic Ca2+. Such amount of Ca2+ allows starting the crossbridge cycle, which is at the basis of the cardiomyocite contraction and continues until Ca2+ is restored to its basal cytosolic concentration. This is mainly done by the SERCA-2 pump, which reabsorbs Ca2+ from the cytosol into the SR. Moreover, cytosolic Ca2+ is extruded into the extracellular space by the Na+/Ca2+ exchanger (INaCa) and by the sarcolemmal Ca2+ pump (IpCa). All these mechanisms make the intracellular Ca2+ concentration change, thus producing Ca2+ transients associated to the APs (Walker and Spinale, 1999).

In 2013 we published the first in silico hiPSC-CMs model (Paci et al., 2013), based on our previous model of cardiomyocytes derived from human embryonic stem cells (Paci et al., 2012) and on the experimental data by Ma et al. (2011). This model has been widely used for computational studies, such as (i) the prediction of drug effects on cardiac electrophysiology (Paci et al., 2015; Lei et al., 2017), (ii) the model extension to multielectrode array simulations (Raphel et al., 2017), and (iii) the assessment of hiPSC-CM electrophysiological variability in control and mutant cells, by means of populations of in silico hiPSC-CMs (Paci et al., 2017). However, one of the Paci2013 model limitations resides in its formulations of the Ca2+ handling system: especially the Ca2+ release from the sarcoplasmic reticulum (SR) is formulated with the functional but quite elementary Ca2+ release from the TenTusscher2004 model (ten Tusscher et al., 2004).

In this work, we propose an updated version of the Paci2013 hiPSC-CM ventricular-like model with a more flexible Ca2+ handling formulation. The Ca2+ transients produced by the model were calibrated on experimental Ca2+ imaging data recorded in our laboratory on hiPSC-CMs. Moreover, the fine tuning of the model parameters was performed by means of an automatic optimization technique (Fabbri et al., 2017), in order to reproduce realistic AP and Ca2+ transient shapes and to speed up the parameter tuning phase. Parameter optimization affected only the parameters representing the Ca2+ SR fluxes and a very limited set of parameters of membrane currents in order to be consistent with the voltage clamp experiment fitting done by Paci et al. (Paci et al., 2013). Finally, the resulting updated model was validated against ion current blocking data.

Materials and Methods

New Ca2+ Handling System Formulation

The main limitation of the original Paci2013 model (Paci et al., 2013) and its following minor updates (Paci et al., 2015, 2017) is the simplified description of the Ca2+ release from SR, i.e., the release current Irel. Irel was formulated as in the TenTusscher2004 model (ten Tusscher et al., 2004) in the following way:

Irel=Krel·(crel+arel·CaSR2brel2+CaSR2)·d·g

where Krel, arel, brel, and crel are constants, CaSR is the Ca2+ concentration in the SR reticulum, d is the L-type Ca2+ current (ICaL) voltage-dependent activation gate, and g is the Irel specific inactivation gate. Therefore, Irel is not activated by the cytosolic Ca2+ concentration Cai sensed by the Ryanodine-sensitive receptor located on the SR membrane, but it is triggered the same way as ICaL. Due to this mechanism, the model was not able to produce proarrhythmic triggers, such as delayed afterdepolarization (DADs) (Fink et al., 2011).

We reformulated Irel according to the formulation used by Koivumäki et al. (2011) for the human atrial myocyte:

Irel=Irel,max·RyRCaSR·RyRo·RyRc·(CaSR-Cai)

where Irel,max represents the maximum Ca2+ release from SR, RyRCaSR is the dependence on CaSR, RyRo is the open (activation) gating variable, RyRc the closed (inactivation) gating variable. A third gating variable, RyRa, was used to modulate the working point (adaptation) of the RyRo and RyRc gates according to the cytosolic Ca2+ concentration as in Koivumäki et al. (2011). The formulation of the single gating variables is the following

RyRCaSR=1-11+eCaSR-0.30.1  RyRa,ss=RyRa1-RyRa21+eCai-RyRa,halfRyRa,k  RyRo,ss=1-11+eCai-(RyRa+RyRo,half)RyRo,k  RyRc,ss=11+eCai-(RyRa+RyRc,half)RyRc,k

where the subscript ss indicates the steady state value of the gating variable and the subscripts half and k indicate the half Ca2+ concentration and the slope, respectively, of the gating variable steady state. The values of the constants RyRa1, RyRa2 and the gating variable half activations RyRa,half, RyRo,half and RyRc,half were then optimized (together with other parameters) as described in section “Parameter Optimization.” We report the full set of equations in section 1 of the Supplementary Material.

The formulation of the other SR fluxes, namely the SERCA pump (Iup) and the leakage current (Ileak) were not changed with respect to the Paci2013 model (Paci et al., 2013). However, their parameters went through the optimization process, as we describe in section Parameter Optimization. In the following we will refer as Paci2013 to the model version presented in Paci et al. (2017).

Parameter Optimization

The parameter optimization process was adapted from Fabbri et al. (2017).

Shortly, parameter optimization was done using the Matlab© function fminsearch, which implements the Nelder-Mead Simplex Method, as reported in Fabbri et al. (2017). Such function minimizes a cost function built on the experimental biomarkers we want the model to simulate. The following equations show the cost function structure (Fabbri et al., 2017):

                            Cost=i=1NbiomarkersCostiCosti={|bi,Expbi,Sim|SDiSDi·wiif|bi,Expbi,Sim|>SD0otherwise

where Costi and SDi represent the cost and the standard deviation for a single biomarker bi.

The structure of the cost function Cost is consistent with the cost function used in Fabbri et al. (2017). The contribution of each biomarker to the overall cost is zero if the biomarker is within its ±SD range, otherwise it grows linearly. Moreover, each contribution Costi was weighted according to the respective weight wi. This results in a non-linear cost function, which is zero if all the biomarkers are within the experimental ranges and greater than zero if at least one biomarker is out of range. The AP experimental biomarkers considered for the cost function were: AP amplitude (APA), maximum diastolic potential (MDP), cycle length (CL), maximum upstroke velocity (Vmax), AP duration at 10, 30, and 90% of repolarization (APD10, APD30, APD90) and AP shape factor (Triangulation). Triangulation is a shape factor used by Ma et al. (2011) to discriminate between atrial-like (Triangulation < 1.5) and ventricular-like (Triangulation>1.5) APs, and it is computed as:

Triangulation=APD30-APD40APD70-APD80

The Ca2+ transient experimental biomarkers were: duration of the Ca2+ transient (DURATION), time to peak (TPEAK), rise time from 10 to 50% (RT1050), rise time from 10 to 90% (RT1090), decay time from 90 to 10% (DT9010) and the Ca2+ transient rate (FREQ). The full list of the parameters optimized by this method and their original values are reported in Table 1. For optimization, we only chose parameters related to Ca2+ handling and additional parameters of INaCa and the Na+/K+ pump (INaK), i.e., those currents for which the Ma et al. dataset (Ma et al., 2011) did not provide experimental data. The parameter values were constrained in a range [−20%, +20%] with respect to their nominal value in the Paci2013 model, in order to avoid non-physiological values, such as negative conductances. The initial parameter values for RyRa1, RyRa2, RyRa,half, RyRo,half, RyRc,half were rescaled to 1/10 of their values in Koivumäki et al. (2011) before the optimization procedure, in order to adapt them to the Ca2+ concentrations in the Paci2013 model, and then optimized as the other parameters. All the other ionic current parameters were kept as in Paci et al. (2013) and Paci et al. (2017) for the late Na+ current (INaL) only.

TABLE 1
www.frontiersin.org

Table 1. Parameters chosen for optimization.

To constrain the model parameters, two different experimental datasets were used. AP data were taken from Ma et al. (2011). This is the same dataset used to calibrate the Paci2013 model and it includes biomarkers computed on spontaneous ventricular-like hiPSC-CM APs recorded at 35–37°C. A second, new, dataset was obtained through Ca2+ transient recordings performed at the BioMediTech Institute (Tampere, Finland).

Calcium Recordings in hiPSC-CMs

This study was carried out in accordance with the recommendations of Guidelines of the Ethics Committee of Pirkanmaa Hospital District (Tampere, Finland). The protocol was approved by the Ethics Committee of Pirkanmaa Hospital District (Aalto-Setälä R08070). All subjects gave written informed consent in accordance with the Declaration of Helsinki. New Ca2+ transient dataset were recorded at the BioMediTech Institute from healthy control hiPSC-CMs at 35–37°C. The generation and characterization of the control hiPSC line and cardiac differentiation were done as described earlier (Ojala et al., 2016). Cardiomyocytes plated on a coverslip were loaded with 4 μM Fluo-4 AM (Thermo Fisher Scientific) for 30 min and de-esterified for 10 min in perfusate medium: (in mM) 137 NaCl, 5 KCl, 0.44 KH2PO4, 20 HEPES, 4.2 NaHCO3, 5 D-glucose, 2 CaCl2, 1.2 MgCl2, and 1 Na-pyruvate dissolved in H2O (all from Sigma Aldrich). pH of the perfusate medium was adjusted to 7.4 with NaOH (Sigma Aldrich). Perfusate was heated with an inline heater SH-27B controlled with a TC-324B controller unit and input into a RC-25 imaging chamber (all Warner instruments Inc., CT, USA). Calcium kinetics of spontaneously beating cardiomyocytes were imaged with an inverted Olympus IX70 microscope using UApo/340 0,75NA 20x air objective (Olympus, Tokyo, Japan) and recorded with ANDOR iXon 885 EM-CCD camera (Andor Technology, Belfast, Northern Ireland) using 2 × 2 binning and synchronized with a Polychrome V light source by a real time DPS control unit. LiveAcquisition software (TILL Photonics, Munich, Germany) was used to control light source and camera during recording. Fluo-4 was excited at 490 nm wavelength and the emission was recorded through Olympus U-MF2 Alexa 488 band-pass filter cube (ex.470–495, em.525/50 nm).

Data were recorded from 15 cells for a total of 218 transients. For each biomarker mean value and standard deviation (SD) were computed for model calibration. The experimental values for the biomarkers and their weights are reported in Table 2.

TABLE 2
www.frontiersin.org

Table 2. Experimental and simulated values of the biomarkers considered for the parameter optimization.

Model Validation With Current Blocker Simulations

In order to validate the new hiPSC-CM model after the introduction of the new Irel formulation and the parameters optimization by means of our Ca2+ transient data, we chose to replicate the current blocker simulations performed in Paci et al. (2013). The model was paced at 1 Hz for 800 s to reach its steady state, then the current blocker was simulated by reducing the maximum conductance of the affected current and finally the AP biomarkers were computed after 400 s. The current blockers considered in the simulations were the same used by Ma et al. (2011): tetrodotoxine (TTX, INa blocker), nifedipine (NIFED, ICaL blocker), E4031 (IKr blocker) and 3R4S-Chromanol 293B (CHR, IKs blocker). For TTX, NIFED and E4031 we considered the following IC50 values: 0.64 μM (Ma et al., 2011), 0.038 μM (Ma et al., 2011) and 100 nM (Sanguinetti and Jurkiewicz, 1990; Gerlach et al., 2010) respectively. In Ma et al. (2011) CHR had small effects on the AP biomarkers, therefore we tested a range of block levels: 30, 50, 70, and 90%. The stimulus current was 550 pA for NIFED, E4031 and CHR, while we chose 750 pA for TTX in order to trigger APs also at the highest blocker concentrations.

Delayed Afterdepolarizations

In order to trigger DADs in spontaneous APs, Ca2+ overload had to be simulated: we chose to increase the superfusate Ca2+ concentration (Volders et al., 2000). To ensure that Ca2+ overload was simulated only once the model was in steady state, we ran a 800 s simulation without external pacing. From these steady state conditions, we increased the extracellular Ca2+ concentration by setting it to 3.945 mM and again we simulated the spontaneous APs. We ran this protocol with both the Paci2013 and the Paci2018 models.

Results

Spontaneous Action Potentials and Ca2+ Transients

The automatic optimization process provided a new set of parameters, which is reported in Table 1, fourth column. The simulated spontaneous APs and Ca2+ transients were in agreement with all the biomarker variability ranges (mean ± SD), as reported in the last column of Table 2. Figure 1 shows a comparison between the Paci2018 model (in solid blue) and the Paci2013 model (in dashed red). Figure 2 illustrates the simulated Ca2+ transients and traces from four illustrative cells: the comparison highlights that the simulated Ca2+ transient contour is fully in agreement with our experiments. In Figure 3 simulated spontaneous AP and Ca2+ traces in steady state are reported, together with the main ionic currents and concentrations. Irel gating variables, together with the Ca2+ transients, are detailed in Figure 4. It shows that the Ca2+ release from SR is not directly dependent on the ICaL activation, but on the cytosolic Ca2+ concentration Cai, which rules the behavior of the RyRo, RyRc, and RyRa gating variables.

FIGURE 1
www.frontiersin.org

Figure 1. Comparison between the APs and Ca2+ transients simulated by the new Paci2018 model (solid blue) and the Paci2013 model (dashed red). (A) membrane potential. (B) Ca2+ transient.

FIGURE 2
www.frontiersin.org

Figure 2. Illustrative experimental Ca2+ transients from four cells (blue, cyan, green, magenta) and the Ca2+ transient simulated by the Paci2018 model (black). The y-axis of the experimental traces (originally ΔF/F) was normalized between the minimum and maximum values of the simulated trace for comparison. The inset contains only the decay part (from 90 to 10% of the transient amplitude) of the experimental and simulated transients reported in the main Figure.

FIGURE 3
www.frontiersin.org

Figure 3. Simulated spontaneous action potentials and ionic currents. (A,G) membrane potential. (B) Fast Na+ current (INa). (C) L-type Ca2+ current (ICaL). (D) Funny current (If) and Late Na+ current (INaL). (E) Transient outward K+ current (Ito) and Na+/K+ pump (INaK). (F) Rapid delayed (IKr), slow delayed (IKs) and inward (IK1) rectifier K+ currents. (H) Na+/Ca2+ exchanger (INaCa). (I) Release current from sarcoplasmic reticulum (Irel). (J) Na+ cytosolic concentration (Nai). (K) Cytosolic Ca2+ concentration (Cai). (L) Sarcoplasmic Ca2+ concentration (CaSR).

FIGURE 4
www.frontiersin.org

Figure 4. Details on the new formulation of the Ca2+ release from sarcoplasmic reticulum. (A) Time-course of the three sarcoplasmic release current gates during spontaneous action potentials. (B) Ca2+ transients during spontaneous action potentials.

Current Blocker Simulations

We compared in Table 3 the current blocker effects on the AP biomarkers (i) experimentally recorded by Ma et al. (2011), (ii) simulated by means of the Paci2013 model (Paci et al., 2013), and (iii) simulated by the Paci2018 model. Figure 5 shows that, despite the changes in the Irel formulation and parameter optimization, the behavior of the Paci2018 model is still consistent with the Paci2013 simulations (Paci et al., 2013) and the experimental data by Ma et al. (2011). TTX affected the upstroke phase, reducing the AP Vmax and delaying the AP Peak. NIFED, by acting on ICaL, reduced the Ca2+ influx through the cell membrane thus reducing the AP duration. Conversely, E4031 reduced the K+ efflux through the cell membrane, thus prolonging the AP. Finally, CHR showed little effect on the AP biomarkers, in agreement with the experiments.

TABLE 3
www.frontiersin.org

Table 3. Quantitative effects of current blockers on stimulated APs.

FIGURE 5
www.frontiersin.org

Figure 5. Simulation of current block effects on hiPSC-CMs paced at 1 Hz. (A) Tetrodotoxine blocks INa, slowing down the upstroke phase. (B) Nifedipine blocks ICaL, shortening APD and triangulating AP profile. (C) E4031 blocks selectively IKr, increasing APD. (D) IKs block by 3R4S-Chromanol 293B does not affect significantly the AP shape.

Delayed Afterdepolarizations

In Figure 6 we compared the behavior of the Paci2013 and the Paci2018 models in conditions of high extracellular Ca2+ concentration (Cao = 3.945 mM instead of 1.8 mM) and no external stimulation, i.e., both models produced spontaneous APs. As shown in Figure 6B, for t = [14, 30] s Irel underwent a small reactivation, slowing down the decay of the Ca2+ transients (Figure 6E), affecting the inward component of INaCa (Figure 6D) and triggering DADs in the membrane potential (Figures 6A,F). Around t = 30 s, Irel underwent a full reactivation, triggering an anticipated Ca2+ transient, which resulted in a full anticipated beat (Volders et al., 2000). We reported the extended AP and Ca2+ transient traces in the Supplementary Figure 1. Conversely, the old Paci2013 model did not show anomalies in the membrane potential (Figure 6G) or Ca2+ transients (Figure 6J) as consequence of the high extracellular Ca2+ concentration, while the Ca2+ in SR grew dramatically. The improvements introduced in the Paci2018 model allowed to simulate DADs not only in case of Ca2+ overload, but offered flexibility to simulate behaviors such as the spontaneous premature Ca2+ releases reported by Kim et al. (2015) in their Figures 2C, 4B and their Supplemental Figure S1B. We tested also the effect of a rate increment on the generation of DADs in case of hypercalcemia, by providing external pacing. In Supplementary Figure 2, we paced the Paci2018 model at its spontaneous basal rate (38 bpm, in solid blue) and at a basal rate increased by 50% (57 bpm, in dashed red). Despite the shorter diastolic depolarization phase, also at 57 bpm we obtained a fully developed DAD.

FIGURE 6
www.frontiersin.org

Figure 6. Comparison of the behavior of the Paci2018 (in blue, A–F) and the Paci2013 model (in red, G–K) in case of spontaneous APs and Ca2+ overload induced by hypercalcemia (Cao = 3.945 mM). In this condition both models can produce spontaneous APs, reported in (A,G). The Paci2018 model shows a saturation of the sarcoplasmic Ca2+ concentration (C), after which reactivation of Irel (B) induces a distortion in the shape of the Ca2+ transients (E) and in the inward component of INaCa (D), which corresponds to DADs in the membrane potential (A,F) in between t = [14, 30] s. At t = 30 s a strong spontaneous release of Ca2+ from SR is shown, which induces a spurious AP. On the contrary, the Paci2013 is not affected in these conditions and does not produce repolarization abnormalities (G).

Figure 7 shows the spontaneous APs, Ca2+ transients and INaCa, together with the release Ca2+ flux (Irel), characterized with premature Ca2+ releases. These traces were obtained with a normal extracellular Ca2+ concentration (Cao = 1.8 mM) but simulating more “immature” RyR machinery (namely, by shifting RyRo,half and RyRc,half by −0.002 and 0.002 mM respectively, doubling RyRo time constant and reducing to half of its nominal value RyRc time constant). In these conditions of more “immature” RyR machinery, we tested also the effects of INaCa block on DADs generation. In Supplementary Figure 3, we tested three INaCa block levels (50, 70, and 90%) while pacing the Paci2018 model at 1 Hz. With no INaCa block (solid blue trace), APs show small DADs comparable to those reported in Figure 7 in t = [10, 25] s. By blocking 50% of INaCa (dashed orange trace), we observed that the amplitude of DADs increases. Finally, for the higher block levels, we did not observe any DAD.

FIGURE 7
www.frontiersin.org

Figure 7. Repolarization abnormalities with standard extracellular Ca2+ concentrations. Such behavior was obtained with a control extracellular Ca2+ concentration Cao = 1.8 mM, shifting RyRo,half and RyRc,half by −0.002 and 0.002 mM respectively, doubling RyRo time constant and reducing to half of its nominal value RyRc time constant. These traces show the ability of the Paci2018 to simulate pathological conditions affecting the Ca2+ from SR. (A) Membrane potential. (B) INaCa. (C) Cytosolic Ca2+ concentration. (D) Irel. A similar morphology of the cytosolic Ca2+ time-course is reported in Kim et al. (2015) in their Figure 4B.

Comparison With Other Experimental Data

In order to compare the Paci2018 model with the experimental data, we challenged our model in the following conditions: (i) If block by ivabradine, (ii) hyperkalemia, (iii) hypocalcemia.

We simulated If block by 3 μM ivabradine as a 41% block of If maximum conductance (Yaniv et al., 2012; Koivumäki et al., 2018). As shown in Figure 8, If block induced only a slight reduction of the frequency of spontaneous APs (−2.3%), in agreement with (Kim et al., 2015), where 3 μM had virtually no effects on the spontaneous activity.

FIGURE 8
www.frontiersin.org

Figure 8. Effect of 3 μM ivabradine on the spontaneous APs. Administration of ivabradine (in dashed red) slightly slows down the rate of the spontaneous APs in comparison to control (in solid blue).

Kim et al. (2015) reported a slowdown of the spontaneous activity in conditions of hyperkalemia (from 4 to 8 mM the median frequency dropped from about 1.2 to 0.2 Hz) and a further increase in extracelluar K+ to 12 mM stopped the spontaneous activity. Our single cell model showed a similar trend (although with less sensitivity on K+): for Ko = 8 mM and Ko = 16 mM the Ca2+ transient spontaneous rate dropped by 11% and 15% respectively. An extracellular K+ concentration equal to 20 mM stopped the hiPSC-CM spontaneous activity as in Kim et al. (2015).

In the previous section, we tested the model in hypercalcemia conditions to simulate Ca2+ overload and the generation of DADs. We also challenged the model with an extracellular Ca2+ concentration of 0.1 mM (hypocalcemia) as in Kim et al. (2015): again we observed a trend qualitatively in agreement with the experimental data, since the rate of the spontaneous Ca2+ transients dropped by 17%.

Discussion

In this paper we present an updated version of the Paci2013 model of ventricular-like hiPSC-CM (Paci et al., 2013), here named Paci2018, developed by exploiting new Ca2+ transients data and an automatic algorithm for the parameter optimization.

The Paci2018 model aims to overcome an important limitation of the Paci2013 model, i.e., the simple formulation of the Ca2+ release from SR, Irel, initially presented by the TenTusscher2004 model of human adult ventricular cardiomyocyte (ten Tusscher et al., 2004). The peculiarity of such Irel is the direct link between the activation of ICaL and the activation of Irel (through the voltage-dependent activation gating variable d). This Irel formulation was then chosen for its simplicity instead of more complex formulations, however it presumes Irel to be directly dependent on the intracellular Ca2+ concentration and it prevents the model to be able to reproduce Ca2+-related anomalies such as DADs. As reported in Fink et al. (2011), models where the open probability of the RyR-sensitive channels is directly dependent on the opening of ICaL are unable to produce DADs. Indeed, DADs start from a repolarized membrane, at potentials where the ICaL gates are closed. Therefore, RyR-sensitive channels have to be able to open even when ICaL is zero, which is not possible in the old Irel formulation. Moreover, the presence of the crel constant in the TenTussher2004 Irel formulation could allow a Ca2+ flux also from an empty SR, which is clearly not possible.

We chose a quite simple Irel formulation, inspired by the Koivumäki2011 model of human adult atrial cardiomyocyte (Koivumäki et al., 2011), characterized by three gating variables. We did not choose more complex formulations, such as the Markov formulation in TenTusscher2006 model (ten Tusscher and Panfilov, 2006), since the latter was implemented mainly to correctly take into account the presence of T-tubules and the associated microdomains in the human adult ventricular cell. However, hiPSC-CMs produced so far with current differentiation protocols have not shown functional T-tubules (Li et al., 2013).

In order to keep the Paci2018 model the most consistent with its predecessor and with the Ma et al. dataset (Ma et al., 2011) of ionic currents and AP biomarkers, we did not change the formulation of INa, ICaL, If, Ito, IKr, IKs, and IK1, which were carefully fitted for the development of the Paci2013 model. In order to tune the Paci2018 model on the new data of Ca2+ transients, we decided to act only on those parameters defining the formulation of currents for which the Ma et al. (2011) dataset did not provide direct experimental data.

The result is a new model where spontaneous APs match the Ma et al. dataset and where Ca2+ transients are in agreement with the new Ca2+ dataset. In this paper, we did not aim to propose only modifications to a previous model in order to simulate a specific phenomenon of interest, but we wanted to develop a new cardiac cell model to substantially increase modeling accuracy. Therefore, it was fundamental to check the Paci2018 model capability to reproduce all the experimental data used for the Paci2013 model validation. In detail, we tested the Paci2018 model's responses to the same prototypical current blockers used to validate the Paci2013 model, showing that the behavior of the two models is consistent and in agreement with the experimental data. We further validated the model to simulate the experimentally observed phenomenon of DADs and other specific results obtained in different experimental conditions, being able to replicate all these experiments. In particular, the Paci2018 model can reproduce Ca2+-related abnormalities such as DAD-like anomalies that were not available in old formalism. We observed that the small or full Irel reactivations increase the cytosolic Ca2+ concentration and this produces an increased activity of the inward component of INaCa, which then translates in DAD-like abnormalities in the membrane potential. This is particularly clear in Supplementary Figure 3, where for high levels of INaCa block (70 and 90%) we do not observe DADs. It is interesting to note that the 50% block is not enough to cancel DADs: in fact, the inward INaCa component is still strong enough and further enhanced by the cytosolic Ca2+ accumulation, due to the very same INaCa block (minimum diastolic Ca2+: 0.028 vs. 0.017 μM; average diastolic Ca2+: 0.060 vs. 0.027 μM, in case of 50% INaCa block or no INaCa block, respectively). This results in even larger DADs than in the no block case. Similar DAD-like abnormalities in the electrical properties of hiPSC-CMs were also observed in mutant hiPSC lines, such as those derived from patients with HCM (Ojala et al., 2016) or CPVT (Kujala et al., 2012). Of note, in this paper we do not aim to simulate these mutations, which would deserve a specific analysis by their own, but to provide a model which can enable there in silico modeling. Notably, the membrane potential and cytosolic Ca2+ time-courses obtained by our model are very similar to the experimental ones, as seen by comparison of our Figures 6, 7 with Figures 2C, 4B and S1 in Kim et al. (2015). In particular, Figure 4B in Kim et al. (2015) shows the effect of isoproterenol, a non-selective β adrenoreceptor agonist, which first induces small delayed spontaneous Ca2+ releases and then triggers repetitive high rate firing. In our Figures 6A,E we show a very similar pattern in conditions of hypercalcemia (see Supplementary Figure 1 for an expanded version of Ca2+ and AP time-courses), where the spontaneous Ca2+ releases, which trigger DADs, finally cause a transition to a higher rate of Ca2+ transients and APs.

To further validate our model, we challenged it to simulate specific conditions such as hyperkalemia and hypocalcemia, to compare its behavior to published experimental data. The model showed qualitative agreement with such experiments, in spite of the fact that the spontaneous rate reduction is stronger in the experiments. However, such differences can emerge from (i) the multicellular embryoid body culturing and measurements, i.e., on multicellular ensembles, and especially (ii) as consequence of the high inter-cell line and inter-laboratory variability showed by hiPSC-CMs (Knollmann, 2013), especially in beating rate and MDP. In fact, the embryoid bodies examined by Kim et al. (2015) exhibited negligible IK1 and a much depolarized MDP (−59.1 ± 3.3 mV), which facilitate the spontaneous electrical activity, as shown also by the high frequency of the spontaneous Ca2+ transients (1 Hz or higher) in control conditions. On the contrary, the experimental datasets we used in this paper showed slower spontaneous activity and more hyperpolarized MDP, as simulated by our model and reported in Table 2 (MDP = −75.8 mV, FREQ = 0.65 Hz).

We tested also the capability of the Paci2018 model to simulate early afterdepolarizations (EADs) in conditions of IKr block: the protocol and the results of the EAD simulations were reported in section 3 of the Supplementary Material. We observed that the Paci2018, as the old Paci2013 model, does not produce EADs. This is not surprising, since we designed the Paci2018 model carefully maintaining the same formulation for the ion currents that we fit on the experimental data by Ma et al. (2011), and that we used for the Paci2013 model (Paci et al., 2013). However, in our previous study (Paci et al., 2016), we demonstrated that IKr block triggered EADs in a population of in silico hiPSC-CM based on the Paci2013 model. Replicating here that study is out of the scope of this paper, however we decided to test if two parameter sets, which in Paci et al. (2016) triggered EADs, and one parameter set, which triggered repolarization failure, could induce repolarization abnormalities also in the Paci2018 model. Supplementary Figure 4 shows that, as expected, two derived models produced EADs and the third one repolarization failure. In particular, as reported in the Supplementary Table 1, these parameter sets showed (compared to baseline): (i) greater ICaL (which corresponds to a greater Ca2+ influx during the repolarization phase in case of ICaL reactivation); (ii) greater INaCa (which is inward during repolarization, thus promoting abnormalities in conditions of compromised repolarization as in case of IKr block), and (iii) smaller IK1 (corresponding to a smaller contribution to the membrane potential stabilization). It is also interesting to notice that the three parameter sets showed greater IKr compared to baseline: this suggests that in models where IKr is highly expressed, a 90% block can have more dramatic effects than in cells expressing smaller IKr. This small test suggests, as we reported among the limitations, that the use of more refined modeling approaches, e.g., populations of in silico models, is advisable to simulate specific phenomena (in this case EADs as consequence of IKr block) in cell types characterized by high electrophysiological variability. In the perspective of offering researchers a tool to investigate the mechanisms of ventricular arrhythmia development in hiPSC-CMs, we ran also two preliminary tests about the occurrence of alternans or alternans-like patterns in the Paci2018 model (Methods and Results are detailed in section 4 of the Supplementary Material). When pacing the Paci2018 model at 120 bpm (2 Hz), the membrane potential exhibited an alternans-like pattern, characterized by a full AP followed by a smaller one (Supplementary Figure 5A). We tested also the effect of very high pacing rate (200 bpm, 3.33 Hz) in ischemia-like conditions (Supplementary Figure 5B): in this case, the membrane potential was strongly depolarized (−60 mV) and exhibited 2:1 alternans. We consider these results valuable also for more advanced studies regarding arrhythmias propagation in case of extension of the simulations to a monodimensional strand or a bidimensional patch.

Overall, we believe that the here presented Paci2018 model substantially increases hiPSC-CM modeling accuracy. Of course, as every in silico model, also the Paci2018 is an approximated description of a real system: consequently, new mechanisms can be explored and improved. This is especially true for such a complex cellular system as hiPSC-CM, which is not fully characterized yet. The main limitation that affects the Paci2018 model, as well the other single cell in silico models available in literature, is the fact that the simulated behavior belongs to a theoretical cell, built on averaged experimental data from many cells. In spite of our efforts to build the Paci2018 in a rigorous way, we do not have the presumption to consider it representative of all the hiPSC-CMs available nowadays, both from commercial lines and laboratory-specific lines. This is particularly true in case of an in vitro model like hiPSC-CMs, since the same experiments that we need to build an in silico model are extremely variable, as reported by multiple studies (Knollmann, 2013; Lu et al., 2015; Paci et al., 2017). As an additional example, in Lu et al. (2015) the authors succeeded in culturing iCell hiPSC-CMs that showed limited biomarker variability. Nevertheless, their cells showed way slower spontaneous beating rate than the cells we used to build our model (~15 bpm vs. ~38 bpm), despite both datasets were recorded at the same 37°C temperature. In order to have a better translation from in silico to in vitro models and finally to (pre)clinical applications, more refined modeling techniques exist, e.g., the population of in silico models (Britton et al., 2013; Passini et al., 2017). In addition to this, in terms of translation of drug test results from hiPSC-CMs to clinical applications, more challenges emerged recently. For example, in Abi-Gerges et al. (2017), testing a panel of 30 drugs on non-paced hiPSC-CMs provided on one hand good predictivity of pro-arrhythmic events, but on the other hand limited predictivity as an early QT screening. Therefore, in spite the high potential of hiPSC-CMs as in vitro models for diseases and drug tests, direct translation of results is still a work in progress, as indicated by the fact that advances in phenotype selection, cell maturation and combined recording platforms are currently being made. One inherent limitation of the Paci2018 model consists in the use of two different datasets of experimental data, one for the AP and ionic current data and the other for Ca2+ data. This is clear, e.g., by comparing the mean values of APD90 and Ca2+ transient DURATION biomarkers: 414.7 and 804.5 ms, respectively. It is known that high affinity Ca2+ indicators, such as Fluo-4, can artificially prolong the Ca2+ transients and that Ca2+ transients last more than APs (Lee et al., 2012); however, the aforementioned Ca2+ transient DURATION value is substantially longer than APD90. Therefore, we tried to find an acceptable tradeoff by means of the parameter optimization algorithm, in order to fulfill all the experimental ranges of both the AP and Ca2+ transient biomarkers. In particular, in the cost function we chose heavier weights (wi = 2) only for two critical biomarkers which characterize hiPSC-CMs, i.e., MDP (depolarized compared to adult cardiomyocytes) and the rate of the spontaneous APs. For all the other biomarkers, the weights were set to one. The tradeoff obtained by the parameter optimization produced a model that simulates Ca2+ transients whose DURATION is 617.8 ms. Despite this value lies within the experimental range 804.5 ± 188.0 ms, it is very close to its lower bound. Conversely, the DT9010 biomarker properly lies within its experimental range, although slightly smaller than the experimental mean DT9010 (409.8 ± 100.1 vs. 367.6 ms). The inset included in Figure 2 shows only the decay part (from 90 to 10% of the transient amplitude) of the experimental and simulated transients reported in the main Figure 2: our simulated decay lies within the illustrative experimental traces, as well as DT9010 of the simulated transient is included in the experimental interval. Of note, the experimental traces reported in Figure 2 are purely illustrative of our experiments: the global quantitative comparison between experiments and simulations is done in Table 2 in terms of biomarkers. It is also known that temperature affects the spontaneous electrical activity of hiPSC-CMs, e.g., in Laurila et al. (2016) the authors showed how the temperature increments and reductions accelerate or slow down hiPSC-CM spontaneous beating, respectively. The fact that the AP and Ca2+ transient experimental datasets, chosen to develop the Paci2018 model, were recorded at the same temperature (35–37°C) further mitigates the differences between the two datasets. Another potential limitation of our model is related to its simple structure, consisting of two compartments only (cytosol and SR). Such compartmentalization does not allow to simulate a spatial inhomogeneous Ca2+ distribution, as conversely recently done by others (Koivumäki et al., 2018). We preferred to keep the model formulation simple, without entering in such detailed description, since it would have required a more complex mathematical formulation, a more detailed knowledge of the subcellular organization of hiPSC-CMs and, consequently, a more challenging parameter identification and higher simulation time. However, our formulation is successful in simulating the most important electrophysiological mechanisms at the whole cell level, including DADs elicited from different kind of stresses, without resorting to random RyR openings, as done in Koivumäki et al. (2018). A final limitation of the Paci2018 model is that it describes only the electrophysiology and the Ca2+ handling, but it does not take into consideration contractility. In literature, a few models of myofilament are available (Negroni and Lascano, 2008; Rice et al., 2008; Negroni et al., 2015). For instance, in Negroni et al. (2015) the authors presented an integrated rabbit ventricular cardiomyocyte model, which included also the mathematical description of force generation and the myofilament Ca2+ kinetics. However, the sarcomeric structure of hiPSC-CMs is immature and the myofibrils in hiPSC-CMs are oriented in multiple directions within the cell: in Bedada et al. (2016) a higher level of sarcomeric organization was reached only by repopulating hiPSC-CMs into a biological cardiac matrix. We acknowledge the potential of an integration of the myofilament description in the Paci2018, nevertheless the translation from a different species model of, e.g., rabbit cardiomyoycte to a hiPSC-CM model is not elementary and out of the scope of this paper. Future developments of this work will include using the Paci2018 model to ran in silico drug tests and to ran monodimensional and bidimensional simulations, as done, e.g., by Raphel et al. with the Paci2013 model (Raphel et al., 2017).

In conclusion, in this work we present an updated and more versatile version of our hiPSC-CM in silico model, based on a new dataset of electrophysiological data. Our model can represents the basis for new in silico studies on the effects of drugs or mutations (e.g., CPVT), which affect the Ca2+ handling in hiPSC-CMs. Due to its relatively light formulation (23 ordinary differential equations), our model is suitable also for very large studies on in silico populations, e.g., to support screening of different drug/compounds at various concentrations.

Author Contributions

All the authors conceived and designed the study. R-PP and KP performed the in vitro measurements and analyzed the in vitro data. MP and DC performed the in silico model development and validation. MP and SS analyzed the in silico data, prepared the figures and drafted the manuscript. All the authors interpreted the results and revised the manuscript.

Funding

MP was financially supported by the Academy of Finland (project CardSiPop, decision number 307967). R-PP was financially supported by Paavo Nurmi Foundation. The Heart Group has been supported by the Finnish Cardiovascular Foundation, the Academy of Finland, TEKES—the Finnish Funding Agency for Innovation and Technology, and Pirkanmaa Hospital District.

Conflict of Interest Statement

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.

The reviewer JB and handling Editor declared their shared affiliation.

Acknowledgments

The authors acknowledge CSC—IT Center for Science, Finland, for generous computational resources, Tampere Facility of Electrophysiological Measurements for their service and Dr. Jussi Koivumäki for sharing his code.

Supplementary Material

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

References

Abi-Gerges, N., Pointon, A., Oldman, K. L., Brown, M. R., Pilling, M. A., Sefton, C. E., et al. (2017). Assessment of extracellular field potential and Ca2+ transient signals for early QT/pro-arrhythmia detection using human induced pluripotent stem cell-derived cardiomyocytes. J. Pharmacol. Toxicol. Methods 83, 1–15. doi: 10.1016/j.vascn.2016.09.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Ahola, A., Pölönen, R. P., Aalto-Setälä, K., and Hyttinen, J. (2018). Simultaneous measurement of contraction and calcium transients in stem cell derived cardiomyocytes. Ann. Biomed. Eng. 46, 148–158. doi: 10.1007/s10439-017-1933-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Bedada, F. B., Wheelwright, M., and Metzger, J. M. (2016). Maturation status of sarcomere structure and function in human iPSC-derived cardiac myocytes. Biochim. Biophys. Acta 1863, 1829–1838. doi: 10.1016/j.bbamcr.2015.11.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Björk, S., Ojala, E. A., Nordström, T., Ahola, A., Liljeström, M., Hyttinen, J., et al. (2017). Evaluation of optogenetic electrophysiology tools in human stem cell-derived cardiomyocytes. Front. Physiol. 8:884. doi: 10.3389/fphys.2017.00884

PubMed Abstract | CrossRef Full Text | Google Scholar

Britton, O. J., Bueno-Orovio, A., Van Ammel, K., Lu, H. R., Towart, R., Gallacher, D. J., et al. (2013). Experimentally calibrated population of models predicts and explains intersubject variability in cardiac cellular electrophysiology. Proc. Nat. Aca. Sci. U.S.A. 110, E2098–E2105. doi: 10.1073/pnas.1304382110

PubMed Abstract | CrossRef Full Text | Google Scholar

Entcheva, E., and Bub, G. (2016). All-optical control of cardiac excitation: combined high-resolution optogenetic actuation and optical mapping. J. Physiol. 9, 2503–2510. doi: 10.1113/JP271559

CrossRef Full Text | Google Scholar

Fabbri, A., Fantini, M., Wilders, R., and Severi, S. (2017). Computational analysis of the human sinus node action potential: model development and effects of mutations. J. Physiol. 7, 2365–2396. doi: 10.1113/JP273259

CrossRef Full Text | Google Scholar

Fink, M., Noble, P. J., Noble, D., Fink, M., Noble, P. J., and Noble, D. (2011). Ca2+-induced delayed afterdepolarizations are triggered by dyadic subspace Ca+ affirming that increasing SERCA reduces aftercontractions. Cardiac Excit. Contract. 301, 921–935. doi: 10.1152/ajpheart.01055.2010

PubMed Abstract | CrossRef Full Text

Gerlach, A. C., Stoehr, S. J., and Castle, N. A. (2010). pharmacological removal of human ether-à-go-go-related gene potassium channel inactivation by 3-nitro-N-(4-phenoxyphenyl) benzamide (ICA-105574). Mol. Pharmacol. 77, 58–68. doi: 10.1124/mol.109.059543

PubMed Abstract | CrossRef Full Text

Kim, J. J., Yang, L., Lin, B., Zhu, X., Sun, B., Kaplan, A. D., et al. (2015). Mechanism of automaticity in cardiomyocytes derived from human induced pluripotent stem cells. J. Mol. Cell. Cardiol. 81, 81–93. doi: 10.1016/j.yjmcc.2015.01.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Klimas, A., Ambrosi, C. M., Yu, J., Williams, J. C., Bien, H., and Entcheva, E. (2016). OptoDyCE as an automated system for high-throughput all-optical dynamic cardiac electrophysiology. Nat. Commun. 7:11542. doi: 10.1038/ncomms11542

PubMed Abstract | CrossRef Full Text | Google Scholar

Knollmann, B. C. (2013). Induced pluripotent stem cell-derived cardiomyocytes: boutique science or valuable arrhythmia model? Circ. Res. 112, 969–976. doi: 10.1161/CIRCRESAHA.112.300567

PubMed Abstract | CrossRef Full Text | Google Scholar

Koivumäki, J. T., Korhonen, T., and Tavi, P. (2011). Impact of sarcoplasmic reticulum calcium release on calcium dynamics and action potential morphology in human atrial myocytes: a computational study. PLoS Comput. Biol. 7:e1001067. doi: 10.1371/journal.pcbi.1001067

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

Kujala, K., Paavola, J., Lahti, A., Larsson, K., Pekkanen-Mattila, M., Viitasalo, M., et al. (2012). Cell model of catecholaminergic polymorphic ventricular tachycardia reveals early and delayed afterdepolarizations. PLoS ONE 7:e44660. doi: 10.1371/journal.pone.0044660

PubMed Abstract | CrossRef Full Text | Google Scholar

Lahti, A. L., Kujala, V. J., Chapman, H., Koivisto, A.-P., Pekkanen-Mattila, M., Kerkela, E., et al. (2012). Model for long QT syndrome type 2 using human iPS cells demonstrates arrhythmogenic characteristics in cell culture. Dis. Model. Mech. 220–230. doi: 10.1242/dmm.008409

PubMed Abstract | CrossRef Full Text

Laurila, E., Ahola, A., Hyttinen, J., and Aalto-Setälä, K. (2016). Methods for in vitro functional analysis of ipsc derived cardiomyocytes—special focus on analyzing the mechanical beating behavior. Biochim. Biophys. Acta 1863, 1864–1872. doi: 10.1016/j.bbamcr.2015.12.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, P., Klos, M., Bollensdorff, C., Hou, L., Ewart, P., Kamp, T. J., et al. (2012). Simultaneous voltage and calcium mapping of genetically purified human induced pluripotent stem cell-derived cardiac myocyte monolayers. Circ. Res. 110, 1556–1563. doi: 10.1161/CIRCRESAHA.111.262535

PubMed Abstract | CrossRef Full Text

Lei, C. L., Wang, K., Clerx, M., Johnstone, R. H., Hortigon-Vinagre, M. P., Zamora, V., et al. (2017). Tailoring mathematical models to stem-cell derived cardiomyocyte lines can improve predictions of drug-induced changes to their electrophysiology. Front. Physiol. 8:986. doi: 10.3389/fphys.2017.00986

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

Lu, H. R., Whittaker, R., Price, J. H., Vega, R., Pfeiffer, E. R., Cerignoli, F., et al. (2015). High throughput measurement of Ca++dynamics in human stem cell-derived cardiomyocytes by kinetic image cytometery: a cardiac risk assessment characterization using a large panel of cardioactive and inactive compounds. Toxicol. Sci. 148, 503–516. doi: 10.1093/toxsci/kfv201

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, D., Wei, H., Zhao, Y., Lu, J., Li, G., Binte, N., et al. (2013). Modeling type 3 Long QT syndrome with cardiomyocytes derived from patient-specific induced pluripotent stem cells. Int. J. Cardiol. 168, 5277–5286. doi: 10.1016/j.ijcard.2013.08.015

PubMed Abstract | 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, H2006–H2017. 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. doi: 10.1056/NEJMoa0908679

PubMed Abstract | CrossRef Full Text | Google Scholar

Negroni, J. A., and Lascano, E. C. (2008). Simulation of steady state and transient cardiac muscle response experiments with a huxley-based contraction model. J. Mol. Cell. Cardiol. 45, 300–312. doi: 10.1016/j.yjmcc.2008.04.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Negroni, J. A., Morotti, S., Lascano, E. C., Gomes, A. V., Grandi, E., Puglisi, J. L., et al. (2015). β-Adrenergic effects on cardiac myofilaments and contraction in an integrated rabbit ventricular myocyte model. J. Mol. Cell. Cardiol. 81, 162–175. doi: 10.1016/j.yjmcc.2015.02.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Ojala, M., Prajapati, C., Pekka Pölönen, R., Rajala, K., Pekkanen-Mattila, M., Rasku, J., et al. (2016). Mutation-Specific Phenotypes in hiPSC-Derived Cardiomyocytes Carrying Either Myosin-Binding Protein C or α-Tropomyosin Mutation for Hypertrophic Cardiomyopathy. Stem Cell Int. 2016:1684792. doi: 10.1155/2016/1684792

CrossRef Full Text | 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. doi: 10.1007/s10439-013-0833-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Paci, M., Hyttinen, J., Rodriguez, B., and Severi, S. (2015). Human induced pluripotent stem cell-derived versus adult cardiomyocytes: an in silico electrophysiological study on effects of ionic current block. Br. J. Pharmacol. 172, 5147–5160. doi: 10.1111/bph.13282

PubMed Abstract | CrossRef Full Text | Google Scholar

Paci, M., Passini, E., Severi, S., Hyttinen, J., and Rodriguez, B. (2016). A population of in silico models to face the variability of human induced pluripotent stem cell-derived cardiomyocytes: the hERG block case study. Comput. Cardiol. 43, 1189–1192. doi: 10.23919/CIC.2016.7868961

CrossRef Full Text | Google Scholar

Paci, M., Passini, E., Severi, S., Hyttinen, J., and Rodriguez, B. (2017). Phenotypic Variability in LQT3 Human induced pluripotent stem cell-derived cardiomyocytes and their response to antiarrhythmic pharmacologic therapy: an in silico approach. Heart Rhythm 14, 1704–1712. doi: 10.1016/j.hrthm.2017.07.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Paci, M., Sartiani, L., Lungo, M., Jaconi, M., Mugelli, A., Cerbai, E., et al. (2012). Mathematical modelling of the action potential of human embryonic stem cell derived cardiomyocytes. Biomed. Eng. Online 11:61. doi: 10.1186/1475-925X-11-61

PubMed Abstract | CrossRef Full Text | Google Scholar

Passini, E., Britton, O. J., Lu, H. R., Rohrbacher, J., Hermans, A. N., Gallacher, D. J., et al. (2017). Human in silico drug trials demonstrate higher accuracy than animal models in predicting clinical pro-arrhythmic cardiotoxicity. Front. Physiol. 8:668. doi: 10.3389/fphys.2017.00668

PubMed Abstract | CrossRef Full Text | Google Scholar

Rajamohan, D., Kalra, S., Hoang, M. D., George, V., Staniforth, A., Russell, H., et al. (2016). Automated electrophysiological and pharmacological evaluation of human pluripotent stem cell-derived cardiomyocytes. Stem Cells Dev. 25, 439–452. doi: 10.1089/scd.2015.0253

PubMed Abstract | CrossRef Full Text | Google Scholar

Raphel, F., Boulakia, M., Zemzemi, N., Coudiere, Y., Guillon, J.-M., Zitoun, P., et al. (2017). identification of ion currents components generating field potential recorded in MEA from hiPSC-CM. IEEE Trans. Biomed. Eng. 9294, 1–1. doi: 10.1109/TBME.2017.2748798

CrossRef Full Text | Google Scholar

Rice, J. J., Wang, F., Bers, D. M., and De Tombe, P. P. (2008). Approximate model of cooperative activation and crossbridge cycling in cardiac muscle using ordinary differential equations. Biophys. J. 95, 2368–2390. doi: 10.1529/biophysj.107.119487

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanguinetti, M. C., and Jurkiewicz, N. K. (1990). Two Components of cardiac delayed rectifier K+ current. differential sensitivity to block by class III antiarrhythmic agents. J. Gen. Physiol. 96, 195–215. doi: 10.1085/jgp.96.1.195

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

ten Tusscher, K. H., and Panfilov, A. V. (2006). Alternans and spiral breakup in a human ventricular tissue model. Am. J. Physiol. Heart Circ. Physiol. 291, H1088–H1100. doi: 10.1152/ajpheart.00109.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Volders, P. G., Vos, M. A., Szabo, B., Sipido, K. R., de Groot, S. H., Gorgels, A. P., et al. (2000). progress in the understanding of cardiac early afterdepolarizations and torsades de pointes: time to revise current concepts. Cardiovas. Res. 46, 376–92. doi: 10.1016/S0008-6363(00)00022-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Walker, C. A., and Spinale, F. G. (1999). The structure and function of the cardiac myocyte: a review of fundamental concepts. J. Thorac. Cardiovasc. Surg. 118, 375–382. doi: 10.1016/S0022-5223(99)70233-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Yaniv, Y., Maltsev, V. A., Ziman, B. D., Spurgeon, H. A., and Lakatta, E. G. (2012). The ‘funny’ current (If) inhibition by ivabradine at membrane potentials encompassing spontaneous depolarization in pacemaker cells. Molecules 17, 8241–8254. doi: 10.3390/molecules17078241

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: human induced pluripotent stem cell-derived cardiomyocyte, action potential, calcium transient, computer simulation, in silico modeling

Citation: Paci M, Pölönen R-P, Cori D, Penttinen K, Aalto-Setälä K, Severi S and Hyttinen J (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

Received: 16 March 2018; Accepted: 22 May 2018;
Published: 26 June 2018.

Edited by:

Jichao Zhao, University of Auckland, New Zealand

Reviewed by:

Jieyun Bai, University of Auckland, New Zealand
Michelle M. Monasky, Policlinico San Donato (IRCCS), Italy
Anuj Agarwal, Signal Solutions LLC, Fairfax Virginia, United States

Copyright © 2018 Paci, Pölönen, Cori, Penttinen, Aalto-Setälä, Severi and Hyttinen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Michelangelo Paci, michelangelo.paci@tut.fi

These authors have contributed equally to this work.