Creating cell-specific computational models of stem cell-derived cardiomyocytes using optical experiments

32 Human induced pluripotent stem cell-derived cardiomyocytes (iPSC-CMs) have gained traction as a powerful 33 model in cardiac disease and therapeutics research, since iPSCs are self-renewing and can be derived from 34 healthy and diseased patients without invasive surgery. However, current iPSC-CM differentiation methods 35 produce cardiomyocytes with immature, fetal-like electrophysiological phenotypes, and the variety of 36 maturation protocols in the literature results in phenotypic differences between labs. Heterogeneity of iPSC 37 donor genetic backgrounds contributes to additional phenotypic variability. Several mathematical models of 38 iPSC-CM electrophysiology have been developed to help understand the ionic underpinnings of, and to 39 simulate, various cell responses, but these models individually do not capture the phenotypic variability 40 observed in iPSC-CMs. Here, we tackle these limitations by developing a computational pipeline to calibrate 41 cell preparation-specific iPSC-CM electrophysiological parameters. 42 We used the genetic algorithm (GA), a heuristic parameter calibration method, to tune ion channel parameters 43 in a mathematical model of iPSC-CM physiology. To systematically optimize an experimental protocol that 44 generates sufficient data for parameter calibration, we created simulated datasets by applying various 45 protocols to a population of in silico cells with known conductance variations, and we fitted to those datasets. 46 We found that calibrating models to voltage and calcium transient data under 3 varied experimental conditions, 47 including electrical pacing combined with ion channel blockade and changing buffer ion concentrations, 48 improved model parameter estimates and model predictions of unseen channel block responses. This 49 observation held regardless of whether the fitted data were normalized, suggesting that normalized 50 fluorescence recordings, which are more accessible and higher throughput than patch clamp recordings, could 51 sufficiently inform conductance parameters. Therefore, this computational pipeline can be applied to different 52 iPSC-CM preparations to determine cell line-specific ion channel properties and understand the mechanisms 53 behind variability in perturbation responses. 54 55 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 3 Author Summary: 56 Many drug treatments or environmental factors can trigger cardiac arrhythmias, which are dangerous and often 57 unpredictable. Human cardiomyocytes derived from donor stem cells have proven to be a promising model for 58 studying these events, but variability in donor genetic background and cell maturation methods, as well as 59 overall immaturity of stem cell-derived cardiomyocytes relative to the adult heart, have hindered reproducibility 60 and reliability of these studies. Mathematical models of these cells can aid in understanding the underlying 61 electrophysiological contributors to this variability, but determining these models’ parameters for multiple cell 62 preparations is challenging. In this study, we tackle these limitations by developing a computational method to 63 simultaneously estimate multiple model parameters using data from imaging-based experiments, which can be 64 easily scaled to rapidly characterize multiple cell lines. This method can generate many personalized models of 65 individual cell preparations, improving drug response predictions and revealing specific differences in 66 electrophysiological properties that contribute to variability in cardiac maturity and arrhythmia susceptibility.67 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 4 ABBREVIATIONS 68  AP: Action potential 69  CaT: Calcium transient 70  GA: Genetic algorithm 71  Gx, Ix, Jx: Maximal conductance (G), current density (I), or flux (J) for x, where x can be: Na (Na+), f 72 (funny Na+), CaL (L-type Ca2+), to (transient outward K+), Ks (slow delayed rectifier K+), Kr (rapid 73 delayed rectifier K+), K1 (inward rectifier K+), PMCA (plasma membrane Ca2+ ATPase), bNa 74 (background Na+), bCa (background Ca2+), Up (SR Ca2+ uptake), rel (SR Ca2+ release), NCX (Na+/Ca2+ 75 exchanger), NaK (Na+/K+ ATPase), SRleak (SR leak), or CaT (T-type Ca2+) 76  iPSC-CM: Induced pluripotent stem cell-derived cardiomyocyte 77  LTCC: L-type calcium channel 78  NCX: Sodium-calcium exchanger 79  SERCA: Sarco/endoplasmic reticulum calcium ATPase 80  SR: Sarcoplasmic reticulum 81 82 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 5 GLOSSARY 83  Model/parameter calibration: tuning one or more parameters in the computational model so that the 84 model output more closely matches experimental data 85  Experiment/protocol optimization: the process of determining what type and amount of data is 86 sufficient but also feasible for our model calibration goals 87 o Protocol conditions – buffer calcium, potassium, or sodium concentrations; addition or removal 88 of stimulus; pacing rates; channel block; etc. 89 o Protocol length(?) – number of protocol conditions 90 o Protocol data type(?) –AP, CaT, or both; normalized or non-normalized data 91  Model prediction: using the calibrated computational model to simulate response to new (unseen) 92 conditions, drugs, or perturbations (in our case, IKr block) 93 o Papers on independent validation/prediction 94 Computational pipeline: the full process of iPSC-CM computational model calibration; Includes iPSC-CM 95 data acquisition/simulation -> data processing -> parameter calibration using genetic algorithm -> validation of 96 calibrated models on an unseen condition (i.e. evaluating model predictions) 97 98 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 6

