{"paper_id":"40a62802-7256-41b4-bca3-8f6230ec29be","body_text":"1\n 1 \n 2 \n 3 \n 4 \nCreating cell-specific computational models of stem cell-derived cardiomyocytes 5 \nusing optical experiments  6 \n 7 \n 8 \nJanice Yang,1 Neil Daily,2 Taylor K. Pullinger,1 Tetsuro Wakatsuki,2 Eric A. Sobie1,*  9 \n 10 \n 11 \n1Department of Pharmacological Sciences & Graduate School of Biomedical Sciences, Icahn School 12 \nof Medicine at Mount Sinai, New York, NY, USA  13 \n2InvivoSciences Inc., Madison, WI 53719, USA. 14 \n 15 \nShort Title: Optimizing parameter identification algorithms  16 \n 17 \n*Corresponding author: eric.sobie@mssm.edu 18 \n 19 \n 20 \nSources of Funding 21 \nSupported by National Institutes of Health grants U01 HL 136297 and R01 HL 167520 (EAS) and 22 \nR44 HL 139248 (TW). JY was partially supported by NIH training grant T32 HD 075735. The authors 23 \nacknowledge computational resources and staff expertise provided by Scientific Computing and Data 24 \nat Mount Sinai, supported by CTSA grant ULTR004419 from NCATS. 25 \n 26 \nAim: To generate predictive cell line-specific computational models of stem cell-derived 27 \ncardiomyocytes 28 \nFocus: Optimizing a computational pipeline for parameterizing models 29 \n 30 \n  31 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 2\nAbstract:  32 \nHuman induced pluripotent stem cell-derived cardiomyocytes (iPSC-CMs) have gained traction as a powerful 33 \nmodel in cardiac disease and therapeutics research, since iPSCs are self-renewing and can be derived from 34 \nhealthy and diseased patients without invasive surgery. However, current iPSC-CM differentiation methods 35 \nproduce cardiomyocytes with immature, fetal-like electrophysiological phenotypes, and the variety of 36 \nmaturation protocols in the literature results in phenotypic differences between labs. Heterogeneity of iPSC 37 \ndonor genetic backgrounds contributes to additional phenotypic variability. Several mathematical models of 38 \niPSC-CM electrophysiology have been developed to help understand the ionic underpinnings of, and to 39 \nsimulate, various cell responses, but these models individually do not capture the phenotypic variability 40 \nobserved in iPSC-CMs. Here, we tackle these limitations by developing a computational pipeline to calibrate 41 \ncell preparation-specific iPSC-CM electrophysiological parameters.  42 \nWe used the genetic algorithm (GA), a heuristic parameter calibration method, to tune ion channel parameters 43 \nin a mathematical model of iPSC-CM physiology. To systematically optimize an experimental protocol that 44 \ngenerates sufficient data for parameter calibration, we created simulated datasets by applying various 45 \nprotocols to a population of in silico cells with known conductance variations, and we fitted to those datasets. 46 \nWe found that calibrating models to voltage and calcium transient data under 3 varied experimental conditions, 47 \nincluding electrical pacing combined with ion channel blockade and changing buffer ion concentrations, 48 \nimproved model parameter estimates and model predictions of unseen channel block responses. This 49 \nobservation held regardless of whether the fitted data were normalized, suggesting that normalized 50 \nfluorescence recordings, which are more accessible and higher throughput than patch clamp recordings, could 51 \nsufficiently inform conductance parameters. Therefore, this computational pipeline can be applied to different 52 \niPSC-CM preparations to determine cell line-specific ion channel properties and understand the mechanisms 53 \nbehind variability in perturbation responses.  54 \n  55 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 3\nAuthor Summary:  56 \nMany drug treatments or environmental factors can trigger cardiac arrhythmias, which are dangerous and often 57 \nunpredictable. Human cardiomyocytes derived from donor stem cells have proven to be a promising model for 58 \nstudying these events, but variability in donor genetic background and cell maturation methods, as well as 59 \noverall immaturity of stem cell-derived cardiomyocytes relative to the adult heart, have hindered reproducibility 60 \nand reliability of these studies. Mathematical models of these cells can aid in understanding the underlying 61 \nelectrophysiological contributors to this variability, but determining these models’ parameters for multiple cell 62 \npreparations is challenging. In this study, we tackle these limitations by developing a computational method to 63 \nsimultaneously estimate multiple model parameters using data from imaging-based experiments, which can be 64 \neasily scaled to rapidly characterize multiple cell lines. This method can generate many personalized models of 65 \nindividual cell preparations, improving drug response predictions and revealing specific differences in 66 \nelectrophysiological properties that contribute to variability in cardiac maturity and arrhythmia susceptibility.67 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 4\nABBREVIATIONS 68 \n AP: Action potential 69 \n CaT: Calcium transient 70 \n GA: Genetic algorithm 71 \n Gx, Ix, Jx: Maximal conductance (G), current density (I), or flux (J) for x, where x can be: Na (Na+), f 72 \n(funny Na+), CaL (L-type Ca2+), to (transient outward K+), Ks (slow delayed rectifier K+), Kr (rapid 73 \ndelayed rectifier K+), K1 (inward rectifier K+), PMCA (plasma membrane Ca2+ ATPase),  bNa 74 \n(background Na+), bCa (background Ca2+), Up (SR Ca2+ uptake), rel (SR Ca2+ release), NCX (Na+/Ca2+ 75 \nexchanger), NaK (Na+/K+ ATPase), SRleak (SR leak), or CaT (T-type Ca2+) 76 \n iPSC-CM: Induced pluripotent stem cell-derived cardiomyocyte 77 \n LTCC: L-type calcium channel 78 \n NCX: Sodium-calcium exchanger 79 \n SERCA: Sarco/endoplasmic reticulum calcium ATPase 80 \n SR: Sarcoplasmic reticulum 81 \n  82 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 5\nGLOSSARY 83 \n Model/parameter calibration: tuning one or more parameters in the computational model so that the 84 \nmodel output more closely matches experimental data 85 \n Experiment/protocol optimization: the process of determining what type and amount of data is 86 \nsufficient but also feasible for our model calibration goals 87 \no Protocol conditions – buffer calcium, potassium, or sodium concentrations; addition or removal 88 \nof stimulus; pacing rates; channel block; etc.  89 \no Protocol length(?) – number of protocol conditions  90 \no Protocol data type(?) –AP, CaT, or both; normalized or non-normalized data 91 \n Model prediction: using the calibrated computational model to simulate response to new (unseen) 92 \nconditions, drugs, or perturbations (in our case, IKr block) 93 \no Papers on independent validation/prediction 94 \nComputational pipeline: the full process of iPSC-CM computational model calibration; Includes iPSC-CM 95 \ndata acquisition/simulation -> data processing -> parameter calibration using genetic algorithm -> validation of 96 \ncalibrated models on an unseen condition (i.e. evaluating model predictions) 97 \n  98 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 6\nINTRODUCTION 99 \nHuman iPSC-derived cardiomyocytes (iPSC-CMs) are derived from patient cells that have been reprogrammed 100 \ninto pluripotency and subsequently manipulated to differentiate into the cardiac cell lineage. The development 101 \nof this in vitro platform marked a technological breakthrough in cardiac research, and iPSC-CMs are now 102 \nwidely-used in cardiac pharmacology and disease research. iPSCs can be sourced through minimally-invasive 103 \nprocedures, such as skin biopsy or blood draw, and maintained or banked for long periods [1,2]. These 104 \nproperties of iPSC-CMs make them an ideal platform for pharmacological studies and personalized disease 105 \nmodeling [3]. However, electrophysiological variability between iPSC-CM preparations from different cell lines 106 \nand differentiation methods limit the potential of this platform [4,5].  107 \nMathematical models of cardiomyocyte electrophysiology can provide valuable insight into iPSC-CM 108 \nelectrophysiological variability. These models contain parameters describing the ionic currents and fluxes that 109 \ncontribute to the properties of cardiac action potential (AP) and calcium transient (CaT) waveforms. Several 110 \nmodels of iPSC-CM physiology exist in the literature [6–8], with the most recent and comprehensive model 111 \npublished by Kernik et al. in 2019, hereon referred to as the “Kernik model” [9]. However, the baseline 112 \nparameter values in these models are not representative of the electrophysiological heterogeneity observed in 113 \niPSC-CMs. Additionally, parameter calibration has often relied on separate patch clamp measurements of each 114 \ntype of current. This process is time-consuming, inaccessible to many research groups, and often captures the 115 \nphysiology of only the average behavior of the population of cells examined. Methods to simultaneously 116 \noptimize several or all model parameters using minimal data are under development [10–14]. Recent work on 117 \nguinea pig ventricular myocyte models [13] and human iPSC-CM models [14] have demonstrated the utility of 118 \nautomated tuning algorithms for improving model parameterization. For example, in Devenyi et al [13], an 119 \nadjusted model, constrained by data, predicted the effect of slow delayed rectifier K+ current (IKs) block more 120 \naccurately compared to baseline model. Analogous work, at the level of individual ion channels, has shown the 121 \nsuperiority of sinusoidal voltage clamp protocols, compared with traditional square-wave protocols, for 122 \ncalibration of channel kinetic parameters [15].  123 \nEven with recent advances in iPSC-CM maturation protocols [16–19], most iPSC-CMs still display embryonic- 124 \nor neonatal-like cardiac phenotypes with different AP and CaT shapes compared with adult myocytes. These 125 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 7\nphenotypic differences, and how they change temporally, can depend on subtle differences in cell culture 126 \nconditions and the genetic background of the cell donor [4,5,20], which likely contributes to the wide variation 127 \nobserved between studies [4,5]. No standardized procedure currently exists for characterizing iPSC-CM 128 \nheterogeneity; prior efforts have generally examined a handful of molecular markers [21], used subjective 129 \nobservations of physiology [22,23], or employed specialized methods such as single cell RNA sequencing 130 \n[24,25]. These factors, and the potential importance of iPSC-CMs as a research platform, highlight the need for 131 \nautomated methods to characterize the cell lines used in each study.  132 \nHere we combined these ideas of simultaneous calibration of multiple parameters and optimization of a 133 \nminimal experimental protocol for generating calibration data. We aimed to create digital twins of iPSC-CM cell 134 \npreparations from the Kernik iPSC-CM model, thereby connecting observed iPSC-CM electrophysiology with 135 \nmolecular function and mechanisms. We hypothesize that fluorescence readouts of iPSC-CM physiology under 136 \nvaried experimental conditions provide enough information to reveal cellular ion channel properties, which can 137 \nthen be incorporated into the iPSC-CM digital twins. To systematically evaluate various data types and 138 \nexperimental protocols in their ability to inform model parameters, we created a simulated, in silico dataset 139 \nfrom Kernik model variations with known parameter values. From this study, we developed a computational 140 \npipeline which includes: 1) an optimized protocol for fluorescence recordings of iPSC-CM preparations; 2) a 141 \ngenetic algorithm process for calibration of ion channel parameters in iPSC-CM computational model, using 142 \nthe experimental recordings; and 3) validation of the resulting calibrated models by evaluating their predictions 143 \non independent, yet physiologically-important, perturbations. We demonstrate the utility of our computational 144 \npipeline in generating iPSC-CM digital twins that capture ionic variability and predict cell-specific 145 \nelectrophysiological phenotypes, including drug-induced arrhythmia susceptibility.   146 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 8\nRESULTS 147 \nOur primary objective was to optimize an experimental protocol for generating sufficient data to calibrate an 148 \niPSC-CM mathematical model, taking experimental feasibility and throughput into account. To systematically 149 \nevaluate how different types of data and experimental conditions impact parameter calibration and model 150 \npredictions of responses to new conditions or drug treatments, we created a simulated in silico dataset, 151 \ngenerated from a population of Kernik models with random variations in their 16 maximal conductance 152 \nparameters (Figure 1A, Supplementary Table S2). We simulated AP and CaT generated by these model cells 153 \nunder various conditions (Supplementary Table S3), and then concatenated these data in various 154 \ncombinations based on the “candidate protocols” we wanted to evaluate. The 16 Kernik model conductance 155 \nparameters were then calibrated, using a genetic algorithm, to best match the steady state AP and CaT traces 156 \nfrom each candidate protocol. These fitted parameters were incorporated into the Kernik model, and the newly-157 \ncalibrated models were evaluated on their ability to predict the original model cell’s responses to new 158 \nconditions (outlined in methods and Figure 1B). Optimizing the experimental protocol using this in silico dataset 159 \nprovided 3 major advantages: 1) knowledge of the ground truth parameter values that generated the dataset, 160 \nso the accuracy of the calibrated parameters can be evaluated, 2) relatively fast and easy simulation of new 161 \nconditions from the same or additional model cells, in comparison to acquiring new in vitro data and cell lines, 162 \nand 3) generating corresponding data types through computational methods and data processing, forgoing the 163 \nneed to set up new experiments to collect each data type.  164 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 9\n 165 \nSimultaneously calibrating to AP and CaT data improves model predictions 166 \nFirst, we examined whether AP or CaT recordings alone provided sufficient information for conductance 167 \nparameter calibration and model predictions. We used simulated recordings from baseline physiological 168 \nconditions without pacing stimuli (spontaneous APs), and supplied to a genetic algorithm (GA) either 1) only 169 \nthe AP traces, 2) only the CaT traces, or 3) both (Figure 2A). The GA calibrates parameters by varying them to 170 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 10\ncreate a population of Kernik model cells, then iteratively selecting and modifying individual models that 171 \ngenerated AP and CaT traces that best matched the supplied data (outlined in Figure 1B ) [26]. We ran the 172 \ngenetic algorithm 10 times per input dataset, with each GA run starting on a different initial population of 173 \nmodels. This allowed us to assess whether each input dataset provided enough information to consistently 174 \nestimate parameter values, even when the initial search space varied. 175 \n 176 \nWe evaluated each candidate protocol on parameter calibration accuracy (calibration error) and consistency 177 \n(calibration spread). Calibration error was represented by the absolute values of calibrated value log-178 \nnormalized to the corresponding ground truth parameter value, and calibration spread was assessed by 179 \ncalculating the standard deviation (spread) in log-normalized calibrated parameter values between the 10 runs 180 \n(Figure 2B). After averaging over all 16 fitted parameters, we found minimal differences in calibration error (p = 181 \n0.11) or spread (p = 0.07) between fitting only voltage, only calcium, or both simultaneously ( Figure 2C). From 182 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 11\nthese initial tests, therefore, there appears to be little difference between experimental protocols in how well (or 183 \npoorly) each can identify model parameters. 184 \nWe next assessed the calibrated models from each candidate protocol on how well each predicted model cell 185 \nresponses to 30% block of rapid delayed rectifier current IKr,as this perturbation is both: 1) independent of the 186 \nconditions used for model calibration and 2) relevant to cardiac pharmacology [27]. With this evaluation (Figure 187 \n3A), we found that models fitted to only AP data or only CaT data differed significantly from the ground truth 188 \nmodel (Figure 3B, Supplementary Figure S1). When both AP and CaT recordings were used simultaneously 189 \nduring parameter calibration, the resulting calibrated models exhibited substantial improvement in their 190 \npredictions of 30% IKr block. Therefore, fitting parameters to both AP and CaT traces simultaneously is 191 \nrequired for accurate prediction of results that were not used in the fitting process. 192 \n 193 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 12\nOverall parameter error and spread show minimal change with varied protocol conditions 194 \nNext, we explored how varying the number and identity of simulated conditions in the candidate protocols 195 \naffected parameter calibration. AP and CaT recordings were simulated under many conditions that could be 196 \nproduced in an electrophysiology laboratory and constructed candidate protocols from various combinations of 197 \nthese conditions, which included: 1) varied buffer calcium, from hypo- to normal to hyper-calcemic conditions, 198 \n2) spontaneous beating and varied pacing rates from 1-2 Hz, 3) varied levels of L-type calcium channel (ICaL) 199 \nblockade, and 4) a mixture of these conditions (see Methods for details). We also separated these data into 200 \nshorter protocols, to see if the experimental setup could be further simplified (Figure 4A). Unexpectedly, we 201 \nfound that parameter errors and spread across all 16 maximal conductance parameters did not differ 202 \nsignificantly between the protocols with varied cell culture conditions (Figure 4B, p = 0.814 and p = 0.673, 203 \nrespectively). Additionally, we observed no significant differences in average calibration errors (p = 0.088) and 204 \nspreads (p = 0.137) when varying the number of conditions in the protocol (Figure 4C). Similar findings were 205 \nobserved when evaluating different combinations of these experimental conditions (Supplementary Figure S2). 206 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 13\n 207 \nA readily-achievable 3-condition protocol improves predictions of IKr block response 208 \nAs with the results shown in Figure 3, we assessed how well calibrated models predicted the response to IKr 209 \nblock, and we found that accuracy depended strongly on the experimental conditions simulated in each 210 \ncandidate protocol (sample AP traces shown in Figure 5A). Out of the protocols using 3 experimental 211 \nconditions in Figure 4A, the mixed-conditions protocol generated calibrated models with the most accurate 212 \npredictions of changes in APs with 30% IKr block, compared with protocols which varied only the buffer 213 \ncalcium, pacing rates, or ICaL blockade (Figure 5B-C). When evaluating whether a shorter version of this 214 \nprotocol could produce equally predictive models, the full protocol with all 3 conditions outperformed protocols 215 \nwith 1 or 2 of the conditions (Figure 5B-C). These data suggest that an optimal experimental protocol for iPSC-216 \nCM model parameter calibration should include voltage and calcium transient fluorescence recordings under at 217 \nleast 3 varieties of cell culture condition changes or perturbations. Of note, this supports our hypothesis that 218 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 14\ndata from iPSC-CMs under a variety of conditions provides sufficient parameter identification to generate 219 \npredictive models for pharmacological applications. 220 \n 221 \nNormalized fluorescence recordings sufficiently inform parameter calibration 222 \nFluorescence recordings generally detect relative changes and do not provide absolute levels of voltage and 223 \ncalcium. Fluorescent dyes can be calibrated [28,29], and patch clamp can detect true transmembrane 224 \npotentials, but these are more challenging and lower-throughput techniques, compared with straightforward 225 \nfluorescence recordings. To assess the potential implications for parameter identification, we compared two 226 \nversions of simulated recordings obtained with our optimized protocol: 1) unscaled (original) in silico data, 227 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 15\nrepresenting patch clamp and calibrated calcium measurements, and 2) normalized versions of the same 228 \nwaveforms, representing simpler fluorescent recordings (Figure 6A).  229 \n 230 \nParameter calibration accuracy and consistency did not differ significantly when comparing calibrated models 231 \nfrom the original data with those from the corresponding normalized data (Figure 6B). Notably, calibrated 232 \nmodels from these two data types also performed similarly when evaluating predictions of 30% I Kr block 233 \nresponse (Figure 6C). This suggests that fluorescence voltage and calcium recordings from our optimized 234 \nprotocol conditions provide sufficient information to calibrate predictive models.  235 \nCalibrated models predict variability in arrhythmia susceptibility in silico 236 \nThe primary goal of optimizing a computational pipeline for model calibration is to create cell preparation-237 \nspecific models that can predict variability in phenotypes, particularly arrhythmia susceptibility, between iPSC-238 \nCM lines. To assess the ability of our optimized protocol to predict cell line-specific arrhythmia susceptibility, 239 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 16\nwe calculated the lowest level of IKr block (i.e. highest % IKr) that induced arrhythmia dynamics 240 \n(afterdepolarizations, Torsades de Pointes, alternans, beating cessation, or tachycardia) in the same 4 Kernik 241 \nmodel cells from our in silico dataset. These simulations defined, for each of the 4 model cells, an “IKr block 242 \ntolerance threshold” (Figure 7A), which ranged from 37% IKr block (highest susceptibility) to 62% IKr block 243 \n(highest tolerance) (Figure 7C) across the 4 model cells. 244 \n 245 \nPredicted IKr block tolerance thresholds were also calculated for the calibrated models generated by our 246 \ncomputational pipeline. These predicted thresholds were assessed against the ground truth threshold of the 247 \nmodel cell. Calibrated models from our computational pipeline predicted IKr block thresholds with high accuracy 248 \nand consistency, with similar error and spread as models calibrated to corresponding non-normalized 249 \nrecordings (Figure 7B). The aggregated calibrated models also predicted relative arrhythmia susceptibility of 250 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 17\nthe 4 model cells, as shown by a significant positive correlation between ground truth and predicted IKr block 251 \nthresholds in Figure 7C. Additionally, there was no correlation between the model cell’s true IKr block threshold 252 \nand its associated calibrated models’ prediction error, suggesting that the accuracy of these predictions would 253 \nnot be affected by variability in arrhythmia susceptibility between iPSC-CM preparations (Figure 7D). 254 \nOptimized model calibration process constrains key ionic parameters  255 \nThroughout the protocol optimization process, we observed that while calibration error and spread largely 256 \nremained constant when averaged over all fitted parameters, certain individual parameters consistently 257 \nshowed low calibration spread (e.g. GKr, GNaK) or high calibration spread (e.g. Grel, GCaT) (Supplementary 258 \nFigure S3). Thus, we hypothesized that our calibration method only needs to tightly constrain a few key 259 \nparameters, tolerating inaccuracy or variation in less important parameters while still generating highly 260 \npredictive calibrated models. To determine whether normalized data from the optimized protocol constrains the 261 \nparameters that play the largest roles in determining AP and CaT morphology, we analyzed the relationship 262 \nbetween parameter calibration accuracy and consistency with the Kernik model’s sensitivity to parameter 263 \nchanges. We used multivariable regression to determine how the 16 calibrated model parameters affect 264 \nrelevant phenotypes such as AP duration, CaT amplitude, and IKr block threshold (Figure 8A) [30]. The 265 \nresulting model coefficients represent the magnitude and direction of the effect of changes in the 266 \ncorresponding parameter on the modeled phenotype (Figure 8A-C). The most important parameters depended 267 \non which model output was considered (Figure 8A-C), but parameters that appeared frequently included  rapid 268 \ndelayed rectifier K+ conductance (GKr), L-type Ca2+ conductance (GCaL), Na+ conductance (GNa), inward rectifier 269 \nK+ conductance (GK1), and Na+/K+ pump conductance (GNaK).   270 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 18\n 271 \nWe then grouped the parameters into the highest and lowest regression coefficient magnitudes for each 272 \nphenotype, and assessed calibration errors and calibration spreads from the optimized calibration protocol 273 \nwithin those groups. The 4 parameters with highest regression coefficient magnitudes for CaT amplitude and 274 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 19\nIKr block threshold showed significantly lower average calibration error and spread compared with the 4 275 \nparameters with the lowest regression coefficient magnitudes (Figure 8E-F). Accordingly, the optimized 276 \ncalibration protocol displayed significant improvement over most other candidate protocols in predictions of 277 \nCaT amplitude and IKr block threshold. We saw similar results when examining calibration errors and spreads 278 \nin parameters grouped by APD90 sensitivity, though the difference in calibration errors was not statistically 279 \nsignificant (Figure 8D). These results suggest that our optimized computational pipeline does constrain the key 280 \nKernik model conductance parameters needed to generate predictive models. 281 \nValidation of optimized computational pipeline on in vitro iPSC-CM recordings 282 \nAfter optimization of our computational pipeline using in silico data, where the ground truth model parameter 283 \nvalues were known, we assessed whether our pipeline could still generate predictive models using in vitro 284 \ndata. We applied the computational pipeline to normalized fluorescence AP and CaT recordings from one 285 \niPSC-CM cell line under different combinations of 3 varied conditions: 1) low buffer calcium (1.0 mM) with 1 Hz 286 \npacing, 2) physiological buffer calcium (1.8 mM) with 1 Hz pacing, and 3) physiological buffer calcium with 1.25 287 \nHz pacing. Figure 9A outlines the data processing and normalization procedures performed on the 288 \nfluorescence recordings prior to model calibration. Recordings from a separate condition (1.0 mM buffer 289 \ncalcium with 2 Hz pacing) were left out from model calibration, so they could be used to independently 290 \nevaluate predictions generated by the calibrated model.  291 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 20\n 292 \nCalibrating the model conductances using data from only one of these conditions resulted in unconstrained 293 \nconductance values and highly variable predictions of response to the left-out validation data (Figure 9B). 294 \nIncluding data from 2 additional, mixed conditions significantly improved consistency of key conductance 295 \nparameter calibrations (Figure 9C). When these calibrated parameters were incorporated into the Kernik model 296 \nand the cell line’s response to a new condition was predicted, we found that models calibrated to recordings 297 \nfrom 3 conditions substantially outperformed models calibrated to only 1 condition. These results further 298 \nsupported our findings from in silico protocol optimization, showing that data from fluorescence voltage and 299 \ncalcium transients under 3 varied conditions can constrain key model parameters and generate predictive, cell 300 \npreparation-specific models.   301 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 21\nDISCUSSION 302 \nWhile the application of iPSC-derived cardiomyocytes to cardiac disease and pharmacology research has 303 \ngrown dramatically in the past decade, intrinsic and extrinsic variability between cell preparations still hinders 304 \nreproducibility and clinical translation of studies that use these cells. Generic out-of-box mathematical models 305 \nof iPSC-CMs do not reflect this variability and potentially generate inaccurate predictions when used to 306 \nsimulate pharmacological ion channel effects or genetic variations [10]. In this study, we aimed to overcome 307 \nthese limitations of the iPSC-CM in vitro platform and their mathematical models by creating a computational 308 \npipeline which can rapidly identify the cell preparation-specific ionic properties contributing to this phenotypic 309 \nvariability. The parameter calibration algorithm we used, a heuristic genetic algorithm, allows for selection of 310 \nspecific parameters of interest or parameter value limits, and is suitable for parallel computing. To determine 311 \nthe feasibility of this approach, we used a model-generated in silico dataset to evaluate candidate protocols 312 \ncontaining varied data types and experimental conditions. We found an optimized protocol for fluorescence AP 313 \nand CaT recordings, consisting of varied buffer calcium concentration, electrical pacing, and ICaL block 314 \nconditions, which provided enough information to calibrate accurate and predictive cell-specific models, while 315 \nremaining short in duration and straightforward to acquire. Subsequent model calibrations with an in vitro 316 \ndataset suggested potential flexibility in which specific conditions are selected for the protocol. In our 317 \ncomputational pipeline development and subsequent validation, these data were able to inform key ionic 318 \ncontributors, generating calibrated models which accurately and consistently predicted cell-specific responses 319 \nto IKr block. 320 \nCalibration of models of cardiac electrophysiology 321 \nParameter values in computational models of cardiac electrophysiology are often determined by manual fits to 322 \nvoltage and current recordings, often collected under one or few experimental conditions. This limits the 323 \nbaseline models’ ability to represent electrophysiological heterogeneity. Additionally, determination of these 324 \nbaseline parameter values may be affected by selection bias for cells with larger currents, inadequate 325 \nseparation of activation and inactivation kinetics, and inter-laboratory differences in current and voltage clamp 326 \nprotocols [31]. These issues have, in many cases, led researchers to recalibrate the model parameters in 327 \nattempts to better reflect the behaviors of their specific cardiomyocyte preparations. For example, Potse et al. 328 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 22\nfitted a reaction-diffusion model of ventricular electrical activity to reproduce phenotypes of patients with heart 329 \nfailure or left branch bundle block [32]. Similarly, Lombardo et al. tuned atrial electrophysiology models for 330 \npatients undergoing ablation therapy [33] whereas Krogh-Madsen et al. used the genetic algorithm to fit a 331 \nhuman ventricular cardiac model to clinical QT interval data [34]. In all three cases, these groups found 332 \nsignificant discrepancies between the respective initial models and their final, calibrated, fit-for-purpose 333 \nmodels. However, to the best of our knowledge, the experimental data used in virtually all previous 334 \ncardiomyocyte model calibration research come from voltage clamp, patch clamp, or clinical data. Currently, 335 \nthese data are still difficult to obtain from a large group of patients or cell preparations. More recent work has 336 \ndemonstrated the utility of genetic algorithms and the resulting recalibrated models in prediction of arrhythmic 337 \nbehaviors and drug mechanisms, but these also used complex voltage step protocols to calibrate model 338 \nparameters [35]. These prior results motivated our search for protocols that could calibrate mathematical 339 \nmodels using voltage- and calcium-sensitive fluorescent dye measurements, which are considerably easier to 340 \nacquire than patch-clamp recordings.  341 \nGuiding experimental design for optimal parameter calibration 342 \nWe considered several factors when optimizing the experimental protocol to generate iPSC-CM data for our 343 \ncomputational pipeline: 1) the information that the data would provide for parameter calibration, 2) the 344 \ncomplexity of the protocol, and 3) the feasibility of the experiment for broad accessibility and high throughput 345 \nsetups. Since the goal of model calibration is usually to improve the accuracy of model parameterizations and 346 \npredictions, previous work has largely focused on designing experiments to maximize the information content 347 \nof the resulting data [36–38]. In particular, for cardiac electrophysiology models, model calibration has almost 348 \nuniversally been performed on data from complex voltage step protocols, in order to accurately determine 349 \nindividual ion channel conductances and kinetics [10,39,40]. Here, we emphasize that it is also important to 350 \nidentify adaptable protocols that can be widely adopted. Therefore, we opted to focus on using fluorescence 351 \nrecordings to calibrate model parameters, instead of the usual techniques such as patch clamp. To evaluate 352 \ndifferent protocols, we created an in silico dataset for testing various data combinations in their ability to inform 353 \nparameter calibration. This approach provided 2 major advantages. First, guided by experimental feasibility 354 \nand relevance, we were able to simulate iPSC-CM electrophysiology under many combinations of conditions 355 \nwithout the time and financial costs associated with in vitro experiments. Second, we knew the ground truth 356 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 23\nionic parameter values that produced these AP and CaT traces, so we could precisely assess which simulated 357 \nconditions informed which parameters. We also ran multiple rounds of genetic algorithms on each tested in 358 \nsilico protocol, with each round starting from a different initial model population, to assess how consistently 359 \neach parameter was identified. We used this as a measure of confidence in each parameter estimate.  360 \nIn their analysis on model parameter constraint using regression and Bayesian approaches, Sarkar and Sobie 361 \nfound that, in both cases, parameter constraint improved when values for more than one output feature were 362 \nprovided [41]. Similarly, we found that calibrating model parameters to fluorescence recordings of AP and CaT 363 \ntraces, as opposed to using only one of the two, resulted in better parameter constraint for some (but not all) 364 \nparameters. Perhaps more importantly, we found that calibrating to both outputs simultaneously resulted in 365 \nsignificant improvement of the calibrated models’ predictions on a new, independent output. This is consistent 366 \nwith prior results which showed that, when translating drug responses across cell types, using multiple 367 \nelectrophysiological features or observations of iPSC-CMs under multiple conditions improved translations 368 \n[42–44]. Similarly, our final optimized protocol included a mixture of 3 experimental conditions. Finally, overall 369 \nparameter calibration accuracy and consistency, as well as independent predictions of IKr block response, were 370 \nsimilar between the models fitted to normalized data (representing fluorescence data) and those fitted to 371 \ncorresponding non-normalized data (representing patch clamp and calibrated Ca2+ transient data). Collectively, 372 \nthese findings support our hypothesis that data from fluorescence recordings under multiple conditions, which 373 \nare more practical compared with microelectrode or patch clamp recordings, can be used to calibrate 374 \npredictive models. Notably, several parameters such as the background Ca2+ and Na+ conductances, GKs, and 375 \nGrel were rarely constrained by any of the datasets we tested. Our subsequent parameter sensitivity analysis 376 \nshowed that model conductances that strongly affect relevant model outputs tend to be well-constrained by our 377 \noptimized protocol.  378 \nPrediction of cell-specific arrhythmia susceptibility 379 \nFinally, we demonstrated that our optimized computational pipeline can generate consistent and predictive 380 \nmathematical models from in vitro iPSC-CM data, using fluorescence voltage and calcium recordings under a 381 \nsimilar set of experimental conditions. Many therapeutics have the potential to induce patient-specific 382 \narrhythmic effects through their intended targets (e.g. quinidine blockade of INa and IKr) or unintended effects 383 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 24\n(e.g. antibiotic off-target blockade of K+ channels) [45,46]. Previous analyses of systems biology models have 384 \nalso argued the importance of a prediction-focused computational modeling approach. In an influential paper, 385 \nGutenkunst et al. posited that a model that accurately predicts relevant phenotypes even with some parameter 386 \nvalue inaccuracies is more useful than a model with high parameter accuracy and consistency, but inaccurate 387 \npredictions [47]. This suggests that attempting to determine parameter values directly through precise 388 \nmeasurements is an inefficient way of optimizing models. Subsequent studies in both cell signaling and cardiac 389 \nelectrophysiology models have demonstrated success in focusing parameter calibrations on optimizing 390 \nprediction accuracy in place of specific parameter accuracy [11,48,49]. These analyses demonstrate the 391 \nimportance of calibrating models to not only fit available experimental data, but also generate accurate 392 \npredictions of responses to novel perturbations. 393 \nTherefore, while optimizing the computational pipeline, we included a prediction metric in addition to parameter 394 \nvalue accuracy and consistency. We created a quantifiable, pharmacologically-relevant measure of cell-395 \nspecific arrhythmia susceptibility by finding the lowest level of IKr block that induced arrhythmic dynamics in 396 \neach Kernik model (in silico) or iPSC-CM preparation (in vitro), which we termed the “IKr block tolerance 397 \nthreshold”. Calibrated models from our optimized computational pipeline were able to predict this threshold 398 \nconsistently and accurately, with similar error as those from calibrating to the original, non-normalized data. 399 \nThe models calibrated to normalized data were also able to correctly rank the IKr block arrhythmia 400 \nsusceptibilities of the fitted Kernik model cells, suggesting that the models generated by this pipeline contain 401 \nenough information to predict cell-specific arrhythmia susceptibilities, despite errors and non-constraint of 402 \nseveral smaller conductance parameters.  403 \nLimitations and potential directions 404 \nThe model calibration process developed and optimized in this study can be applied broadly to research that 405 \ninvestigates the ionic mechanisms behind iPSC-CM phenotype variability or predicts differential effects of 406 \ntherapeutics on various iPSC-CM cell lines and preparations. However, phenotypic variability also exists within 407 \niPSC-CMs from the same preparation. While it is possible to use fluorescent dyes to record electrophysiology 408 \nof individual cells, we currently focused our efforts on calibrating models to represent multicellular preparations. 409 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 25\nFuture work could include adapting this computational pipeline to calibrate models representing single 410 \ncardiomyocytes within an iPSC-CM preparation.  411 \nSeveral conductance parameters fitted in our current pipeline did not contribute significantly to 412 \nelectrophysiological phenotype and responses (e.g. Gf, GSRleak, and Grel), according to our parameter sensitivity 413 \nanalyses. Many of these same parameters were also less accurate and less constrained by our computational 414 \npipeline. While we showed that the resulting calibrated models could accurately predict several major 415 \nelectrophysiology phenotypes and channel block responses, we do not know whether these models could still 416 \npredict responses that are strongly affected by one or more of the unconstrained parameters. Calibrated 417 \nmodels generated by our computational pipeline should be evaluated on their ability to predict cell-specific 418 \nsusceptibility to other pro- or anti-arrhythmia triggers of interest, such as ICaL or buffer [K+] changes. If 419 \nimproving parameter calibration accuracy and constraint of particular ionic current(s) is found to be necessary 420 \nfor model predictive power, one option would be to replace the conductances with lowest parameter sensitivity 421 \nregression coefficient magnitudes with representative kinetic parameters, such as time constants or activation 422 \ngates. If these kinetic parameters are found to impact electrophysiology in a manner that is detectable by the 423 \ngenetic algorithm, their calibrations will likely be more accurate and constrained, generating more precise cell 424 \npreparation-specific iPSC-CM models.  425 \nConclusions 426 \nThe optimized model calibration computational pipeline we developed in this study can be applied to ongoing 427 \npharmacological and disease studies to understand cell response variability in iPSC-CMs and advance 428 \nprecision medicine. Our pipeline simultaneously fits multiple model parameters to fluorescence voltage and 429 \ncalcium recordings, which are quicker and simpler to acquire than the patch clamp or other more involved 430 \nelectrophysiological techniques used in prior model calibration work. The ability to use fluorescence recordings 431 \nincreases this pipeline’s accessibility and throughput compared to prior methods. This unique characteristic 432 \nmakes our computational pipeline suitable for rapid generation of cell-specific models for therapeutic, disease, 433 \nor population studies. Comparing cell-specific model parameters between iPSC-CM cell lines or to a 434 \nbenchmark, such as an adult cardiac model, would reveal ionic mechanisms behind this variability in iPSC-CM 435 \nmaturation and phenotypes. In addition, the personalized models created by this computational pipeline can 436 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 26\nsimulate patient-specific drug and perturbation responses to quickly predict therapeutic efficacy or drug 437 \ncardiotoxicity. In future studies, these calibrated iPSC-CM models can inform translation of ionic properties and 438 \nphysiological responses into adult cardiac models, which could, in turn, generate even more accurate, 439 \nclinically-relevant cardiac electrophysiology phenotype and drug response predictions.   440 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 27\nMETHODS 441 \nSimulations of iPSC-CM electrophysiology 442 \nIn silico simulations of iPSC-CM electrophysiology (e.g. membrane potential, calcium transients, changes in 443 \nmaximal conductances and ionic currents) were carried out in MATLAB R2020b by solving the ordinary 444 \ndifferential equations defined in the Kernik et al. (2019) mathematical model of human iPSC-CMs [9], using the 445 \ninitial conditions listed in Supplementary Table S1. In cases where a protocol included cardiomyocyte pacing, 446 \nelectrical stimuli of 60 pA/pF were simulated for a duration of 1 ms at the specified intervals. All simulations 447 \nwere run for 5 minutes, which was sufficient to achieve a steady state using the baseline model and 448 \nphysiological, unperturbed conditions (Supplementary Table S3, protocol #19). The voltage and intracellular 449 \ncalcium concentration from the final 5 seconds of each simulation were stored and included in the in silico 450 \ndataset.  451 \nIn silico dataset for optimization of iPSC-CM experimental protocol 452 \nTo simulate a heterogeneous set of iPSC-CM cell lines, a “population” of Kernik model cells was created by 453 \ndrawing multiplier factors for each maximal conductance parameter from a log-normal distribution with mean 454 \nmultiplier = 1, spread = 0.2. This population was then simulated under physiological, unperturbed conditions 455 \n(Supplementary Table S3, protocol 19), and filtered for cells which showed a spontaneous beating frequency 456 \nbetween 0.3 Hz and 1.0 Hz, nearly matching the automaticity rate range for human iPSC-CMs reported in [50]. 457 \nOut of the remaining model cells, 4 were randomly selected for simulation under the 19 different conditions, to 458 \ngenerate the final in silico dataset. The Kernik model’s baseline conductance parameter values and the scaled 459 \nconductances for these 4 model cells are listed in Supplementary Table S2. A list of the various conditions 460 \nthese 4 cells were simulated under can be found in Supplementary Table S3. All simulations were run for a 5-461 \nminute time span (steady state under physiological buffer conditions, without pacing), and membrane potential 462 \nand intracellular calcium recordings from the final 5 seconds of these simulations were extracted and stored for 463 \nthe in silico dataset. The built-in interp1 function in MATLAB was used to interpolate these data in equal time 464 \nsteps (0.1 ms). To mimic processed data from fluorescence recordings, we also created “normalized” versions 465 \nof each dataset by scaling the data from a minimum of 0 to a maximum of 1. This dataset allowed for 466 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 28\ncomparison of conductance parameter estimates produced during model calibration against known ground 467 \ntruth conductance values.  468 \nParameter calibration using the genetic algorithm 469 \nThe data used for calibrating model parameters included the steady state AP and CaT waves from the 470 \npreviously-generated in silico dataset. Individual simulation runs are either included or excluded from the fitting 471 \nprocedure to compare different experimental “candidate protocols”. Rapid delayed rectifier potassium current 472 \n(IKr) block conditions were left out from protocol optimization to use for conductance estimate validation. If the 473 \ncandidate protocol consists of more than 1 condition, the final values of membrane potential, intracellular 474 \ncalcium concentration, and other steady state parameters are used as the starting values for the subsequent 475 \nprotocol condition. To determine the optimal protocol for accurate and consistent parameter identification, 476 \nKernik model maximal conductance parameters were fitted to in silico from each candidate protocol simulated 477 \nfor 4 of the dataset’s model cells, using the genetic algorithm (GA) [10,12]. Briefly, the GA creates a new 478 \npopulation of Kernik model cells, each randomly assigned a set of conductance parameter scale factors, and 479 \nsimulates the candidate protocol for each of these model cells. Then, the mean squared error between 480 \ncorresponding data points in the GA-generated trace and the trace from the in silico dataset is calculated. The 481 \nmodel cells with the lowest errors are retained, while higher-scoring cells have their parameter values altered 482 \nin the next algorithm iteration. Details about specific GA settings, such as population size and parameter 483 \nretention/alteration criteria, can be found in Supplementary Table S4. This process repeats for 20 iterations, 484 \nwhere each new population has a lower average error than the previous. The parameters of the model cell with 485 \nthe lowest error from the final population are selected as the final calibrated parameter values from that run. 486 \nParameter calibrations are conducted 10 times per candidate protocol per model cell, with each of the 10 GA 487 \nruns starting with a different set of initial model populations. The 10 sets of conductance estimates from each 488 \nGA run are evaluated on 1) accuracy to the ground truth parameter values of the model cells and 2) ability to 489 \npredict the model cell’s response to a perturbation that was not seen during parameter calibration, such as IKr 490 \nblock.  491 \nIn vitro iPSC-CM tissue culture and optical recordings 492 \nHuman induced pluripotent stem cells (iPSCs) derived from a single cell line were differentiated into 493 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 29\ncardiomyocytes (iPSC-CMs) by modulating canonical Wnt signaling [51]. The cardiomyocytes were enriched, 494 \nthen combined with cardiac fibroblasts in a ratio of 90% iPSC-CMs to 10% fibroblasts within a collagen-fibrin 495 \nhydrogel solution (InvivoSciences, Inc., Madison, WI, USA). The engineered heart tissues (EHTs) thus formed 496 \nwere maintained in a serum-free cardiac maintenance medium, supplemented with penicillin and streptomycin 497 \n(Thermo Fisher Scientific, Waltham, MA, USA), within 96-well micro culture plates (MC-96, InvivoSciences). 498 \nEHT maintenance was carried out as previously described in relevant literature [52]. Following a five-day 499 \nremodeling phase, the EHTs underwent further maturation with biphasic constant current electrical stimulations 500 \nat a frequency of 1Hz for 8 days.  501 \nFor optical voltage and calcium transient recordings, the EHTs were loaded with Fluovolt (at a dilution of 1:500, 502 \nThermo Fisher Scientific), or Cal-520 AM (also at 1:500, AAT Bioquest), along with PowerLoad (Thermo Fisher 503 \nScientific), by incubating them for one hour in Tyrode’s solution. This solution was prepared with either 1.0 mM 504 \nor 1.8 mM [Ca2+]. Following the removal of the dyes, the solution was pH adjusted to 7.4 and warmed, and 505 \nTyrode’s solution with the corresponding calcium concentration was reintroduced to aid in recovery from dye 506 \nloading stress over a 30-minute period. The EHTs were then paced at frequencies of 1.0 Hz, 1.25 Hz, or 2.0 507 \nHz for five minutes. The steady-state membrane potential and intracellular calcium transients during this period 508 \nwere recorded using a high-throughput fluorescence plate imager, FDSS/µCell (Hamamatsu Photonics K.K., 509 \nJapan), utilizing 470nm excitation and 540nm emission at a rate of 125 data points per second. The collected 510 \ndata were analyzed with the iVSurfer™ software (InvivoSciences), specifically designed for high-throughput 511 \nwaveform data analysis. 512 \nProcessing of in vitro iPSC-CM recordings for the computational pipeline 513 \nFluorescence voltage and calcium recordings were processed in MATLAB with baseline drift subtraction 514 \n(imerode function), median filtering to remove noise (medfilt1 function), and the same normalization steps as 515 \nthe pseudo-data. Kernik model maximal conductance parameters were fitted to these data using the same 516 \nworkflow as prior fittings to pseudo-data, with the 1.0mM buffer [Ca2+], 2.0Hz pacing data left out for validation.517 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 30\nREFERENCES 518 \n1.  Karakikes I, Ameen M, Termglinchan V, Wu JC. Human Induced Pluripotent Stem Cell-Derived 519 \nCardiomyocytes: Insights into Molecular, Cellular, and Functional Phenotypes. Circ Res. 2015;117: 80–88. 520 \ndoi:10.1161/CIRCRESAHA.117.305365 521 \n2.  Huang C-Y, Liu C-L, Ting C-Y, Chiu Y-T, Cheng Y-C, Nicholson MW, et al. Human iPSC banking: barriers 522 \nand opportunities. J Biomed Sci. 2019;26: 87. doi:10.1186/s12929-019-0578-x 523 \n3.  Hnatiuk AP, Briganti F, Staudt DW, Mercola M. Human iPSC modeling of heart disease for drug 524 \ndevelopment. Cell Chem Biol. 2021;28: 271–282. doi:10.1016/j.chembiol.2021.02.016 525 \n4.  Yang J, Argenziano MA, Burgos Angulo M, Bertalovitz A, Beidokhti MN, McDonald TV. Phenotypic 526 \nVariability in iPSC-Induced Cardiomyocytes and Cardiac Fibroblasts Carrying Diverse LMNA Mutations. 527 \nFront Physiol. 2021;12: 778982. doi:10.3389/fphys.2021.778982 528 \n5.  Nakagawa M, Ooie T, Ou B, Ichinose M, Takahashi N, Hara M, et al. Gender differences in autonomic 529 \nmodulation of ventricular repolarization in humans. J Cardiovasc Electrophysiol. 2005;16: 278–284. 530 \ndoi:10.1046/j.1540-8167.2005.40455.x 531 \n6.  Paci M, Hyttinen J, Aalto-Setälä K, Severi S. Computational models of ventricular- and atrial-like human 532 \ninduced pluripotent stem cell derived cardiomyocytes. Ann Biomed Eng. 2013;41: 2334–2348. 533 \ndoi:10.1007/s10439-013-0833-3 534 \n7.  Koivumäki JT, Naumenko N, Tuomainen T, Takalo J, Oksanen M, Puttonen KA, et al. Structural Immaturity 535 \nof Human iPSC-Derived Cardiomyocytes: In Silico Investigation of Effects on Function and Disease 536 \nModeling. Front Physiol. 2018;9: 80. doi:10.3389/fphys.2018.00080 537 \n8.  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 \nSilico Model of Human iPSC Derived Cardiomyocytes Recapitulating Calcium Handling Abnormalities. 539 \nFront Physiol. 2018;9: 709. doi:10.3389/fphys.2018.00709 540 \n9.  Kernik DC, Morotti S, Wu H, Garg P, Duff HJ, Kurokawa J, et al. A computational model of induced 541 \npluripotent stem-cell derived cardiomyocytes incorporating experimental variability from multiple data 542 \nsources. J Physiol. 2019;597: 4533–4564. doi:10.1113/JP277724 543 \n10.  Groenendaal W, Ortega FA, Kherlopian AR, Zygmunt AC, Krogh-Madsen T, Christini DJ. Cell-Specific 544 \nCardiac Electrophysiology Models. PLoS Comput Biol. 2015;11: e1004242. 545 \ndoi:10.1371/journal.pcbi.1004242 546 \n11.  Dokos S, Lovell NH. Parameter estimation in cardiac ionic models. Prog Biophys Mol Biol. 2004;85: 407–547 \n431. doi:10.1016/j.pbiomolbio.2004.02.002 548 \n12.  Syed Z, Vigmond E, Nattel S, Leon LJ. Atrial cell action potential parameter fitting using genetic algorithms. 549 \nMed Biol Eng Comput. 2005;43: 561–571. doi:10.1007/BF02351029 550 \n13.  Devenyi RA, Ortega FA, Groenendaal W, Krogh-Madsen T, Christini DJ, Sobie EA. Differential roles of two 551 \ndelayed rectifier potassium currents in regulation of ventricular action potential duration and arrhythmia 552 \nsusceptibility. J Physiol. 2017;595: 2301–2317. doi:10.1113/JP273191 553 \n14.  Akwaboah AD, Tsevi B, Yamlome P, Treat JA, Brucal-Hallare M, Cordeiro JM, et al. An in silico hiPSC-554 \nDerived Cardiomyocyte Model Built With Genetic Algorithm. Front Physiol. 2021;12: 675867. 555 \ndoi:10.3389/fphys.2021.675867 556 \n15.  Beattie KA, Hill AP, Bardenet R, Cui Y, Vandenberg JI, Gavaghan DJ, et al. Sinusoidal voltage protocols 557 \nfor rapid characterisation of ion channel kinetics. J Physiol. 2018;596: 1813–1828. doi:10.1113/JP275733 558 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 31\n16.  Ahmed RE, Anzai T, Chanthra N, Uosaki H. A Brief Review of Current Maturation Methods for Human 559 \nInduced Pluripotent Stem Cells-Derived Cardiomyocytes. Front Cell Dev Biol. 2020;8. Available: 560 \nhttps://www.frontiersin.org/article/10.3389/fcell.2020.00178 561 \n17.  Karbassi E, Fenix A, Marchiano S, Muraoka N, Nakamura K, Yang X, et al. Cardiomyocyte maturation: 562 \nadvances in knowledge and implications for regenerative medicine. Nat Rev Cardiol. 2020;17: 341–359. 563 \ndoi:10.1038/s41569-019-0331-x 564 \n18.  James EC, Tomaskovic-Crook E, Crook JM. Bioengineering Clinically Relevant Cardiomyocytes and 565 \nCardiac Tissues from Pluripotent Stem Cells. Int J Mol Sci. 2021;22: 3005. doi:10.3390/ijms22063005 566 \n19.  Guo Y, Pu W. Cardiomyocyte Maturation: New Phase in Development. Circ Res. 2020;126: 1086–1106. 567 \ndoi:10.1161/CIRCRESAHA.119.315862 568 \n20.  Mannhardt I, Saleem U, Mosqueira D, Loos MF, Ulmer BM, Lemoine MD, et al. Comparison of 10 Control 569 \nhPSC Lines for Drug Screening in an Engineered Heart Tissue Format. Stem Cell Rep. 2020;15: 983–998. 570 \ndoi:10.1016/j.stemcr.2020.09.002 571 \n21.  Lauschke K, Volpini L, Liu Y, Vinggaard AM, Hall VJ. A Comparative Assessment of Marker Expression 572 \nBetween Cardiomyocyte Differentiation of Human Induced Pluripotent Stem Cells and the Developing Pig 573 \nHeart. Stem Cells Dev. 2021;30: 374–385. doi:10.1089/scd.2020.0184 574 \n22.  Yang X, Pabon L, Murry CE. Engineering Adolescence: Maturation of Human Pluripotent Stem Cell-575 \nderived Cardiomyocytes. Circ Res. 2014;114: 511–523. doi:10.1161/CIRCRESAHA.114.300558 576 \n23.  Yoshida Y, Yamanaka S. Induced Pluripotent Stem Cells 10 Years Later. Circ Res. 2017;120: 1958–1968. 577 \ndoi:10.1161/CIRCRESAHA.117.311080 578 \n24.  Kannan S, Farid M, Lin BL, Miyamoto M, Kwon C. Transcriptomic entropy benchmarks stem cell-derived 579 \ncardiomyocyte maturation against endogenous tissue at single cell level. PLoS Comput Biol. 2021;17: 580 \ne1009305. doi:10.1371/journal.pcbi.1009305 581 \n25.  Grancharova T, Gerbin KA, Rosenberg AB, Roco CM, Arakaki JE, DeLizo CM, et al. A comprehensive 582 \nanalysis of gene expression changes in a high replicate and open-source dataset of differentiating hiPSC-583 \nderived cardiomyocytes. Sci Rep. 2021;11: 15845. doi:10.1038/s41598-021-94732-1 584 \n26.  Wang F-S, Chen L-H. Heuristic Optimization. In: Dubitzky W, Wolkenhauer O, Cho K-H, Yokota H, editors. 585 \nEncyclopedia of Systems Biology. New York, NY: Springer; 2013. pp. 885–885. doi:10.1007/978-1-4419-586 \n9863-7_411 587 \n27.  Carusi A, Burrage K, Rodríguez B. Bridging experiments, models and simulations: an integrative approach 588 \nto validation in computational cardiac electrophysiology. Am J Physiol-Heart Circ Physiol. 2012;303: H144–589 \nH155. doi:10.1152/ajpheart.01151.2011 590 \n28.  Grynkiewicz G, Poenie M, Tsien RY. A new generation of Ca2+ indicators with greatly improved 591 \nfluorescence properties. J Biol Chem. 1985;260: 3440–3450. doi:10.1016/S0021-9258(19)83641-4 592 \n29.  Bräuner T, Hülser DF, Strasser RJ. Comparative measurements of membrane potentials with 593 \nmicroelectrodes and voltage-sensitive dyes. Biochim Biophys Acta BBA - Biomembr. 1984;771: 208–216. 594 \ndoi:10.1016/0005-2736(84)90535-2 595 \n30.  Sobie EA. Parameter Sensitivity Analysis in Electrophysiological Models Using Multivariable Regression. 596 \nBiophys J. 2009;96: 1264–1274. doi:10.1016/j.bpj.2008.10.056 597 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 32\n31.  Krogh-Madsen T, Sobie EA, Christini DJ. Improving cardiomyocyte model fidelity and utility via dynamic 598 \nelectrophysiology protocols and optimization algorithms. J Physiol. 2016;594: 2525–2536. 599 \ndoi:10.1113/JP270618 600 \n32.  Potse M, Krause D, Kroon W, Murzilli R, Muzzarelli S, Regoli F, et al. Patient-specific modelling of cardiac 601 \nelectrophysiology in heart-failure patients. Europace. 2014;16: iv56–iv61. doi:10.1093/europace/euu257 602 \n33.  Lombardo DM, Fenton FH, Narayan SM, Rappel W-J. Comparison of Detailed and Simplified Models of 603 \nHuman Atrial Myocytes to Recapitulate Patient Specific Properties. PLOS Comput Biol. 2016;12: 604 \ne1005060. doi:10.1371/journal.pcbi.1005060 605 \n34.  Krogh-Madsen T, Jacobson AF, Ortega FA, Christini DJ. Global Optimization of Ventricular Myocyte Model 606 \nto Multi-Variable Objective Improves Predictions of Drug-Induced Torsades de Pointes. Front Physiol. 607 \n2017;8. Available: https://www.frontiersin.org/articles/10.3389/fphys.2017.01059 608 \n35.  Clark AP, Wei S, Kalola D, Krogh-Madsen T, Christini DJ. An in silico–in vitro pipeline for drug 609 \ncardiotoxicity screening identifies ionic pro-arrhythmia mechanisms. Br J Pharmacol. 2022;179: 4829–610 \n4843. doi:10.1111/bph.15915 611 \n36.  Liepe J, Filippi S, Komorowski M, Stumpf MPH. Maximizing the Information Content of Experiments in 612 \nSystems Biology. PLOS Comput Biol. 2013;9: e1002888. doi:10.1371/journal.pcbi.1002888 613 \n37.  Chaloner K, Verdinelli I. Bayesian Experimental Design: A Review. Stat Sci. 1995;10: 273–304. 614 \ndoi:10.1214/ss/1177009939 615 \n38.  Ryan EG, Drovandi CC, McGree JM, Pettitt AN. A Review of Modern Computational Algorithms for 616 \nBayesian Optimal Design. Int Stat Rev. 2016;84: 128–154. doi:10.1111/insr.12107 617 \n39.  Whittaker DG, Clerx M, Lei CL, Christini DJ, Mirams GR. Calibration of ionic and cellular cardiac 618 \nelectrophysiology models. Wiley Interdiscip Rev Syst Biol Med. 2020;12: e1482. doi:10.1002/wsbm.1482 619 \n40.  Lei CL, Clerx M, Beattie KA, Melgari D, Hancox JC, Gavaghan DJ, et al. Rapid Characterization of hERG 620 \nChannel Kinetics II: Temperature Dependence. Biophys J. 2019;117: 2455–2470. 621 \ndoi:10.1016/j.bpj.2019.07.030 622 \n41.  Sarkar AX, Sobie EA. Regression Analysis for Constraining Free Parameters in Electrophysiological 623 \nModels of Cardiac Cells. PLOS Comput Biol. 2010;6: e1000914. doi:10.1371/journal.pcbi.1000914 624 \n42.  Gong JQX, Sobie EA. Population-based mechanistic modeling allows for quantitative predictions of drug 625 \nresponses across cell types. Npj Syst Biol Appl. 2018;4: 1–11. doi:10.1038/s41540-018-0047-2 626 \n43.  Tveito A, Jæger KH, Huebsch N, Charrez B, Edwards AG, Wall S, et al. Inversion and computational 627 \nmaturation of drug response using human stem cell derived cardiomyocytes in microphysiological systems. 628 \nSci Rep. 2018;8: 17626. doi:10.1038/s41598-018-35858-7 629 \n44.  Morotti S, Liu C, Hegyi B, Ni H, Fogli Iseppe A, Wang L, et al. Quantitative cross-species translators of 630 \ncardiac myocyte electrophysiology: Model training, experimental validation, and applications. Sci Adv. 631 \n2021;7: eabg0927. doi:10.1126/sciadv.abg0927 632 \n45.  Kannankeril PJ, Roden DM, Norris KJ, Whalen SP, George AL, Murray KT. Genetic susceptibility to 633 \nacquired long QT syndrome: Pharmacologic challenge in first-degree relatives. Heart Rhythm. 2005;2: 634 \n134–140. doi:10.1016/j.hrthm.2004.10.039 635 \n46.  Zequn Z, Yujia W, Dingding Q, Jiangfang L. Off-label use of chloroquine, hydroxychloroquine, azithromycin 636 \nand lopinavir/ritonavir in COVID-19 risks prolonging the QT interval by targeting the hERG channel. Eur J 637 \nPharmacol. 2021;893: 173813. doi:10.1016/j.ejphar.2020.173813 638 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint \n\n \n 33\n47.  Gutenkunst RN, Waterfall JJ, Casey FP, Brown KS, Myers CR, Sethna JP. Universally Sloppy Parameter 639 \nSensitivities in Systems Biology Models. PLoS Comput Biol. 2007;3: e189. 640 \ndoi:10.1371/journal.pcbi.0030189 641 \n48.  Brown KS, Hill CC, Calero GA, Myers CR, Lee KH, Sethna JP, et al. The statistical mechanics of complex 642 \nsignaling networks: nerve growth factor signaling. Phys Biol. 2004;1: 184. doi:10.1088/1478-3967/1/3/006 643 \n49.  Casey FP, Baird D, Feng Q, Gutenkunst RN, Waterfall JJ, Myers CR, et al. Optimal experimental design in 644 \nan epidermal growth factor receptor signalling and down-regulation model. IET Syst Biol. 2007;1: 190–202. 645 \ndoi:10.1049/iet-syb:20060065 646 \n50.  Forny C, Sube R, Ertel EA. Contractions of Human-iPSC-derived Cardiomyocyte Syncytia Measured with a 647 \nCa-sensitive Fluorescent Dye in Temperature-controlled 384-well Plates. J Vis Exp JoVE. 2018; 58290. 648 \ndoi:10.3791/58290 649 \n51.  Lian X, Hsiao C, Wilson G, Zhu K, Hazeltine LB, Azarin SM, et al. Robust cardiomyocyte differentiation 650 \nfrom human pluripotent stem cells via temporal modulation of canonical Wnt signaling. Proc Natl Acad Sci 651 \nU S A. 2012;109: E1848–E1857. doi:10.1073/pnas.1200250109 652 \n52.  Greenberg MJ, Daily NJ, Wang A, Conway MK, Wakatsuki T. Genetic and Tissue Engineering Approaches 653 \nto Modeling the Mechanics of Human Heart Failure for Drug Discovery. Front Cardiovasc Med. 2018;5. 654 \nAvailable: https://www.frontiersin.org/articles/10.3389/fcvm.2018.00120 655 \n 656 \n.CC-BY 4.0 International licenseavailable under a \n(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 \nThe copyright holder for this preprintthis version posted January 8, 2024. ; https://doi.org/10.1101/2024.01.07.574577doi: bioRxiv preprint","source_license":"CC-BY-4.0","license_restricted":false}