A computational model of induced pluripotent stem-cell derived cardiomyocytes incorporating experimental variability from multiple data sources

Abstract

Induced pluripotent stem cell-derived cardiomyocytes (iPSC-CMs) capture patient-specific genotype-phenotype relationships, as well as cell-to-cell variability of cardiac electrical activity Computational modelling and simulation provide a high throughput approach to reconcile multiple datasets describing physiological variability, and also identify vulnerable parameter regimes We have developed a whole-cell model of iPSC-CMs, composed of single exponential voltage-dependent gating variable rate constants, parameterized to fit experimental iPSC-CM outputs We have utilized experimental data across multiple laboratories to model experimental variability and investigate subcellular phenotypic mechanisms in iPSC-CMs This framework links molecular mechanisms to cellular-level outputs by revealing unique subsets of model parameters linked to known iPSC-CM phenotypes ABSTRACT: There is a profound need to develop a strategy for predicting patient-to-patient vulnerability in the emergence of cardiac arrhythmia. A promising in vitro method to address patient-specific proclivity to cardiac disease utilizes induced pluripotent stem cell-derived cardiomyocytes (iPSC-CMs). A major strength of this approach is that iPSC-CMs contain donor genetic information and therefore capture patient-specific genotype-phenotype relationships. A cited detriment of iPSC-CMs is the cell-to-cell variability observed in electrical activity. We postulated, however, that cell-to-cell variability may constitute a strength when appropriately utilized in a computational framework to build cell populations that can be employed to identify phenotypic mechanisms and pinpoint key sensitive parameters. Thus, we have exploited variation in experimental data across multiple laboratories to develop a computational framework for investigating subcellular phenotypic mechanisms. We have developed a whole-cell model of iPSC-CMs composed of simple model components comprising ion channel models with single exponential voltage-dependent gating variable rate constants, parameterized to fit experimental iPSC-CM data for all major ionic currents. By optimizing ionic current model parameters to multiple experimental datasets, we incorporate experimentally-observed variability in the ionic currents. The resulting population of cellular models predicts robust inter-subject variability in iPSC-CMs. This approach links molecular mechanisms to known cellular-level iPSC-CM phenotypes, as shown by comparing immature and mature subpopulations of models to analyse the contributing factors underlying each phenotype. In the future, the presented models can be readily expanded to include genetic mutations and pharmacological interventions for studying the mechanisms of rare events, such as arrhythmia triggers.

Keywords: computer modelling; iPSC-CMs; variability.

Figure 1. A schematic of the computational iPSC‐CM model

Red stars indicate individual currents (*I _Na, I _CaL, I _Kr, I _Ks, I _K1, I _to and I _f), formulated using single‐exponential voltage‐dependent rate constants. Parameters were optimized to experimental iPSC‐CM kinetic data. The mathematical formulation for an example current, I _Na, is shown on the right. All gating variables in the starred currents were modelled using the example formula for gating variable x. Additional calcium‐dependent currents (I _NCX, I _PMCA, I _CaT and SR currents: I _SERCA, I _RyR and I _leak) were modelled using previously published model formulations, optimized to calcium transient data from iPSC‐CMs. Remaining currents (I _NaK, I _bCa, I _bNa) were modelled using ten Tusscher 2004 formulations optimized to recapitulate iPSC‐CM AP data.

Figure 2. Flow chart

Flow chart showing the methodology for building the iPSC‐CM model populations.

Figure 3. Sodium current (I _Na) model optimization

A, steady‐state inactivation and activation curves. Dataset‐specific model fits (lines) optimized to experimental data (points). The sodium current model used in the baseline whole‐cell model is shown in black. Colours represent distinct experimental iPSC‐CM data from Ma et al. (2011) and from immature and mature cell preparations from the Jalife lab (Herron et al. 2016). B, I–V curves for I _Na. Dataset‐specific models were simulated using the experimental conditions of the corresponding experimental data. C, I _Na activation (m‐gate) time constants. D, I _Na fast‐inactivation (h‐gate) time constants. E, I_Na slow‐inactivation (j‐gate) time constants. j‐gate time constant parameters for all I _Na models were optimized to experimental iPSC‐CM data from the Kurokawa lab (Li et al. 2017).

Figure 4. Calcium current model optimization

