Equations

Let Ω be the domain representing a MEA's well. The thickness of the layer of cells being small compared to the size of the well, the problem is assumed to be two-dimensional. We denote by A_m, C_m the surface area of membrane per unit volume of tissue, the membrane capacitance, and the thickness of the cell layer, respectively. The intra and extra-cellular conductivity tensors σ_i and σ_e are assumed to be scalar. The parameters values are reported in Table ​Table1.1. The propagation of the transmembrane potential V_m and the extracellular potential ϕ_e are modeled in Ω with the bidomain model (Tung, 1978):

Bidomain model parameters.

In the second equation, Ielk is the electric current which goes through the electrode located at e_k, |e_k| is the electrode surface and χ_ek is the characteristic function of e_k (which takes the value 1 on the electrode and 0 elsewhere). An imperfect model for the electrode is used to compute Ielk and described in the Supplementary Material. The activation is assumed to be triggered by a current I_app that is applied in an arbitrary region of the well with a cycle length of 1,200 ms. The locations of the stimulations are randomized to model the uncertainties of the spontaneous stimulus locations in in vitro measurements. This is further explained in the Heterogeneity modeling subsection. The computational domain Ω corresponds to one well of the MEA device as shown in Figure ​Figure1.1. Let n be the outward normal to the boundary of the domain Ω. Equations (1) are completed with the following boundary conditions: σ_i ∇ ϕ_i · n = 0 (where ϕ_i = V_m + ϕ_e), and either ϕ_e = 0 on the region connected to the ground or σ_e ∇ ϕ_e · n = 0 elsewhere. The ground location is indicated in Figure ​Figure11.

(A) Schematic of one well of the nine-electrode MEA device. The bidomain equations are solved in the domain Ω with homogeneous Neumann boundary conditions on Γ_N: ∇ϕe·n→=0 and homogeneous Dirichlet boundary conditions on Γ_D: ϕ_e = 0 where the ground is located. (B) Corresponding finite element mesh.

The transmembrane ionic current I_ion is described with the Minimal Ventricular (MV) model (Bueno-Orovio et al., 2008) which includes three currents: fast inward (fi), slow inward (si) and outward (so) currents. The reader is referred to the original publication for more details. Schematically, I_ion depends on V_m and on gating variables w = (w_j)1≤j≤3, solution of a system of three non-linear ordinary differential equations. A conductance coefficient g_s, with s = fi, si or so, controls the activity of the idealized channels associated with each of three currents of the model.

The partial differential equations are discretized in space by means of P1 finite elements, and in time by using backward differentiation formula (BDF) schemes with adaptive time steps and order provided by Sundials' CVODE library (Hindmarsh et al., 2005). The quantity of interest is the extra-cellular potential, also referred to as field potential (FP). It is a function of time and recorded at the electrodes locations.

In the present work, the computational model is used to generate synthetic MEA measurements. The main idea is to enrich the experimental data set with in silico measurements to make the classification more robust, in particular by exploring regions of the parametric space that are not covered by the experience. For a given set of conductances, the model is evaluated and the electrodes FPs are recorded. The conductances are chosen as to represent meaningful scenarios, as explained later in the Results section. To mimic experimental measurements, a zero-mean Gaussian noise of standard deviation 10 μV is added to the FPs (see Figure ​Figure2).2). A heterogeneity model of some ionic parameters is also considered to replicate the variability exhibited by the experimental measurements. This model is described later in this section. The stimulation location is also varied to model the uncertainty of the spontaneous stimulus location in the experiments. Figure ​Figure33 shows examples of synthetic recordings generated using the aforementioned in silico model. The FPs are simulated for three different scenarios. The scenarios consist in simulating the effects of sodium, calcium and potassium antagonist drugs, in each case with five different concentrations. In Supplementary Figure 1 a simulated FP recorded on an electrode is shown with the simulated action potential recorded on the same electrode.

Experimental recording of MEA field potential. Eight biomarkers are extracted from the time series: DA, depolarization amplitude; DW, depolarization width; RA, repolarization amplitude; FPD, field potential duration; AUCr, area under repolarization curve; RC, repolarizarion center; RW, repolarization width; FPN, field potential notch.

Comparison between in vitro and in silico MEA FP recordings. In each case, the FPs are recorded in control conditions and for five different drug concentrations. For the in silico measurements, the drugs effects are modeled using Equation (2), which amounts to reducing g_fi, g_si or g_so depending on the ion channel affected by the drug. (A) Effect of flecainide (sodium antagonist drug) on experimental recordings (left) and effect of a virtual sodium antagonist drug on simulated MEA FPs (right). (B) Effect of diltiazem (assumed to be mainly calcium antagonist in this study) on experimental recordings (left) and effect of a virtual calcium antagonist drug on simulated MEA FPs (right). (C) Effect of moxifloxacin (potassium antagonist drug) on experimental recordings (left) and effect of a virtual potassium antagonist drug on simulated MEA FPs (right).

