Identifying Drug Response by Combining Measurements of the Membrane Potential, the Cytosolic Calcium Concentration, and the Extracellular Potential in Microphysiological Systems

Cardiomyocytes derived from human induced pluripotent stem cells (hiPSC-CMs) offer a new means to study and understand the human cardiac action potential, and can give key insight into how compounds may interact with important molecular pathways to destabilize the electrical function of the heart. Important features of the action potential can be readily measured using standard experimental techniques, such as the use of voltage sensitive dyes and fluorescent genetic reporters to estimate transmembrane potentials and cytosolic calcium concentrations. Using previously introduced computational procedures, such measurements can be used to estimate the current density of major ion channels present in hiPSC-CMs, and how compounds may alter their behavior. However, due to the limitations of optical recordings, resolving the sodium current remains difficult from these data. Here we show that if these optical measurements are complemented with observations of the extracellular potential using multi electrode arrays (MEAs), we can accurately estimate the current density of the sodium channels. This inversion of the sodium current relies on observation of the conduction velocity which turns out to be straightforwardly computed using measurements of extracellular waves across the electrodes. The combined data including the membrane potential, the cytosolic calcium concentration and the extracellular potential further opens up for the possibility of accurately estimating the effect of novel drugs applied to hiPSC-CMs.


INTRODUCTION
In recent reports (Tveito et al., 2018;Jaeger et al., 2020a) we have demonstrated how microphysiological systems utilizing human induced pluripotent stem cell derived cardiomyocytes (hiPSC-CMs) (Mathur et al., 2015;Mathur et al., 2016) can be used to estimate drug induced changes to the cardiac action potential using computational approaches. These methods use optical measurements of the membrane potential and the cytosolic calcium concentration to quantitate changes in underlying ion channel conductances and calcium handling pathways using a mathematical model of the hiPSC-CM dynamics. We have further shown how these estimates, at least in principle, carry over from immature cells to adult cardiomyocytes. This methodology provides information on a number of the major ion channels and when compared to data presented in (Mohammad et al., 1997;Di Stilo et al., 1998;Zhang et al., 1999;Zhabyeyev et al., 2000;Mirams et al., 2011;Kramer et al., 2013;Crumb et al., 2016), the method is able to provide reasonable estimates of the IC50 values of well-known drugs like Nifedipine, Lidocaine, Cisapride, Flecainide and Verapamil; see Table 2 of (Jaeger et al., 2020a). These drug affects the I CaL , I NaL or I Kr currents and the effect is well estimated by our methodology.
However, difficulties remain in the characterization of the fast sodium current, I Na . This is a major issue since this current more or less completely governs the rapid upstroke of the action potential and thus also the conduction velocity. Therefore, it is of great importance to characterize the effect of drugs on this current. The reason for this deficiency in our methodology is the time resolution of the data obtained by fluorescence; the data used in the inversions are provided with a resolution of 10 ms and this is far too coarse to be able to estimate the strength of I Na . Time resolution can be improved but at the cost of less accurate data and therefore another experimental technique is needed to pin down the channel density of and drug effects on I Na . It is well known that the extracellular potential can be measured in microphysical systems using multielectrode arrays; see, e.g., (Zwi et al., 2009;Clements and Thomas, 2014;Asakura et al., 2015;Bouyssier and Zemzemi, 2017;Tixier et al., 2018). In this report we will show that the extracellular data can be used to determine the sodium current. And therefore, by combining imaging data for the membrane potential (V) and the cytosolic calcium concentration (Ca) with data for the extracellular potential (U), we are able to identify both the fast sodium current and other major currents characterizing the action potential of the hiPSC-CMs.
The main challenge in combining V, Ca and U data is that a spatial problem needs to be resolved. When the data are given by V and Ca only, we have simply used a data trace obtained by taking the average over the whole chip (see Mathur et al., 2015;Tveito et al., 2018) and the inversion of the data has amounted to estimating parameters describing a system of ordinary differential equations. But when U is added to the data, the extracellular potential needs to be calculated. In our present implementation, we use the bidomain model (see, e.g., Franzone et al., 2014) for this purpose. The bidomain model has already been used for inversion of U data (but not V and Ca) by several authors; see (Bouyssier and Zemzemi, 2017;Raphel et al., 2017;Abbate et al., 2018;Tixier et al., 2018;Raphel et al., 2020). However, it is demonstrated in (Abbate et al., 2018) that the bidomain model does not provide an extracellular repolarization wave. Such a wave is clearly present in the experimental data and inhomogeneities have to be introduced in the bidomain model in order to enforce a repolarization wave. These inhomogeneities are difficult to obtain from measurements, and therefore we choose to use U only to estimate the currents involved in generating the upstroke and not the whole action potential. This turns out to determine I Na accurately, -at least in data generated by simulations.

METHODS
In this section, we describe the methods applied to identify drug response by combining measurements of the membrane potential, the cytosolic calcium concentration and the extracellular potential in microphysiological systems of hiPSC-CMs.

