Inversion and computational maturation of drug response using human stem cell derived cardiomyocytes in microphysiological systems

While cardiomyocytes differentiated from human induced pluripotent stems cells (hiPSCs) hold great promise for drug screening, the electrophysiological properties of these cells can be variable and immature, producing results that are significantly different from their human adult counterparts. Here, we describe a computational framework to address this limitation, and show how in silico methods, applied to measurements on immature cardiomyocytes, can be used to both identify drug action and to predict its effect in mature cells. Our synthetic and experimental results indicate that optically obtained waveforms of voltage and calcium from microphysiological systems can be inverted into information on drug ion channel blockage, and then, through assuming functional invariance of proteins during maturation, this data can be used to predict drug induced changes in mature ventricular cells. Together, this pipeline of measurements and computational analysis could significantly improve the ability of hiPSC derived cardiomycocytes to predict dangerous drug side effects.


Modification of the action potential model
In the process of adjusting the Paci et al. [1] model to the data obtained from an MPS (microphysiological system, [2]), we have to run the model many thousand times with varying choices of parameters. One difficulty encountered in this process is drift of ion concentrations. This is a well-known problem of mathematical models of electrophysiology; see e.g. [3,4,5]. In Figure S1, we illustrate this problem for the original Paci et al. model. One approach to solve this problem is to decompose stimulus currents into ion concentrations and thereby retain conservation of the ion concentrations, see e.g. [3]. A problem with this approach is that drift is observed also when no stimulus is applied (see Figure S1). Another approach relies on the fact that some ion concentrations vary little and can therefore be kept constant. Here, we follow this latter approach and freeze the intracellular sodium concentration and the SR calcium concentration at their initial value. In Figure S2, we show the properties of this approximation. In the right panel, we note that the cytosolic calcium concentration no longer drifts even for very long simulations. In the left panel, we show that the effect of this approximation Figure S1: Example of drift of ionic concentrations in the Paci et al. model [1] with no stimulus current applied. First, we compute the steady state solution of the original Paci et al. model. Then, we reduce the I Kr current by 50% and run a simulation of this adjusted model for 200 seconds (corresponding to approximately 120 AP cycles). The plots show how the cytosolic calcium concentration (left panel), the SR calcium concentration (center panel), and the intracellular sodium concentration (right panel) change with time during this long simulation.   on the transmembrane potential and the cytosolic calcium concentration is very small. With this approximation, convergence towards the steady state solution (a steady periodic solution) is rapid and the solutions appears to be stable. This is demonstrated in Figure S3 where convergence to steady state is illustrated. First, we compute the steady state solution of the modified model using the original parameters of the Paci et al. model. Then, we reduce I Kr by 50% and note that the solution rapidly reaches equilibrium.

The maturation map of the Calcium dynamics
We consider how the Ca-dynamics change under maturation. As for the membrane ion channel case, we do this by illustrating the maturation process for a very simple model.
We consider an intracellular volume consisting of the cytosol (c) surrounding the sarcoplasmic reticulum (SR (s)); see Figure S4.
We let N c denote the number of Ca 2+ -ions in the cytosol and N s denote the number of Ca 2+ -ions in the sarcoplasmic reticulum; both given in mmol. The associated volumes are given by V c and V s , both given in L. Let J c,s and J s,c denote the flux (in mmol/ms) of Ca 2+ -ions from the cytosol to the SR, and from the SR to the cytosol, -respectively. Conservation of Ca 2+ -ions yields the model The fluxes are models of proteins carrying ions from one volume to the other. Let g 0 c,s (in mmol/ms) be the flux representing one single protein transporting Ca 2+ions from the cytosol to the SR. Similarly, g 0 s,c (in mmol/ms) is the flux representing one single protein releasing Ca 2+ -ions from the SR to the cytosol. The number of such proteins are given by N c,s and N s,c . Then, the system (1) and (2) takes the form By defining the fluxes (in mM/ms) the system takes the form where C c and C s are the concentrations (in mM) of Ca 2+ -ions in the cytosol and SR, respectively; For maturation, we can now follow the same steps as for the membrane proteins. During maturation, the properties of the single proteins will remain constant, but the number of proteins and the volumes will increase. Therefore, we introduce constants q Vc , q Vs , q Nc,s and q Ns,c such that With we get Consequently, we have the IM model and the associated M model Again, we observe that the M model is obtained simply by multiplication by a set of maturation factors.