99 Human iPSC-derived cardiomyocytes (iPSC-CMs) are derived from patient cells that have been reprogrammed 100 into pluripotency and subsequently manipulated to differentiate into the cardiac cell lineage. The development 101 of this in vitro platform marked a technological breakthrough in cardiac research, and iPSC-CMs are now 102 widely-used in cardiac pharmacology and disease research. iPSCs can be sourced through minimally-invasive 103 procedures, such as skin biopsy or blood draw, and maintained or banked for long periods [1,2]. These 104 properties of iPSC-CMs make them an ideal platform for pharmacological studies and personalized disease 105 modeling [3]. However, electrophysiological variability between iPSC-CM preparations from different cell lines 106 and differentiation methods limit the potential of this platform [4,5]. 107 Mathematical models of cardiomyocyte electrophysiology can provide valuable insight into iPSC-CM 108 electrophysiological variability. These models contain parameters describing the ionic currents and fluxes that 109 contribute to the properties of cardiac action potential (AP) and calcium transient (CaT) waveforms. Several 110 models of iPSC-CM physiology exist in the literature [6–8], with the most recent and comprehensive model 111 published by Kernik et al. in 2019, hereon referred to as the “Kernik model” [9]. However, the baseline 112 parameter values in these models are not representative of the electrophysiological heterogeneity observed in 113 iPSC-CMs. Additionally, parameter calibration has often relied on separate patch clamp measurements of each 114 type of current. This process is time-consuming, inaccessible to many research groups, and often captures the 115 physiology of only the average behavior of the population of cells examined. Methods to simultaneously 116 optimize several or all model parameters using minimal data are under development [10–14]. Recent work on 117 guinea pig ventricular myocyte models [13] and human iPSC-CM models [14] have demonstrated the utility of 118 automated tuning algorithms for improving model parameterization. For example, in Devenyi et al [13], an 119 adjusted model, constrained by data, predicted the effect of slow delayed rectifier K+ current (IKs) block more 120 accurately compared to baseline model. Analogous work, at the level of individual ion channels, has shown the 121 superiority of sinusoidal voltage clamp protocols, compared with traditional square-wave protocols, for 122 calibration of channel kinetic parameters [15]. 123 Even with recent advances in iPSC-CM maturation protocols [16–19], most iPSC-CMs still display embryonic- 124 or neonatal-like cardiac phenotypes with different AP and CaT shapes compared with adult myocytes. These 125 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 7 phenotypic differences, and how they change temporally, can depend on subtle differences in cell culture 126 conditions and the genetic background of the cell donor [4,5,20], which likely contributes to the wide variation 127 observed between studies [4,5]. No standardized procedure currently exists for characterizing iPSC-CM 128 heterogeneity; prior efforts have generally examined a handful of molecular markers [21], used subjective 129 observations of physiology [22,23], or employed specialized methods such as single cell RNA sequencing 130 [24,25]. These factors, and the potential importance of iPSC-CMs as a research platform, highlight the need for 131 automated methods to characterize the cell lines used in each study. 132 Here we combined these ideas of simultaneous calibration of multiple parameters and optimization of a 133 minimal experimental protocol for generating calibration data. We aimed to create digital twins of iPSC-CM cell 134 preparations from the Kernik iPSC-CM model, thereby connecting observed iPSC-CM electrophysiology with 135 molecular function and mechanisms. We hypothesize that fluorescence readouts of iPSC-CM physiology under 136 varied experimental conditions provide enough information to reveal cellular ion channel properties, which can 137 then be incorporated into the iPSC-CM digital twins. To systematically evaluate various data types and 138 experimental protocols in their ability to inform model parameters, we created a simulated, in silico dataset 139 from Kernik model variations with known parameter values. From this study, we developed a computational 140 pipeline which includes: 1) an optimized protocol for fluorescence recordings of iPSC-CM preparations; 2) a 141 genetic algorithm process for calibration of ion channel parameters in iPSC-CM computational model, using 142 the experimental recordings; and 3) validation of the resulting calibrated models by evaluating their predictions 143 on independent, yet physiologically-important, perturbations. We demonstrate the utility of our computational 144 pipeline in generating iPSC-CM digital twins that capture ionic variability and predict cell-specific 145 electrophysiological phenotypes, including drug-induced arrhythmia susceptibility. 146 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 8

