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.

1 Introduction

In recent reports (Tveito et al., 2018; Jæger 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 (Jæger 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.

2 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.

2.1 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

$\chi C_m\left(\frac{\partial v}{\partial t}+I_{\text{ion}}\left(v,s\right)+I_{\text{stim}}\right)=\nabla \cdot \left({\sigma }_i\nabla v\right)+\nabla \cdot \left({\sigma }_i\nabla u_e\right),(1)$

$0=\nabla \cdot \left({\sigma }_i\nabla v\right)+\nabla \cdot \left[\left({\sigma }_i+{\sigma }_e\right)\nabla u_e\right],(2)$

$\frac{\partial s}{\partial t}=F\left(v,s\right)(3)$

(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.

TABLE 1

TABLE 1. Parameter values used in the bidomain-base model simulations.

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 (Jæger et al., 2020a) to represent these currents. The current density is then given on the form

$I_{\text{ion}}\left(v,s\right)={\displaystyle \sum }_i{I}_i\left(v,s\right),(4)$

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²⁺ 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.

TABLE 2

TABLE 2. Parameter values used in the hiPSC-CM base model. The remaining parameter values are the same as specified in (Jæger et al., 2020a).

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.

FIGURE 1

FIGURE 1. Geometry used in the bidomain-base model simulations. The domain contains 8 × 8 evenly spaced electrodes, and the propagating wave is initiated by stimulating the lower left corner of the domain. The conduction velocity is computed between the electrodes marked as a and b, and when we plot the extracellular potential of just a single electrode, we consider the electrode marked as c.

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

$\frac{dv}{dt}=-\frac{1}{C_m}\left(I_{\text{ion}}+I_{\text{stim}}\right),(5)$

$\frac{ds}{dt}=F\left(v,s\right).(6)$

Below, we refer to the system Eqs 5 and 6 as the pure ODE system.

2.1.1 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

$I_{\text{ion}}\left(v,s\right)={\displaystyle \sum }_i\left(1+{\lambda }_i\right)I_i\left(v,s\right).(7)$

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.

2.1.2 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.

FIGURE 2

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.

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).

2.1.3 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 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.

2.1.4 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² 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². 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.

2.1.5 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.

2.2 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).

2.2.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

$H\left(\lambda \right)={\displaystyle \sum }_j\left(w_j{H}_j{\left(\lambda \right)}^2\right)+\delta {\left({\displaystyle \sum }_i\left|{\lambda }_i\right|\right)}^2,(8)$

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 ${\left({\displaystyle {\sum }_i\left|{\lambda }_i\right|}\right)}^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.

2.2.1.1 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

$H_{\text{APD}p}\left(\lambda \right)=\frac{|\text{APD}p\left(\lambda \right)-\text{APD}{p}^{\text{*}}|}{\left|\text{APD}{p}^{\text{*}}\right|},(9)$

$H_{\text{CaD}p}\left(\lambda \right)=\frac{|\text{CaD}p\left(\lambda \right)-\text{CaD}{p}^{\text{*}}|}{\left|\text{CaD}{p}^{\text{*}}\right|},(10)$

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).

FIGURE 3

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).

2.2.1.2 Membrane Potential Integral

The cost function also includes a term of the form

$H_{\text{Int}30}\left(\lambda \right)=\frac{|\text{Int}30\left(\lambda \right)-\text{Int}{30}^{\text{*}}|}{\left|\text{Int}{30}^{\text{*}}\right|},(11)$

where Int30 is defined as

$\text{Int}30={{\displaystyle \int }}_{t_1}^{t_2}[v-v\left(t_1\right)]dt(12)$

and v is the membrane potential. The values t₁ and t₂ 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).

2.2.1.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

$H_{t,\text{Ca}}\left(\lambda \right)=\frac{|\text{CaP}\left(\lambda \right)-{\text{CaP}}^{\text{*}}|}{\left|{\text{CaP}}^{\text{*}}\right|},(13)$

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*.

2.2.1.4 Extracellular Potenial

In order to detect changes in U measured at the electrodes, we include a cost function term of the form