Bidomain-Base Model Simulations
In order to represent the electrical properties of a microphysiological system, we conduct simulations of the bidomain model of the form (see, e.g., Tung, 1978;Sundnes et al., 2007;Keener and James, 2009;Franzone et al., 2014). Here, v and u e are the membrane potential and the extracellular potential, respectively. In addition, σ i and σ e are the bidomain conductivities of the extracellular and intracellular spaces, respectively, C m is the specific membrane capacitance, and χ is the surface-to-volume ratio of the cell membrane. The values chosen for these parameters are given in Table 1. Furthermore, I ion represents the density of currents through different types of ion channels, pumps and exchangers on the cell membrane. We use an adjusted version of the hiPSC-CM base model introduced in (Jaeger et al., 2020a) to represent these currents. The current density is then given on the form where each I i represents the current through a specific type of ion channel, pump or exchanger. The hiPSC-CM base model includes a number of additional state variables representing the gating of ion channels and intracellular Ca 2+ concentrations. These variables are represented by s in the model above, and their dynamics are modeled by a set of ordinary differential equations (ODEs) given by F(v, s). A number of parameters in the hiPSC-CM base model have been adjusted to make the size of eight currents of particular interest (I Kr , I NaL , I CaL , I Na , I to , I Ks , I K1 , and I f ) close to the size of the currents in the model described in (Kernik et al., 2019), which is fitted to recordings of hiPSC-CMs from several different studies. The adjusted parameter values of the hiPSC-CM base model are given in Table 2, and the base model is, for completeness, given in the Supplementary Material. The geometry of the domain used to represent a microphysiological system is described in Figure 1 and Table 1. We consider a two-dimensional domain, and initiate a traveling wave by stimulating the lower left corner. Furthermore, 8 × 8 electrodes are distributed in the center of the domain, and we record the extracellular potential in these electrodes. On the boundary of the domain, we apply the Dirichlet boundary condition u e 0.
In addition to the bidomain simulations, we in some cases conduct simulations in which the spatial variation of the variables are ignored. In that case, the system Eqs 1-3 can be written as a pure ODE system of form Below, we refer to the system Eqs 5 and 6 as the pure ODE system.

Adjustment Factors
In order to investigate whether we are able to identify drug responses from combined measurements of V, Ca and U data, we perform numerical simulations representing a number of different drugs reported in (Crumb et al., 2016). More specifically, we simulate ion channel blockers by introducing adjustment factors −1 ≤ λ i ≤ 0, such that I ion is given by In our simulation, we consider the drug effect on three major ionic currents, I Kr , I CaL and I Na , believed to be of importance in evaluation of drug safety (see, e.g., Crumb et al., 2016), and we therefore introduce the four adjustment factors λ Kr , λ CaL and λ Na .

Bidomain-Base Model Simulations Used to Generate Data
We wish to investigate how V, Ca and U data from a microphysiological system can be used to identify the effect of drugs. In order to generate data representing this type of recordings, we perform bidomain-base model simulations. The simulation procedure used to generate the data is illustrated in Figure 2.
As illustrated in the upper panel of Figure 2, for a given combination of λ-values, we first perform bidomain-base model simulations for the duration of an entire action potential and store the extracellular potential, the membrane potential, and the cytosolic calcium concentration in each mesh point for each time step. We then wish to convert these recorded solutions into corresponding measurements that may be performed in microphysiological systems.
The first considered type of measurement is optical measurements of voltage-sensitive dyes. This type of measurement is mimicked in the computations by extracting the mean membrane potential over the entire domain for each time step. Because the exact conversion factor between the pixel intensity of the optical measurements and the associated membrane potential (in mV) is not known for this type of optical measurement, we normalize the mean V-trace to values between 0 and 1. Moreover, because the time resolution of the measurement equipment for the optical measurements typically is quite limited (see, e.g., Tveito et al., 2018), we only store the membrane potential at relatively large time steps (every 10 ms). Traces representing optical measurements of the cytosolic calcium concentration is generated in exactly the same manner.
In addition to the optical measurements of V and Ca, the extracellular potential may also be measured in electrodes located in the microphysiological system. In the simulations, we extract this type of measurements by storing the mean extracellular potential in the grid points overlapping the location of the electrodes. We are primarily interested in the extracellular potential during the depolarization wave, and we therefore only store the extracellular potential in the beginning of the simulation (typically the first 50 ms).