Because the initial conditions of the ionic model do not correspond to those of a steady-state regime, several beats may need to be simulated before reaching a regime where there is negligible beat-to-beat variations. A numerical experiment was carried out to determine when this regime is reached. Figure ​Figure44 shows superimposed consecutive simulated FPs and the normalized beat-to-beat variations in the FP. When considering noisy synthetic measurements as described above, the steady-state is assumed to be reached when the beat-to-beat variations are comparable to variations induced by noise only. The beat-to-beat variability observed after this beat may be imputed to the coarseness of the mesh, the time discretization errors and the fluctuations of the ionic model itself. In the present work, the steady-state is assumed to be reached at the second beat. Therefore, the simulations are run for two cardiac cycles and the second beat is recorded to be used as a synthetic measurement.

Steady-state analysis: the Bidomain equations are solved for 100 consecutive beats. Qualitatively, a satisfactory steady state is reached at the second beat (left). The beat-to-beat relative difference of the FP is monitored (right) and is to be compared to the relative difference between two identical solutions, each polluted by an independent noise (right).

We chose to model the action of drugs on the ion channels by the conductance-block formulation of the pore block model (Bottino et al., 2006; Mirams et al., 2011; Zemzemi et al., 2013). This simple approach, which relies on a small number of parameters, was shown in Abbate et al. (2017) to be able to reproduce the expected effects of several drugs on MEA signals. The conductance of a given channel s is given by:

where g_control,s is the drug-free maximal conductance, [D] is the drug concentration, IC₅₀ is the value of the drug concentration at which the peak current is reduced of 50%, n is the Hill coefficient. In this work, n will be assumed to be equal to 1.

A typical experimental MEA FP measurement exhibits both a depolarization spike and a repolarization wave (see Figure ​Figure2).2). Using the computational model described above, the repolarization wave is usually too small compared to what is observed in experiments. As noted in Abbate et al. (2017), the repolarization wave provided by this model is larger when the domain includes cells with different APDs. In Abbate et al. (2017), the cell heterogeneity was defined on a checkerboard arbitrarily chosen in the MEA's well. We propose here a different approach, based on a probabilistic description of the heterogeneity. The tissue is supposed to be a continuous mixture of two cell types: A and B. We make the assumption that the transition between these two types can be described by a single space dependent parameter c(x, y) as follows:

where c is a random process with values in [0, 1] and p^(A) (resp. p^(B)) the set of 19 parameters of the MV model corresponding to cell type A (resp. B). The values of p^(A) and p^(B) are given in Supplementary Table 1. The APs corresponding to different realizations of c are shown in Figure ​Figure5.5. We make the hypothesis that the spatial variations of c are structured by a normal correlation function f_c:

Heterogeneity modeling: different APs obtained by simulating the MV model with different values of the heterogeneity parameter c. The heterogeneity parameter is a function of space and its pattern differs from one well to another (see Figure ​Figure66).

where l_c is the correlation length, set to l_c = 0.25 mm in the present work. To discretize the random process c, we compute the correlation matrix on the finite element mesh used for the discretization of the bidomain equations. The correlation matrix C = [C_i,j] ∈ ℝ^{N_mesh×N_mesh} reads:

where N_mesh is the total number of mesh nodes and (x^i,ŷi) are the coordinates of the ith node. The eigenpairs of C are denoted by (λ_i, Φ_i), and ordered by decreasing order of the eigenvalues λ_i. By a convenient abuse of notation, we denote by (x^,ŷ)→Φi(x^,ŷ) the function of the finite element space associated with the eigenmode Φ_i. Finally, the discretized heterogeneity field is approximated by the following truncated expansion:

where ξ = (ξ_i)i = 1…n_c is a random vector and n_c a truncation index chosen so that the truncation explains at least 99% of the variance. In other words, n_c is the smallest index n such that the following criterion is verified:

In our case, the choice of l_c and the domain geometry yields n_c = 14. Heterogeneity fields can now be generated simply by sampling the random variable ξ. In the present work, N_h = 128 heterogeneity fields were generated by sampling ξ from an uncorrelated uniform distribution over [-1,1]^n_c, and each sample is rescaled to range between 0 and 1. An example of heterogeneity field is presented in Figure ​Figure66.

One sample of cell heterogeneity field c(x, y) generated using the correlation matrix method. As c ranges from 0 to 1, the cell action potential varies from that of cell type “A” to cell type “B” (see Figure ​Figure55).

The observed variations in the experimental MEA FP recordings are also attributable to fluctuations in the stimulation location. In practice, the hiPSC-CM are not electrically stimulated: a stimulus arises spontaneously in the medium, probably due to the presence of pacemaker cells. The location of the spontaneous stimulation is not known to the experimentalist. We make the assumption that the location is random and therefore model it with a random uniform law over the square [0.15, 0.85]² where Ω = [0, 1]² is the complete domain.

To conclude, in a given experimental setting, we know neither the stimulation position nor the cell distribution inside the well and we would like the classification method to be robust with respect to all these unknown, random elements. This is why, when generating synthetic MEA FPs using our in silico model, we introduce two sources of uncertainty: the heterogeneous CM field and the stimulation location.