Sodium, calcium, and lumped outward currents were measured in Cor.4U cells in the whole-cell patch clamp configuration using the Nanion SyncroPatch 96 platform (Nanion Technologies GmbH, Germany). Sodium and lumped outward currents were measured using an intracellular solution containing (in mM) 50 KCl, 60 KF, 10 NaCl, 10 HEPES, and 20 EGTA (pH: 7.2), and a bath solution containing (in mM) 150 NaCl, 4 KCl, 1 MgCl₂, 1.2 CaCl₂, 10 HEPES, and 5 glucose (pH: 7.4). Calcium current recordings were made using an intracellular solution containing (in mM) 50 CsCl, 60 CsF, 10 TEA-Cl (a potassium current blocker), 5 HEPES, 10 EGTA, 4 Na2-ATP, 0.1 Na-GTP, and 0.1 cAMP (pH: 7.2) and a bath solution containing (in mM) 130 NMDG, 10 BaCl₂, 4 CsCl₂, 1 MgCl₂, 2 CaCl₂, 10 HEPES, and 5 glucose (pH: 7.4). All currents were recorded at room temperature.

For the sodium current measurements, cells were held at −80 mV and then stepped to potentials ranging from −60 to 60 mV with 10 mV increments, before returning to the holding potential. The step duration was 20 ms and the interval between steps was 5 s. The calcium current experiments used a similar protocol, but with 200 ms steps from −40 to 40 mV. Outward current was measured with 500 ms steps from −40 to 50 mV, with a 10 s interval between steps. All three protocols are shown in Supplementary Figure S1.

For the outward current experiments, we fitted directly to the experimental current traces (see section 2.5), and so leak correction was applied using I_leak = V/R_leak where R_leak was the leak resistance estimated at the holding potential. Capacitance artifacts were filtered out by omitting the first 10 ms after each change in potential (see e.g., Ogden and Stanfield, 1994).

Action potentials in iCell iPSC-CMs (CDI, USA) plated on coverslips were measured in whole-cell patch clamp configuration using a HEKA amplifier (EPC 10 USB Triple, HEKA Elektronik, Germany). Recordings were made using a pipette solution containing (in mM) 10 NaCl, 125 KCl, 1 MgCl₂, 10 HEPES, 0.1 Na₃GTP, 5 Mg-ATP, 5 EGTA (pH 7.2) and a bath solution containing (in mM) 150 NaCl, 4 KCl, 1.2 CaCl₂, 1 MgCl₂, 10 HEPES (pH 7.4). Cells were stimulated at a frequency of 1.0 Hz, for at least 50 cycles before recording.

Action potentials were recorded from Cor.4U cultures with optical mapping using the CellOPTIQ electrophysiology platform (Clyde Biosciences Ltd). Cells were incubated in serum-free media at 35 ± 2°C, and transiently loaded with voltage sensitive fluorescent dye di-4-ANEPPS (20 μL of stock solution 27 mM in ethanol; University of Connecticut Health Center). The loaded dye was then excited with a peak wavelength 470 nm LED, and the emitted fluorescence from the Cor.4U iPSC-CMs was recorded at a sample frequency of 10 kHz. Measurements were performed before and after addition of Dofetilide, Quinidine, Sotalol and Verapamil at the concentrations shown in Table 1. Paracetamol was applied as a negative control.

A semi-automatic data analysis method based in Wang et al. (2015) was employed to normalize the data. In short, heuristics were used to form an initial estimate of the start and end time of the AP. The region just before the estimated upstroke was used to determine V_normalized = 0, while the 95th percentile of the data during the (estimated) AP was used as V_normalized = 1. We then calculated the final APD₉₀ and APD₅₀ from this normalized signal.

Simulations of the patch clamp protocols were carried out using the model by Paci et al. (2013). Initial intracellular and extracellular ion concentrations were set to the values used in the experiments. For the I_Na, I_CaL, and I_outward voltage clamp experiments, concentrations were clamped (corresponding to the buffering effects of the pipette), but for AP simulations concentrations were allowed to vary following model equations. The temperature parameter in the model, which affects reversal potentials as well as I_CaL permeability and I_Kr, I_NaK, and I_NaCa kinetics, was set to 25°C (298 K) to match the experimental temperature. Simulations were run using Myokit (Clerx et al., 2016), with CVODE (Hindmarsh et al., 2005) set to the default tolerance settings of abs_tol = 10⁻⁶ and rel_tol = 10⁻⁴. Model code was imported from a CellML (Cuellar et al., 2003) file downloaded from the Physiome model repository (Yu et al., 2011). Numerical integration was carried out using NumPy/SciPy (Jones et al., 2001). All codes and data are freely available from https://gitlab.com/MichaelClerx/tailored-ipsc-models.