Pure ODE Simplification Used in the Inversion Procedure
During the inversion procedure, we wish to identify the λ-values associated with some considered data, generated as explained in the previous section. In this inversion procedure, we need to conduct simulations using a large number of different λ-values FIGURE 2 | Simulation procedure used to generate data for the inversions. In Step 1, we conduct bidomain-base model simulations, recording V, Ca and U for a number of time steps. In Step 2, we extract V and Ca data for inversion by recording the mean V and Ca over the domain at each time step and normalizing the traces to values between 0 and 1. In addition, U is recorded (in mV) in 8 × 8 electrodes at each time step during the depolarization wave.
Frontiers in Pharmacology | www.frontiersin.org February 2021 | Volume 11 | Article 569489 for comparison to the data. As mentioned above, the repolarization wave is known to be poorly represented by the bidomain model when applied to small collections of hiPSC-CMs. Therefore, we perform spatial simulations only for the start of the action potential. The data from this part of the simulation is used to determine the conduction velocity, the Ca transient time to peak and the terms of the cost function related to the extracellular potential, U. By performing spatial simulations only for the beginning of the action potential, we save time by avoiding a full bidomain simulation for all of the different parameter combinations. To generate the full V and Ca traces, we instead run a simple pure ODE simulation of the form Eqs 5 and 6 of the full action potential. The solution of this ODE simulation is used to compute all the terms in the cost function (see Section 2.2.1), except for the ones related to U, the Ca transient time to peak and the conduction velocity.

Stimulation Protocol and Update of Initial Conditions
In the pure ODE simulations, the cell is stimulated by a constant stimulus current of −5 μA/cm 2 until v is above −40 mV. In the bidomain simulations, we stimulate the domain in the lower left corner as illustrated in Figure 1 by a 5-ms-long constant membrane stimulus current of −40 μA/cm 2 . Note that after each parameter change, and before any bidomain or pure ODE simulation is performed, we conduct an ODE simulation for 10 AP cycles, stimulating at 1 Hz, to update the initial conditions.

Numerical Methods
The numerical simulations of the bidomain model are performed using an operator splitting procedure (see, e.g., Sundnes et al., 2005;Sundnes et al., 2007;Schroll et al., 2007). For each time step, we first solve the non-linear part of the system (i.e., Eqs 1 and 3 with the left-hand side of Eq. 1 set to zero), using the GRL2 method (Sundnes et al., 2009). Then, we solve the remaining linear part of the system (i.e., Eqs 1 and 2 with I ion 0) using a finite difference discretization in space and a backward Euler discretization in time. The discretization parameters used in the numerical simulations are given in Table 1. In the pure ODE simulations used in the inversion procedure we apply the same GRL2 method as in the non-linear part of the bidomain simulations, and in the pure ODE simulations used to update the initial conditions after a parameter change, we use the ode15s ODE-solver in Matlab.

Inversion Procedure
In this section, we will describe the inversion procedure applied to identify the effect of drugs based on V, Ca and U data in the form of adjustment factors, λ (see Section 2.1.1).

Cost Function
In the inversion procedure, we wish to find appropriate adjustment factors, λ, such that the solution of the model specified by λ is as similar as possible to the considered data.
To this end, we define a cost function H(λ), measuring the difference between the data and the model solution. This cost function is defined as where each H j (λ) represent various differences between the data and the model solution specified by λ, and w j are weights for each of these terms. In addition, δ is the weight of the regularization term ( i |λ i |) 2 , which is included so that small adjustments, λ, are preferred if several choices of λ result in almost identical solutions. The individual cost function terms, H j (λ), are defined below.

Action Potential and Calcium Transient Durations
A number of terms in the cost function measure differences in the action potential and calcium transient durations. These terms take the form for p 20, 30, . . ., 70, 80 The APDp value is measured as the time from V is p% below its maximum value during the upstroke of the action potential to the time at which V reaches a value below p% of the maximum during the repolarization phase; see Figure 3. APDp(λ) is the value obtained from the solution of the model given by the parameter vector λ, whereas APDp* is the value obtained from the data. The calcium transient durations, CaDp, are defined just like the action potential durations. Note also that the notation of a * marking the data values is used for all the terms in the cost function (see below).

Membrane Potential Integral
The cost function also includes a term of the form where Int30 is defined as and v is the membrane potential. The values t 1 and t 2 are here time points corresponding to when V is 30% below the maximum value during the depolarization and repolarization phases of the action potential (see Figure 3).

Calcium Transient Time to Peak
The typical time resolution of the V and Ca data (10 ms) is too coarse to be able to detect changes in the upstroke velocity of the action potential. However, the upstroke of the calcium transient is slower, and therefore, changes in the upstroke of the calcium transient may be detected. To measure this upstroke velocity, we include a term of the form where CaP is the time from Ca is 10% above its minimum value until it reaches its maximum (see Figure 3).
As mentioned above, we use a bidomain simulation instead of a pure ODE simulation to compute an estimate of CaP(λ) (and CaP* in the case of simulated data). We conduct this bidomain simulation for a rather short time interval (typically 50 ms) and record Ca in all the grid points. Then, we extract the grid points that have reached their peak Ca concentration during this simulation, and set up a normalized mean Ca-trace for these grid points. The value of CaP is then computed from this normalized mean Ca-trace. When the data is generated, we use the same procedure of considering the solution for only the first part of the simulation (e.g., 50 ms) when we compute CaP*.