147 Our primary objective was to optimize an experimental protocol for generating sufficient data to calibrate an 148 iPSC-CM mathematical model, taking experimental feasibility and throughput into account. To systematically 149 evaluate how different types of data and experimental conditions impact parameter calibration and model 150 predictions of responses to new conditions or drug treatments, we created a simulated in silico dataset, 151 generated from a population of Kernik models with random variations in their 16 maximal conductance 152 parameters (Figure 1A, Supplementary Table S2). We simulated AP and CaT generated by these model cells 153 under various conditions (Supplementary Table S3), and then concatenated these data in various 154 combinations based on the “candidate protocols” we wanted to evaluate. The 16 Kernik model conductance 155 parameters were then calibrated, using a genetic algorithm, to best match the steady state AP and CaT traces 156 from each candidate protocol. These fitted parameters were incorporated into the Kernik model, and the newly-157 calibrated models were evaluated on their ability to predict the original model cell’s responses to new 158 conditions (outlined in methods and Figure 1B). Optimizing the experimental protocol using this in silico dataset 159 provided 3 major advantages: 1) knowledge of the ground truth parameter values that generated the dataset, 160 so the accuracy of the calibrated parameters can be evaluated, 2) relatively fast and easy simulation of new 161 conditions from the same or additional model cells, in comparison to acquiring new in vitro data and cell lines, 162 and 3) generating corresponding data types through computational methods and data processing, forgoing the 163 need to set up new experiments to collect each data type. 164 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 9 165 Simultaneously calibrating to AP and CaT data improves model predictions 166 First, we examined whether AP or CaT recordings alone provided sufficient information for conductance 167 parameter calibration and model predictions. We used simulated recordings from baseline physiological 168 conditions without pacing stimuli (spontaneous APs), and supplied to a genetic algorithm (GA) either 1) only 169 the AP traces, 2) only the CaT traces, or 3) both (Figure 2A). The GA calibrates parameters by varying them to 170 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 10 create a population of Kernik model cells, then iteratively selecting and modifying individual models that 171 generated AP and CaT traces that best matched the supplied data (outlined in Figure 1B ) [26]. We ran the 172 genetic algorithm 10 times per input dataset, with each GA run starting on a different initial population of 173 models. This allowed us to assess whether each input dataset provided enough information to consistently 174 estimate parameter values, even when the initial search space varied. 175 176 We evaluated each candidate protocol on parameter calibration accuracy (calibration error) and consistency 177 (calibration spread). Calibration error was represented by the absolute values of calibrated value log-178 normalized to the corresponding ground truth parameter value, and calibration spread was assessed by 179 calculating the standard deviation (spread) in log-normalized calibrated parameter values between the 10 runs 180 (Figure 2B). After averaging over all 16 fitted parameters, we found minimal differences in calibration error (p = 181 0.11) or spread (p = 0.07) between fitting only voltage, only calcium, or both simultaneously ( Figure 2C). From 182 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 11 these initial tests, therefore, there appears to be little difference between experimental protocols in how well (or 183 poorly) each can identify model parameters. 184 We next assessed the calibrated models from each candidate protocol on how well each predicted model cell 185 responses to 30% block of rapid delayed rectifier current IKr,as this perturbation is both: 1) independent of the 186 conditions used for model calibration and 2) relevant to cardiac pharmacology [27]. With this evaluation (Figure 187 3A), we found that models fitted to only AP data or only CaT data differed significantly from the ground truth 188 model (Figure 3B, Supplementary Figure S1). When both AP and CaT recordings were used simultaneously 189 during parameter calibration, the resulting calibrated models exhibited substantial improvement in their 190 predictions of 30% IKr block. Therefore, fitting parameters to both AP and CaT traces simultaneously is 191 required for accurate prediction of results that were not used in the fitting process. 192 193 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 12 Overall parameter error and spread show minimal change with varied protocol conditions 194 Next, we explored how varying the number and identity of simulated conditions in the candidate protocols 195 affected parameter calibration. AP and CaT recordings were simulated under many conditions that could be 196 produced in an electrophysiology laboratory and constructed candidate protocols from various combinations of 197 these conditions, which included: 1) varied buffer calcium, from hypo- to normal to hyper-calcemic conditions, 198 2) spontaneous beating and varied pacing rates from 1-2 Hz, 3) varied levels of L-type calcium channel (ICaL) 199 blockade, and 4) a mixture of these conditions (see Methods for details). We also separated these data into 200 shorter protocols, to see if the experimental setup could be further simplified (Figure 4A). Unexpectedly, we 201 found that parameter errors and spread across all 16 maximal conductance parameters did not differ 202 significantly between the protocols with varied cell culture conditions (Figure 4B, p = 0.814 and p = 0.673, 203 respectively). Additionally, we observed no significant differences in average calibration errors (p = 0.088) and 204 spreads (p = 0.137) when varying the number of conditions in the protocol (Figure 4C). Similar findings were 205 observed when evaluating different combinations of these experimental conditions (Supplementary Figure S2). 206 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 13 207 A readily-achievable 3-condition protocol improves predictions of IKr block response 208 As with the results shown in Figure 3, we assessed how well calibrated models predicted the response to IKr 209 block, and we found that accuracy depended strongly on the experimental conditions simulated in each 210 candidate protocol (sample AP traces shown in Figure 5A). Out of the protocols using 3 experimental 211 conditions in Figure 4A, the mixed-conditions protocol generated calibrated models with the most accurate 212 predictions of changes in APs with 30% IKr block, compared with protocols which varied only the buffer 213 calcium, pacing rates, or ICaL blockade (Figure 5B-C). When evaluating whether a shorter version of this 214 protocol could produce equally predictive models, the full protocol with all 3 conditions outperformed protocols 215 with 1 or 2 of the conditions (Figure 5B-C). These data suggest that an optimal experimental protocol for iPSC-216 CM model parameter calibration should include voltage and calcium transient fluorescence recordings under at 217 least 3 varieties of cell culture condition changes or perturbations. Of note, this supports our hypothesis that 218 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 14 data from iPSC-CMs under a variety of conditions provides sufficient parameter identification to generate 219 predictive models for pharmacological applications. 220 221 Normalized fluorescence recordings sufficiently inform parameter calibration 222 Fluorescence recordings generally detect relative changes and do not provide absolute levels of voltage and 223 calcium. Fluorescent dyes can be calibrated [28,29], and patch clamp can detect true transmembrane 224 potentials, but these are more challenging and lower-throughput techniques, compared with straightforward 225 fluorescence recordings. To assess the potential implications for parameter identification, we compared two 226 versions of simulated recordings obtained with our optimized protocol: 1) unscaled (original) in silico data, 227 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 15 representing patch clamp and calibrated calcium measurements, and 2) normalized versions of the same 228 waveforms, representing simpler fluorescent recordings (Figure 6A). 229 230 Parameter calibration accuracy and consistency did not differ significantly when comparing calibrated models 231 from the original data with those from the corresponding normalized data (Figure 6B). Notably, calibrated 232 models from these two data types also performed similarly when evaluating predictions of 30% I Kr block 233 response (Figure 6C). This suggests that fluorescence voltage and calcium recordings from our optimized 234 protocol conditions provide sufficient information to calibrate predictive models. 235 Calibrated models predict variability in arrhythmia susceptibility in silico 236 The primary goal of optimizing a computational pipeline for model calibration is to create cell preparation-237 specific models that can predict variability in phenotypes, particularly arrhythmia susceptibility, between iPSC-238 CM lines. To assess the ability of our optimized protocol to predict cell line-specific arrhythmia susceptibility, 239 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 16 we calculated the lowest level of IKr block (i.e. highest % IKr) that induced arrhythmia dynamics 240 (afterdepolarizations, Torsades de Pointes, alternans, beating cessation, or tachycardia) in the same 4 Kernik 241 model cells from our in silico dataset. These simulations defined, for each of the 4 model cells, an “IKr block 242 tolerance threshold” (Figure 7A), which ranged from 37% IKr block (highest susceptibility) to 62% IKr block 243 (highest tolerance) (Figure 7C) across the 4 model cells. 244 245 Predicted IKr block tolerance thresholds were also calculated for the calibrated models generated by our 246 computational pipeline. These predicted thresholds were assessed against the ground truth threshold of the 247 model cell. Calibrated models from our computational pipeline predicted IKr block thresholds with high accuracy 248 and consistency, with similar error and spread as models calibrated to corresponding non-normalized 249 recordings (Figure 7B). The aggregated calibrated models also predicted relative arrhythmia susceptibility of 250 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 17 the 4 model cells, as shown by a significant positive correlation between ground truth and predicted IKr block 251 thresholds in Figure 7C. Additionally, there was no correlation between the model cell’s true IKr block threshold 252 and its associated calibrated models’ prediction error, suggesting that the accuracy of these predictions would 253 not be affected by variability in arrhythmia susceptibility between iPSC-CM preparations (Figure 7D). 254 Optimized model calibration process constrains key ionic parameters 255 Throughout the protocol optimization process, we observed that while calibration error and spread largely 256 remained constant when averaged over all fitted parameters, certain individual parameters consistently 257 showed low calibration spread (e.g. GKr, GNaK) or high calibration spread (e.g. Grel, GCaT) (Supplementary 258 Figure S3). Thus, we hypothesized that our calibration method only needs to tightly constrain a few key 259 parameters, tolerating inaccuracy or variation in less important parameters while still generating highly 260 predictive calibrated models. To determine whether normalized data from the optimized protocol constrains the 261 parameters that play the largest roles in determining AP and CaT morphology, we analyzed the relationship 262 between parameter calibration accuracy and consistency with the Kernik model’s sensitivity to parameter 263 changes. We used multivariable regression to determine how the 16 calibrated model parameters affect 264 relevant phenotypes such as AP duration, CaT amplitude, and IKr block threshold (Figure 8A) [30]. The 265 resulting model coefficients represent the magnitude and direction of the effect of changes in the 266 corresponding parameter on the modeled phenotype (Figure 8A-C). The most important parameters depended 267 on which model output was considered (Figure 8A-C), but parameters that appeared frequently included rapid 268 delayed rectifier K+ conductance (GKr), L-type Ca2+ conductance (GCaL), Na+ conductance (GNa), inward rectifier 269 K+ conductance (GK1), and Na+/K+ pump conductance (GNaK). 270 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 18 271 We then grouped the parameters into the highest and lowest regression coefficient magnitudes for each 272 phenotype, and assessed calibration errors and calibration spreads from the optimized calibration protocol 273 within those groups. The 4 parameters with highest regression coefficient magnitudes for CaT amplitude and 274 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 19 IKr block threshold showed significantly lower average calibration error and spread compared with the 4 275 parameters with the lowest regression coefficient magnitudes (Figure 8E-F). Accordingly, the optimized 276 calibration protocol displayed significant improvement over most other candidate protocols in predictions of 277 CaT amplitude and IKr block threshold. We saw similar results when examining calibration errors and spreads 278 in parameters grouped by APD90 sensitivity, though the difference in calibration errors was not statistically 279 significant (Figure 8D). These results suggest that our optimized computational pipeline does constrain the key 280 Kernik model conductance parameters needed to generate predictive models. 281 Validation of optimized computational pipeline on in vitro iPSC-CM recordings 282 After optimization of our computational pipeline using in silico data, where the ground truth model parameter 283 values were known, we assessed whether our pipeline could still generate predictive models using in vitro 284 data. We applied the computational pipeline to normalized fluorescence AP and CaT recordings from one 285 iPSC-CM cell line under different combinations of 3 varied conditions: 1) low buffer calcium (1.0 mM) with 1 Hz 286 pacing, 2) physiological buffer calcium (1.8 mM) with 1 Hz pacing, and 3) physiological buffer calcium with 1.25 287 Hz pacing. Figure 9A outlines the data processing and normalization procedures performed on the 288 fluorescence recordings prior to model calibration. Recordings from a separate condition (1.0 mM buffer 289 calcium with 2 Hz pacing) were left out from model calibration, so they could be used to independently 290 evaluate predictions generated by the calibrated model. 291 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 20 292 Calibrating the model conductances using data from only one of these conditions resulted in unconstrained 293 conductance values and highly variable predictions of response to the left-out validation data (Figure 9B). 294 Including data from 2 additional, mixed conditions significantly improved consistency of key conductance 295 parameter calibrations (Figure 9C). When these calibrated parameters were incorporated into the Kernik model 296 and the cell line’s response to a new condition was predicted, we found that models calibrated to recordings 297 from 3 conditions substantially outperformed models calibrated to only 1 condition. These results further 298 supported our findings from in silico protocol optimization, showing that data from fluorescence voltage and 299 calcium transients under 3 varied conditions can constrain key model parameters and generate predictive, cell 300 preparation-specific models. 301 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 21