The maximum conductance of I_Na was estimated by scaling the I_Na conductance in the Paci model to match the peak current recorded experimentally with the sodium protocol (n = 35 cells), based on the assumption that the peak current is composed of I_Na alone. We tested this assumption by running a simulated experiment, where we observed that I_Na alone would reach 1.01 × the initial inward deflection after each voltage step. So the peak is almost entirely due to sodium and only decreased slightly by the presence of other currents. Similarly, the recordings made with the calcium protocol (n = 25 cells) were used to directly infer the maximum conductance of I_CaL.

To estimate the conductances of the remaining major currents, we used the recordings made with the outward-current protocol (n = 22 cells). Using the iPSC model by Paci et al. (2013) we simulated the response to this protocol of I_Na, I_CaL, I_K1, I_Kr, I_Ks, I_to, I_f, and I_NaCa (see Supplementary Figure S2). We then tried to find a weighted sum of these simulated currents that could replicate the measured signal. This was done by minimizing the sum of square errors between measured and simulated current during the voltage steps, using the optimization method CMA-ES (Hansen, 2006). The procedure was repeated for each of the 22 measured cells, resulting in a unique set of scaling factors per cell.

While this protocol was intended to find values for only the outward currents (such as I_Kr, I_Ks, and I_to) we chose to vary the inward currents I_Na and I_CaL in the optimization to reduce the risk that any inward currents in the signal would erroneously be attributed to the outward currents (see e.g., Sarkar and Sobie, 2010). Note that we do not use these fitted I_Na and I_CaL conductances because we fit these from other dedicated experiments; they were included here just to yield more accurate outward current fits. After finding that the most prominent outward currents were I_Ks and I_NaCa (see section 3.2), we ran a second optimization with only I_Ks and I_NaCa: little change was observed, but we expect that fitting outward currents together with I_Na and I_CaL is likely to yield slightly more accurate results.

Convergence of the optimization results was verified by repeating the process 10 times, using different random seeds for each run. We found that the L² norm of the difference between the first and repeated scaling factor vectors was smaller than 10⁻⁵ for all 10 random starting points. To further verify the identifiability of the problem, we performed the same analysis on synthetic data (with synthetic noise), and were able to successfully infer the conductance scaling factors (see Supplementary Figure S4 and Supplementary Table S1). We note that this analysis assumes the kinetics of the currents have low model discrepancy, i.e., reflect the kinetics of the real currents well.

Finally, to quantify the contribution of each current to the total outward current, we defined a contribution score c_i for each current I_i as:

where I_{final, i} was defined as the current measured at the end of the final step of the outward-current protocol. This measure simply gives us a sense of the proportion of outward current that is contributed by each individual component during the end of the 50 mV step, but it is not used in tailoring the models.

Next, the estimated maximum conductances were used to tailor the Paci et al. (2013) model to individual cells from the Cor.4U cell line. A total of 22 model variants were parameterized, corresponding to the 22 cells for which the outward current was measured. Since we found many currents were not discernible in the recorded outward current (see Table 2), we only applied the cell-specific scaling factors for I_Ks and I_NaCa. All tailored models used the same I_Na and I_CaL scaling factors, found in the inward current experiments which were measured in different cells and hence any covariance could not be accounted for. The remaining currents were left unchanged, as they are necessary for other cellular behavior, such as homeostasis, even though they might not contribute strongly to the recorded outward current. Note that we have used only linear scaling of the conductances, and the current kinetics of the original currents were not altered. Finally, Na⁺, K⁺, and Ca²⁺ evolve in time according to the Paci et al. (2013) model, to mimic the intact cell conditions of our optical mapping experiments.

These tailored models were then used to simulate baseline APs, as well as APs with drug perturbation. The effects of drugs on ion current maximum conductances were modeled using the Hill equation (Hill, 1910; Weiss, 1997).

where x denotes the concentration of the applied drug, IC50 is the inhibitory concentration 50% value, h is the Hill coefficient, and f(x) is a scaling factor for the maximum conductance that varies from 0 (full block) to 1 (no block).

For each cell and each drug, a model was created where the maximum conductances of the ion currents were scaled according to Equation (2) using the IC50 values from Table 1 and a Hill coefficient of 1.0. For comparison, the same scaling factors were applied to an original model with the untailored conductance values from Paci et al. (2013).