Extracellular Potenial
In order to detect changes in U measured at the electrodes, we include a cost function term of the form where u k e is a vector denoting U (in mV) in electrode k for each of the recorded time steps.

Average Conduction Velocity
In addition, the cost function includes a term for the average conduction velocity of the form where cv denotes the average conduction velocity (in cm/s) computed from U recorded in the electrodes. More specifically, the average conduction velocity is computed as the distance between the center of the two electrodes marked as a and b in Figure 1 divided by the time between U crosses 0 mV after the peak in these two electrodes (see Figure 3).

Cost Function Weights
In the applications of the inversion procedure reported below, we have used the weight w j 1 for each of the cost function terms, except for H APD80 and H CaD80 , which are given the weight of 5. In addition, we use the value δ 0.01 for the regularization term. All terms of the cost function are scaled and thus have no unit, and therefore also w and δ are unit free.

Continuation-Based Optimization
In order to find the optimal λ-values minimizing the cost function H defined in Eq. 8, we apply a continuation-based minimization procedure. This approach is described in detail in (Jaeger et al., 2020a). In short, we attempt to find the optimal λ-values by gradually moving from a known solution λ 0 0 to the final λ fitting the data as best possible. To this end, we introduce a parameter θ that is gradually increased from 0 to 1, and define an intermediate cost function where H j (λ, θ) are adjusted versions of each of the cost function terms defined above and ξ is some small number (e.g., 10 −10 ). In general, H j (λ, θ) is defined as where R j is some property of the solution, for example an action potential duration. Moreover, R j (λ) is the property computed in the solution of the model defined by λ, and where R j (λ 0 ) is the property computed from the model solution defined by λ 0 0, and R * j is the property computed from the data we are trying to invert.
Considering the cost function terms, H j , defined above, all the terms may straightforwardly be defined on the form Eq. 17 except for the term H U defined in Eq. 14. For this term, we instead define FIGURE 3 | Illustration of some of the properties included in the cost function Eq. 8. From the V-trace, we include a number of APD-values (see Eq. 9), and the integral of V above APD30 (see Eq. 11). From the Ca-trace, we include a number of CaD-values (see Eq. 10) in addition to the Ca transient time to peak, CaP (see Eq. 13). From the U data, we include the value in each time point and each electrode (see Eq. 14) and the conduction velocity computed using the time points in which U crosses 0 mV after the peak in different electrodes (see Eq. 15 and Figure 1).

Frontiers in
where H 0 U is defined by Eq. 14 with u k,p e replaced by the model solution u k e (λ 0 ), and H U (λ) is defined by Eq. 14. From the definitions Eqs 16-19, we observe that for θ 0, R * j R j (λ 0 ) and H j (λ 0 , 0) 0, implying that the known solution λ λ 0 minimize Eq. 16 for θ 0. Furthermore, for θ 1, R * j (1) R * j defined by the data, H j (λ, 1) H j (λ), and H(λ, 1) H(λ). Consequently, for θ 1, the minimum of H defined in Eq. 16 is the same as the minimum of H defined in Eq. 8.
The advantage of the continuation algorithm is that we can move from the known optimal value λ 0 to the final optimal value λ gradually in a number of iterations, and that we for each iteration can assume that the new estimated λ is in the vicinity of the optimal λ from the previous iteration. In the applications of the inversion procedure reported below, we use four θ-iterations in the continuation algorithm (θ 1/4, θ 1/2, θ 3/4 and θ 1). In the first three iterations, m, we draw 63 random initial guesses of λ(θ m ) in the vicinity of λ(θ m−1 ), and in the last iteration, we draw 126 initial guesses. More specifically, the initial guesses for λ i (θ m ) are drawn from [λ i (θ m−1 ) − 0.2, λ i (θ m−1 ) + 0.2], where λ i (θ 0 ) is defined to be 0. From these initial guesses we minimize the cost function H(λ, θ) using the Nelder-Mead algorithm (Nelder and Mead, 1965). In the first three θ-iterations, we use 10 iterations of the Nelder-Mead algorithm for each initial guess, and in the last θ-iteration, we use 25 iterations.

Adjustment of Extracellular Concentrations
In order to better distinguish between different ion channel blockers, we consider the drug effects under different extracellular conditions. In particular, we vary the extracellular calcium concentration by introducing a number of known adjustment factors c Ca such that where [Ca 2+ ] e is the extracellular calcium concentration and [Ca 2+ ] * e is the default extracellular calcium concentration reported in Table 2.
In the inversions reported below, we use the two values c Ca,1 0 and c Ca,2 0.25. We generate V, Ca and U data using the approach described in Section 2.1.2 for both of these values of c Ca . Furthermore, for a given choice of λ in the inversion procedure, we compute a version of the solutions for each of these c Ca -values, and define the cost function where q counts the different extracellular conditions and H j (λ, c Ca,q ) represents the cost function terms computed using the extracellular calcium concentration defined by c Ca,q .