302 While the application of iPSC-derived cardiomyocytes to cardiac disease and pharmacology research has 303 grown dramatically in the past decade, intrinsic and extrinsic variability between cell preparations still hinders 304 reproducibility and clinical translation of studies that use these cells. Generic out-of-box mathematical models 305 of iPSC-CMs do not reflect this variability and potentially generate inaccurate predictions when used to 306 simulate pharmacological ion channel effects or genetic variations [10]. In this study, we aimed to overcome 307 these limitations of the iPSC-CM in vitro platform and their mathematical models by creating a computational 308 pipeline which can rapidly identify the cell preparation-specific ionic properties contributing to this phenotypic 309 variability. The parameter calibration algorithm we used, a heuristic genetic algorithm, allows for selection of 310 specific parameters of interest or parameter value limits, and is suitable for parallel computing. To determine 311 the feasibility of this approach, we used a model-generated in silico dataset to evaluate candidate protocols 312 containing varied data types and experimental conditions. We found an optimized protocol for fluorescence AP 313 and CaT recordings, consisting of varied buffer calcium concentration, electrical pacing, and ICaL block 314 conditions, which provided enough information to calibrate accurate and predictive cell-specific models, while 315 remaining short in duration and straightforward to acquire. Subsequent model calibrations with an in vitro 316 dataset suggested potential flexibility in which specific conditions are selected for the protocol. In our 317 computational pipeline development and subsequent validation, these data were able to inform key ionic 318 contributors, generating calibrated models which accurately and consistently predicted cell-specific responses 319 to IKr block. 320 Calibration of models of cardiac electrophysiology 321 Parameter values in computational models of cardiac electrophysiology are often determined by manual fits to 322 voltage and current recordings, often collected under one or few experimental conditions. This limits the 323 baseline models’ ability to represent electrophysiological heterogeneity. Additionally, determination of these 324 baseline parameter values may be affected by selection bias for cells with larger currents, inadequate 325 separation of activation and inactivation kinetics, and inter-laboratory differences in current and voltage clamp 326 protocols [31]. These issues have, in many cases, led researchers to recalibrate the model parameters in 327 attempts to better reflect the behaviors of their specific cardiomyocyte preparations. For example, Potse et al. 328 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 22 fitted a reaction-diffusion model of ventricular electrical activity to reproduce phenotypes of patients with heart 329 failure or left branch bundle block [32]. Similarly, Lombardo et al. tuned atrial electrophysiology models for 330 patients undergoing ablation therapy [33] whereas Krogh-Madsen et al. used the genetic algorithm to fit a 331 human ventricular cardiac model to clinical QT interval data [34]. In all three cases, these groups found 332 significant discrepancies between the respective initial models and their final, calibrated, fit-for-purpose 333 models. However, to the best of our knowledge, the experimental data used in virtually all previous 334 cardiomyocyte model calibration research come from voltage clamp, patch clamp, or clinical data. Currently, 335 these data are still difficult to obtain from a large group of patients or cell preparations. More recent work has 336 demonstrated the utility of genetic algorithms and the resulting recalibrated models in prediction of arrhythmic 337 behaviors and drug mechanisms, but these also used complex voltage step protocols to calibrate model 338 parameters [35]. These prior results motivated our search for protocols that could calibrate mathematical 339 models using voltage- and calcium-sensitive fluorescent dye measurements, which are considerably easier to 340 acquire than patch-clamp recordings. 341 Guiding experimental design for optimal parameter calibration 342 We considered several factors when optimizing the experimental protocol to generate iPSC-CM data for our 343 computational pipeline: 1) the information that the data would provide for parameter calibration, 2) the 344 complexity of the protocol, and 3) the feasibility of the experiment for broad accessibility and high throughput 345 setups. Since the goal of model calibration is usually to improve the accuracy of model parameterizations and 346 predictions, previous work has largely focused on designing experiments to maximize the information content 347 of the resulting data [36–38]. In particular, for cardiac electrophysiology models, model calibration has almost 348 universally been performed on data from complex voltage step protocols, in order to accurately determine 349 individual ion channel conductances and kinetics [10,39,40]. Here, we emphasize that it is also important to 350 identify adaptable protocols that can be widely adopted. Therefore, we opted to focus on using fluorescence 351 recordings to calibrate model parameters, instead of the usual techniques such as patch clamp. To evaluate 352 different protocols, we created an in silico dataset for testing various data combinations in their ability to inform 353 parameter calibration. This approach provided 2 major advantages. First, guided by experimental feasibility 354 and relevance, we were able to simulate iPSC-CM electrophysiology under many combinations of conditions 355 without the time and financial costs associated with in vitro experiments. Second, we knew the ground truth 356 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 23 ionic parameter values that produced these AP and CaT traces, so we could precisely assess which simulated 357 conditions informed which parameters. We also ran multiple rounds of genetic algorithms on each tested in 358 silico protocol, with each round starting from a different initial model population, to assess how consistently 359 each parameter was identified. We used this as a measure of confidence in each parameter estimate. 360 In their analysis on model parameter constraint using regression and Bayesian approaches, Sarkar and Sobie 361 found that, in both cases, parameter constraint improved when values for more than one output feature were 362 provided [41]. Similarly, we found that calibrating model parameters to fluorescence recordings of AP and CaT 363 traces, as opposed to using only one of the two, resulted in better parameter constraint for some (but not all) 364 parameters. Perhaps more importantly, we found that calibrating to both outputs simultaneously resulted in 365 significant improvement of the calibrated models’ predictions on a new, independent output. This is consistent 366 with prior results which showed that, when translating drug responses across cell types, using multiple 367 electrophysiological features or observations of iPSC-CMs under multiple conditions improved translations 368 [42–44]. Similarly, our final optimized protocol included a mixture of 3 experimental conditions. Finally, overall 369 parameter calibration accuracy and consistency, as well as independent predictions of IKr block response, were 370 similar between the models fitted to normalized data (representing fluorescence data) and those fitted to 371 corresponding non-normalized data (representing patch clamp and calibrated Ca2+ transient data). Collectively, 372 these findings support our hypothesis that data from fluorescence recordings under multiple conditions, which 373 are more practical compared with microelectrode or patch clamp recordings, can be used to calibrate 374 predictive models. Notably, several parameters such as the background Ca2+ and Na+ conductances, GKs, and 375 Grel were rarely constrained by any of the datasets we tested. Our subsequent parameter sensitivity analysis 376 showed that model conductances that strongly affect relevant model outputs tend to be well-constrained by our 377 optimized protocol. 378 Prediction of cell-specific arrhythmia susceptibility 379 Finally, we demonstrated that our optimized computational pipeline can generate consistent and predictive 380 mathematical models from in vitro iPSC-CM data, using fluorescence voltage and calcium recordings under a 381 similar set of experimental conditions. Many therapeutics have the potential to induce patient-specific 382 arrhythmic effects through their intended targets (e.g. quinidine blockade of INa and IKr) or unintended effects 383 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 24 (e.g. antibiotic off-target blockade of K+ channels) [45,46]. Previous analyses of systems biology models have 384 also argued the importance of a prediction-focused computational modeling approach. In an influential paper, 385 Gutenkunst et al. posited that a model that accurately predicts relevant phenotypes even with some parameter 386 value inaccuracies is more useful than a model with high parameter accuracy and consistency, but inaccurate 387 predictions [47]. This suggests that attempting to determine parameter values directly through precise 388 measurements is an inefficient way of optimizing models. Subsequent studies in both cell signaling and cardiac 389 electrophysiology models have demonstrated success in focusing parameter calibrations on optimizing 390 prediction accuracy in place of specific parameter accuracy [11,48,49]. These analyses demonstrate the 391 importance of calibrating models to not only fit available experimental data, but also generate accurate 392 predictions of responses to novel perturbations. 393 Therefore, while optimizing the computational pipeline, we included a prediction metric in addition to parameter 394 value accuracy and consistency. We created a quantifiable, pharmacologically-relevant measure of cell-395 specific arrhythmia susceptibility by finding the lowest level of IKr block that induced arrhythmic dynamics in 396 each Kernik model (in silico) or iPSC-CM preparation (in vitro), which we termed the “IKr block tolerance 397 threshold”. Calibrated models from our optimized computational pipeline were able to predict this threshold 398 consistently and accurately, with similar error as those from calibrating to the original, non-normalized data. 399 The models calibrated to normalized data were also able to correctly rank the IKr block arrhythmia 400 susceptibilities of the fitted Kernik model cells, suggesting that the models generated by this pipeline contain 401 enough information to predict cell-specific arrhythmia susceptibilities, despite errors and non-constraint of 402 several smaller conductance parameters. 403