$H_{\text{U}}\left(\lambda \right)=\frac{{{\displaystyle \sum }}_k{\Vert u_e^k\left(\lambda \right)-u_e^{k,\ast }\Vert }_2}{{{\displaystyle \sum }}_k{\Vert u_e^{k,\ast }\Vert }_2},(14)$

where $u_e^k$ is a vector denoting U (in mV) in electrode k for each of the recorded time steps.

2.2.1.5 Average Conduction Velocity

In addition, the cost function includes a term for the average conduction velocity of the form

$H_{\text{cv}}\left(\lambda \right)=\frac{|\text{cv}\left(\lambda \right)-{\text{cv}}^{\text{*}}|}{\left|{\text{cv}}^{\text{*}}\right|},(15)$

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).

2.2.1.6 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.

2.2.2 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 (Jæger et al., 2020a). In short, we attempt to find the optimal λ-values by gradually moving from a known solution λ⁰ = 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

$\overline{H}\left(\lambda ,\theta \right)={\displaystyle \sum }_j\left(w_j{\overline{H}}_j{\left(\lambda ,\theta \right)}^2\right)+\frac{\delta }{\text{max}\left(\theta ,\xi \right)}{\left({\displaystyle \sum }_i\left|{\lambda }_i\right|\right)}^2,(16)$

where ${\overline{H}}_j\left(\lambda ,\theta \right)$ are adjusted versions of each of the cost function terms defined above and ξ is some small number (e.g., 10⁻¹⁰). In general, ${\overline{H}}_j\left(\lambda ,\theta \right)$ is defined as

${\overline{H}}_j\left(\lambda ,\theta \right)=\frac{|R_j\left(\lambda \right)-R_j^{\text{*}}\left(\theta \right)|}{\left|R_j^{\text{*}}\left(\theta \right)\right|},(17)$

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 $R_j^{\text{*}}\left(\theta \right)$ is defined as

$R_j^{\text{*}}\left(\theta \right)=\left(1-\theta \right)R_j\left({\lambda }_0\right)+\theta R_j^{\text{*}},(18)$

where R_j(λ₀) is the property computed from the model solution defined by λ⁰ = 0, and $R_j^{\text{*}}$ 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

${\overline{H}}_{\text{U}}\left(\lambda ,\theta \right)=\left(1-\theta \right)H_{\text{U}}^0+\theta H_{\text{U}}\left(\lambda \right),(19)$

where $H_{\text{U}}^0$ is defined by Eq. 14 with $u_e^{k,\ast }$ replaced by the model solution $u_e^k\left({\lambda }_0\right)$, and H_U(λ) is defined by Eq. 14.