In our optical mapping experiments, cells formed a spontaneously-beating and electrotonically-coupled monolayer. However, in this preparation not all cells beat at their spontaneous rates. Most cells will fire an AP when triggered by an activation wave from their neighbors rather than spontaneously, and a relatively small region of (by definition) faster spontaneously-beating cells sets the pacing rate for the entire monolayer. Therefore, to mimic this effect, we paced the cells at the mean rates observed in the optical mapping experiments, for a given compound, to account for any AP rate dependency. We used the cycle lengths of 1.375 s for Dofetilide, 1.176 s for Quinidine, 0.933 s for Sotalol, 0.905 s for Verapamil and 1.0 s for Paracetamol and all other experiments. To allow direct comparison with the optical mapping data, the simulated AP was normalized using the same algorithm (see section 2.3).

Figure 1 shows the peak current-voltage relationships for I_Na and I_CaL, measured in Cor.4U cells. The mean peak current in 35 cells (I_Na) and 25 cells (I_CaL) is plotted, as are the 25th and 75th percentiles. Compared to the prediction of the original Paci et al. (2013) model (created from iCell iPSC-CM data), the experimental data show a lower amplitude of both currents in Cor.4U cells. We had to scale by a factor 0.69 to match the mean peak I_Na, and 0.80 to match the mean peak I_CaL recordings. The simulated I_Na peaked at the same potential as the experimental data, suggesting the activation kinetics of I_Na in iCell and Cor.4U cell lines are similar. The simulated I_CaL kinetics followed the Paci et al. (2013) model, and were left-shifted relative to experimental data. Further experiments established that this shift was due to a right-shift in the experimental IV curve due to Ba²⁺ being present in the I_CaL voltage clamp experiment bath solution (see Supplementary Figure S3), hence we do not adjust the kinetic terms and tailor only the maximum conductance.

Figure 2 (left panel) shows the current measured with the outward-current protocol in a single Cor.4U cell. To analyse the composition of this current, we simulated the same protocol, and looked for a sum of scaled transmembrane currents from the Paci et al. (2013) model that gave a similar result (see section 2.5). We repeated this process for each of the 22 cells with measured outward current, and obtained the scaling factors s for each cell and current shown in Table 2. Note that the scaling factors s are relative to the original Paci et al. (2013) model. In many cases, the optimization routine indicated that the kinetic profile of certain currents was not discernible in the measured outward current. This is indicated in the table with a dash (—) for any scaling factor smaller than 10^-10. After seeing these results, as a comparison, we also tried fitting by varying only I_Ks and I_NaCa (and using the scaling factors for I_Na and I_CaL determined previously), and the results are similar (see Supplementary Table S2).

For most cells, we found that the measured responses differed greatly, in both the shape and size of the currents, from the original model predictions (see Supplementary Figure S4), leading to a poor quality of fit (see Supplementary Figure S5). As a result, the best reconstructions of the simulated current relied almost entirely on a greatly amplified I_Ks current, along with strong I_NaCa, while other currents such as I_Kr and I_f were notably absent. Based on this, we might assume that I_Ks and I_NaCa are more strongly expressed in Cor.4U cells than in the iCell cells the Paci et al. (2013) model was based on. As an initial verification of these findings, we repeated some outward current measurements in the presence of Chromanol (an I_Ks blocker), see Supplementary Figure S6 for an example where I_Ks is indeed significant. The near-zero contributions of other currents (e.g., I_Kr) does not imply that these currents are completely absent in Cor.4U cells, but instead suggests that the currents as simulated from the model could not be found in our recordings using the specified patch clamp protocol. This is a strong hint that changes to the kinetics of the currents will be required to accurately simulate the ion currents in Cor.4U cells at this temperature using the model by Paci et al. (2013). Such a mismatch in kinetics would also explain the large remaining errors between measurements and fit seen in Supplementary Figure S5, causing other currents, such as I_Kr, to be fitted as absent. This is discussed further in section 4.4.

We then created tailored models by modifying the original Paci et al. (2013) model in two ways: First, we scaled the maximum conductances of I_Na and I_CaL by a factor 0.69 and 0.8 respectively, to match the averaged data from the inward current experiments. We then further modified this model to create 22 tailored models based on the 22 cells in which outward current was measured, by applying the I_Ks and I_NaCa scaling factors from Table 2.

Significant variability in the outward currents was observed among the Cor.4U cells. This can be seen from the scaling factors in Table 2, but it is also evident when directly inspecting the currents measured from different cells (see Supplementary Figure S5) or when looking at peak I_outward (see Supplementary Figure S7).

Figure 3 (left panel) shows APs simulated with the tailored models. A wide variety of APs could be seen, with some models showing a spike-and-dome waveform, some showing a more triangular waveform, and with a varying slope in resting potential (leading to different degrees of auto-excitation). Some models also show beat-to-beat alternans, or fail to completely depolarize. The corresponding contribution of the major currents throughout the APs are also shown in Supplementary Figure S8.

Recordings of APs in single iPSC-CMs show a similar variety of AP waveforms. Figure 2 (right panel) shows APs measured in 7 different iCell iPSC-CMs. Again, various waveform morphologies (roughly corresponding to atrial, ventricular and sinoatrial node APs) and differing levels of auto-excitability can be distinguished. Whilst these recordings are for a different cell line than our tailored models, the inter-cell variability in channel expression within a batch of iPSC-CMs has not been observed to be markedly different between cell lines (see e.g., the relative size of the “error bars” in Figure 2 of Blinova et al., 2017).

Figure 4 (left panel) shows the median of all simulated traces as shown in Figure 3, along with the 25th and 75th percentiles. The optically recorded APs from the Cor.4U cells were plotted on the same graph (the median shown as black line and the 25th and 75th percentiles shown as gray shading). Due to the increased outward current, the tailored models exhibit a shorter APD than the original model, that matches the measured APDs more closely in the early and late repolarization phase. A histogram of APDs in measured and simulated cells is shown in Figure 4 (right panel), with the blue line representing the result from the original model. A similar histogram for APD₅₀ is shown in Supplementary Figure S9.

Figure 5 shows the dose-response curves of the APD₉₀ of four drugs, measured experimentally and simulated using the original and tailored models. Equivalent results using the of the APD₅₀ are shown in Supplementary Figure S9. Results for the control drug paracetamol are shown in Supplementary Figure S10. For all four drugs tested, although not fitting the experimental data exactly, the tailored models match the measured data more closely than the original model. For Dofetilide in particular, the tailored models show a realistically smaller increase in APD than the original model, which shows alternans and then repolarization failure at higher drug concentrations. For Quinidine, although the tailored models do not fit better at the highest concentration, we improve the predictions at lower concentrations.

Predictions made with the adult-CM model by O'Hara et al. (2011) are shown for comparison. Note how the adult CM model predicts APD prolongation with Verapamil, whereas both the tailored and original iPSC-CM models accurately predict the shortening that is observed in iPSC-CM optical mapping. Such qualitative differences highlight the need for models specific to iPSC-CMs to interpret experimental findings in these cells.

We obtained maximum conductance values for the inward currents I_Na and I_CaL that are lower than suggested by the Paci et al. (2013) model based on iCell cells, while I_Ks and I_NaCa were increased in most Cor.4U cells. The slight reduction in I_CaL and increase in I_Ks suggests a decrease in APD. This was borne out by the AP simulations, and was consistent with our optical mapping measurements which showed shorter APDs compared to the AP simulated by the original Paci model. Our simulations displayed a similar degree of AP variability to the experimental iPSC recordings, but larger variability than the optical mapping measurements. Both findings are consistent given that electrotonic coupling of cells (present in the optical mapping experiments) reduces variability.

A potential explanation of the large I_Ks current is suggested by Lei et al. (2017). It shows that both the KCNQ1 and KCNE1 (subunits of the channel carrying I_Ks) were present in our Cor.4U cells, however, KCNE1 was not as well expressed in iCells. The difference in KCNE1 expression could lead to the observed larger I_Ks currents in the Cor.4U cells compared to iCells, and hence a shorter APD and less prolongation under I_Kr blockers, which is in agreement with Blinova et al. (2017). Our observation is supported by Silva and Rudy (2005) who found that native I_Ks (from channels comprised of both KCNQ1 and KCNE1) activates more than with KCNQ1 only.

iPSC cardiomyocytes, from the same donor and differentiated/matured in the same way, can display vastly different AP waveforms, reminiscent of those of ventricular, atrial, and sinoatrial-node cells. Our tailored models, created by varying the maximum conductances of I_NaCa and I_Ks, showed a similar model-to-model (cell-to-cell) variety in style of generated APs. This shows that variation in genetic expression, which correlates directly with maximum conductance (Schulz et al., 2006), could be enough to explain the different AP waveforms observed in iPSC-CMs. However, it does not preclude other explanations, and it is possible the APs could take on a more distinct shape if differences in ion channel kinetics were also included. As discussed in a recent white paper, the inclusion of cell-cell variability, as well as variability between cell lines is an important research area (Johnstone et al., 2016).