A, L‐type calcium current (I _CaL) steady‐state inactivation and activation curves with dataset‐specific model fits (lines) optimized to experimental data (points). The L‐type calcium model used in the baseline cellular model is shown in black. Coloured symbols represent experimental iPSC‐CM data from Ma et al. (2011), Veerman et al. (2016) and Es‐Salah‐Lamoureux et al. (2016). B, I–V curves for I _CaL. Calcium‐dependent gating model formulation retained from ten Tusscher 2004 adult cardiomyocyte model with parameter optimization to fit whole cell iPSC‐CM outputs. C, time constants of I _CaL activation gate. Time constant parameters for all I _CaL models were optimized to experimental iPSC‐CM data from Ma et al. (2011). D, time constants of I _CaL inactivation gate. E, optimization of peak T‐type calcium current (I _CaT) to experimental iPSC‐CM data from the Kurokawa lab (Li et al. 2017). Model formulation of I _CaT was retained from the Maltsev and Lakatta sinoatrial node model.

Figure 5. Rapid delayed rectifier potassium current (I _Kr) model optimization

A, steady‐state activation with dataset‐specific model fits (lines) optimized to experimental data (points). The I _Kr model used in the baseline cellular model is shown in black. Coloured symbols represent experimental iPSC‐CM data from Ma et al. (2011), the Wu lab (Garg et al. 2018), Es‐Salah‐Lamoureux et al. (2016) and Bellin et al. (2013). For I _Kr inactivation gating, existing ten Tusscher 2004 model components were reformulated to single exponential forms. B, I–V curves for I _Kr.C, time constants of the I _Kr activation gate. Activation time constants for the models of Ma et al. (2011), the Wu lab (Garg et al. 2018) and Es‐Salah‐Lamoureux et al. (2016) were optimized to experimental iPSC‐CM data from Ma et al. (2011). D, time constants of I _Kr inactivation gate using the ten Tusscher 2004 model reformulated to single exponential forms.

Figure 6. Transient outward potassium current (I _to) model optimization

A, steady‐state activation and inactivation curves with dataset‐specific model fits (lines) optimized to experimental data (points). The I _to model used in the baseline cellular model is shown in black. Coloured symbols represent experimental iPSC‐CM data from Veerman et al. (2016), Ma et al. (2011) and Cordeiro et al. (2013). B, I–V curves for I _to. C, time constants of I _to activation gate. For activation time constants in all I_to models, the ten Tusscher 2004 I_to activation time constants were reformulated to single exponential forms. D, time constants of I_to inactivation gate. Model time constant parameters of Veerman et al. (2016) were optimized to iPSC‐CM experimental data from Ma et al. (2011).

Figure 7. Slow delayed rectifier potassium current (I _Ks) model optimization

A, steady‐state activation with dataset‐specific model fits (lines) optimized to experimental data (points). The I _Ks model used in the baseline cellular model is shown in black. Coloured symbols represent experimental iPSC‐CM data from Ma et al. (2011) and two separate iPSC‐CM cell‐line datasets in Ma et al. (2015). B, time constants of the I _Ks activation gate. Time constants for all I_Ks models were optimized to experimental iPSC‐CM data from Ma et al. (2011). C, I–V curves for I _Ks.

Figure 8. Pacemaker/funny current (I _f) model optimization

A, steady‐state activation with dataset‐specific model fits (lines) optimized to experimental data (points). The I _f model used in the baseline cellular model is shown in black. Coloured symbols represent experimental iPSC‐CM data from the Kurokawa lab (Li et al. 2017) and Ma et al. (2011). B, time constants of the I _f inactivation gate. C, I–V curves for I _f.

Figure 9. Inward rectifier potassium current (I _K1) model optimization

I–V curves for I _K1 with dataset‐specific model fits (lines) optimized to experimental data (points). The I _K1 model used in baseline cellular model is shown in black. Coloured symbols represent experimental iPSC‐CM data from Ma et al. (2011), the Kurokawa lab (Li et al. 2017) and immature and mature cell preparations from the Jalife lab (Herron et al. 2016).

Figure 10. Optimization of calcium handling in the iPSC‐CM baseline model

A, experimental iPSC‐CM CaT traces from the Wu lab (grey) (Garg et al. 2018) with baseline model CaT (red). Experimental data were reported as the normalized Ca²⁺ florescence (F _ratio). Separately, average iPSC‐CM peak and diastolic Ca²⁺ concentrations were measured by the Wu lab. The two y‐axes are plotted so that the average F_ratio peak and diastolic values of the experimental dataset shown correspond to the average experimental concentration of peak and diastolic Ca²⁺ (B). The baseline model CaT output is nm. B, comparison of baseline model CaT morphology markers with experimental iPSC‐CM data from the Wu lab. C, relative contribution of calcium from I _SERCA, I _NCX and I _PMCA to the CaT during a single AP in the baseline model. D, comparison of experimental (black and white) and baseline model (coloured) relative contribution of calcium flux from I _SERCA, I _NCX and I _PMCA during the CaT. Experimental data from Hwang et al. (2015).