Technical specifications of the model formulation and inversion procedure
In this section, technical specifications regarding the model formulation used in the simulations and the inversion procedure will be provided.

Intracellular concentrations
In almost all of our computations, we use the modified version of the Paci et al. model described above with fixed intracellular sodium and SR calcium concentrations. The only exception is that we also run some simulations of ten Tusscher et al. model [6] in Figure 9 of the paper. In these simulations, the intracellular potassium, sodium and SR calcium concentrations are also fixed at constant values. The intracellular concentrations used in the IM and M formulations of the modified Paci et al. model and the similarly modified ten Tusscher et al. model are given in Table S1.

Numerical stimulation protocol
In all simulations, the cells are stimulated every 1000 ms by a 5 ms long stimulus current of 8 µA/µF. The simulations are run for five AP cycles before recording the action potential and calcium transient for each new parameter combination. [

Technical specifications for the drug inversions
When the inversion procedure is used to fit simulated or measured drug and control data, we only consider adjustments of the q Na , q CaL , q Kr and q K1 factors, unless otherwise specified. Note, however, that for the inversion of the Verapamil data in Figures 6 and 7 of the paper, the I Na current was reduced by 50%, the I NaK current was reduced by 60%, the I CaL current was increased by 60%, and the I up and I rel fluxes were increased by 30% before running the inversion of the q Na , q CaL , q Kr and q K1 factors, in order to make the base model used in the inversion more similar to the control data.

Technical specifications for the construction of the maturation map
In the construction of the maturation maps demonstrated in Figure 9 of the paper, we use the inversion procedure to fit an immature model (Paci et al. [1]) to a mature model (ten Tusscher et al. [6]) and to fit the mature model to the immature model. In these inversions, we consider adjustments of the q Na , q CaL , q to , q Ks , q Kr , q K1 , q NaCa , q NaK , q pCa , q f , q bNa , q bCa , q leak , q up , and q rel factors (in addition to the q pK -factor for the ten Tusscher et al. model). Note that the I f current is added to the ten Tusscher et al. model in these simulations using the same formulation as in the default Paci et al. model, but with a conductance reduced by a factor of 10 for the mature ten Tusscher model, i.e. g f = 0.003 mS/µF. Because of the large number of free parameters, we conducted a more detailed inversion procedure in this case with twelve iterations and 15000 randomly chosen adjustment factors in each iteration. In addition, we included some additional terms in the cost function containing information that is not available from the optical measurements, but may be obtained from the mathematical models. More specifically, we used a cost function of the form where H 1 − H 8 are the same as in the remaining applications of the inversion procedure, that is Here, v max and c max denote the maximum value of the membrane potential and the calcium concentration, respectively. Similarly, v rest and c rest denote the resting membrane potential and calcium concentration, respectively, defined as the values obtained 10 ms before stimulation. Moreover, I 2 and I ∞ are defined as where n runs over all the time steps of an action potential. The currents I Na , I CaL , I Kr , and I K1 are chosen to be included in the cost function because we are especially interested in obtaining realistic behaviors for these currents since these are the currents considered in the drug inversions.

Identification of simulated single-channel block using H V andH Ca
In Figure 5 of the paper we showed the value of H V +Ca for pairwise perturbations of the maximum conductance of four major currents for simulated single channel block of each of the currents. Figures S5 and S6 show corresponding plots for the cost functions H V and H Ca , respectively. In these figures, we observe that the terms of H V seem to contain the main part of H V +Ca observed in Figure 5 of the paper. Figure S5: The cost function H V with ε = 0.2 evaluated for pairwise perturbations of the maximum conductances of four major currents for simulated single-channel block of each of the currents. In the upper panel, I Na is blocked by 50%, and in the next panels, I CaL , I Kr and I K1 are similarly blocked by 50%. Figure S6: Like Figure S5, except that we cosnider the cost function H Ca instead of H V .