and potential directions 404 The model calibration process developed and optimized in this study can be applied broadly to research that 405 investigates the ionic mechanisms behind iPSC-CM phenotype variability or predicts differential effects of 406 therapeutics on various iPSC-CM cell lines and preparations. However, phenotypic variability also exists within 407 iPSC-CMs from the same preparation. While it is possible to use fluorescent dyes to record electrophysiology 408 of individual cells, we currently focused our efforts on calibrating models to represent multicellular preparations. 409 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 25 Future work could include adapting this computational pipeline to calibrate models representing single 410 cardiomyocytes within an iPSC-CM preparation. 411 Several conductance parameters fitted in our current pipeline did not contribute significantly to 412 electrophysiological phenotype and responses (e.g. Gf, GSRleak, and Grel), according to our parameter sensitivity 413 analyses. Many of these same parameters were also less accurate and less constrained by our computational 414 pipeline. While we showed that the resulting calibrated models could accurately predict several major 415 electrophysiology phenotypes and channel block responses, we do not know whether these models could still 416 predict responses that are strongly affected by one or more of the unconstrained parameters. Calibrated 417 models generated by our computational pipeline should be evaluated on their ability to predict cell-specific 418 susceptibility to other pro- or anti-arrhythmia triggers of interest, such as ICaL or buffer [K+] changes. If 419 improving parameter calibration accuracy and constraint of particular ionic current(s) is found to be necessary 420 for model predictive power, one option would be to replace the conductances with lowest parameter sensitivity 421 regression coefficient magnitudes with representative kinetic parameters, such as time constants or activation 422 gates. If these kinetic parameters are found to impact electrophysiology in a manner that is detectable by the 423 genetic algorithm, their calibrations will likely be more accurate and constrained, generating more precise cell 424 preparation-specific iPSC-CM models. 425

426 The optimized model calibration computational pipeline we developed in this study can be applied to ongoing 427 pharmacological and disease studies to understand cell response variability in iPSC-CMs and advance 428 precision medicine. Our pipeline simultaneously fits multiple model parameters to fluorescence voltage and 429 calcium recordings, which are quicker and simpler to acquire than the patch clamp or other more involved 430 electrophysiological techniques used in prior model calibration work. The ability to use fluorescence recordings 431 increases this pipeline’s accessibility and throughput compared to prior methods. This unique characteristic 432 makes our computational pipeline suitable for rapid generation of cell-specific models for therapeutic, disease, 433 or population studies. Comparing cell-specific model parameters between iPSC-CM cell lines or to a 434 benchmark, such as an adult cardiac model, would reveal ionic mechanisms behind this variability in iPSC-CM 435 maturation and phenotypes. In addition, the personalized models created by this computational pipeline can 436 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 26 simulate patient-specific drug and perturbation responses to quickly predict therapeutic efficacy or drug 437 cardiotoxicity. In future studies, these calibrated iPSC-CM models can inform translation of ionic properties and 438 physiological responses into adult cardiac models, which could, in turn, generate even more accurate, 439 clinically-relevant cardiac electrophysiology phenotype and drug response predictions. 440 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 27