Figure 11. Characterization of the baseline model AP

A, time course of the spontaneously beating APs in the baseline model. B–D, comparison of AP morphology in the baseline model (red) and experimental iPSC‐CM data (black). Experimental data from the Wu Lab (Garg et al. 2018), the Jalife Lab (Herron et al. 2016), Ma et al. (2011), Doss et al. (2012), Cordeiro et al. (2013), Es‐Salah‐Lamoureux et al. (2016) and Ma et al. (2015). E, sensitivity analysis using multivariable regression in the baseline model. Only parameters with regression coefficients >0.3 are shown.

Figure 12. Kinetic variability generated by varying individual current model parameters

Steady‐state and time constant curves for each gate in (A) I _Na, (B) I _CaL, (C) I _Kr and (D) I _f. E, I–V curves for I _K1. Dataset‐specific model fits (black lines, also shown in Figs 3, 4, 5, 6, 7, 8, 9, 10), randomly‐parameterized models resulting in spontaneous AP generation in the cell models (coloured lines) and randomly‐parameterized models resulting in non‐spontaneous or non‐AP generating model cells (grey lines) are all shown.

Figure 13. Variation of action potential morphology in model iPSC‐CM populations

A, APs of spontaneously beating cells (n = 25,434) generated by varying one current at a time (I _Na, I _CaL, I _Kr, I _f and I _K1). B, APs of spontaneously beating cells (n = 17,139) generated by varying the same five currents simultaneously. C, representative AP time courses of spontaneously beating cells at various pacing frequencies. D–F, comparison of AP morphology in the populations of models (colour) and experimental iPSC‐CM data (black). Each coloured point represents a spontaneously beating cell created by varying a single current (A), or by varying all five currents simultaneously (B). G, mean ± SD of AP morphology measures for each population, normalized to the baseline model AP.

Figure 14. Sample APs showing the effect of ion channel blockers within the model population

Showing the same cellular models in the control (solid lines) and drugged (dashed lines) conditions. Three cells are shown for each drug, representing a cell with a change in the given AP parameter near the population mean (cell 2) (± 1 SD). The mean ± SD for the full population are described in Table 4. Drug effects are shown for (A) TTX, (B) nifedipine and (C) E‐4031.

Figure 15. Comparison of immature and mature cellular models

A, AP for the immature (baseline) cellular model compared to AP for a representative mature cellular model. Immature and mature cellular models determined by scaling g _K1 and g _Na based on peak currents reported in iPSC‐CMs with control (immature) and maturation‐promoting cell preparations (Herron et al. 2016). B, comparison of models and experimental AP morphology for mature and immature cell‐types. Experimental data are from the Jalife Lab (Herron et al. 2016). All of the model and experimental data were normalized to the Jalife Lab mature experimental dataset average. C, comparison of sensitivity analysis immature and mature models using multivariable regression. Only parameters with regression coefficients >0.3 are shown.

Figure 16. Comparison of mature and immature iPSC‐CM model subpopulations

A, division of model population into mature and immature phenotypes (using five‐current variation population with simultaneous variation in I _Na, I _CaL, I _Kr, I _F and I _K1 parameters). Experimental data are from the Jalife Lab (Herron et al. 2016). Model subpopulation shown in red (n = 325) represents mature phenotypes with MDP <−75 mV and maximal upstroke velocity >85 V s⁻¹. The model subpopulation shown in blue (n = 13 759) represents immature phenotypes with MDP >−75 mV and maximal upstroke velocity <85 V s⁻¹. The model subpopulation shown in grey was not analysed in this comparison. B, the four model parameters with the largest difference between the mature and immature model subpopulations. For each subpopulation, parameter averages and SDs are shown as the percentage change from the baseline model parameter value. C, steady‐state inactivation for I _Na in the mature and immature model subpopulations. Individual cells (light colours) and subpopulation average parameter values (darker colored lines) are shown. D, peak I _K1 and I _Na for the I–V relationship of each cell in the model subpopulations were compared with data reported in Herron et al. (2016). Model and experimental values are shown as the percentage change from the immature to mature phenotype.