The sharp increase in I_Ks seen in our Cor.4U tailored models suggests Cor.4U cells have a stronger reliance on I_Ks as a repolarizing force, and will therefore be less likely to show AP prolongation when treated with I_Kr blocking drugs (see, e.g., the Figure 5 of Blinova et al., 2017, which shows, for 8 out of 12 drugs with comparable concentrations and prolongation in iCells, the Cor.4U cells have a smaller APD prolongation than the iCells). Consistent with this suggestion, simulations of treatment with the potent I_Kr blocker Dofetilide showed only a modest increase in APD at concentrations that caused the iCell-cell based model to display excessive AP prolongation resulting in alternans. Treatment with Quinidine, a less potent I_Kr blocker, showed similar results. The modest APD increase predicted by the tailored models underestimated the APD prolongation observed in the data, suggesting the role of I_Kr as a repolarizing force was underestimated in these models. More refined experiments will need to be conducted to separate the outward currents and to better estimate I_Kr conductance. Application of Verapamil, which blocks I_CaL as well as I_Kr, had a smaller effect in our tailored models than in the original model, which is consistent with the lowered levels of I_Kr and I_CaL.

The strong I_Ks-reliance we observed may be problematic when using these iPSC-CMs as models for ventricular myocytes, where I_Ks only plays a major part when other repolarizing currents are blocked or in the presence of sympathetic stimulation (Jost et al., 2005).

This study showed the need to build cell-line or even cell-specific models for iPSC-CMs, and this work serves as a pilot attempt for such an approach. However, a refined study with additional experiments will be needed to improve the tailored models further.

As might be expected, the ion current profiles during voltage steps could not be recreated well using this approach. It is likely that ion channel kinetics also vary between cell lines due to (e.g.,) differences in subunit expression (Lei et al., 2017), although this could also be partly due to the difference of temperature, and that the model we used does not accurately capture the kinetics of ion currents in Cor.4U cells. Channel kinetics play an important role in the contribution of a current to the different phases of the AP. Modifying the kinetic parameters which characterize the voltage-current relationship for the activation, inactivation, deactivation, etc. of a channel could change both the current and the AP, and would influence responses to drugs. Varying the kinetic parameters would also alter the conductances we estimated by fitting the outward current. Further tailoring the models to include refitted kinetic parameters may lead to further improvements in predictive power. However, since models of ion channel kinetics contain many parameters, specifically designed experiments (e.g., with channel blockers and/or specialized voltage protocols) will be required to refine these tailored models.

The method of fitting multiple currents to a single experimental recording is a highly useful approach, as it reduces the number of experiments needed to tailor a model. However, due to the limitation of experiments being performed in different cells, we were not able to examine the covariance between the inward and outward currents. Also, since it depends on the number of current conductances to be fitted and the experimental data (e.g., the quality of the data and the actual current shape), one may run into problems of practical identifiability if one tries to refit kinetic parameters here (e.g., multiple combinations of conductance and kinetic parameters that can provide an equally good fit, as in Fink and Noble, 2009). Additional experimental data with refined experimental designs will be needed to identify all parameters; for example, to perform experiments with channel blockers to isolate the contribution of particular currents, or to iteratively refine the models using a dynamic clamp approach (Devenyi et al., 2017).

Variability/noise on drug-ion channel interaction parameters (IC50s) from different labs or repeats of experiments will also impact our simulation predictions. A probabilistic uncertainty quantification framework using the techniques proposed in Elkins et al. (2013), Johnstone et al. (2017) could be used in future to address this.

Our Cor.4U-tailored predictions of both baseline AP and drug responses matched the optical mapping data more closely than the un-tailored model. However, these optical mapping data were gathered from cultures of spontaneously beating electrotonically coupled cells, while our simulations are of paced single iPSC-CMs. Another avenue for future work would be to combine (a representative distribution of) tailored cell-specific models into heterogeneous tissue models (Bowler et al., 2016).

iPSC-CMs have gained significant popularity as an in vitro model for drug screening and, as one pillar of the CiPA strategy, are anticipated to become a routine part of the cardiac safety pipeline. It is therefore critical to understand how to interpret the iPSC-CM data variability (intra- and inter-cell line variability) and to translate these data to the adult human situation. Mathematical models are a promising tool to integrate data, gain mechanistic insights and perform this translation.

Our results show that differences between iPSC-CM cell lines can be analyzed and understood using tailored computational models. Furthermore, even models based on relatively simple methods (e.g., scaling maximum conductances) and a limited set of measurements (two inward current and one outward current experiments) can lead to improved predictions of baseline and drug-blocked electrophysiology parameters.

LP and KW were employed by company F. Hoffmann—La Roche. MH-V, VZ, AA, and GS were employed by Clyde Biosciences. The reviewer DK and handling Editor declared their shared affiliation, and the handling Editor states that the process met the standards of a fair and objective review. The other 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.