441 Simulations of iPSC-CM electrophysiology 442 In silico simulations of iPSC-CM electrophysiology (e.g. membrane potential, calcium transients, changes in 443 maximal conductances and ionic currents) were carried out in MATLAB R2020b by solving the ordinary 444 differential equations defined in the Kernik et al. (2019) mathematical model of human iPSC-CMs [9], using the 445 initial conditions listed in Supplementary Table S1. In cases where a protocol included cardiomyocyte pacing, 446 electrical stimuli of 60 pA/pF were simulated for a duration of 1 ms at the specified intervals. All simulations 447 were run for 5 minutes, which was sufficient to achieve a steady state using the baseline model and 448 physiological, unperturbed conditions (Supplementary Table S3, protocol #19). The voltage and intracellular 449 calcium concentration from the final 5 seconds of each simulation were stored and included in the in silico 450 dataset. 451 In silico dataset for optimization of iPSC-CM experimental protocol 452 To simulate a heterogeneous set of iPSC-CM cell lines, a “population” of Kernik model cells was created by 453 drawing multiplier factors for each maximal conductance parameter from a log-normal distribution with mean 454 multiplier = 1, spread = 0.2. This population was then simulated under physiological, unperturbed conditions 455 (Supplementary Table S3, protocol 19), and filtered for cells which showed a spontaneous beating frequency 456 between 0.3 Hz and 1.0 Hz, nearly matching the automaticity rate range for human iPSC-CMs reported in [50]. 457 Out of the remaining model cells, 4 were randomly selected for simulation under the 19 different conditions, to 458 generate the final in silico dataset. The Kernik model’s baseline conductance parameter values and the scaled 459 conductances for these 4 model cells are listed in Supplementary Table S2. A list of the various conditions 460 these 4 cells were simulated under can be found in Supplementary Table S3. All simulations were run for a 5-461 minute time span (steady state under physiological buffer conditions, without pacing), and membrane potential 462 and intracellular calcium recordings from the final 5 seconds of these simulations were extracted and stored for 463 the in silico dataset. The built-in interp1 function in MATLAB was used to interpolate these data in equal time 464 steps (0.1 ms). To mimic processed data from fluorescence recordings, we also created “normalized” versions 465 of each dataset by scaling the data from a minimum of 0 to a maximum of 1. This dataset allowed for 466 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 28 comparison of conductance parameter estimates produced during model calibration against known ground 467 truth conductance values. 468 Parameter calibration using the genetic algorithm 469 The data used for calibrating model parameters included the steady state AP and CaT waves from the 470 previously-generated in silico dataset. Individual simulation runs are either included or excluded from the fitting 471 procedure to compare different experimental “candidate protocols”. Rapid delayed rectifier potassium current 472 (IKr) block conditions were left out from protocol optimization to use for conductance estimate validation. If the 473 candidate protocol consists of more than 1 condition, the final values of membrane potential, intracellular 474 calcium concentration, and other steady state parameters are used as the starting values for the subsequent 475 protocol condition. To determine the optimal protocol for accurate and consistent parameter identification, 476 Kernik model maximal conductance parameters were fitted to in silico from each candidate protocol simulated 477 for 4 of the dataset’s model cells, using the genetic algorithm (GA) [10,12]. Briefly, the GA creates a new 478 population of Kernik model cells, each randomly assigned a set of conductance parameter scale factors, and 479 simulates the candidate protocol for each of these model cells. Then, the mean squared error between 480 corresponding data points in the GA-generated trace and the trace from the in silico dataset is calculated. The 481 model cells with the lowest errors are retained, while higher-scoring cells have their parameter values altered 482 in the next algorithm iteration. Details about specific GA settings, such as population size and parameter 483 retention/alteration criteria, can be found in Supplementary Table S4. This process repeats for 20 iterations, 484 where each new population has a lower average error than the previous. The parameters of the model cell with 485 the lowest error from the final population are selected as the final calibrated parameter values from that run. 486 Parameter calibrations are conducted 10 times per candidate protocol per model cell, with each of the 10 GA 487 runs starting with a different set of initial model populations. The 10 sets of conductance estimates from each 488 GA run are evaluated on 1) accuracy to the ground truth parameter values of the model cells and 2) ability to 489 predict the model cell’s response to a perturbation that was not seen during parameter calibration, such as IKr 490 block. 491 In vitro iPSC-CM tissue culture and optical recordings 492 Human induced pluripotent stem cells (iPSCs) derived from a single cell line were differentiated into 493 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 29 cardiomyocytes (iPSC-CMs) by modulating canonical Wnt signaling [51]. The cardiomyocytes were enriched, 494 then combined with cardiac fibroblasts in a ratio of 90% iPSC-CMs to 10% fibroblasts within a collagen-fibrin 495 hydrogel solution (InvivoSciences, Inc., Madison, WI, USA). The engineered heart tissues (EHTs) thus formed 496 were maintained in a serum-free cardiac maintenance medium, supplemented with penicillin and streptomycin 497 (Thermo Fisher Scientific, Waltham, MA, USA), within 96-well micro culture plates (MC-96, InvivoSciences). 498 EHT maintenance was carried out as previously described in relevant literature [52]. Following a five-day 499 remodeling phase, the EHTs underwent further maturation with biphasic constant current electrical stimulations 500 at a frequency of 1Hz for 8 days. 501 For optical voltage and calcium transient recordings, the EHTs were loaded with Fluovolt (at a dilution of 1:500, 502 Thermo Fisher Scientific), or Cal-520 AM (also at 1:500, AAT Bioquest), along with PowerLoad (Thermo Fisher 503 Scientific), by incubating them for one hour in Tyrode’s solution. This solution was prepared with either 1.0 mM 504 or 1.8 mM [Ca2+]. Following the removal of the dyes, the solution was pH adjusted to 7.4 and warmed, and 505 Tyrode’s solution with the corresponding calcium concentration was reintroduced to aid in recovery from dye 506 loading stress over a 30-minute period. The EHTs were then paced at frequencies of 1.0 Hz, 1.25 Hz, or 2.0 507 Hz for five minutes. The steady-state membrane potential and intracellular calcium transients during this period 508 were recorded using a high-throughput fluorescence plate imager, FDSS/µCell (Hamamatsu Photonics K.K., 509 Japan), utilizing 470nm excitation and 540nm emission at a rate of 125 data points per second. The collected 510 data were analyzed with the iVSurfer™ software (InvivoSciences), specifically designed for high-throughput 511 waveform data analysis. 512 Processing of in vitro iPSC-CM recordings for the computational pipeline 513 Fluorescence voltage and calcium recordings were processed in MATLAB with baseline drift subtraction 514 (imerode function), median filtering to remove noise (medfilt1 function), and the same normalization steps as 515 the pseudo-data. Kernik model maximal conductance parameters were fitted to these data using the same 516 workflow as prior fittings to pseudo-data, with the 1.0mM buffer [Ca2+], 2.0Hz pacing data left out for validation.517 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 30

518 1. Karakikes I, Ameen M, Termglinchan V, Wu JC. Human Induced Pluripotent Stem Cell-Derived 519 Cardiomyocytes: Insights into Molecular, Cellular, and Functional Phenotypes. Circ Res. 2015;117: 80–88. 520 doi:10.1161/CIRCRESAHA.117.305365 521 2. Huang C-Y, Liu C-L, Ting C-Y, Chiu Y-T, Cheng Y-C, Nicholson MW, et al. Human iPSC banking: barriers 522 and opportunities. J Biomed Sci. 2019;26: 87. doi:10.1186/s12929-019-0578-x 523 3. Hnatiuk AP, Briganti F, Staudt DW, Mercola M. Human iPSC modeling of heart disease for drug 524 development. Cell Chem Biol. 2021;28: 271–282. doi:10.1016/j.chembiol.2021.02.016 525 4. Yang J, Argenziano MA, Burgos Angulo M, Bertalovitz A, Beidokhti MN, McDonald TV. Phenotypic 526 Variability in iPSC-Induced Cardiomyocytes and Cardiac Fibroblasts Carrying Diverse LMNA Mutations. 527 Front Physiol. 2021;12: 778982. doi:10.3389/fphys.2021.778982 528 5. Nakagawa M, Ooie T, Ou B, Ichinose M, Takahashi N, Hara M, et al. Gender differences in autonomic 529 modulation of ventricular repolarization in humans. J Cardiovasc Electrophysiol. 2005;16: 278–284. 530 doi:10.1046/j.1540-8167.2005.40455.x 531 6. Paci M, Hyttinen J, Aalto-Setälä K, Severi S. Computational models of ventricular- and atrial-like human 532 induced pluripotent stem cell derived cardiomyocytes. Ann Biomed Eng. 2013;41: 2334–2348. 533 doi:10.1007/s10439-013-0833-3 534 7. Koivumäki JT, Naumenko N, Tuomainen T, Takalo J, Oksanen M, Puttonen KA, et al. Structural Immaturity 535 of Human iPSC-Derived Cardiomyocytes: In Silico Investigation of Effects on Function and Disease 536 Modeling. Front Physiol. 2018;9: 80. doi:10.3389/fphys.2018.00080 537 8. Paci M, Pölönen R-P, Cori D, Penttinen K, Aalto-Setälä K, Severi S, et al. Automatic Optimization of an in 538 Silico Model of Human iPSC Derived Cardiomyocytes Recapitulating Calcium Handling Abnormalities. 539 Front Physiol. 2018;9: 709. doi:10.3389/fphys.2018.00709 540 9. Kernik DC, Morotti S, Wu H, Garg P, Duff HJ, Kurokawa J, et al. A computational model of induced 541 pluripotent stem-cell derived cardiomyocytes incorporating experimental variability from multiple data 542 sources. J Physiol. 2019;597: 4533–4564. doi:10.1113/JP277724 543 10. Groenendaal W, Ortega FA, Kherlopian AR, Zygmunt AC, Krogh-Madsen T, Christini DJ. Cell-Specific 544 Cardiac Electrophysiology Models. PLoS Comput Biol. 2015;11: e1004242. 545 doi:10.1371/journal.pcbi.1004242 546 11. Dokos S, Lovell NH. Parameter estimation in cardiac ionic models. Prog Biophys Mol Biol. 2004;85: 407–547 431. doi:10.1016/j.pbiomolbio.2004.02.002 548 12. Syed Z, Vigmond E, Nattel S, Leon LJ. Atrial cell action potential parameter fitting using genetic algorithms. 549 Med Biol Eng Comput. 2005;43: 561–571. doi:10.1007/BF02351029 550 13. Devenyi RA, Ortega FA, Groenendaal W, Krogh-Madsen T, Christini DJ, Sobie EA. Differential roles of two 551 delayed rectifier potassium currents in regulation of ventricular action potential duration and arrhythmia 552 susceptibility. J Physiol. 2017;595: 2301–2317. doi:10.1113/JP273191 553 14. Akwaboah AD, Tsevi B, Yamlome P, Treat JA, Brucal-Hallare M, Cordeiro JM, et al. An in silico hiPSC-554 Derived Cardiomyocyte Model Built With Genetic Algorithm. Front Physiol. 2021;12: 675867. 555 doi:10.3389/fphys.2021.675867 556 15. Beattie KA, Hill AP, Bardenet R, Cui Y, Vandenberg JI, Gavaghan DJ, et al. Sinusoidal voltage protocols 557 for rapid characterisation of ion channel kinetics. J Physiol. 2018;596: 1813–1828. doi:10.1113/JP275733 558 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 31 16. Ahmed RE, Anzai T, Chanthra N, Uosaki H. A Brief Review of Current Maturation Methods for Human 559 Induced Pluripotent Stem Cells-Derived Cardiomyocytes. Front Cell Dev Biol. 2020;8. Available: 560 https://www.frontiersin.org/article/10.3389/fcell.2020.00178 561 17. Karbassi E, Fenix A, Marchiano S, Muraoka N, Nakamura K, Yang X, et al. Cardiomyocyte maturation: 562 advances in knowledge and implications for regenerative medicine. Nat Rev Cardiol. 2020;17: 341–359. 563 doi:10.1038/s41569-019-0331-x 564 18. James EC, Tomaskovic-Crook E, Crook JM. Bioengineering Clinically Relevant Cardiomyocytes and 565 Cardiac Tissues from Pluripotent Stem Cells. Int J Mol Sci. 2021;22: 3005. doi:10.3390/ijms22063005 566 19. Guo Y, Pu W. Cardiomyocyte Maturation: New Phase in Development. Circ Res. 2020;126: 1086–1106. 567 doi:10.1161/CIRCRESAHA.119.315862 568 20. Mannhardt I, Saleem U, Mosqueira D, Loos MF, Ulmer BM, Lemoine MD, et al. Comparison of 10 Control 569 hPSC Lines for Drug Screening in an Engineered Heart Tissue Format. Stem Cell Rep. 2020;15: 983–998. 570 doi:10.1016/j.stemcr.2020.09.002 571 21. Lauschke K, Volpini L, Liu Y, Vinggaard AM, Hall VJ. A Comparative Assessment of Marker Expression 572 Between Cardiomyocyte Differentiation of Human Induced Pluripotent Stem Cells and the Developing Pig 573 Heart. Stem Cells Dev. 2021;30: 374–385. doi:10.1089/scd.2020.0184 574 22. Yang X, Pabon L, Murry CE. Engineering Adolescence: Maturation of Human Pluripotent Stem Cell-575 derived Cardiomyocytes. Circ Res. 2014;114: 511–523. doi:10.1161/CIRCRESAHA.114.300558 576 23. Yoshida Y, Yamanaka S. Induced Pluripotent Stem Cells 10 Years Later. Circ Res. 2017;120: 1958–1968. 577 doi:10.1161/CIRCRESAHA.117.311080 578 24. Kannan S, Farid M, Lin BL, Miyamoto M, Kwon C. Transcriptomic entropy benchmarks stem cell-derived 579 cardiomyocyte maturation against endogenous tissue at single cell level. PLoS Comput Biol. 2021;17: 580 e1009305. doi:10.1371/journal.pcbi.1009305 581 25. Grancharova T, Gerbin KA, Rosenberg AB, Roco CM, Arakaki JE, DeLizo CM, et al. A comprehensive 582 analysis of gene expression changes in a high replicate and open-source dataset of differentiating hiPSC-583 derived cardiomyocytes. Sci Rep. 2021;11: 15845. doi:10.1038/s41598-021-94732-1 584 26. Wang F-S, Chen L-H. Heuristic Optimization. In: Dubitzky W, Wolkenhauer O, Cho K-H, Yokota H, editors. 585 Encyclopedia of Systems Biology. New York, NY: Springer; 2013. pp. 885–885. doi:10.1007/978-1-4419-586 9863-7_411 587 27. Carusi A, Burrage K, Rodríguez B. Bridging experiments, models and simulations: an integrative approach 588 to validation in computational cardiac electrophysiology. Am J Physiol-Heart Circ Physiol. 2012;303: H144–589 H155. doi:10.1152/ajpheart.01151.2011 590 28. Grynkiewicz G, Poenie M, Tsien RY. A new generation of Ca2+ indicators with greatly improved 591 fluorescence properties. J Biol Chem. 1985;260: 3440–3450. doi:10.1016/S0021-9258(19)83641-4 592 29. Bräuner T, Hülser DF, Strasser RJ. Comparative measurements of membrane potentials with 593 microelectrodes and voltage-sensitive dyes. Biochim Biophys Acta BBA - Biomembr. 1984;771: 208–216. 594 doi:10.1016/0005-2736(84)90535-2 595 30. Sobie EA. Parameter Sensitivity Analysis in Electrophysiological Models Using Multivariable Regression. 596 Biophys J. 2009;96: 1264–1274. doi:10.1016/j.bpj.2008.10.056 597 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 32 31. Krogh-Madsen T, Sobie EA, Christini DJ. Improving cardiomyocyte model fidelity and utility via dynamic 598 electrophysiology protocols and optimization algorithms. J Physiol. 2016;594: 2525–2536. 599 doi:10.1113/JP270618 600 32. Potse M, Krause D, Kroon W, Murzilli R, Muzzarelli S, Regoli F, et al. Patient-specific modelling of cardiac 601 electrophysiology in heart-failure patients. Europace. 2014;16: iv56–iv61. doi:10.1093/europace/euu257 602 33. Lombardo DM, Fenton FH, Narayan SM, Rappel W-J. Comparison of Detailed and Simplified Models of 603 Human Atrial Myocytes to Recapitulate Patient Specific Properties. PLOS Comput Biol. 2016;12: 604 e1005060. doi:10.1371/journal.pcbi.1005060 605 34. Krogh-Madsen T, Jacobson AF, Ortega FA, Christini DJ. Global Optimization of Ventricular Myocyte Model 606 to Multi-Variable Objective Improves Predictions of Drug-Induced Torsades de Pointes. Front Physiol. 607 2017;8. Available: https://www.frontiersin.org/articles/10.3389/fphys.2017.01059 608 35. Clark AP, Wei S, Kalola D, Krogh-Madsen T, Christini DJ. An in silico–in vitro pipeline for drug 609 cardiotoxicity screening identifies ionic pro-arrhythmia mechanisms. Br J Pharmacol. 2022;179: 4829–610 4843. doi:10.1111/bph.15915 611 36. Liepe J, Filippi S, Komorowski M, Stumpf MPH. Maximizing the Information Content of Experiments in 612 Systems Biology. PLOS Comput Biol. 2013;9: e1002888. doi:10.1371/journal.pcbi.1002888 613 37. Chaloner K, Verdinelli I. Bayesian Experimental Design: A Review. Stat Sci. 1995;10: 273–304. 614 doi:10.1214/ss/1177009939 615 38. Ryan EG, Drovandi CC, McGree JM, Pettitt AN. A Review of Modern Computational Algorithms for 616 Bayesian Optimal Design. Int Stat Rev. 2016;84: 128–154. doi:10.1111/insr.12107 617 39. Whittaker DG, Clerx M, Lei CL, Christini DJ, Mirams GR. Calibration of ionic and cellular cardiac 618 electrophysiology models. Wiley Interdiscip Rev Syst Biol Med. 2020;12: e1482. doi:10.1002/wsbm.1482 619 40. Lei CL, Clerx M, Beattie KA, Melgari D, Hancox JC, Gavaghan DJ, et al. Rapid Characterization of hERG 620 Channel Kinetics II: Temperature Dependence. Biophys J. 2019;117: 2455–2470. 621 doi:10.1016/j.bpj.2019.07.030 622 41. Sarkar AX, Sobie EA. Regression Analysis for Constraining Free Parameters in Electrophysiological 623 Models of Cardiac Cells. PLOS Comput Biol. 2010;6: e1000914. doi:10.1371/journal.pcbi.1000914 624 42. Gong JQX, Sobie EA. Population-based mechanistic modeling allows for quantitative predictions of drug 625 responses across cell types. Npj Syst Biol Appl. 2018;4: 1–11. doi:10.1038/s41540-018-0047-2 626 43. Tveito A, Jæger KH, Huebsch N, Charrez B, Edwards AG, Wall S, et al. Inversion and computational 627 maturation of drug response using human stem cell derived cardiomyocytes in microphysiological systems. 628 Sci Rep. 2018;8: 17626. doi:10.1038/s41598-018-35858-7 629 44. Morotti S, Liu C, Hegyi B, Ni H, Fogli Iseppe A, Wang L, et al. Quantitative cross-species translators of 630 cardiac myocyte electrophysiology: Model training, experimental validation, and applications. Sci Adv. 631 2021;7: eabg0927. doi:10.1126/sciadv.abg0927 632 45. Kannankeril PJ, Roden DM, Norris KJ, Whalen SP, George AL, Murray KT. Genetic susceptibility to 633 acquired long QT syndrome: Pharmacologic challenge in first-degree relatives. Heart Rhythm. 2005;2: 634 134–140. doi:10.1016/j.hrthm.2004.10.039 635 46. Zequn Z, Yujia W, Dingding Q, Jiangfang L. Off-label use of chloroquine, hydroxychloroquine, azithromycin 636 and lopinavir/ritonavir in COVID-19 risks prolonging the QT interval by targeting the hERG channel. Eur J 637 Pharmacol. 2021;893: 173813. doi:10.1016/j.ejphar.2020.173813 638 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint 33 47. Gutenkunst RN, Waterfall JJ, Casey FP, Brown KS, Myers CR, Sethna JP. Universally Sloppy Parameter 639 Sensitivities in Systems Biology Models. PLoS Comput Biol. 2007;3: e189. 640 doi:10.1371/journal.pcbi.0030189 641 48. Brown KS, Hill CC, Calero GA, Myers CR, Lee KH, Sethna JP, et al. The statistical mechanics of complex 642 signaling networks: nerve growth factor signaling. Phys Biol. 2004;1: 184. doi:10.1088/1478-3967/1/3/006 643 49. Casey FP, Baird D, Feng Q, Gutenkunst RN, Waterfall JJ, Myers CR, et al. Optimal experimental design in 644 an epidermal growth factor receptor signalling and down-regulation model. IET Syst Biol. 2007;1: 190–202. 645 doi:10.1049/iet-syb:20060065 646 50. Forny C, Sube R, Ertel EA. Contractions of Human-iPSC-derived Cardiomyocyte Syncytia Measured with a 647 Ca-sensitive Fluorescent Dye in Temperature-controlled 384-well Plates. J Vis Exp JoVE. 2018; 58290. 648 doi:10.3791/58290 649 51. Lian X, Hsiao C, Wilson G, Zhu K, Hazeltine LB, Azarin SM, et al. Robust cardiomyocyte differentiation 650 from human pluripotent stem cells via temporal modulation of canonical Wnt signaling. Proc Natl Acad Sci 651 U S A. 2012;109: E1848–E1857. doi:10.1073/pnas.1200250109 652 52. Greenberg MJ, Daily NJ, Wang A, Conway MK, Wakatsuki T. Genetic and Tissue Engineering Approaches 653 to Modeling the Mechanics of Human Heart Failure for Drug Discovery. Front Cardiovasc Med. 2018;5. 654 Available: https://www.frontiersin.org/articles/10.3389/fcvm.2018.00120 655 656 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint