The simplified Kirchhoff network model (SKNM): a cell-based reaction–diffusion model of excitable tissue

Cell-based models of excitable tissues offer the advantage of cell-level precision, which cannot be achieved using traditional homogenized electrophysiological models. However, this enhanced accuracy comes at the cost of increased computational demands, necessitating the development of efficient cell-based models. The widely-accepted bidomain model serves as the standard in computational cardiac electrophysiology, and under certain anisotropy ratio conditions, it is well known that it can be reduced to the simpler monodomain model. Recently, the Kirchhoff Network Model (KNM) was developed as a cell-based counterpart to the bidomain model. In this paper, we aim to demonstrate that KNM can be simplified using the same steps employed to derive the monodomain model from the bidomain model. We present the cell-based Simplified Kirchhoff Network Model (SKNM), which produces results closely aligned with those of KNM while requiring significantly less computational resources.


S1 Derivation of the KNM parameters
In this section, we give an explanation for the definition of the KNM conductance parameters, G j,k e and G j,k i .
S1. 1 The extracellular conductance, G j,k e The parameter G j,k e represents the conductance of the electrical currents between the centers of the extracellular compartments surrounding cells j and k.In general, the conductance through a medium of length l (in cm), with conductivity σ (in mS/cm) and with cross-sectional area A (in cm 2 ) is given by [1] In the special case of the current between the centers of two extracellular compartments, j and k, the distance, l is given by l j,k and the conductivity, σ, is given by the extracellular conductivity, σ e .Furthermore, the cross-sectional area, A, of the extracellular compartments can be approximated by δ j,k e A j,k , where δ j,k e is the mean extracellular volume fraction of cell compartments j and k, and A j,k is the mean total cross-sectional area between the centers of compartments j and k including both the intracellular and extracellular parts of the compartments.Summing up, this gives G j,k e = δ j,k e A j,k σ e l j,k . (S2) S1. 2 The conductance between cells, G j,k i In order to model the conductance for the currents between connected cells, it is convenient to first consider the reciprocal of the conductance, i.e., the resistance.We assume that this resistance is given as the sum of the purely intracellular resistance (R j,k c ) and the resistance across the gap junctions (R j,k g ): Here, the gap junction resistance, R j,k g is defined as where G j,k g (in mS) is the conductance through the gap junctions connecting cells j and k.Furthermore, the purely intracellular resistance, R j,k c , is defined as 1/G j,k c , where G j,k c is given by (S1) for the purely intracellular space between the centers of cells j and k.In this case, the distance, l is given by l j,k , the conductivity, σ, is given by the intracellular conductivity, σ i , and the cross-sectional area, A, can be approximated by δ j,k i A j,k , where δ j,k i is the mean intracellular volume fraction of cell compartments j and k, and A j,k is, as above, the mean total cross-sectional area between the centers of compartments j and k, including both the intracellular and extracellular parts.Consequently, we get Inserting (S4) and (S5) into (S3), the conductance between cells j and k is defined as the reciprocal of R j,k i , i.e., by

S2 Convergence of solutions
In order to investigate the numerical resolution required for the simulations, we perform simulations using a number of different resolutions and compare the conduction velocity (CV) and maximal upstroke velocity ( dv dt max ) for each resolution, h (∆t or ∆x), to those computed for the finest considered resolution, h .More specifically, we compute the error where CV h and dv dt max h are the CV and maximal upstroke velocity, respectively, computed for the resolution h.The computation of the CV is as explained in the paper and the computation of the maximal upstroke velocity is as described below.

S2.1 Computation of the maximal upstroke velocity
In order to compute the maximal upstroke velocity, we record the membrane potential for each time step in the center cell for KNM and SKNM, and in a point in the center of the domain for BD and MD.From this recorded solution, the maximal upstroke velocity is computed as where v n is the recorded membrane potential solution in time step n, and ∆t is the size of the time step.

S2.2 Convergence of KNM and SKNM solutions for hiPSC-CMs
In Table S1, we investigate how the CV and maximal upstroke velocity computed using KNM and SKNM changes as the time step, ∆t, used in the simulation is decreased.We report the CV and maximal upstroke velocity computed for a number of different values of ∆t as well as the error (in percent) defined by (S7) with h = 0.001 ms.We observe that a time step of ∆t = 0.02 ms appears to be sufficient to achieve an error below 2%.In the KNM and SKNM simulations of hiPSC-CMs we therefore use a time step of ∆t = 0.02 ms.

S2.3 Convergence of BD and MD solutions for hiPSC-CMs
In Table S2 the CV and maximal upstroke velocity computed using BD and MD with a number of different spatial discretization parameters (∆x) are reported and compared to those computed for ∆x = 5 µm.We observe that a spatial resolution of 10 µm is sufficient for the error to be less than 2%.Furthermore, in Table S3, we report similar results for a number of different time steps, ∆t.Here, we observe that a time step of ∆t = 0.01 ms is sufficient for the error to be below 2%.We therefore use a time step of ∆t = 0.01 ms and spatial discretization of ∆x = 10 µm for the BD and MD simulations.

S2.4 Convergence of KNM and SKNM for pancreatic β cells
In Table S4, we report the CV, maximal upstroke velocity and the error of SKNM and KNM simulations for pancreatic β cells for different values of ∆t.Again, we observe that a time step of ∆t = 0.02 ms seems to be sufficient to ensure an error below 2%.We therefore use ∆t = 0.02 ms for the KNM and SKNM β cell simulations in this study.S3 SKNM based on negligible extracellular potential, SKNM(u e = 0) In the derivation of KNM (see [2]), Kirchhoff's current law is first applied for each cell k in the tissue.This results in the equation where and To derive KNM, we additionally need to apply Kirchhoff's current law for the extracellular compartment associated with each cell.However, in some cases it could be appropriate to instead assume that the extracellular potential, u k e , surrounding the cells is negligible (close to 0).In that case, (S11) gives v k = u k i , which can be inserted into (S10) and (S9) to yield where we have also included the equation (S13) for the additional state variables of the membrane model.We observe that the model (S12)-(S13) is the same as SKNM, except that the factor λ 1+λ is excluded from (S12) (i.e., the factor is set to 1).Note that we would also arrive at the same model if we assumed that G j,k e G j,k i .In that case, λ would be very large and the factor λ 1+λ in SKNM would be very close to 1.We here refer to the alternative model (S12)-(S13) resembling SKNM as SKNM(u e = 0).In the next subsections, we will investigate how well this further simplification of KNM approximates KNM in the two applications investigated in the main manuscript.

S3.1 SKNM(u e = 0) for a collection of β-cells
In Figure S1, we have computed the conduction velocity for several choices of the gap junction variation, γ, and the extracellular volume fraction for both KNM, SKNM and SKNM(u e = 0) in the case of a collection of β cells.We observe that KNM, SKNM and SKNM(u e = 0) all provide identical solutions in terms of conduction velocity for these examples.Figure S1: Comparison of the conduction velocity (CV) computed using KNM, SKNM and SKNM assuming a negligible extracellular potential, SKNM(u e = 0), for a collection of pancreatic β cells with an increasing gap junction conductance variation, γ.We consider the cases of 50%, 20%, 10%, and 2% extracellular space as indicated in the panel titles.

S3.2 SKNM(u e = 0) for a collection of hiPSC-CMs
In Figure S1, we observed that both SKNM and the version of SKNM where the extracellular potential is assumed to be negligible, i.e., SKNM(u e = 0), provided very good approximations of KNM for a collection of pancreatic β cells.In Figure S2, however, we observe that the assumption of negligible extracellular potential does not provide an equally good approximation for a collection of hiPSC-CMs.For a large extracellular volume percentage (50%), the conduction velocities computed for SKNM(u e = 0) seem quite close to those computed for KNM and SKNM, but as the extracellular volume percentage is decreased, SKNM(u e = 0) gradually approximate KNM less accurately.We consider the cases of 50%, 20%, 10%, and 2% extracellular space as indicated in the panel titles.
S4 A special case where the discretized BD and MD are equal to KNM and SKNM In this section, we will show that in the special case where a collection of cells is organized in a structured grid and BD and MD are spatially discretized using a finite difference approximation with resolution equal to the cell size, then the spatially discretized BD and MD equal the KNM and SKNM, respectively.Furthermore, if the SKNM assumption holds as well, then the algebraic equations defining all four models coincide.

S4.1 Modeling setup
For simplicity, we consider the 2D versions of the BD and MD models, corresponding to a 2D collection of cells in KNM and SKNM.We assume that each cell has the maximal length of l x in the x-direction, l y in the y-direction and l z in the z-direction, and that the cells are organized in a structured grid like, e.g., illustrated in Figures 1-2 and Figures 10-11 in the paper.Note that we do not assume that a cell makes up the entire volume l x l y l z .This volume is also assumed to contain some extracellular space.Rather, l x , l y and l z are the maximal cell lengths in each spatial direction, corresponding to the distance between the points where the cell connections in each direction are located.

S4.2 Finite difference approximation of BD equal to KNM
Let us first recall the definition of the bidomain model: with and See the paper for definitions and units for each of the model parameters and variables.
We wish to spatially discretize this model with a resolution equal to the cell size.In other words, we use the spatial discretization parameters ∆x = l x and ∆y = l y in the x-and y-directions, respectively.In addition, we let the grid points correspond to the center of the cells.We let p denote the cell number in the x-direction, which equals the grid point number in the xdirection and q denote the cell (and grid point) number in the y-direction.A finite difference discretization of (S14)-(S16) for points not on the boundary then reads x,p−1/2,q e (u p,q e −u p−1,q e ) l 2 x + M y,p,q+1/2 e (u p,q+1 e −u p,q e )−M y,p,q−1/2 e (u p,q e −u p,q−1 e ) l 2 y − I ion (v q,p , s q,p ) 0 = M x,p+1/2,q i (v p+1,q −v p,q )−M x,p−1/2,q i (v p,q −v p−1,q ) l 2 x + M y,p,q+1/2 i (v p,q+1 −v p,q )−M y,p,q−1/2 i (v p,q −v p,q−1 ) l 2 y (S21) x,p+1/2,q i +M x,p+1/2,q e )(u p+1,q e −u p,q e )−(M x,p−1/2,q i +M x,p−1/2,q e )(u p,q e −u p−1,q e ) l 2 x + (M y,p,q+1/2 i +M y,p,q+1/2 e )(u p,q+1 e −u p,q e )−(M y,p,q−1/2 i +M y,p,q−1/2 e )(u p,q e −u p,q−1 e ) l 2 y ds p,q dt = F (s p,q , v p,q ), (S22) where v p,q and u p,q e are the finite difference approximations of v and u e , respectively, in the grid point corresponding to the center of cell (p, q).Moreover, as an example, M x,p+1/2,q i is the value of M x i between the centers of cell (p, q) and cell (p + 1, q).