Measuring the Extracellular Potential
Experimentally, the extracellular potential was measured in a monolayer of cardiomyocytes using the microelectrode array system MED64-Basic with P515 electrode dishes. Cardiomyocytes (CM) were differentiated from the human induced-pluripotent stem cell (hiPSC) line Wild Type C (WTC, Coriell Repository, # GM25256) expressing the fluorescent calcium reporter GCaMP6f . Cells were differentiated into CM using a modified version of the Palecek protocol (Lian et al., 2012) applying 6 μm CHIR, 5 μM IWP4 and minus insulin media for 48 h each. PCR analysis of these cells, see , revealed MYL2 expression, which is indicative of an atrial phenotype, however at lower levels than found in adult atria. Similarly, the levels of MYL7 were higher than in adult atria but lower than in adult ventricles.
Overall, MLY2 and MLY7 expression together with the ventricular-like action potential waveform indicate these are immature ventricular cardiomyocytes (Veevers et al., 2018). MED64 electrode dishes were coated with polyethyleneimine followed by Matrigel and handled as suggested in the vendor application notes (MED64 Application Note, 2015; MED64 Application Note, 2016). CM were seeded to form a confluent cell layer and allowed to stabilize for 30 days with media changes of RPMI 1640 + B27 supplement (Gibco) every 2-3 days. Flecainide (Abcam, ab120504) doses were prepared in RPMI 1640 + B27 from a 25 mM stock in DMSO. Each drug dose (0, 1, 2.5 and 10 μM flecainide) was incubated for at least 30 min before measurements. Extracellular potential recordings were performed using the Mobius software (Version Win 7 0.5.1) with the template for spontaneous QT recording. Traces of spontaneous beating activity were recorded on all 64 electrodes and directly exported as raw data without any prefiltering or peak extraction. Recordings were performed on a heated stage (37°C sample temperature) in ambient atmosphere.

RESULTS
In this section, we report the results of some applications of the inversion procedure defined above. First, we investigate the sensitivity of the V, Ca and U data in response to perturbations of the I Kr , I CaL and I Na currents, and how this sensitivity may be increased when data for several different extracellular calcium concentrations are included. Next, we show some examples of extracellular repolarization waves, and explain why we have chosen to only consider the extracellular depolarization waves in the inversion procedure. Afterward, we show some examples of how bidomain-base model simulations are able to reproduce measured drug induced effects on the average conduction velocity, illustrating that including U data in the inversion procedure improves the identifiability of drug effects on I Na . Finally, we test the inversion procedure by using it to identify drug effects for a number of simulated drugs and investigate how the accuracy of the inversion is affected when noise is included in the simulated data.

Sensitivity of Currents
In Figure 4, we investigate the effect on the simulated V, Ca and U data of perturbing the I Kr , I CaL and I Na currents. We compute the data for the default base model (λ 0) and for two perturbations λ i −0.2 and λ i −0.4 for each of the currents i Kr, CaL, Na, corresponding to 20% and 40% block of the currents, respectively (see Section 2.1.1). Note that in the upper panel (investigating block of I Kr ), we have used a pacing frequency of 0.5 Hz instead of 1 Hz in order to allow for increased action potential durations resulting from block of I Kr .
We observe that block of the I Kr current results in increased action potential and calcium transient durations, whereas block of I CaL results in decreased action potential and calcium transient durations. Moreover, no effects are visible in the U data resulting from block of I Kr or I CaL (recall that the U data is only considered for the first 20 ms). For block of the I Na current, on the other hand, no effects on the action potential and calcium transient durations are visible, but there are clearly visible effects on the U data. Specifically, the amplitude of U is decreased and the timing of the peaks is delayed in response to block of the I Na current. Furthermore, in Figure 5, we report the conduction velocities computed from the U data in response to perturbations of the three currents. We observe that for perturbations of I Kr and I CaL , the effects on the conduction velocity are quite small, whereas block of I Na results in a considerably decreased conduction velocity.

Adjustment of Extracellular Concentrations
Because the effect on the V, Ca and U data of blocking I Kr is almost exactly the opposite of the effect of blocking I CaL (see Figure 4), we expect that determining the correct combination of block for these two currents will be difficult. Therefore, it could be useful to consider drug effects under different extracellular conditions in order to increase the chance of identifying the correct block for the two currents. Figure 6 shows the effect of block of the I Kr , I CaL and I Na currents for three different extracellular calcium concentrations. We consider the default extracellular calcium concentration specified in Table 2 (c Ca 0), a 10% increased concentration   In Figure 7, we show an example of a case where including an additional extracellular calcium concentration could help identify the correct block of I Kr and I CaL . In the upper panel, we use the default extracellular calcium concentration of Table 2. The solid line shows V and Ca for the default base model, and the dotted line shows the solution for a case with 20% block of I Kr and 33% block of I CaL . We observe that V and Ca look very similar in these two cases. As a result, the inversion procedure might mistake the  case of (λ Kr −0.2, λ CaL −0.33) to the case with no block. In the lower panel, however, we compare the two solutions for the case with an increased extracellular calcium concentration, and in this case, the difference between the solutions is more prominent. Consequently, the inversion procedure can more easily distinguish between the two cases when the solutions for both extracellular calcium concentrations are included. Figure 8A shows measurements of the extracellular potential in 64 electrodes recorded for collections of hiPSC-CMs. In the upper panel, we observe that there is some early activity, corresponding to the time of the upstroke of the action potential, in addition to a weaker signal after some hundred milliseconds. This weaker signal occurs at the time of repolarization of the action potential and is therefore referred  Table 1, except for σ i and σ e which are reported in the plot titles.

Estimating the Action Potential Duration
Frontiers in Pharmacology | www.frontiersin.org February 2021 | Volume 11 | Article 569489 10 to as the repolarization wave. In the lower panel, we zoom in on this repolarization wave, and we observe that the extracellular potential reaches a magnitude of up to about 0.1-0.2 mV in this period.
It has been demonstrated (see, e.g., Abbate et al., 2018), that bidomain simulations of small collections of hiPSC-CMs tend to give rise to very weak or even non-exsisting repolarization waves. Indeed, in the leftmost panel of Figure 8B, we plot the extracellular potential in a bidomain-base model simulation using the default parameter values of Table 1, and we observe that the magnitude of the wave in this case is very small, and completely invisible in the upper plot over the entire action potential duration. If we decrease the intracellular conductivity, σ i , on the other hand, a repolarization wave is visible, but the entire extracellular signal is quite weak (see the center panel of Figure 8B). However, if the extracellular conductivity, σ e , is decreased as well, the size of both the depolarization wave and the repolarization wave are quite similar to the recorded data (see the rightmost panel of Figure 8B).
In theory, the repolarization waves observed in the data and simulations could be used to estimate the action potential duration in the inversion procedure. However, as observed in the lower panels of Figure 8, the repolarization waves are quite smooth, and it is not clear which time points represent different degrees of repolarization. For the V-traces on the other hand (see, e.g., Figure 2), it is straightforward to define accurate measures of different degrees of repolarization in the form of APD-values (see Section 2.2.1 and Figure 3). Therefore, we use the optical measurements of V to define the action potential durations and only use the U data for information regarding the depolarization wave.

Estimating Drug Induced Effects on the Average Conduction Velocity
One of the advantages of including measurements of the extracellular potential in addition to optical measurements of V and Ca in the inversion procedure is that the extracellular measurements can be used to estimate the average conduction velocity of the cell collection. This information could be useful for determining drug effects on the I Na current (see Figure 5). The left column of Table 3 reports the conduction velocity computed from measurements of the extracellular potential in a collection of hiPSC-CMs exposed to different doses of the drug Flecainide (measured data, not simulated). We observe that as the drug dose is increased, the conduction velocity is decreased. This could indicate that the drug blocks the I Na current, which has also been found in previous studies (Kramer et al., 2013;Crumb et al., 2016).
In the paper (Jaeger et al., 2020a), we estimated IC50-values for I CaL , I NaL and I Kr based on optical measurements of V and Ca of hiPSC-CMs exposed to Flecainide, but we were not able to estimate the effect on I Na due to the low time resolution for the optical measurements. The IC50-values were estimated to IC50 CaL 9 μM, IC50 NaL 47 μM, and IC50 Kr 1.9 μM, where the conductance of each current was scaled according to where g d i is the conductance of current i in presence of the drug dose D and g c i is the conductance in the control case with no drug present. Furthermore, h i is the so-called Hill coefficient, assumed to be one in (Jaeger et al., 2020a). Incorporating these IC50-values in addition to some estimated IC50-values for I Na in bidomainbase model simulations, we obtain the conduction velocities reported in the center and right columns of Table 3. In the center column we use IC50 Na 2.2 μM and h Na 1, and in the right column we use IC50 Na 0.2 μM and h Na 0.4. The IC50 Na value in the case of h Na 1 compares relatively well with literature values; 6.7, 6.2 and 4.4 μM from Crumb et al., (2016), Kramer et al., (2013) and Qu and Vargas, (2015), respectively.
Note that since high doses of Flecainide result in increased action potential durations (see, e.g., Jaeger et al., 2020a), we have used a pacing frequency of 0.5 Hz instead of 1 Hz when updating the initial conditions for these simulations. In addition, to match the conduction velocity in the control case, the default value of g Na value of Table 2 is increased by 45%. We observe that the bidomain-base model simulations are able to roughly reproduce the drug induced reduction in conduction velocity observed in the measurements, indicating that a comparison of measured and simulated conduction velocities could help identify drug effects on I Na .
The data and model solutions are further compared in Figure 9. Here, we show the U solutions and the measured U data in the 64 electrodes for the control case and for the 10 μM dose case at some different points in time. In the simulations, we use IC50 Na 0.2 μM and h Na 0.4, and we have adjusted the stimulation location to correspond to the propagation direction observed in the measured data. In the control case, the wave moves in the y-direction from the lower part of the domain to the upper part, and in the drug case, the wave moves from the upper left corner to the lower right corner of the domain. In both the data and the simulations, we observe that the extracellular depolarization wave moves more slowly across the domain for 10 μM Flecainide than in the control case, consistent with the reduced conduction velocities observed in Table 3. The drug effect on I Na is set to IC50 Na 2.2 μM, h Na 1 in the center column and IC50 Na 0.2 μM, h Na 0.4 in the right column. The parameters used in the simulations are given in Tables 1 and 2, and the default g Na value is increased by 45% to match the conduction velocity in the control case.

Inversion of Simulated Drugs
In order to investigate how well the inversion procedure outlined above is able to identify the effect of drugs on I Kr , I CaL and I Na , we generate simulated data for twelve drugs whose effect on I Kr , I CaL and I Na was investigated in (Crumb et al., 2016). The block percentages used to generate the data are based on the block percentages for drug concentrations corresponding to three times the free plasma C max reported in (Crumb et al., 2016). However, large block percentages for I Kr have been reduced in order to obtain reasonable V and Ca traces. In the inversions, we use two extracellular calcium concentrations, as explained in Section 2.3. Furthermore, we save U data and information about the Ca peak time extracted from the first 50 ms of the bidomain simulations for all the drugs, except for the drug Diltiazem. For Diltiazem we had to increase the bidomain simulation time to 100 ms in order to ensure that some of the grid points in the domain had reached their peak calcium concentration to compute the Ca peak time (see Section 2.2.1).
In Figure 10, we compare the block percentages estimated by the inversion procedure to those used to generate the data for the twelve drugs. We observe that for all the considered drugs, the inversion procedure is able to identify the block of the three currents quite accurately.

Effect of Noise
In Figure 10, we considered how well the inversion procedure was able to identify drugs based on data generated by bidomanbase model simulations, and we observed that the inversion procedure was able to identify the correct channel blocking quite accurately. However, when the V, Ca and U data are recorded from real measurements from microphysiological systems, the data will include some noise. In order to investigate how well the inversion procedure is able to identify the effect of drugs from data including noise, we include 5% noise in the V data, 3% noise in the Ca data and 1% noise in the U data and repeat the inversions shown in Figure 10. The noise is added by drawing a random number between −p·A and p·A for each point in time, where A is the difference between the maximum and minimum value of the considered data and p is the noise percentage. This random number is then added to the V, Ca and U traces.
The result of the inversions with noise included in the data is given in Figure 11. We observe that the inversion procedure is able to estimate the block percentage quite accurately in most cases, but that the accuracy is reduced compared to the case with no noise in Figure 10.

Sensitivity of Parameters Is Necessary for Identification
A generic problem in mathematical models of physiology is to determine the parameters included in the model. For models of the action potential of excitable cells, there is a large number of parameters that needs to be determined in order to use the model; see, e.g., (Sobie, 2009;Otte et al., 2016;Mora et al., 2017;Jaeger et al., 2019b). If a model is insensitive to changes in a specific parameter, that parameter is impossible to determine by comparing experimental results and results of numerical simulations. Contrarily, if the model is sensitive to a parameter, numerical experiments can be used to match the model to the data by fine-tuning the parameter. In Figure 4 above, we saw that the membrane potential is clearly sensitive to changes in the I Kr and I CaL currents. The upstroke is also sensitive to changes in the I Na current, but since the upstroke is very fast, the sensitivity is difficult to see in the plot of the entire AP. Since optical measurements have a relatively coarse time resolution, we have been unable to identify the strength of the sodium current based solely on voltage and calcium traces. However, we note from Figure 4 that the extracellular potential is indeed sensitive to changes in the sodium current. This sensitivity also carries over to sensitivity of the average conduction velocity (CV) with respect to changes in the sodium current, and this sensitivity is utilized in the cost function; see Eq. 15.

The Conduction Velocity Is Governed by the Sodium Current
The conduction of the electrochemical signal through the cardiac muscle is essential for the functioning of the heart. Surprisingly, the conduction process is still not completely understood; see, e.g., (Veeraraghavan et al., 2014), for a review of the development of the understanding of cardiac conduction. Globally (mm scale), cardiac conduction is usually modeled using homogenized models like the bidomain or monodomain models (see, e.g., FIGURE 10 | Result of the inversion procedure applied to simulated data of twelve drugs, based on the block percentages corresponding to three times the free plasma C max reported in (Crumb et al., 2016).
FIGURE 11 | Result of the inversion procedure applied to simulated data of the twelve drugs of Figure 10.We have here included 5% noise in the V data, 3% noise in the Ca data and 1% noise in the U data.
Frontiers in Pharmacology | www.frontiersin.org February 2021 | Volume 11 | Article 569489 Sundnes et al., 2007;Keener and James, 2009;Franzone et al., 2014), but locally (μm scale) more detailed models are used; see, e.g., (Kucera et al., 2002;Kucera et al., 2017;Jaeger et al., 2019a). The microphysiological system is somewhere in between these scales (the size is about 1.2 mm × 1.2 mm) and we would have liked to use the detailed EMI model introduced for cardiac conduction in (Tveito et al., 2017). However, the EMI approach is challenging from a computational point of view and we have therefore used the much simpler bidomain model to consider how the CV depends on changing the different ion currents.
In Figure 5 we observe that the CV is very sensitive to changes of the I Na current, but almost insensitive to changes in the I Kr and I CaL currents. This means that we can use the measurements of U to identify the I Na current first and then only look for the I Kr and I CaL currents based on the V-and Ca-traces.

Improving Visibility of Parameters by Changing the Extracellular Concentrations
In using experimental data to determine parameters of a computational model, it is desirable to use protocols designed to highlight the effect of different parameters. For instance, it is argued in (Groenendaal et al., 2015) that random stimulation protocols can improve visibility of the parameters, and this is confirmed in (Jaeger et al., 2019b). Another experimental parameter that can be changed is the ionic concentrations of the extracellular environment. In Figures 6 and 7 we show that using two different extracellular calcium concentration can improve identifiability of the model parameters. One particular important effect of this is that multiple extracellular calcium concentration can aid in distinguishing changes of the I Kr and I CaL currents. Since blocking I Kr has more or less the same effect as increasing I CaL , it is difficult to distinguish these effects in measurements of the AP. But it is apparent from Figure 7 that the level of blocking changes with the extracellular calcium concentration, and this can be used to distinguish between different blocking combinations for these two currents.

Inversion of Simulated Drugs
Simulated data is often used as a proxy for real data when the real data are cumbersome to obtain. Here, we have used real data for U both for the case of no drug and when various doses of Flecainide have been applied. Furthermore, the base model used to represent the membrane dynamics has been parameterized using real data; see (Tveito et al., 2018;Jaeger et al., 2020a). However, in order to study inversion for the range of data provided in the CIPA report (Crumb et al., 2016), we have used simulated data. One advantage of this is that we can get data to any desired accuracy and that we can control the level of noise inevitably introduced in the measurements. The disadvantage is clearly the reduced realism of the data and it is a priority of future work to use combined V, Ca, and U data to do inversion of both well characterized and novel drugs with hitherto unknown properties.
For simulated data, we notice that the inversion procedure provides quite accurate estimates. Systematically, we observe a tendency to overestimate the block of both the I Kr and I CaL currents. As alluded to above, it is notoriously difficult to distinguish reduced I Kr from increased I CaL and vice versa. This can be improved (Figures 6 and 7), by using several values of the extracellular calcium concentration, but the problem is not completely removed.

Repolarization Waves in Bidomain Simulations
As mentioned above, the bidoman model has been used by many authors (see, e.g., Abbate et al., 2018;Bouyssier and Zemzemi, 2017;Raphel et al., 2017;Raphel et al., 2020;Tixier et al., 2018) to simulate the electrophysiology of collections of hiPSC-CMs. One problem pointed out by several authors is the lack of a repolarization wave in the simulated results although the repolarization wave is clearly present in measurements of the extracellular potential; see, e.g., Figures 3 and 4 of (Abbate et al., 2018). This feature of the bidomain solution is repaired by introducing heterogeneities in the tissue which lead to a repolarization wave. In Figure 8A, we show that there is indeed a repolarization wave present in the extracellular data obtained from collections of hiPSC-CMs. The repolarization wave is also present in our bidomain simulations, but the strength of the wave depends critically on the intracellular conductivity, σ i . The repolarization wave becomes stronger as σ i is reduced. Since the intracellular conductivity represents the geometrical average of the intercellular conductivity (regulated by gap junctions) and the cytosolic conductivity, it is reasonable to use a reduced value of σ i for immature cells, since the gap junctions are most likely less developed in collections of hiPSC-CMs than for collections of adult cardiomyocytes. Experimentally, this has been observed in these cells, with indications of incomplete maturity and underdeveloped gap junctions by connexin-43 staining that localized in small speckles  rather than distinct disks as found in the adult human heart. Furthermore, immature tissue geometry, in terms of cell alignment, shape, and orientation of connectivity may also change these dynamics, and while these cells have demonstrated good agreement of geometric cell parameters, such as alignment within a microphysiological system compared to the adult human heart , more study is needed on the relative effects of these features on conduction properties. In (Abbate et al., 2018), it is argued that it is reasonable to introduce heterogeneities in the tissue and we agree. However, our results indicate that it is not necessary in order to see a repolarization wave in the simulation results.

CONCLUSION
We have shown that by using data traces of the membrane potential, the intracellular calcium concentration and the extracellular potential, we can estimate the major sodium, calcium and potassium currents in the base model. It remains to enable concurrent observation of all three modalities, but when such data become available, the methodology described in the present report may be used to invert the data and thus obtain channel densities and estimate drug effects on the channels.