From the definitions Eqs 16–19, we observe that for θ = 0, $R_j^{\text{*}}=R_j\left({\lambda }_0\right)$ and ${\overline{H}}_j\left({\lambda }_0,0\right)=0$, implying that the known solution λ = λ₀ minimize Eq. 16 for θ = 0. Furthermore, for θ = 1, $R_j^{\text{*}}\left(1\right)=R_j^{\text{*}}$ defined by the data, ${\overline{H}}_j\left(\lambda ,1\right)=H_j\left(\lambda \right)$, and $\overline{H}\left(\lambda ,1\right)=H\left(\lambda \right)$. Consequently, for θ = 1, the minimum of $\overline{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 λ₀ 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(θ₀) is defined to be 0. From these initial guesses we minimize the cost function $\overline{H}\left(\lambda ,\theta \right)$ 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.

2.3 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 γ_Ca such that

${\left[{\text{Ca}}^{2+}\right]}_e=\left(1+{\gamma }_{\text{Ca}}\right){\left[{\text{Ca}}^{2+}\right]}_e^{\text{*}},(20)$

where [Ca²⁺]_e is the extracellular calcium concentration and ${\left[{\text{Ca}}^{2+}\right]}_e^{\text{*}}$ is the default extracellular calcium concentration reported in Table 2.

In the inversions reported below, we use the two values γ_Ca,1 = 0 and γ_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 γ_Ca. Furthermore, for a given choice of λ in the inversion procedure, we compute a version of the solutions for each of these γ_Ca-values, and define the cost function

$H(\lambda )={\displaystyle \sum _q{\displaystyle \sum _j\left(w_j{H}_j{(\lambda ,{\gamma }_{\text{Ca},q})}^2\right)}}+\delta {\left({\displaystyle \sum _i}\left|{\lambda }_i\right|\right)}^2,(21)$

where q counts the different extracellular conditions and H_j(λ, γ_Ca,q) represents the cost function terms computed using the extracellular calcium concentration defined by γ_Ca,q.

2.4 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 (Huebsch et al., 2015). 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 (Huebsch et al., 2018), 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 pre-filtering or peak extraction. Recordings were performed on a heated stage (37°C sample temperature) in ambient atmosphere.

3 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.

3.1 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.

FIGURE 4

FIGURE 4. Effect of perturbing I_Kr, I_CaL and I_Na on the V, Ca and U data. For the U data, we show the extracellular potential in the electrode marked as c in Figure 1 and the time scale is zoomed in on the first 20 ms of the simulation.

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.

FIGURE 5

FIGURE 5. Effect of perturbing I_Kr, I_CaL and I_Na on the computed conduction velocity. The conduction velocities are computed from the U data as explained in Section 2.2.1.

3.2 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 (γ_Ca = 0), a 10% increased concentration (γ_Ca = 0.1), and a 25% increased concentration (γ_Ca = 0.25). The solid lines show the default solutions for each of these extracellular environments, and we observe that the action potential and calcium transient durations are increased as the extracellular calcium concentration is increased. Furthermore, the dotted lines show the solutions corresponding to 20% block of the considered currents. For all the considered U data and for the V and Ca data for block of I_Na, the different extracellular calcium concentrations do not seem to have a significant effect. However, for the V and Ca data for block of I_Kr and I_CaL, we observe that the effect of the block varies for the different concentrations. For example, the effect of block of I_Kr is more prominent for an increased extracellular calcium concentration.

FIGURE 6

FIGURE 6. Effect of perturbing I_Kr, I_CaL and I_Na on the V, Ca and U data for three different adjustments, γ_Ca, of the extracellular calcium concentration (see Section 2.3). For the U data, we show the extracellular potential in the electrode marked as c in Figure 1 and the time scale is zoomed in on the first 20 ms of the simulation.

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 7

FIGURE 7. Illustration of how including an extra extracellular calcium concentration adjustment may improve the identifiability of I_Kr and I_CaL. In the upper panel, we compare V and Ca for two different block combinations for the default extracellular calcium concentration (γ_Ca = 0), and in the lower panel we consider the solutions for a 25% increased extracellular calcium concentration (γ_Ca = 0.25). For γ_Ca = 0, the two solutions are very similar, but the difference is increased for γ_Ca = 0.25.

3.3 Estimating the Action Potential Duration

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 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.

FIGURE 8

FIGURE 8. Repolarization wave in hiPSC-CM data and bidomain-base model simulations. The upper panel shows measurements of U during an action potential for three different data sets. The traces for the 64 electrodes are overlaid in the plots. In the second panel, we zoom in on the repolarization wave. The third panel shows U in three bidomain-base model simulations, and the bottom panel shows the corresponding repolarization waves. The parameters used in the simulations are given in Table 1, except for σ_i and σ_e which are reported in the plot titles.

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.

3.4 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).

TABLE 3

TABLE 3. Effect of the drug Flecainide on the conduction velocity computed from the extracellular potential.

In the paper (Jæger 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

$g_i^d=\left[\frac{1}{1+{\left(\frac{D}{\text{IC}{50}_i}\right)}^{h_i}}\right]\cdot g_i^c,\text{for }i=\text{ CaL},\text{ NaL},\text{ Kr},(22)$

where $g_i^d$ is the conductance of current i in presence of the drug dose D and $g_i^c$ 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 (Jæger et al., 2020a). Incorporating these IC50-values in addition to some estimated IC50-values for I_Na in bidomain-base 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., Jæger 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.

FIGURE 9

FIGURE 9. Measured and simulated U data in 64 electrodes for the control case and for 10 μM Flecainide, modeled using IC50_CaL = 9 μM, IC50_NaL = 47 μM, IC50_Kr = 1.9 μM, IC50_Na = 0.2 μM, h_i = 1 for i = CaL, NaL, and Kr, and h_Na = 0.4. The parameters used in the simulations are given in Tables 1 and 2, and the default g_Na value is increased by 45%. The stimulation location is adjusted to correspond to the propagation direction in the data.

3.5 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.

FIGURE 10

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).

