Improved Computational Identification of Drug Response Using Optical Measurements of Human Stem Cell Derived Cardiomyocytes in Microphysiological Systems

Cardiomyocytes derived from human induced pluripotent stem cells (hiPSC-CMs) hold great potential for drug screening applications. However, their usefulness is limited by the relative immaturity of the cells’ electrophysiological properties as compared to native cardiomyocytes in the adult human heart. In this work, we extend and improve on methodology to address this limitation, building on previously introduced computational procedures which predict drug effects for adult cells based on changes in optical measurements of action potentials and Ca2+ transients made in stem cell derived cardiac microtissues. This methodology quantifies ion channel changes through the inversion of data into a mathematical model, and maps this response to an adult phenotype through the assumption of functional invariance of fundamental intracellular and membrane channels during maturation. Here, we utilize an updated action potential model to represent both hiPSC-CMs and adult cardiomyocytes, apply an IC50-based model of dose-dependent drug effects, and introduce a continuation-based optimization algorithm for analysis of dose escalation measurements using five drugs with known effects. The improved methodology can identify drug induced changes more efficiently, and quantitate important metrics such as IC50 in line with published values. Consequently, the updated methodology is a step towards employing computational procedures to elucidate drug effects in adult cardiomyocytes for new drugs using stem cell-derived experimental tissues.


INTRODUCTION
The development of human induced pluripotent stem cells (hiPSCs) opens promising avenues of investigation into a wide variety of fundamental questions in cell physiology and beyond [for recent reviews, see, e.g., (Yoshida and Yamanaka, 2017;Di Baldassarre et al., 2018;Ye et al., 2018)]. One of the more immediately tractable applications of hiPSCs is the creation of specific human cell and tissue samples to augment drug discovery and development pipelines. These pipelines have traditionally relied on animal models in key areas of testing, but are limited by significant physiological differences between animal and human cells [see, e.g., (Mathur et al., 2016;Fine and Vunjak-Novakovic, 2017;Yoshida and Yamanaka, 2017;Ye et al., 2018)]. These differences, both at the genetic and proteomic levels, give rise to distinctly non-human system dynamics, for example, a mouse's heart rate is much more rapid than a human's (∼600 bpm vs. ∼60 bpm), such that it is often difficult to translate drug effects from one species to another [see, e.g., (Mathur et al., 2016;Fine and Vunjak-Novakovic, 2017;Yoshida and Yamanaka, 2017;Ye et al., 2018)].
By using hiPSC-derived cells, it is possible to measure drug effects directly in human-based systems, and therapeutics can eventually, in principle, be tested and adjusted at the level of the individual patient. This hiPSC-based, patient-centric approach opens up great possibilities for drug development, both in terms of the scope of illnesses approachable, including disorders caused by rare mutations, as well as improved safety by the early identification of drug side effects in human cells. Nevertheless, hiPSCs are also associated with a variety of scientific challenges that must be resolved to realize the full potential of the technology [see, e.g., (Mathur et al., 2015;Mathur et al., 2016;Mora et al., 2017;Christensen et al., 2018;Ronaldson-Bouchard et al., 2018;Zhao et al., 2018)].
Maturity of generated cells and tissues is one of these key challenges, a prominent example being the maturation of hiPSCderived cardiomyocytes (hiPSC-CMs) (Chen et al., 2016). Human cardiomyocytes develop over many years [see (Hille, 2001), ch. 21], and during this period the density of specific ion channels changes significantly, due both to increased area of the cell membrane and proliferation of membrane channels [see, e.g., (Sontheimer et al., 1992;Moody and Bosma, 2005;Bedada et al., 2016)]. Therefore, the physiological response of immature hiPSC-CMs to a drug cannot necessarily be used to infer the properties of the drug, nor the response of adult human cardiomyocytes. Even if it is known exactly how a drug affects an hiPSC-CM, it is difficult to deduce its effect on adult cells; direct interpretation may in fact lead to both false positives and false negatives [see (Liang et al., 2013;Mathur et al., 2015)].
In , we used mathematical modeling of cardiac cell dynamics to address these challenges associated with the application of hiPSC-CMs. Such mathematical modeling of the cardiac action potential (AP) remains an active area of research, and sophisticated models have been developed in order to accurately simulate both single cells and cardiac tissue dynamics [see e.g., (CellML Model Repository; Rudy and Silva, 2006;Grandi et al., 2010;O'Hara et al., 2011;Rudy, 2012;Franzone et al., 2014;Qu et al., 2014;Edwards and Louch, 2017;Quarteroni et al., 2017;Tveito et al., 2017)]. We presented an algorithm for inverting experimental measurements of the membrane potential and the cytosolic calcium (Ca 2+ ) concentration in order to obtain parameters for a mathematical model of the hiPSC-CM AP. We then demonstrated how this model of hiPSC-CMs can be mapped to an AP model representing adult cells. We were able to estimate the effect of a drug on essential ion currents for hiPSC-CMs as based on measurements from a microphysiological system (Mathur et al., 2015), and then to map this effect onto the adult cardiomyocyte AP model. The combination of these two methods permitted, in principle, to deduce drug effects on adult human cardiomyocytes as based on measurements of hiPSC-CMs in a microphysiological system.
The overall method developed in ) is illustrated in Figure 1. In this procedure, we take optical measurements using fluorescent dyes in a microphysiological system to define relative traces of the membrane potential and cytosolic Ca 2+ concentration for cells under normal media conditions and in the presence of drugs. We then define a mathematical model for the control (undrugged) cases by identifying parameters denoted by p hiPSC,C (hiPSC is for hiPSC-derived, C is for control) in an AP model that matches the experimental waveforms. Using this model of hiPSC-CMs, we then define a maturation matrix Q such that Qp hiPSC,C = p A,B , where p A,B (A is for adult, B is for base) are known parameters representing a generic AP model of an adult human cardiomyocyte. Here, the matrix Q represents the developmental change in ion channel density and geometry from hiPSC-CMs to adult cardiomyocytes, independent of drug effects on individual channels.
Next, experimental traces of the membrane potential and cytosolic Ca 2+ concentration are taken for the same cells in the presence of drugs, and these traces are used to define a new parameter-vector p hiPSC,D (hiPSC-CM, Drug) that matches the data. This new parameterization gives us information about what modeled channels have been altered by the drug. Then, by assuming that the drug affects every individual ion channel in the same manner for the hiPSC-CMs and the adult cells, the parameter vector for the adult case is given by p A,D = Qp hiPSC,D . Hence, we can find an AP model for adult human cardiomyocytes under the influence of the drug, even though only the effect for hiPSC-CMs has been measured.
The present report aims to present a number of modifications to improve the accuracy and reliability of these methods. First, using the AP models of Grandi et al. (2010), O'Hara et al. (2011), Paci et al. (2013), Paci et al. (2015), Paci et al. (2017), and Paci et al. (2018) as a basis, we have derived a new AP model to improve representation of experimental data. As our inversion algorithm is based on conducting a huge number of simulations with varying parameter values, it is essential to have a model that is stable with respect to perturbations of the parameters.
Therefore, the new model is designed for improved stability. In particular, the model of the intracellular Ca 2+ dynamics has been modified to avoid instabilities in the balance between the influx and efflux of Ca 2+ to the sarcoplasmic reticulum (SR).
In addition, our aim has been to create models that can be mapped back and forth between hiPSC-CMs and adult cardiomyocytes. A vital modeling assumption has been that the individual channels are the same in these two cases, and that only channel density should change. However, existing AP models are not derived with such a mapping in mind, and models of identical single channel dynamics vary significantly among models. Therefore, we have derived a new AP model which strictly adheres to the principle that every current (and flux) should be written as a product of the ion channel density and the dynamics of a single channel; identical ion channels are represented by identical mathematical models. Consequently, the mathematical representation of a single channel is the same for the hiPSC-CMs and adult cardiomyocytes in the novel AP model presented here.
Finally, we present a new method for inverting experimental data into parameters for the AP model by introducing a continuation-based approach, searching for optimal parameters by gradually moving from known parameters to the parameters we want to identify. Continuation methods are well developed in scientific computing [see e.g., (Keller, 1987;Allgower and Georg, 2012)] and offer significant computational savings to find optimal solutions.
In this manuscript, we first motivate and describe the approaches outlined above. We then evaluate these methods with respect to accuracy using simulated data. Subsequently, the new methods are used to identify the effect of five wellcharacterized drugs based on optical measurements of hiPSC-CMs. In all cases considered, the predicted effects are consistent with known drug effects, lending credence to the principle that novel drug effects on adult cardiomyocytes could reliably be estimated using measurements of hiPSC-CMs and the described methodology.

METHODS
Here, we offer a detailed presentation of all steps illustrated in Figure 1. First, we present the derivation of a new AP model. FIGURE 1 | The effect of a drug on adult cardiomyocytes can be identified by the process illustrated. The cytosolic calcium concentration (Ca 2+ ) and the membrane potential (V) are measured in a microphysiological system (Mathur et al., 2015;Mathur et al., 2016) of hiPSC-CMs. Data are collected when no drug has been applied (control, C) and when a drug has been applied (D). The data are used to parameterize a model for both cases, represented by the parameter vectors p hiPSC,C and p hiPSC,D for the control and drugged cases, respectively. The control parameter vector p hiPSC,C is used to define the maturation matrix Q such that Qp hiPSC,C = p A,B , where p A,B is the parameter vector of a generic base model of adult human cardiomyocytes. By comparing the adult parameter vectors for the control and drugged cases, the effect of the drug can be identified.
Next, we describe the inversion method used in our computations. Finally, we discuss how to characterize the identifiability of the parameters involved in the inversion as based on singular value decomposition (SVD) of model currents.

The Base Model
As noted above, we aim to define an AP model that can be scaled from very early stages of human development (days) to fully developed adult cardiomyocytes. To review, for one specific membrane current, we assume that the only difference between the hiPSC-CM and adult cases is that the number of channels and the membrane area has changed; thus, the density of the specific ion channel carrying the current has changed, but the properties of every individual channel remains the same. The same principle holds for the intracellular Ca 2+ machinery; the individual channels and buffers remain the same, but both the intracellular volumes and the number of channels change from hiPSC-CMs to adult cardiomyocytes. Our model will, therefore, be based on models of single ion channel dynamics and only the density of these single channels will change. When a drug is involved, we assume that the effect of the drug on a single channel is the same in the hiPSC-CM and adult cases, and therefore one can use the effect in the hiPSC-CM case to estimate the effect for the adult case.

Modeling the Membrane Currents
The standard model [see, e.g., (Izhikevich, 2007;Plonsey and Barr, 2007;Ermentrout and Terman, 2010;Sterratt et al., 2011)] of the membrane potential of an excitable cell is given by the equation dv dt where v is the membrane potential (in mV), and I x are the membrane currents through ion channels of different types, as well as pumps and exchangers located on the cell membrane. These currents are all given in units of A/F, and may be written on the form where N x is the number of channels of type x on the cell membrane, A is the area of the cell membrane (in mm 2 ) and C m is the specific capacitance of the cell membrane (in pF/mm 2 ). Furthermore, i x represents the average current through a single channel of type x (in pA). For voltage-gated ion channels, this average single-channel current is given on the form where g x 0 is the conductance of a single open channel (in nS), E x is the equilibrium potential of the channel (in mV), and o x is the unitless open probability of the channel. Note that in models given on this form, it is common to consider a lumped parameter g x , given by and parameters of this type are given for each of the ion channels considered in the base model in the Supplementary Information. For membrane pumps and exchangers, the single-channel current is given on a similar form. The specific currents included in the model will be described below.

Scaling of the Membrane Currents
As mentioned above, we assume that the specific membrane capacitance and the ion channels responsible for each of the membrane currents are the same during different stages of development for the cell, but that the number of ion channels, N x , and the membrane area, A, may differ. Therefore, currents can be mapped from one stage of development, S 1 , to another stage of development, S 2 , simply by adjusting the channel density of the currents. More specifically, for the formulation (1)-(2), this means that we assume that the parameter C m and the expressions for the single-channel currents, i x , are the same for S 1 and S 2 , but that the channel density N x A can be different. Let A S 1 x , A S 2 x and N S 1 x , N S 2 x denote the membrane area and number of channels of type x for the S 1 and S 2 cases, respectively. Furthermore, let l x represent the change of channel density in the sense that Now, the S 1 and S 2 currents are related according to for each of the currents x.

The Base Model is the Generic Adult Model
Based on these considerations, it is convenient to define one default base model from which all other models are derived to simplify a mapping procedure between different development stages. Defining a base model as representing hiPSC-CMs, from which adult cardiomyocytes subsequently develop, may seem to be a natural choice. However, in the scheme illustrated in Figure 1, there is only one fixed model-the generic adult model-while all other models will change depending on the experimental measurements. For simplicity, we, therefore, define the generic adult model to be the default base model, and scale all other models relative to this model.

Main Currents Present in Human Cardiomyocytes
Modern models of human cardiomyocytes are complex and the models for the individual currents are based on years of experience using patch-clamp measurements. In the formulation (1), our aim has been to include the main currents present in human cardiomyocytes, but to keep the number of currents as low as feasible in order to keep the base model relatively simple. The experimental inputs in the present methodology are optically-derived, and data based on sensitive dyes are not expected to be able to uncover equally fine details of the dynamics as compared to traditional electrophysiological measurements derived via patch clamp. It is, therefore, reasonable to represent the data using simpler models. Our choice of currents is based on the O'Hara et al. (2011) model and the Grandi et al. (2010) model for human adult ventricular cardiomyocytes, in addition to the Paci et al. (2013) model for hiPSC-CMs. Furthermore, we have focused on including currents considered to be critical for depolarization and repolarization of the AP and, therefore, those typically investigated for response to drugs [see, e.g., (Crumb et al., 2016)].
In (Crumb et al., 2016), the fast sodium current, I Na , the late sodium current, I NaL , the L-type Ca 2+ current, I CaL , the transient outward potassium current, I to , the rapid and slow delayed rectifier potassium currents, I Kr and I Ks , and the inward rectifier potassium current, I K1 , were investigated for their drug responses, and we have included each of these currents in our model. In addition, we have included the sodium-potassium pump, I NaK , the sodium-calcium exchanger, I NaCa , the Ca 2+ pump, I pCa , the background Ca 2+ current, I bCa , and the background chloride current, I bCl , as they all appear to have a significant effect on the computed AP and Ca 2+ transient of the Grandi et al. (2010) model. Furthermore, we have included the hyperpolarization-activated cyclic nucleotide-gated funny current, I f . While this current is very small for adult ventricular cardiomyocytes, it is substantial for hiPSC-CMs (Garg et al., 2018). The formulation used for each of the currents is given in the Supplementary Information. The formulations are based on those of the currents in the Paci et al. (2013), the Grandi et al. (2010), and the O'Hara et al. (2011) models, and we have chosen formulations that seems to work well for both the hiPSC-CM and adult cases and that are able to provide good fits for our considered data of hiPSC-CMs.

Modeling Intracellular Ca 2+ Dynamics
In addition to the membrane potential, we also want the base model to represent the dynamics of the intracellular Ca 2+ concentration. We consider the following five intracellular compartments [based on (Grandi et al., 2010)]: 1. The dyad, representing the small cytosolic subspace between the L-type Ca 2+ channels and the ryanodine receptors (RyRs), 2. The subsarcolemmal space, representing the remaining part of the cytosolic space that is located close to the membrane, 3. The bulk cytosolic space, 4. The junctional sarcoplasmic reticulum (jSR), representing the part of the SR that is close to the RyR-channels, 5. The network sarcoplasmic reticulum (nSR), representing the remaining part of the SR.
The Ca 2+ concentrations and volume fractions defined for each of these compartments are given in Figure 2. In all compartments except the nSR, we consider both the concentration of free Ca 2+ and the concentration of Ca 2+ bound to a buffer. The Ca 2+ concentration in the extracellular space is assumed to remain constant. The intracellular Ca 2+ fluxes between compartments are illustrated in Figure 2, and the model takes the form See Section S2.1 in the Supplementary Information for a derivation of these equations.

Definition of Ca 2+ Fluxes
Every model flux J x representing the flux through a type of channel can be written on the form where N x is the number of channels of type x, V cell is the cell volume (in L) and j x is the average flux through a single channel of type x (in mmol/ms). Below, we will describe the formulation chosen for the flux through the RyR channels.
Definitions of each of the remaining fluxes are specified in the Supplementary Information.

Modeling Release From the SR
In our model of Ca 2+ dynamics, we deviate from previous modeling approaches in two specific ways: 1. Ca 2+ is released through RyR channels from the SR directly to the subsarcolemmal space (SL) and not to the dyad. 2. Release of Ca 2+ through the RyR channels is a product of two factors; one factor models the open probability of the RyR channels, whereas the other models the availability of channels that can be opened. We assume that each channel can only process a certain amount of Ca 2+ before it deactivates.
We will see below that these two modeling assumptions lead to a model that exhibits two key physiological features of Ca 2+ release from the SR of cardiomyocytes, so-called high gain and graded release (see Section S2.2 in the Supplementary Information for explanations of these terms).

Flux through RyRs (J sl s )
As we will employ the base model for several different parameter combinations, the model for the RyR flux must be stable, in the sense that careful tuning of the model is not requisite to ensure reasonable activation and deactivation of the RyRs. As outlined above, we let the Ca 2+ released from the SR enter the SL space rather than the dyad. This is done in order to achieve graded release (see the Supplementary Information), in the sense that the amount of Ca 2+ released from the SR through the RyRs should depend directly upon the amount of Ca 2+ entering the cell through L-type Ca 2+ channels. If the model were to be formulated such that Ca 2+ released from the jSR instead entered the dyad, it would be difficult to distinguish the increase in dyadic Ca 2+ concentration resulting from L-type Ca 2+ channel flux as opposed to release via RyRs. Directing the RyR flux into the SL, the concentration change in the dyad is almost exclusively due to the influx through L-type Ca 2+ channels, and by letting the flux through the RyRs depend on the Ca 2+ concentration in the dyad, we achieve graded release.
Furthermore, a common modeling approach for the RyRs is to govern inactivation by a decreased jSR concentration [see e.g., (Sobie et al., 2010)]. However, for large variations in parameter values, this may lead to model instabilities, because the jSR concentration depends upon the balance between the flux through the SERCA pumps and the RyRs, which depend upon the balance between the Ca 2+ fluxes in and out of the cell. In order to avoid an RyR model whose inactivation mechanism depends on the jSR concentration, we instead introduce a new assumption that some RyRs are only able to carry a given amount of Ca 2+ ions during each AP.
We then assume that a small portion of the RyR channels are always open (type 0), while the remaining channels (type 1) are activated by an increased dyadic Ca 2+ concentration and are inactivated after they have transported a given amount of Ca 2+ ions. Therefore, the total flux through the RyRs may be expressed as where J RyR represents the flux through the RyR channels of type 1 and J leak represents the flux through the RyR channels of type 0. We assume that the flux through the two types of RyR channels are given by expressions of the form where j RyR denotes the flux through a single open RyR channel (in mmol/ms) and V cell denotes the total cell volume (in L). In addition, M 0 denotes the number of RyR channels that are always open (type 0), M denotes the number of available RyR channels of type 1, and p is the open probability of the channels of type 1. The single channel flux through the RyRs is given by where a RyR,0 (in L/ms) represents the rate of release.
Furthermore, the open probability of the RyR channels of type 1 is modeled by a simple function that increases sigmoidally with the dyadic Ca 2+ concentration, c d , based on the model in (Friel, 1995): We let the total number of RyR channels of type 1 be given by N RyR and the total number of RyR channels of type 0 be given by In other words, the total number of RyR channels (of both types) is given by (1+g RyR )N RyR.
We assume that every RyR of type 1 is able to transport a fixed amount of Z Ca 2+ ions during an AP. After Z ions have been transported, the channel becomes inactivated. However, we assume that as the dyadic Ca 2+ concentration, c d , returns to rest and the open probability, p, consequently decreases, the inactivated channels gradually recover from inactivation. We let the number of available channels of type 1 be governed by Here, the first term dominates for large values of p, driving M towards zero as more Ca 2+ is transported through the RyR channels of type 1. Furthermore, for small values of p (i.e., at rest), the second term dominates and drives M towards the maximum value N RyR .
In order to reduce the number of free parameters in the model, we define a scaled variable r, defined as r = M N RyR , and divide both sides of equation (12) by N RyR . The equation then reads Inserting M = rN RyR into (7) and defining we get the following expression for active RyR flux where we recall that Moreover, inserting (11) and (15) into (8), we obtain

Scaling the RyR flux
When considering cells of different levels of maturity, we assume that the number of RyRs and the cell volume may be different, but that the function of a single RyR channel is the same for different levels of maturity. We also assume that the ratio between RyR channels of type 0 and 1, g RyR , and the number of Ca 2+ ions that each RyR channel of type 1 can transport, Z, is the same for the different maturity levels. Considering the model (13)-(18), this means that the only adjustment necessary between two maturity levels S 1 and S 2 is an adjustment of the density N RyR V cell in the definition of a RyR and b RyR . We, therefore, introduce a scaling factor l RyR such that and represent this adjustment of the RyR density in the model by scaling a RyR and b RyR by where superscript S 1 and S 2 denote the S 1 and S 2 versions of the parameters, respectively.

Inversion of Optical Measurements
The inversion procedure, used to construct base model representations of data obtained from optical measurements of the AP and Ca 2+ transient of hiPSC-CMs, is described below. First, in Section Optical Measurements, we describe how optical measurements of hiPSC-CMs are obtained. Next, in Section Definition of Adjustment Factors, we describe how adjustment factors l are set up to represent control (non-drugged) cells from different data sets. In Section IC50 Modeling of Drug Effects, we describe how the effect of a drug is modeled using IC50 values and corresponding factors, denoted by e. The aim of the inversion procedure is to find optimal parameter vectors l and e so that the model parameterized by l and e aligns to the measured data as best possible. This is explained in more detail in Section Coupled Inversion of Data From Several Doses. In Section Properties of the Cost Function, we describe the cost function constructed to measure the difference between the model and the data. Finally, in Section A Continuation-Based Minimization Method, we describe the continuation-based minimization method used to minimize the cost function in our computations.

Optical Measurements
Using previously developed techniques (Mathur et al., 2015), cardiac microphysiological systems derived from a single line of hiPSCs were loaded and matured prior to drug exposure. The resulting tissues consisted of approximately 90% cardiomyoctyes, with a small population of stromal support cells. On the day upon which studies were performed, freshly measured drugs (Nifedipine, Lidocaine, Cisapride, Flecainide, and Verapamil) were dissolved into DMSO or media and serially diluted. Each concentration of the drug to be tested was preheated for 15-20 min in a water bath at 37°C and subsequently sequentially injected in the device. At each dose, after 20 min of exposure, the drug's response on the microtissue was recorded using a Nikon Eclipse TE300 microscope fitted with a QImaging camera. Fluorescent images were acquired at 200 frames per second using filters to capture GCaMP and BeRST-1 fluorescence as previously described. Images were obtained across the entire chip for 6-8 sec at a resolution of 511 x 222 square 1.3 micron pixels. Excitation was paced at 1 Hz, to capture multiple beats for processing. Fluorescence videos were analyzed using custom Python software utilizing the open source Bio-Formats tool to produce characteristic AP and Ca 2+ waveforms for each chip and tested drug dose. Briefly, for each analysis, the fluorescent signal was averaged over the entire microtissue. The signal was then smoothed using a 3-point median filter, and five to seven individual action potentials or calcium transients were overlayed by aligning the maximum dF/dt and then averaged into a single transient. For each drug escalation study, we chose the single series with the most continuity between control cases and subsequent drug doses for both AP and Ca 2+ transient for inversion and mapping analysis.

Definition of Adjustment Factors
In order to make base model representations of control cells from different data sets, we must define adjustment factors l for a base model parameter set. These adjustment factors represent alterations of the channel densities and geometry of the cells under consideration, as explained above. For example, for each membrane channel type x, the adjustment factor l x is defined as where N x A is the channel density on the cell membrane for the fitted model and N b x A b is the channel density in the default base model. We generally consider adjustment factors for the membrane channel densities for all the currents of the model, i.e., l Na , l NaL , l CaL , l to , l Kr , l Ks , l K1 , l NaCa , l NaK , l pCa , l bCl , l bCa , and l f , although some of the factors are fixed in some cases (see Section Identifiability of the Currents in the hiPSC-CM Base Model).
For the density of an intracellular channel type x, the adjustment factor l x is similarly defined as where N x V cell and N b x V b cell are the channel densities for the fitted model and the default base model, respectively. We consider the adjustment factors l RyR and l SERCA for intracellular channel densities, and the factors l d B , l sl B , l c B and l s B for intracellular calcium buffers (see the Supplementary Information). In addition, we consider adjustments to the intracellular diffusion coefficients, l c d , l c sl , and l c n (see the Supplementary Information). In order to reduce the number offree parameters to be determined in the inversion procedure for different control data, we assume that the buffer concentrations change at the same rate in all intracellular compartments, so that we only consider a single adjustment factor Similarly, we assume that the intracellular diffusion coefficients change at the same rate, so that l c d = l c sl = l s n : = l a : Furthermore, because we wish to avoid ending up with unrealistic values of the surface-to-volume ratio, c, we assume that the scaling factor for the cell surface-to-volume ratio varies little between data sets and only employ two different values of c in the computations. We use the value c = 0.6 mm ˗1 for adult cells and the value c = 0.9 mm ˗1 for hiPSC-CMs, based on the values used in the Grandi et al. AP model for adult cardiomyocytes (Grandi et al., 2010) and the Paci et al. AP model for hiPSC-CMs (Paci et al., 2013). Note here that t-tubules (i.e., invaginations of the cell membrane extending into the center of the cell) are present for adult ventricular cardiomyocytes (Orchard et al., 2009), and this is incorporated into the adult version of c by increasing the cell surface area by a factor of about two compared to the geometrical surface of the cylinder shape of the cell [see e.g., (Luo and Rudy, 1994)]. For hiPSC-CMs, t-tubules are believed to be absent or underdeveloped [see e.g., (Di Baldassarre et al., 2018;Jiang et al., 2018)], and in our choice of c for hiPSC-CMs, we have assumed that t-tubules are not present for hiPSC-CMs.

IC50 Modeling of Drug Effects
Following previous modeling of channel blockers [see, e.g., (Brennan et al., 2009;Davies et al., 2012;Zemzemi et al., 2013;Paci et al., 2015)], we model the dose-dependent effect of a drug by scaling the channel conductances according to where g D i is the conductance of channel i in the presence of a drug with concentration D, IC50 i is the drug concentration that leads to 50% block of channel i, and g C i is the channel conductance in the control case (i.e., in the absence of drugs). Specifically, this means that if the drug concentration D equals the IC50 value, we have g D i = 1 2 g C i . It should be mentioned that a drug may certainly affect a channel in a more complex manner than is assumed here. The effect of drugs can realistically be represented by introducing new states in Markov models of the ion channel. In such models, the transition rates between different model states are able to represent the properties of drugs [see e.g., (Clancy et al., 2007;Tveito et al., 2011;Tveito and Lines, 2016;]. Although Markov model representations of drug effects are more versatile and realistic than the simple blocking assumption employed here (Tveito et al., 2017), it would greatly increase the complexity of the inversion process, as many more parameters would have to be computed.
From (26), we see that for a given drug dose D > 0, the effect of the drug would increase if the IC50 value were decreased, and the effect of the drug would be very small if the IC50 value were much larger than the considered dose. In the continuation-based minimization method applied in our computations (see the section A Continuation-Based Minimization Method below), it is most practical to deal with parameters that are small when no change occurs and large when large changes occur. Therefore, we introduce the parameters Here, a small value of e i represents small effects of a drug while a large value of e i represents large effects, and channel blocking is given by In our computations, we assume that the considered drugs block either I CaL , I NaL , I Kr , or a combination of these currents, and we therefore only consider the e-parameters e CaL , e NaL , e Kr .

Coupled Inversion of Data From Several Doses
The control data obtained from different optical experiments tend to vary significantly, and in order to be able to accurately estimate the drug effect from these measurements, the lparameters must be tuned so that the control model fits the control data as best possible. In addition, we want the l parameters to be constructed such that that the scaling (28) for e CaL , e NaL , and e Kr is sufficient to fit the model to the drug doses under consideration. In order to increase the chance of obtaining such a control model, we fit the control parameters, l, and the drug parameters, e, simultaneously, instead of first finding the optimal control parameters, l, by fitting the base model to the control data, and then subsequently finding appropriate drug parameters, e, for each dose. In addition, all doses are included in the inversion, so that the estimated values of e are based on all the drug doses included in the data set.
In order to illustrate the role of the land e-parameters more clearly, consider a simplified model consisting of just two currents, and assume that the base model is given by (see the section Modeling the Membrane Currents) Assume further that we have data from cells with both no drug present and with different doses of a drug (e.g., one low dose and one high dose). We assume that the drug may block any of the two model currents. In the inversion procedure, we try to find optimal values of the four parameters l 1 , l 2 , e 1 , and e 2 so that the adjusted model of the form fits the data both for the control case (D = 0) and for the considered drug doses. In other words, for a given parameter set l 1 , l 2 , e 1 , and e 2 , we need to compute the solution of the model (30) both for the control case (D = 0) and for the considered drug doses and compare the obtained solutions to the corresponding experimental data. The more general case considered in our computations is conceptually identical; however, as we also consider scaling of parameters that are not assumed to possibly be affected by the drug, we also have some parameters simply scaled by a factor (1+l i ) instead of by 1+l i 1+De i .

Properties of the Cost Function
In order to find the optimal parameters for fitting the model to data, we need to define a cost function that measures the difference between a given model solution and the data. This cost function is defined as Here, d runs over each of the considered drug doses, D d , including the control case (D 0 = 0), and j runs over each cost function term, H j , representing various differences between the data and the model solution. The parameters w d,j represent weights for each of the cost function terms for each of the doses. Each of the cost function terms, H j , are defined in Section S3.1 of the Supplementary Information.

A Continuation-Based Minimization Method
As outlined above, we wish to adjust the base model to data by finding land e-parameters that minimize a cost function of the form (31), measuring the difference between the input data and the model solution. In order to search for the optimal values of l and e, we apply a continuation-based optimization method [see e.g., (Keller, 1987;Allgower and Georg, 2012)]. Briefly, continuation is used to simplify the solution of equations or of optimization problems by introducing a q-parameterization such that the solution is known for one value of q. Suppose, for instance, that the parameterization is defined such that the solution is known for q = 0 and the problem we want to solve is defined by q = 1. Then the solution at q = 1 can be found by starting at q = 0 and carefully step towards the solution at q = 1. One advantage with this method is that we can start at a solution that we know is correct (at q = 0) and then take small steps towards the goal at q = 1. For the problem of inverting membrane potential and Ca 2+ traces, this method has proven to be useful.

Cost function in the continuation case
More specifically, we assume that, for each drug dose, D d , (including the control case) the data we are trying to invert are given by some vector pair [v 1 (D d ), c 1 (D d )], where v 1 (D d ) is the membrane potential and c 1 (D d ) is the Ca 2+ concentration. In addition, from the default base model specified by l = e = 0, we can compute a vector pair (v 0 , c 0 ) for the membrane potential and Ca 2+ concentration as the starting point of the inversion.
The goal of the continuation method is to compute a path for l and e from l = e = 0, which fit (v 0 , c 0 ) perfectly, to some l and e that fit the final data [v 1 (D d ), c 1 (D d )] for each of the drug doses, D d , as best as possible. This is done by defining a cost function of the form for the intermediate steps in the algorithm. Here, q is a parameter that is gradually increased from zero to one. In the definition (32), the terms H j (q, l, e, D d ) correspond to each of the terms H j (l, e, D d ) defined in Section S3.1 of the Supplementary Information. Specifically, the terms take the form where R j (v,c) represent different characteristics of the AP or Ca 2+ transient, e.g., the AP duration at some percentage or the upstroke velocity (see Section S3.1 of the Supplementary  Information 1 ). In the case q = 0, R q j (D d ) is equal to the terms defined by the default model (l = e = 0) for all the doses D d . Therefore, H(0, 0, 0, D d ) = 0, 2 so the optimal solution for q = 0 is l = e = 0. In the case q = 1, the terms R q j (D d ) are equal to the characteristics computed for the data we wish to invert. In other words, H(1, l, e) = H(l, e), where H(l, e) is defined in (31). For the intermediate values of q, the characteristics R q j (D d ) represent weighted averages of characteristics of the model used as a staring point for the inversion (l = e = 0) and the data we are trying to invert. Therefore, we expect the optimal values of l and e to gradually move from zero to the optimal values for the data as q is increased from zero to one.

The minimization algorithm
In the minimization algorithm, we find the optimal solution in M iterations. We define q m = Dq·m for m = 0,…,M, where Dq = 1/M. For m = 1,…,M, we assume that the optimal values l(q m˗1 ) and e(q m˗1 ) have been computed, and we want to find l(q m ) and e(q m ) by finding the minimum of H(q m , l, e). Since the step in q is small, we assume that the changes in l and e are also relatively small. We use the Nelder-Mead algorithm (Nelder and Mead, 1965) to minimize H(q m , l, e), and we use l(q m˗1 ) and e(q m-1 ) as suggestions for the starting vectors to find l(q m ) and e(q m ). However, in order to increase the chance of finding the true optimal value in every iteration, we start the Nelder-Mead algorithm from several randomly chosen starting vectors in the vicinity of l(q m˗1 ) and e(q m˗1 ). Figure S3 in the Supplementary Information illustrates the development of the e-values in an inversion aiming to characterize a drug.

Technical specifications
In the applications presented below, we use M = 20, and in each iteration m, we draw 63 guesses (as the specific computer used for these simulations has 64 cores) for the starting vectors for the Nelder-Mead algorithm from [l(q m-1 )−0.2, l(q m-1 )+0.2] and ½ e(q m−1 ) 5 , 5 Á e(q m−1 ) for l and e, respectively. In the first 15 iterations, we use five iterations of the Nelder-Mead algorithm for each guess, and for the last five iterations we use 25 iterations of the Nelder-Mead algorithm. For each new parameter set, we generally run the simulation for 15 AP cycles using 1 Hz pacing before measuring the AP and Ca 2+ transient, unless otherwise specified. The choice of 15 AP cycles is selected as a compromise between the desire of minimizing the computational efforts required for each cost function evaluation and the desire of reaching new stable steady state values for the state variables following a parameter change. Presently, each cost function evaluation requires about 14 sec of computing time for a data set including four doses in addition to the control case.

Identifiability of the Base Model Based on Singular Value Decomposition of Currents
In the inversion procedure outlined above, we try to find the optimal adjustment factors l and e for the model so that the AP and cytosolic Ca 2+ transient in the model solution match measurements of the AP and Ca 2+ transient as best possible. An important element to consider in this process is whether the identified adjustment factors found by the inversion procedure are the only combination of adjustment factors that fit the data, or whether other adjustment factors might exist which fit the data equally well.
In order to investigate the identifiability of the adjustment factors for the currents in the base model, we apply a method based on singular value decomposition [see, e.g., (Liesen and Mehrmann, 2015;Lyche, 2017)] of the currents. This approach is described in detail in (Jaeger et al., 2019a). In short, the identifiability of the currents is investigated by collecting the model currents at time points t n = nDt, for n = 1, …, N t into a matrix A ∈ R N t ,N c , where N c is the number of model currents. Then, the singular value decomposition of the matrix A = USV T is computed. Here, the matrices U ∈ R N t ,N t and V ∈ R N c ,N c are unitary matrices, and the matrix S ∈ R N t ,N c is a diagonal matrix with singular values s i along the diagonal. The columns u i and v i of U and V, respectively, are the associated singular vectors.
From the properties of the singular value decomposition it can be shown that perturbations of the adjustment factors along singular vectors v i associated with large singular values s i are expected to result in significant changes in the AP, whereas perturbations of the adjustment factors along singular vectors v i associated with small singular values are, accordingly, expected to result in small changes in the AP.
In (Jaeger et al., 2019a) it was shown that this expected result seemed to hold in the case of three well-known AP models of adult ventricular cardiomyocytes. In addition, it was demonstrated how this analysis could be used to define an identifiability index for individual model currents. This index was defined for each current j=1,...,N c as k(e j ) = jje j − P N e j jj 2 , where e j ∈ R N c is the vector that is one in element number j and zero elsewhere. Moreover, P N e j ∈ R N c is the projection of e j onto the unidentifiable space spanned by the singular vectors v i associated with small singular values (or small perturbation effects). In other words, if k(e j ) is close to zero, almost the entire current I j is in the unidentifiable space, and we cannot be sure that the value of the associated adjustment factors l j or e j are the only values that fit the data (i.e., result in the same AP).
On the other hand, if k(e j ) is close to one, we expect that other values of l j or e j would not fit the data as well as the currently assumed values, as perturbations of the adjustment factors would result in large changes in the AP. Note that this approach only aims to determine the identifiability of the adjustment factors for the membrane currents. The analysis could be extended to include other state variables than just the membrane potential (e.g., the Ca 2+ concentrations). In this case, the identifiability of the remaining adjustment factors might also be suggested. However, at this stage primary focus is on identifying drug effects on membrane ion channels, so we are principally interested in ensuring that the adjustment factors for the currents are unique.

RESULTS
Below, we demonstrate a few applications of the method outlined above. First, in the section The Base Model, we define the default hiPSC-CM and adult parameterizations of the general base model formulation. We also demonstrate that these models exhibit high gain and graded release of the Ca 2+ fluxes. In addition, we illustrate the identifiability of the model currents using SVD analysis, as described above. This analysis is used to determine which model currents should be fixed in the applications of the inversion procedure. Next, in the section Identification of Drug Effects on hiPSC-CMs Based on Simulated Data, we use the inversion procedure to identify drug effects for data generated by simulations. Finally, in the section Identification of Drug Effects on hiPSC-CMs Based on Optical Measurements, we apply the inversion procedure to identify drug effects from data obtained from optical measurements of hiPSC-CMs.

The Base Model
Here, we set up the default adult and hiPSC-CM base model formulations used in the inversion procedure in the following sections.

Base Model Approximation of the Grandi Model
The adult base model is fitted to approximate the Grandi et al. model using the inversion procedure described above. The upper right panel of Figure 3 shows the AP and Ca 2+ transient of the Grandi et al. (2010)  In the inversion, the cost function includes all the terms specified in Section S3.1 of the Supplementary Information, except for the regularization terms. For the cost function terms involving information about currents and fluxes, we have included I Na , I CaL , I to , I Kr , I Ks , I K1 , I NaCa , I pCa , and I bCa , as well as the fluxes J RyR and J SERCA (see Section S3.1.7 of the Supplementary Information). All terms measuring the difference in membrane potential or Ca 2+ concentration are given the weight w j = 1 and the terms measuring differences in the currents are given the weight w j = 0.5. The initial conditions for the parameters included in the inversion are specified in Table S9 of the Supplementary Information.
As mentioned above, we define the default base model as the adult base model because this model will remain fixed, whereas the hiPSC-CM models will change depending on experimental data. The parameter values obtained in the inversion procedure therefore define the default base model and are specified in Section S1 of the Supplementary Information.

Default Base Model for hiPSC-CMs
The left panel of Figure 3 shows the solution of the default base model for hiPSC-CMs fitted to optical measurements of the AP and Ca 2+ transient.  Supplementary Information), where the information about the currents is obtained from the Paci et al. model (Paci et al., 2013) which is based on patch-clamp recordings of the ionic currents of hiPSC-CMs from (Ma et al., 2011 The mapping between hiPSC-CMs and adult cardiomyocytes returned by the inversion procedure are reported in Table 1. Note that these factors represent the default hiPSC-CM base model to be used as a starting point for the inversion of the remaining control data sets. In other words, the specific adjustment factors between the hiPSC-CM and adult models will differ for each new data set. Note also that in our applications of the inversion procedure, we consider data from stimulated hiPSC-CMs, and therefore, the hiPSC-CMs are not required to be spontaneously beating. The default hiPSC-CM version of the base model in Figure 3 is not spontaneously beating, and whether the model fitted to a specific data set is spontaneously beating or not will depend on the fitted parameter values. In Figure S4

High Gain and Graded Release of the Base Model
As mentioned above, the base model formulation of Ca 2+ release is designed to exhibit both high gain and graded release. This has proved impossible to achieve using common pool models [see, e.g., (Rice et al., 1999;Dupont et al., 2016)], as discussed in more detail in the Supplementary Information. The Ca 2+ release   model we have designed differs from the classical common pool models in two ways: first, release of Ca 2+ from the SR is not directed into the dyad (d), but rather directly to the subsarcolemmal (SL) space (see Figure 2), and, second; the release mechanism is formulated in terms of both an availability rate and open probability for the RyRs [see (16)].
In Figure 4, we show that this model exhibits high gain and graded release both when the hiPSC-CM and adult parameters are applied. In the figure, we report the peak of the J CaL and J RyR fluxes as well as the integrated fluxes for simulations in which the membrane potential is fixed at specific values. The remaining state variables of the model start at the default initial conditions corresponding to the default resting membrane potential of the model, and the simulations record the J CaL and J RyR fluxes resulting from the clamped membrane potential.
We observe that for most values of v, the J RyR flux is considerably larger than the J CaL flux, indicating high gain. Furthermore, a small J CaL flux seems to be associated with a small J RyR flux, whereas a large J CaL flux is associated with a large J RyR flux, indicating graded release.

Identifiability of the Currents in the hiPSC-CM Base Model
In order to investigate the identifiability of the individual model currents, we apply the singular value decomposition analysis from (Jaeger et al., 2019a) described in the section Identifiability of the Base Model Based on Singular Value Decomposition of Currents.
In Figure 5, titles above each plot indicate the value of each of the singular values of the current matrix, A. The upper plots below the singular values show the singular vectors corresponding to each of the singular values. Here, each letter corresponds to a single current specified in the table on the righthand side. The below left plots show the values of the cost function (31) evaluated using the default base model for hiPSC-CMs as data and a perturbed model as the model solution. In the perturbed model, the maximum conductances are perturbed with l-values [see (5)] equal to w •v i , where v i is the considered singular vector and w is varied between zero and one. The cost function includes the terms H APD30 , H APD50 , H APD80 , and H Int30 with weight 1 for all terms except H APD80 , which is given the weight 5 (see the Supplementary Information for definitions). The maximum values of H are given in the top of the plots. The right plots show the solutions resulting from the perturbations for a few selections of w.
In (Jaeger et al., 2019a) it was shown that perturbations of the maximum conductances along singular vectors corresponding to large singular values generally resulted in large perturbation effects, whereas perturbations along singular vectors corresponding to small singular values generally resulted in small perturbation effects for the Ten Tusscher and Panfilov,  Figure 5 shows correspondence to this result for the hiPSC-CM base model; the main discrepancy is observed for s 2 , which corresponds to a singular vector consisting almost exclusively of the fast sodium current, I Na . The perturbation FIGURE 4 | Graded release for the hiPSC-CM (left) and adult (right) versions of the base model. In the upper panel, we report the peak of the J CaL and J RyR fluxes for simulations in which the membrane potential is fixed at specific values between −50 mV and 80 mV. In the lower panel, we show the fluxes integrated with respect to time from t = 0 ms to t = 100 ms. After 100 ms, both J CaL and J RyR have roughly returned to their resting levels. effects may be very small for this singular value because the upstroke velocity, physiologically governed by I Na , is not included in the cost function [cf. (Jaeger et al., 2019a)].
In order to quantify the identifiability of the individual currents, we compute the identifiability index, k, defined in (35). The unidentifiable space is defined as the space spanned by the singular vectors v i whose maximum value of H(w·v i ) for 0≤w≤1 is smaller than 0.05. The computed values of the identifiability index for each of the model currents are given in the orange box on the right-hand side of Figure 5. A value of k close to 1 indicates a high degree of identifiability, while a value of k close to 0 indicates an unidentifiable current.
From the indices in Figure 5, we see that I CaL , I Kr , and I NaCa are highly identifiable, but that the currents I NaL , I Na , I bCa , I Ks , and I bCl has an identifiability index below 0.5. As a consequence, we fix the conductance of I Na , I bCa , I Ks , and I bCl in the applications of the inversion procedure presented below. In addition, we are aware that the I NaL current might be hard to identify, and that estimated drug effects for this current are associated with a level of uncertainty [see also (Poulet et al., 2015)].

Identification of Drug Effects on hiPSC-CMs Based on Simulated Data
Our first application of the inversion procedure is to identify drug effects as based on simulated data. To generate these data, we set l CaL = l NaL = l Kr = 0.1 in the default hiPSC-CM base model. In addition, we apply a set of e-values to represent five specific drugs-Nifedipine, Lidocaine, Cisapride, Flecainide, and Verapamil. We assume that Nifedipine is a pure I CaL -blocker with an IC50 value of 10 nM, that Lidocaine is a pure I NaLblocker with an IC50 value of 10 mM, and that Cisapride is a pure I Kr -blocker with an IC50 value of 10 nM. Furthermore, Flecainide is assumed to block a combination of all three currents with IC50 values of 25 mM, 20 mM and 10 mM for I CaL , I NaL , and I Kr , respectively. Verapamil is assumed to block I CaL with an IC50 value of 200 nM and I Kr with an IC50 value of 500 nM. Both when the data is generated and in the inversion procedure, we record the sixth generated AP following each parameter change. Figure 6 shows the result of the inversion procedure using the l-values l CaL , l NaL , and l Kr and the e-values e CaL , e NaL , and e Kr as free parameters in the inversion procedure. The left panel shows the e-values used to generate the data (yellow) and the corresponding e-values returned by the inversion procedure (pink). The center and right panels show the AP and Ca 2+ transient, respectively, for the control case and for each of the drug doses included in the data sets. The solid lines show the simulated data and the dotted lines show the solutions generated by the model using the land e-values returned by the inversion procedure. Note that to clearly see differences in the Ca 2+ transient amplitude, the Ca 2+ transients are adjusted so that the Ca 2+ transient amplitude is preserved, but the minimum Ca 2+ concentration is set to zero. We observe that the inversion FIGURE 6 | Identification of drug effects for five drugs based on simulated data. The l-values l CaL , l NaL , and l Kr and the e-values e CaL , e NaL , and e Kr are allowed to vary in the inversion. The left panel shows the e-values used to generate the simulated drug data (yellow) and the corresponding e-values estimated by the inversion procedure (pink). The center and right panels show the AP and Ca 2+ transients, respectively, for each of the drug doses included in the data sets. Solid lines represent the simulated data and dotted lines show the fitted model solutions returned by the inversion procedure. Note that to clearly see changes in the Ca 2+ transient amplitude, the Ca 2+ transients are adjusted so that the Ca 2+ transient amplitude is preserved, but the minimum value is set to zero in all cases. procedure is able to identify the correct e-values accurately, excepting the e-value for Lidocaine, which is predicted to be considerably lower than the value used to generate the data.

Identification of Drug Effects on hiPSC-CMs Based on Optical Measurements
Below, we present use of the inversion procedure to identify drug effects on hiPSC-CMs from optical measurements of the AP and Ca 2+ transient. Figure 7 shows the result of the inversion procedure applied to data from optical measurements of hiPSC-CMs exposed to the drug Nifedipine. The data includes the control case with no drug present and four different drug doses (3 nM, 30 nM, 300 nM, and 3,000 nM). The left panel of Figure 7A shows the membrane potential and Ca 2+ traces obtained from optical measurements, and the center panel shows the corresponding solutions of the hiPSC-CM version of the base model fitted to the optical measurements. Note that the values of the data are mapped so that the maximum and minimum values of the membrane potential and Ca 2+ concentration match those of the fitted hiPSC-CM model. Panel C of Figure 7 compares the experimentally measured data and the fitted model for each of the doses. We observe that the model seems to fit the data quite well for most of the doses, but that the Ca 2+ transient appears to last a bit longer in the model than in the data for the highest considered drug doses.

Nifedipine
The dose-dependent effect of the drug on the I CaL , I NaL , and I Kr currents are modeled using IC50 values (see Section IC50 Modeling of Drug Effects). The values of e i = 1 IC50 i for i = CaL, NaL, and Kr are given in Figure 7B. A large value of e i corresponds to a large drug effect on the current i, and a small value of e i corresponds to a small drug effect on the current i. From Figure 7B, we observe that the inversion procedure predicts that Nifedipine primarily blocks I CaL .
The IC50 values corresponding to the estimated e-values for I CaL , I NaL and I Kr are given and compared to literature values in Table 2 for all the five considered drugs of this section (see the Discussion section for a discussion of these results). Figure 8 shows similar results for inversion of measurements of hiPSC-CMs exposed to the drug Lidocaine. In panel A, we observe that the AP duration is reduced by the drug, and in panel B, we observe that the inversion procedure predicts that the drug primarily blocks the I NaL current. Figure 9 shows the result of the inversion procedure applied to a data set for hiPSC-CMs exposed to the drug Cisapride. In panel A, we observe that the drug increases the AP duration. In panel B, we observe that the inversion procedure predicts that Cisapride primarily blocks the I Kr current. Figure 10 shows the result for the inversion procedure applied to optical measurements of hiPSC-CMs exposed to the drug Flecainide. In panel A, we observe that the drug causes increased AP duration. In panel C, we observe that the fitted model fits the data quite well, excepting that the AP duration at high percentages of repolarization is longer for the data than for the model for the highest considered dose. In addition, the shape of the Ca 2+ transient for the low doses is not entirely captured in the model. In panel B, we observe that the inversion procedure estimates that the drug primarily blocks I Kr and, to some degree, I CaL .

Flecainide
Verapamil Figure 11 shows the result of the inversion procedure applied to measurements of hiPSC-CMs exposed to the drug Verapamil. In panel A, we observe that the drug leads to decreased AP duration.
Panel B shows that the inversion procedure predicts that Verapamil primarily blocks I CaL and, to some extent, I Kr .

Mapping of Drug Effects From hiPSC-CMs to Adult Cells
The rightmost plots of panel A of Figures 7-11 show the predicted drug effects for adult cells for each of the drugs considered. More specifically, the plots show the solution of the adult base model exposed to each drug's effect (e-values) as estimated by the inversion procedure for each of the drug doses included in the data set. To review, this represents the predicted drug response for an adult AP and Ca 2+ transient exposed to each of the drugs, based on the optical measurements of the AP and Ca 2+ transient as obtained in a microphysiological system of hiPSC-CMs. The predictions are made by first using the inversion procedure to estimate the effect of the drug on the I CaL , I NaL , and I Kr currents in the hiPSC-CM case and then mapping the corresponding drug effects to an adult cell using the determined maturation map based on the assumptions of differences in the channel densities and geometry between hiPSC-CMs and adult cardiomyocytes (see the Methods section).

DISCUSSION
Here, we have presented an improved version of the methods presented in  for estimating drug effects for adult human cardiomyocytes based on optical measurements of the AP and Ca 2+ transient of hiPSC-CMs in a microphysiological system. First, we introduce a new base model formulation for representing both adult cells and hiPSC-CMs via different parameter sets. A model for intracellular Ca 2+ dynamics is updated to a formulation constructed for stability with respect to parameter changes. In addition, we use IC50-based modeling of dose-dependent drug effects and find optimal parameters by running a coupled inversion of both the control data and the drug data for several different doses. The cost function measuring  (Kramer et al., 2013) 11 800 nM (Kramer et al., 2013) 26 mM (Crumb et al., 2016) 202 nM (Crumb et al., 2016) 60 nM (Di Stilo et al., 1998 2 7 mM (Kramer et al., 2013) 200 nM (Kramer et al., 2013) (Zhabyeyev et al., 2013) 20 nM (Kramer et al., 2013) 1.5 mM (Kramer et al., 2013) 250 nM (Kramer et al., 2013) 6.5 nM (Mohammad et al., 1997) 143 nM (Zhang et al., 1999) the difference between the data and the model has also been redefined, and we now apply a continuation-based minimization method to minimize the cost function.

Summary of Method Performance and Main Results: Identification of Drug Effects Based on Simulated Data and Optical
Measurements of hiPSC-CMs Figure 6 shows the result of the inversion procedure applied to simulated data. As noted above, we observe that the inversion procedure is able to identify the correct e-values accurately, excepting the e-value for Lidocaine, which is predicted to be considerably lower than the value used to generate the data. This suggests that it might be difficult to obtain correct values of e NaL , as also supported by the low identifiability index for I NaL reported in Figure 5. In addition, we observe that the inversion procedure predicts some block of I NaL for the drug Cisapride, even though only I Kr was blocked when the data was generated.
We additionally presented the use of the inversion procedure to identify drug effects on hiPSC-CMs from optical measurements of the AP and Ca 2+ transients.
Panel C of Figure 7 compares the experimentally measured data and the fitted model for each of the doses of Nifedipine applied to the microphysiological system. The model fits both the membrane potential and the Ca 2+ data well for most doses applied, although the Ca 2+ transient duration is longer in the model than in the data for the highest drug doses considered. Furthermore, in panel B, we observe that the inversion procedure predicts that Nifedipine primarily blocks I CaL . In Table 2, we observe that the IC50 value for I CaL is estimated to be 38 nM, in agreement with values found in literature (12 nM-60 nM (Di Stilo et al., 1998;Kramer et al., 2013)]. The IC50 value for I NaL and I Kr are estimated to be 23 600 nM and 40 200 nM, respectively-considerably larger than the doses considered in the data set. We have not found an IC50 value for I NaL for comparison in literature, but the IC50 values found for I Kr support the claim that the IC50 value is much larger than the drug doses included in the data set, although the literature values [275,000-440,000 nM (Zhabyeyev et al., 2000;Kramer et al., 2013)] are higher than the value predicted by the inversion procedure. Figure 8 shows similar results for inversion of measurements of hiPSC-CMs exposed to the drug Lidocaine. The AP duration is reduced by the drug and the inversion procedure predicts that the drug primarily blocks the I NaL current. The IC50 value estimate for I NaL is 4.3 mM (see Table 2), in rough agreement with values found in literature [11 mM (Crumb et al., 2016)]. We observe that the model fits the data quite well, but that the AP duration for the drug dose of 10 mM is longer in the model than in the data. Figure 9 shows the result of the inversion procedure applied to a data set for hiPSC-CMs exposed to the drug Cisapride. Considering the leftmost and center panels of Figure 9A, we observe that the prolongation of the AP duration is much more prominent in the data as compared to the fitted hiPSC-CM model for a drug dose of 1 nM. This is also confirmed in Figure  9C, where we observe that the model does not fit the membrane potential data for the control case and the 1 nM dose case well.
The fit for the largest dose, however, is quite good. In Figure 9B, we observe that the inversion procedure predicts that Cisapride primarily blocks I Kr . In Table 2, we observe that the IC50 value for I Kr is estimated to be 13 nM, in good agreement with values found in literature [6.5 nM-20 nM (Mohammad et al., 1997;Kramer et al., 2013;Crumb et al., 2016)]. Figure 10 shows the result for the method as applied to measurements of hiPSC-CMs exposed to the drug Flecainide, know to prolong the AP duration. In Table 2, we observe that the IC50 value for I Kr predicted by the inversion procedure (1.9 mM) is in quite good agreement with literature values [0.7-1.5 mM (Kramer et al., 2013;Crumb et al., 2016)], but that the predicted IC50 value for I CaL (9 mM) is too low compared to the reported literature values [26-27 mM (Kramer et al., 2013;Crumb et al., 2016)]. In addition, the estimated IC50 value for I NaL (47 mM) is larger than the literature value of 19 mM (Crumb et al., 2016).
In Figure 11, the method is applied to measurements of hiPSC-CMs exposed to Verapamil. In panel A, the effect on the AP duration for the smallest dose (100 nM) appears to be more prominent in the data than in the fitted model. This is confirmed in panel C, where we observe that the fitted model seems to fit the Ca 2+ data considerably better than the membrane potential data. In particular, the AP duration is too short in the control case and too long for the smallest dose of 100 nM. Panel B shows that the inversion procedure predicts that Verapamil primarily blocks I CaL and to some extent I Kr . The predicted IC50 values from  (Mirams et al., 2011;Kramer et al., 2013;Crumb et al., 2016)] and 143-499 nM for I Kr (Zhang et al., 1999;Kramer et al., 2013;Crumb et al., 2016)].
The rightmost plots of panel A of Figures 7-11 show the predicted drug effects for adult cells for each of the drugs considered. We observe that for some drugs (e.g., Nifedipine and Lidocaine), the drug effects for adult cells are predicted to be approximately as severe as for the hiPSC-CMs. For other drugs, on the other hand, (e.g., Flecainide), the drug effect is predicted to be less severe for the adult cells than for hiPSC-CMs, highlighting the critical importance of maturation phenotype in predictive biophysical models of hiPSC-CMs in pharmacological studies.

Note on Ongoing Complementary Studies
Other recent work has made strong progress in terms of enabling biophysical modeling approaches to assimilate and otherwise make best use of tailored experimental measurements of hiPSC-CMs. For example, in (Gong and Sobie, 2018), the authors address the need to bridge the gap between the effect of drugs on human adult ventricular cardiomyocytes and the effect on animal or hiPSC experimental models often used in drug screening. This work also successfully generated accurate predictions of the effect of ion channel blocking drugs on human adult ventricular cardiomyocytes as based on simulations of hiPSC-CMs via a regression strategy.
Additional recent studies have advanced specific models and methodological approaches for hiPSC-CMs which incorporate experimental variability from multiple data sources [see, e.g., (Kernik et al., 2019)] with the goal of identifying phenotypic mechanisms and identify key parameter sensitivity. The authors introduce a computational whole-cell electrophysiological model of hiPSC-CMs, composed of single exponential voltagedependent gating variable rate constants, which are then parameterized to fit experimental measurements of hiPSC-CMs from multiple laboratories (and thus incorporate variability in the single-cell measurements of ionic currents observed experimentally). The authors compare hiPSC-CM and adult cell models to elucidate the primary properties underpinning the phenotype, a mechanistically central goal that was not an aim of our present study.

Limitations and Notes on Future Work
The model for the intracellular Ca 2+ dynamics in the base model exhibits both high gain and graded release for the hiPSC-CM and adult parameter sets (see Figure 4). However, the assumptions underlying the release model are introduced to obtain a stable model, and not necessarily to represent the underlying physiological mechanisms accurately. Future work necessitates assessment of the Ca 2+ machinery in the base model and potential redevelopment to more accurately represent physiological Ca 2+ release from the SR, as relevant.
In addition, the intracellular Na + and K + concentrations are assumed to be constants in the base model formulation. This was done in order to avoid problems with drift of the concentrations following from parameter changes [see, e.g., (Hund et al., 2001;Wilders, 2007;]. Moreover, freezing the intracellular Na + and K + concentrations during an action potential had very limited effects of the computed membrane potential and cytosolic Ca 2+ concentration. However, in some cases, drugs are believed to lead to significant changes in, e.g., the intracellular Na + concentration, which could affect the AP shape [see, e.g., (Brill and Andrew Wasserstrom, 1986;Faber and Rudy, 2000)]. Freezing the Na + and K + concentrations could therefore potentially make the base model less suitable for detecting such drug effects.
Note also that the terms included in the cost function (31) applied in the inversion procedure may in future work be adapted as required by the specific application under consideration. For example, the cost function could be extended to include information about the frequency dependence of important action potential and Ca 2+ transient features (see Figure S4 in the Supplementary Information) or to include information about the spontaneous activity of the hiPSC-CMs used in the experiments.
In the construction of the maturation map, we currently assume that the function of a single channel is the same for different levels of maturity and that only the geometry of the cell and the number of channels change with maturation. It is, however, perhaps possible for the function of some channels to change during maturation as well. If the conductance of a channel changes during maturation, the same adjustment factors as we have considered may still be applied, but if the dynamics governing the open probability of the channel change, additional adjustment factors would have to be included in the models for the channel open probability.
In addition, we have only looked at a single dose escalation study for each of the drugs investigated, as well as only tested the method on tissues derived from a single stem cell line. Future work will assess the variability of the inversion methodology in combination with noisy and incomplete experimental data obtained through these systems across a range of biological maturation approaches and starting stem cell lines.
Furthermore, we have only considered data of the drug response for hiPSC-CMs in microphysiological systems, and we have not been able to obtain corresponding data for adult human cardiomyocytes. Drug response data for isolated healthy adult human cardiomyocytes are generally very limited for both technical and ethical reasons Sala et al., 2017). We have therefore not been able to validate the predicted adult drug effects against experimental data.
Further validation of the methodology will primarily be based on further analysis of optical data. Hopefully, we will be able to perform analysis of many drugs and also include blind testing for drugs with well-known properties. When more data from patchclamp measurements become available, that will also be very useful for improvements and validation of our methods.
In future work, we will seek to find ways to estimate the sodium current which is not possible to estimate today because of low time resolution in the optical data. Furthermore, we will combine the base model with the bidomain model [see e.g., (Franzone et al., 2014)] to study extracellular waveforms in the chips, and also the more detailed EMI model [see e.g., (Tveito et al., 2017;Jaeger et al., 2019b)] where individual cells can be represented. Hopefully, the spatially resolved models can provide improved accuracy of the inversion process as well as test important considerations such as the effect of tissue composition.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/ Supplementary Material.

AUTHOR CONTRIBUTIONS
KJ and AT are responsible for the development of the mathematical framework and computer modeling. VC, BC, and KH are responsible for the generation of data provided from microphysiological systems. HF and SW are responsible for the analysis of data from microphysiological systems. KJ, AT, SW, MM, and HF wrote the manuscript text and created the figures. All authors reviewed the manuscript.