S4.2.1 Rewriting the discretization in terms of general neighbors
In order to compare this scheme to KNM, we now let the cell (p, q) be denoted as cell number k.Furthermore, if we, for instance, let cell (p+1, q) be denoted by j, then according to (S18), M x,p+1/2,q i is given by Following similar arguments for each instance of M i and M e appearing in (S20)-(S21) and defining the discrete system (S20)-(S22) can be written as where N k is the collection of all the neighbors of cell k.For the notation (p, q) for cell k, the neighbors are cells (p + 1, q), (p − 1, q), (p, q + 1) and (p, q − 1).We assume that the geometry of all cells is the same.Therefore, the membrane surface to volume ratio is given by where A m is the membrane area of a cell and l x l y l z is the total volume (both intracellular and extracellular) of each compartment.We also note that for any cell neighbours, j and k, the volume l x l y l z = A j,k l j,k , where A j,k is the total cross-sectional area between the centers of cell j and k including both the intracellular and extracellular spaces, and l j,k is the distance between the cell centers.Inserting this in (S26), we get Using (S24) and (S25), we observe that A j,k M j,k e l j,k = A j,k l j,k δ j,k e σ e = G j,k e , (S31) which equals the KNM parameters G j,k e and G j,k i defined in (S2) and (S6).Now, multiplying both sides of (S27) by l x l y l z and inserting (S30) and (S31) into (S26)-(S28), we get the following discrete version of BD, which is exactly the same as KNM with a constant value of A m (see the paper).

S4.2.2 Boundary conditions
So far, we have only considered the discretized BD in the grid points that are not located on the boundary of the domain.However, we may also discretize the boundary conditions of BD such that the resulting system of equations equals KNM also at the boundary.We consider homogenous Neumann boundary conditions for both v and u e on the entire boundary except for a homogeneous Dirichlet boundary condition for u e in a corner.
To make sure that the discrete version of BD equals KNM, the boundary conditions for u e should be discretized in the same manner as what we use for the boundary conditions for u e in KNM, for example a first order finite difference approximation like we have done in the computations in the paper.We should also use a first order difference for the boundary conditions for v.For example, for the grid point in the upper right corner, the homogenoeus Neumann boundary conditions state that the normal derivative of v should be zero in the positive x-and y-directions.This can be approximated by the

Figure S2 :
FigureS2: Comparison of the conduction velocity (CV) computed using KNM, SKNM, and SKNM with the assumption of a negligible extracellular potential, SKNM(u e = 0), for a collection of hiPSC-CMs with an increasing gap junction conductance variation, γ.We consider the cases of 50%, 20%, 10%, and 2% extracellular space as indicated in the panel titles.

Table S1 :
Conduction velocity (CV) and maximal upstroke velocity, dv dt max , computed using different values of ∆t in KNM and SKNM for a collection of hiPSC-CMs.The errors are computed by comparing the conduction velocity and maximal upstroke velocity to those computed using KNM or SKNM with h = 0.001 ms (see (S7)).

Table S2 :
Conduction velocity (CV) and maximal upstroke velocity, dv dt max , computed using different values of the spatial discretization parameter ∆x in BD and MD for a collection of hiPSC-CMs.The errors are computed by comparing the conduction velocity and maximal upstroke velocity to those computed using BD or MD with h = 5 µm (see (S7)).All the simulations in this table are performed using a time step of ∆t = 0.001 ms.

Table S3 :
Conduction velocity (CV) and maximal upstroke velocity, dv dt max , computed using different values of ∆t in BD and MD for a collection of hiPSC-CMs.The errors are computed by comparing the conduction velocity and maximal upstroke velocity to those computed using BD or MD with h = 0.001 ms (see (S7)).All the simulations in this table are performed using a spatial discretization of ∆x = 5 µm.

Table S4 :
Conduction velocity (CV) and maximal upstroke velocity, dv dt max , computed using different values of ∆t in KNM and SKNM for a collection of β cells.The errors are computed by comparing the conduction velocity and maximal upstroke velocity to those computed using KNM or SKNM with h = 0.005 ms (see (S7)).