3.6 Effect of Noise

In Figure 10, we considered how well the inversion procedure was able to identify drugs based on data generated by bidoman-base 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.

FIGURE 11

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.

4 Discussion

4.1 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; Jæger 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.

4.2 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., 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; Jæger 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.

4.3 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 (Jæger 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.

4.4 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; Jæger 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.

4.5 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 (Huebsch et al., 2018) 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 (Huebsch et al., 2018), 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.

5 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.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

KJ and AT developed the mathematical framework, performed the numerical experiments and wrote the paper. SW carefully checked all parts of the paper and improved the wording of the paper. VC and KH generated data of hiPSC-CMs from microphysiological systems. All authors read and approved the final manuscript.

Funding

We would like to acknowledge the following funding sources: The Research Council of Norway funded INTPART Project 249885, the SUURPh program funded by the Norwegian Ministry of Education and Research, the Peder Sather Center for Advanced Study, NIH-NCATS UH3TR000487, NIH-NHLBI 2485 HL130417, and the California Institute for Regenerative Medicine DISC2-10090. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Conflict of Interest

KJ, AT, SW, and KH have financial relationships with Organos Inc., and the company may benefit from commercialization of the results of this research.

The remaining author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This manuscript has been released as a pre-print at bioRxiv, (Jæger et al., 2020b).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2020.569489/full#supplementary-material

References

Abbate, E., Boulakia, M., Coudière, Y., Gerbeau, J. F., Zitoun, P., and Zemzemi, N. (2018). In silico assessment of the effects of various compounds in mea/hipsc-cm assays: modeling and numerical simulations. J. Pharmacol. Toxicol. Methods 89, 59–72. doi:10.1016/j.vascn.2017.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Asakura, K., Hayashi, S., Ojima, A., Taniguchi, T., Miyamoto, N., Nakamori, C., et al. (2015). Improvement of acquisition and analysis methods in multi-electrode array experiments with ips cell-derived cardiomyocytes. J. Pharmacol. Toxicol. Methods 75, 17–26. doi:10.1016/j.vascn.2015.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Bouyssier, J., and Zemzemi, N. (2017). “Parameters estimation approach for the mea/hipsc-cm asaays,” in 2017 computing in cardiology (CinC), Rennes, France, September 24–27, 2017 (Lille, France: IEEE), 1–4.

CrossRef Full Text | Google Scholar

Clements, M., and Thomas, N. (2014). High-throughput multi-parameter profiling of electrophysiological drug effects in human embryonic stem cell derived cardiomyocytes using multi-electrode arrays. Toxicol. Sci 140 (2), 445–461. doi:10.1093/toxsci/kfu084

PubMed Abstract | CrossRef Full Text | Google Scholar

Crumb, W. J., Vicente, J., Johannesen, L., and Strauss, D. G. (2016). An evaluation of 30 clinical drugs against the comprehensive in vitro proarrhythmia assay (CiPA) proposed ion channel panel. J. Pharmacol. Toxicol. Methods 81, 251–262. doi:10.1016/j.vascn.2016.03.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Di Stilo, A., Visentin, S., Cena, C., Gasco, A. M., Ermondi, G., and Gasco, A. (1998). New 1,4-dihydropyridines conjugated to furoxanyl moieties, endowed with both nitric oxide-like and calcium channel antagonist vasodilator activities. J. Med. Chem 41 (27), 5393–5401. doi:10.1021/jm9803267

PubMed Abstract | CrossRef Full Text | Google Scholar

Franzone, P. C., Pavarino, L. F., and Scacchi, S. (2014). Mathematical cardiac electrophysiology. Cham, Switzerland: Springer.

CrossRef Full Text | Google Scholar

Groenendaal, W., Ortega, F. A., Kherlopian, A. R., Zygmunt, A. C., Krogh-Madsen, T., and Christini, D. J. (2015). Cell-specific cardiac electrophysiology models. PLoS Comput. Biol 11 (4), e1004242. doi:10.1371/journal.pcbi.1004242

PubMed Abstract | CrossRef Full Text | Google Scholar

Huebsch, N., Loskill, P., Mandegar, M. A., Marks, N. C., Sheehan, A. S., Ma, Z., et al. (2015). Automated video-based analysis of contractility and calcium flux in human-induced pluripotent stem cell-derived cardiomyocytes cultured over different spatial scales. Tissue Eng. C Methods 21 (5), 467–479. doi:10.1089/ten.TEC.2014.0283

PubMed Abstract | CrossRef Full Text | Google Scholar

Huebsch, N., Charrez, B., Siemons, B., Boggess, S. C., Wall, S., Charwat, V., et al. (2018). Metabolically-driven maturation of hipsc-cell derived heart-on-a-chip. Available at: https://www.biorxiv.org/content/10.1101/485169v3.

CrossRef Full Text | Google Scholar

Jæger, K. H., Edwards, A. G., McCulloch, A., and Tveito, A. (2019a). Properties of cardiac conduction in a cell-based computational model. PLoS Comput. Biol 15 (5), e1007042. doi:10.1371/journal.pcbi.1007042

PubMed Abstract | CrossRef Full Text | Google Scholar

Jæger, K. H., Wall, S., and Tveito, A. (2019b). Detecting undetectables: can conductances of action potential models be changed without appreciable change in the transmembrane potential? Chaos 29 (7), 073102. doi:10.1063/1.5087629

PubMed Abstract | CrossRef Full Text | Google Scholar

Jæger, K. H., Charwat, V., Charrez, B., Finsberg, H., Maleckar, M. M., and Wall, S. (2020a). Improved computational identification of drug response using optical measurements of human stem cell derived cardiomyocytes in microphysiological systems. Front. Pharmacol. 10, 1648. doi:10.3389/fphar.2019.01648

CrossRef Full Text | Google Scholar

Jæger, K. H., Charwat, V., Wall, S., Healy, K. E., and Tveito, A. (2020b). Identifying drug response by combining measurements of the membrane potential, the cytosolic calcium concentration, and the extracellular potential in microphysiological systems. Available at: https://www.biorxiv.org/content/10.1101/2020.05.29.122747v1.

CrossRef Full Text | Google Scholar

Keener, J. P., and James, S. (2009). Mathematical physiology. New York, NY: Springer.

CrossRef Full Text | Google Scholar

Kernik, D. C., Morotti, S., Wu, H. D., Garg, P., Henry, J. D., Kurokawa, J., et al. (2019). A computational model of induced pluripotent stem-cell derived cardiomyocytes incorporating experimental variability from multiple data sources. J. Physiol 597 (17), 4533–4564. doi:10.1113/JP277724

PubMed Abstract | CrossRef Full Text | Google Scholar

Kramer, J., Obejero-Paz, C. A., Myatt, G., Kuryshev, Y. A., Bruening-Wright, A., et al. (2013). MICE models: superior to the HERG model in predicting Torsade de Pointes. Sci. Rep 3, 2100. doi:10.1038/srep02100

PubMed Abstract | CrossRef Full Text | Google Scholar

Kucera, J. P., Rohr, S., and Kleber, A. G. (2017). Microstructure, cell-to-cell coupling, and ion currents as determinants of electrical propagation and arrhythmogenesis. Circ. Arrhythm. Electrophysiol 10 (9), e004665. doi:10.1161/CIRCEP.117.004665

PubMed Abstract | CrossRef Full Text | Google Scholar

Kucera, J. P., Rohr, S., and Rudy, Y. (2002). Localization of sodium channels in intercalated disks modulates cardiac conduction. Circ. Res 91 (12), 1176–1182. doi:10.1161/01.res.0000046237.54156.0a

PubMed Abstract | CrossRef Full Text | Google Scholar

Lian, X., Hsiao, C., Wilson, G., Zhu, K., Hazeltine, L. B., Azarin, S. M., et al. (2012). Robust cardiomyocyte differentiation from human pluripotent stem cells via temporal modulation of canonical wnt signaling. Proc. Natl. Acad. Sci. U.S.A 109 (27), E1848–E1857. doi:10.1073/pnas.1200250109

PubMed Abstract | CrossRef Full Text | Google Scholar

Mathur, A., Loskill, P., Shao, K., Huebsch, N., Hong, S., Marcus, S. G., et al. (2015). Human iPSC-based cardiac microphysiological system for drug screening applications. Sci. Rep 5, 8883. doi:10.1038/srep08883

PubMed Abstract | CrossRef Full Text | Google Scholar

Mathur, A., Ma, Z., Loskill, P., Jeeawoody, S., and Healy, K. E. (2016). In vitro cardiac tissue models: current status and future prospects. Adv. Drug Deliv. Rev 96, 203–213. doi:10.1016/j.addr.2015.09.011

PubMed Abstract | CrossRef Full Text | Google Scholar

MED64 Application Note (2015). Cellartis® cardiomyocytes. Available at: http://www.med64.com/support/Application_Note_Cellartis_hPSC-CM.pdf.

CrossRef Full Text | Google Scholar

MED64 Application Note (2016). Human stem cell-derived cardiomyocytes. Available at: http://www.med64.com/wp-content/uploads/2016/04/Application_Note_hSC-CMs.pdf.

CrossRef Full Text | Google Scholar

Mirams, G. R., Cui, Y., Sher, A., Fink, M., Cooper, J., Heath, B. M., et al. (2011). Simulation of multiple ion channel block provides improved early prediction of compounds’ clinical torsadogenic risk. Cardiovasc. Res 91 (1), 53–61. doi:10.1093/cvr/cvr044

PubMed Abstract | CrossRef Full Text | Google Scholar

Mohammad, S., Zhou, Z., Gong, Q., and January, C. T. (1997). Blockage of the HERG human cardiac K+ channel by the gastrointestinal prokinetic agent cisapride. Am. J. Physiol 273 (5), H2534–H2538. doi:10.1152/ajpheart.1997.273.5.H2534

PubMed Abstract | CrossRef Full Text | Google Scholar

Mora, M. T., Ferrero, J. M., Romero, L., and Trenor, B. (2017). Sensitivity analysis revealing the effect of modulating ionic mechanisms on calcium dynamics in simulated human heart failure. PloS One 12 (11), e0187739. doi:10.1371/journal.pone.0187739

PubMed Abstract | CrossRef Full Text | Google Scholar

Nelder, J. A., and Mead, R. (1965). A simplex method for function minimization. Comput. J 7 (4), 308–313. doi:10.1093/comjnl/7.4.308

CrossRef Full Text | Google Scholar

Otte, S., Berg, S., Luther, S., and Ulrich, P. (2016). Bifurcations, chaos, and sensitivity to parameter variations in the sato cardiac cell model. Commun. Nonlinear Sci. Numer. Simulat 37, 265–281. doi:10.1016/j.cnsns.2016.01.014

CrossRef Full Text | Google Scholar

Qu, Y., and Vargas, H. M. (2015). Proarrhythmia risk assessment in human induced pluripotent stem cell-derived cardiomyocytes using the maestro mea platform. Toxicol. Sci 147 (1), 286–295. doi:10.1093/toxsci/kfv128

PubMed Abstract | CrossRef Full Text | Google Scholar

Raphel, F., Boulakia, M., Zemzemi, N., Coudière, Y., Guillon, J. M., Zitoun, P., et al. (2017). Identification of ion currents components generating field potential recorded in mea from hipsc-cm. IEEE Trans. Biomed. Eng 65 (6), 1311–1319. doi:10.1109/TBME.2017.2748798

PubMed Abstract | CrossRef Full Text | Google Scholar

Raphel, F., de Korte, T., Lombardi, D., Braam, S., and Gerbeau, J.-F. (2020). A greedy classifier optimisation strategy to assess ion channel blocking activity and pro-arrhythmia in hipsc-cardiomyocytes. PLoS Comput. Biol 16, e1008203. doi:10.1371/journal.pcbi.1008203

PubMed Abstract | CrossRef Full Text | Google Scholar

Schroll, H. J., Lines, G. T., and Tveito, A. (2007). On the accuracy of operator splitting for the monodomain model of electrophysiology. Int. J. Comput. Math 84 (6), 871–885. doi:10.1080/00207160701458724

CrossRef Full Text | Google Scholar

Sobie, E. A. (2009). Parameter sensitivity analysis in electrophysiological models using multivariable regression. Biophys. J 96 (4), 1264–1274. doi:10.1016/j.bpj.2008.10.056

PubMed Abstract | CrossRef Full Text | Google Scholar

Sundnes, J., Artebrant, R., Skavhaug, O., and Tveito, A. (2009). A second-order algorithm for solving dynamic cell membrane equations. IEEE Trans. Biomed. Eng 56 (10), 2546–2548. doi:10.1109/TBME.2009.2014739

PubMed Abstract | CrossRef Full Text | Google Scholar

Sundnes, J., Lines, G. T., and Tveito, A. (2005). An operator splitting method for solving the bidomain equations coupled to a volume conductor model for the torso. Math. Biosci 194 (2), 233–248. doi:10.1016/j.mbs.2005.01.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Sundnes, J., Lines, G. T., Cai, X., Nielsen, B. F., Mardal, K.-A., and Tveito, A. (2007). Computing the electrical activity in the heart. Cham, Switzerland: Springer Science and Business Media.

CrossRef Full Text | Google Scholar

Tixier, E., Raphel, F., Lombardi, D., and Gerbeau, J. F. (2017). Composite biomarkers derived from micro-electrode array measurements and computer simulations improve the classification of drug-induced channel block. Front. Physiol 8, 1096. doi:10.3389/fphys.2017.01096

PubMed Abstract | CrossRef Full Text | Google Scholar

Tung, L. (1978). A bi-domain model for describing ischemic myocardial dc potentials. PhD thesis. Cambridge, MA: Massachusetts Institute of Technology.

CrossRef Full Text | Google Scholar

Tveito, A., Jæger, K. H., Huebsch, N., Charrez, B., Edwards, A. G., Wall, S., et al. (2018). Inversion and computational maturation of drug response using human stem cell derived cardiomyocytes in microphysiological systems. Sci. Rep 8 (1), 17626. doi:10.1038/s41598-018-35858-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Tveito, A., Karoline, H., Kuchta, M., Mardal, K.-A., and Rognes, M. E. (2017). A cell-based framework for numerical modeling of electrical conduction in cardiac tissue. Front. Phys 5, 48. doi:10.3389/fphy.2017.00048

CrossRef Full Text | Google Scholar

Veeraraghavan, R., Gourdie, R. G., and Poelzing, S. (2014). Mechanisms of cardiac conduction: a history of revisions. Am. J. Physiol. Heart Circ. Physiol 306 (5), H619–H627. doi:10.1152/ajpheart.00760.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Veevers, J., Farah, E. N., Corselli, M., Witty, A. D., Palomares, K., Vidal, J. G., et al. (2018). Cell-surface marker signature for enrichment of ventricular cardiomyocytes derived from human embryonic stem cells. Stem Cell Rep 11 (3), 828–841. doi:10.1016/j.stemcr.2018.07.007

CrossRef Full Text | Google Scholar

Zhabyeyev, P., Missan, S., Jones, S. E., and McDonald, T. F. (2000). Low-affinity block of cardiac K(+) currents by nifedipine. Eur. J. Pharmacol 401 (2), 137–143. doi:10.1016/s0014-2999(00)00413-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, S., Zhou, Z., Gong, Q., Makielski, J. C., and January, C. T. (1999). Mechanism of block and identification of the verapamil binding domain to HERG potassium channels. Circ. Res 84 (9), 989–998. doi:10.1161/01.res.84.9.989

PubMed Abstract | CrossRef Full Text | Google Scholar

Zwi, L., Caspi, O., Arbel, G., Huber, I., Gepstein, A., Park, I. H., et al. (2009). Cardiomyocyte differentiation of human induced pluripotent stem cells. Circulation 120 (15), 1513. doi:10.1161/CIRCULATIONAHA.109.868885

PubMed Abstract | CrossRef Full Text | Google Scholar