# Non-ohmic tissue conduction in cardiac electrophysiology: Upscaling the non-linear voltage-dependent conductance of gap junctions

Daniel E. Hurtado, Javiera Jilberto, Grigory Panasenko

► To cite this version:

Daniel E. Hurtado, Javiera Jilberto, Grigory Panasenko. Non-ohmic tissue conduction in cardiac electrophysiology: Upscaling the non-linear voltage-dependent conductance of gap junctions. PLoS Computational Biology, 2020, 16 (2), pp.e1007232. 10.1371/journal.pcbi.1007232. hal-04135130

HAL Id: hal-04135130

https://hal.science/hal-04135130v1

Submitted on 24 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

RESEARCH ARTICLE

# Non-ohmic tissue conduction in cardiac electrophysiology: Upscaling the non-linear voltage-dependent conductance of gap junctions

Daniel E. Hurtado<sup>1,2,3*</sup>, Javiera Jilberto<sup>1,3</sup>, Grigory Panasenko<sup>4,5,6</sup>

1 Institute for Biological and Medical Engineering, Schools of Engineering, Medicine and Biological Sciences, Pontificia Universidad Católica de Chile, Santiago, Chile, 2 Department of Structural and Geotechnical Engineering, School of Engineering, Pontificia Universidad Católica de Chile, Santiago, Chile, 3 Millennium Nucleus for Cardiovascular Magnetic Resonance, Chile, 4 Institute Camille Jordan, Université Jean Monnet, Université de Lyon, Saint-Etienne, France, 5 Institute of Applied Mathematics, Vilnius University, Vilnius, Lithuania, 6 National Research University "Moscow Power Engineering Institute", Moscow, Russia

* dhurtado@ing.puc.cl

<figure><img src="image_3.png" /></figure>

OPEN ACCESS

**Citation:** Hurtado DE, Jilberto J, Panasenko G (2020) Non-ohmic tissue conduction in cardiac electrophysiology: Upscaling the non-linear voltage-dependent conductance of gap junctions. PLoS Comput Biol 16(2): e1007232. https://doi.org/10.1371/journal.pcbi.1007232

**Editor:** Aslak Tveito, Simula Research Laboratory, NORWAY

**Received:** June 28, 2019

**Accepted:** January 15, 2020

**Published:** February 25, 2020

**Peer Review History:** PLOS recognizes the benefits of transparency in the peer review process; therefore, we enable the publication of all of the content of peer review and author responses alongside final, published articles. The editorial history of this article is available here: https://doi.org/10.1371/journal.pcbi.1007232

**Copyright:** © 2020 Hurtado et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All codes that generated the data included in the paper are freely available at the GitHub repository https://github.com/dehurtado/NonOhmicConduction.

## Abstract

Gap junctions are key mediators of intercellular communication in cardiac tissue, and their function is vital to sustaining normal cardiac electrical activity. Conduction through gap junctions strongly depends on the hemichannel arrangement and transjunctional voltage, rendering the intercellular conductance highly non-Ohmic, particularly under steady-state regimes of conduction. Despite this marked non-linear behavior, current tissue-level models of cardiac conduction are rooted in the assumption that gap-junctions conductance is constant (Ohmic), which results in inaccurate predictions of electrical propagation, particularly in the low junctional-coupling regime observed under pathological conditions. In this work, we present a novel non-Ohmic homogenization model (NOHM) of cardiac conduction that is suitable to tissue-scale simulations. Using non-linear homogenization theory, we develop a conductivity model that seamlessly upscales the voltage-dependent conductance of gap junctions, without the need of explicitly modeling gap junctions. The NOHM model allows for the simulation of electrical propagation in tissue-level cardiac domains that accurately resemble that of cell-based microscopic models for a wide range of junctional coupling scenarios, recovering key conduction features at a fraction of the computational complexity. A unique feature of the NOHM model is the possibility of upscaling the response of non-symmetric gap-junction conductance distributions, which result in conduction velocities that strongly depend on the direction of propagation, thus allowing to model the normal and retrograde conduction observed in certain regions of the heart. We envision that the NOHM model will enable organ-level simulations that are informed by sub- and inter-cellular mechanisms, delivering an accurate and predictive in-silico tool for understanding the heart function. Codes are available for download at https://github.com/dehurtado/NonOhmicConduction.
**Funding:** DEH acknowledges the support from CONICYT thorough grant FONDECYT Regular 1180832. DEH and JJ received funding from Millennium Science Initiative of the Ministry of Economy, Development and Tourism of Chile, grant Nucleus for Cardiovascular Magnetic Resonance. GP received funding from Russian Science Foundation Grant #19-11-00033. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing interests:** The authors have declared that no competing interests exist.

## Author summary

The heart relies on the propagation of electrical impulses that are mediated gap junctions, whose conduction properties vary depending on the transjunctional voltage. Despite this non-linear feature, current mathematical models assume that cardiac tissue behaves like an Ohmic (linear) material, thus delivering inaccurate results when simulated in a computer. Here we present a novel mathematical multiscale model that explicitly includes the non-Ohmic response of gap junctions in its predictions. Our results show that the proposed model recovers important conduction features modulated by gap junctions at a fraction of the computational complexity. This contribution represents an important step towards constructing computer models of a whole heart that can predict organ-level behavior in reasonable computing times.

## Introduction

The conduction of electrical waves in cardiac tissue is key to human life, as the synchronized contraction of the cardiac muscle is controlled by electrical impulses that travel in a coordinated manner throughout the heart chambers. Under pathological conditions cardiac conduction can be severely reduced, potentially leading to reentrant arrhythmias and ultimately death if normal propagation is not restored properly [1]. At a subcellular level, electrical communication in cardiac tissue occurs by means of a rapid flow of ions moving through the cytoplasm of cardiac cells, and a slower intercellular flow mediated by gap junctions embedded in the intercalated discs. Gap junctions are intercellular channels composed by hemichannels of specialized proteins, known as connexins, that control the passage of ions between neighboring cells [2]. The regulation of ionic flow through gap junctions has been established for a variety of connexin types and hexameric arrangements, which under dynamic conditions result in a markedly non-linear relation between the electric conductance and the transjunctional voltage [3], revealing a non-ohmic electrical behavior. Further, it has been shown that ionic flow through cell junctions can take up to 50% of the total conduction time in cultured strands of myocytes with normal coupling levels [4], and that conduction velocity (CV) is largely controlled by the level of gap-junctional communication [1, 5], which highlights the key physiological relevance of gap-junction conductivity and coupling in tissue electrical conduction.

Cardiac modeling and simulation has strongly motivated the development of tissue-level mathematical models of electrophysiology, as they have the ability to connect subcellular mechanisms to whole-organ behavior [6]. To date, the vast majority of continuum models assume a linear conduction model of spatial communication, based on the assumption that electrical current in cardiac tissue follows Ohm's law, i.e, that current is linearly proportional to gradients in the intra-cellular potential [7, 8]. From a mathematical perspective, the assumption that conduction in cardiac tissue follows Ohm's law is conveniently represented by a linear diffusion term, where gradients are modulated by a conductivity tensor that is independent of the local electrical activity. The most complete electrophysiology formulation is given by the bidomain model [9] where both the intra-cellular and extra-cellular potential fields are considered. Further, by assuming that the intra- and extra-cellular conductivity tensors have the same anisotropy ratio, the bidomain equations can be conveniently represented by a non-linear reaction-diffusion partial differential equation known as the monodomain (cable) model [10].
Using two-scale asymptotic homogenization techniques, analytic expressions have been obtained for the effective conductivity tensor, which is then used to model the electrical current in an average macroscopic sense [11–13]. To this end, periodicity at the microstructural level of cardiac tissue is assumed, and a representative tissue unit is partitioned in regions of high and low conductivity that represent the cytoplasm and intercalated discs with gap junctions, respectively. While this approach allows for the explicit consideration of regions with decreased conductivity, e.g. membranes where flow is mediated by gap junctions, Ohm’s law is still assumed to hold throughout the microstructural domain [14]. As a result, the non-Ohmic behavior of gap junctions and their impact on tissue-level conduction continues to be neglected [15]. In particular, it has been shown that continuum models that consider effective conductivity tensors described above fail to capture the slow conduction of electrical impulses in cases of low gap-junctional coupling [13, 16], limiting their applicability to the simulation of pathological conditions in excitable tissue. Non-linear diffusion models that replace the conduction term in the monodomain equation either by a fractional laplacian [17, 18] or a porous-medium-like diffusive term [8, 19] have been recently proposed. While these formulations have shown to modulate the shape of propagating waves and other restitution properties, they remain largely phenomenological, and have the disadvantage of not being able to upscale microscopic physical information, neither have been assessed for cases of low junctional coupling.

In this work, we present a multiscale continuum model of cardiac tissue conduction that accounts for the nonlinear communication between adjacent cells. We argue that the explicit consideration of the non-ohmic behavior of gap junctions can be seamlessly embedded into continuum tissue-scale models of electrophysiology using an asymptotic homogenization approach, which delivers nonlinear continuum equations for characterizing the electrical conduction in excitable media.

# Methods

## Multiscale tissue model for non-ohmic conduction

In the following we consider the microscopic problem of non-linear conduction in a strand of cardiac cells with domain $\Omega = (0, L)$, see Fig 1. We let $\varepsilon$ be the cell length, $\delta\varepsilon$ be the length of gap junctions, and assume that $\delta\varepsilon \ll \varepsilon \ll L$. Further, we let $u_{\varepsilon,\delta}$ be the microscopic transmembrane potential field, and $j_{\varepsilon,\delta}$ be the microscopic current density. The time-independent problem of conduction resulting from current balance reads

$$ -\frac{\partial}{\partial x} j_{\varepsilon,\delta}(u_{\varepsilon,\delta}) = 0, \quad x \in \Omega, \qquad (1) $$

and we note that in writing (1) it has been assumed that the extra-cellular potential is constant along the cardiac strand. We denote the space occupied by the cytoplasm by

$B_{\varepsilon,\delta}^{cyt} = \bigcup_{k=-\infty}^{\infty} ((k + \frac{\delta}{2})\varepsilon, (k + 1 - \frac{\delta}{2})\varepsilon)$, and the space occupied by gap junctions by $B_{\varepsilon,\delta}^{gap} = \bigcup_{k=-\infty}^{\infty} ((k - \frac{\delta}{2})\varepsilon, (k + \frac{\delta}{2})\varepsilon)$. Further, we assume that current is governed by Ohm’s law inside the cytoplasm with conductivity $\sigma_c$, but is non-linearly regulated at the gap junctions, which we express by the following microscopic constitutive law

$$ j_{\varepsilon,\delta}(u_{\varepsilon,\delta}) = -\sigma(x, \{u_{\varepsilon,\delta}\}) \frac{\partial u_{\varepsilon,\delta}}{\partial x} \qquad (2) $$
<figure><img src="image_4.png" /><figcaption>Fig 1. Schematic of the multiscale model of cardiac conduction. Ionic currents are linearly proportional to gradients of transmembrane potential inside the cytoplasm, but are non-linearly mediated by gap junctions located at the intercalated discs.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1007232.g001

where the conductivity is described by the following relation

$$ \sigma(x, u) = \begin{cases} \sigma_c, & x \in B_{\varepsilon, \delta}^{cyt}, \\ \delta\sigma_g(1 + \mu a(S[u]_{j,\varepsilon})), & x \in B_{\varepsilon, \delta}^{gap}, \end{cases} \quad (3) $$

where $\delta\sigma_g$ is a representative conductivity for the intercalated disc with gap junctions, $\mu$ is a positive constant, and $a$ is a smooth bounded function that depends on the transjunctional voltage jump defined as

$$ [u]_{j,\varepsilon} = u\left(k + \frac{\delta}{2}\varepsilon\right) - u\left(k - \frac{\delta}{2}\right), \quad (4) $$

where $S$ is a scaling parameter ($S \le \varepsilon^{-1}$). From (3) and the relationship between electrical conductance and conductivity for a cylindrical domain we write

$$ \delta\sigma_g(1 + \mu a(S[u]_{j,\varepsilon})) = \beta \frac{g_{jo}\delta\varepsilon}{A_{cell}} g_j(S[u]_{j,\varepsilon}), \quad (5) $$

where $g_{jo}$ is a representative conductance for the intercalated disc with gap junctions, $A_{cell}$ the cross-sectional area of the cell, and the parameter $\beta$ represents the level of gap-junctional coupling (GJc), with $\beta \in [0, 1]$. From (5), we set $\delta\sigma_g = \frac{g_{jo}\delta\varepsilon}{A_{cell}}$. As a result, we get

$$ \mu a(S[u]_{j,\varepsilon}) = \beta g_j(S[u]_{j,\varepsilon}) - 1. \quad (6) $$

Using asymptotic analysis (see S1 Appendix for technical details and proofs) we show that the macroscopic (tissue-level) current conservation for the steady-state problem is governed
by the homogenized equation

$$ \frac{\partial}{\partial x}(\hat{\sigma}(\frac{\partial v}{dx})\frac{\partial v}{\partial x}) = 0, \quad x \in \Omega, \qquad (7) $$

where $v$ is the macroscopic transmembrane potential, and the effective conductivity modulating conduction at the macroscopic scale takes the form

$$ \hat{\sigma}(y) = \sigma_c \left\{ \frac{1 + \mu a(S\varepsilon[N](y))}{\frac{\sigma_c}{\sigma_g} + (1 - \delta)(1 + \mu a(S\varepsilon[N](y)))} \right\} \qquad (8) $$

where

$$ [N](y) = -(1 - \delta)\left(\frac{\hat{\sigma}(y)}{\sigma_c} - 1\right)y, \qquad (9) $$

with $[N] = [u]_{j,\varepsilon}/\varepsilon$, and we note that for a given transmembrane potential gradient $y$, the effective conductivity $\hat{\sigma}(y)$ is implicitly solved from (8) and (9). Further, we show that under reasonable assumptions, the following error estimate for the macroscopic transmembrane potential holds

$$ \|u_{\varepsilon,\delta} - v\|_{L^{\infty}((0,1))} = O(\varepsilon + \delta^2). \qquad (10) $$

We now focus on the time-dependent macroscopic model of cardiac electrophysiology for the time interval $(0, T)$. The homogenized electrical flux described in the right-hand side of (7) is then balanced by the transmembrane current, leading to the non-Ohmic cable equation

$$ \frac{\partial}{\partial x}(\hat{\sigma}(\frac{\partial v}{dx})\frac{\partial v}{\partial x}) = A_m \left\{ C_m \frac{\partial v}{\partial t} + I_{\text{ion}} \right\} \quad \text{in } \Omega \times (0, T), \qquad (11) $$

where $I_{\text{ion}}: \mathbb{R} \times \mathbb{R}^M \to \mathbb{R}$ represents the transmembrane ionic current, $C_m$ is the membrane capacity and $A_m$ is the surface-to-volume ratio, and we note that the right-hand side of (11) accounts for the amount of charge that leaves the intra-cellular domain and enters the extra-cellular domain. Further, we will assume that the transmembrane ionic current $I_{\text{ion}}$ is governed by $v$ and by gating variables $w: \Omega \times (0, T) \to \mathbb{R}^M$ that modulate the conductance of ion channels, pumps and exchangers, i.e., $I_{\text{ion}} = I_{\text{ion}}(v, w)$, where the exact functional form of $I_{\text{ion}}$ will depend on the choice of ionic model. The evolution of gating variables is determined by kinetic equations of the form

$$ \frac{\partial w}{\partial t} = g(v, w), \qquad (12) $$

where the form of $g: \mathbb{R} \times \mathbb{R}^M \to \mathbb{R}^M$ will also depend on chosen the ionic model. The Eqs (11) and (12) are supplemented with initial and boundary conditions for the transmembrane potential and gating variables to form an initial boundary value problem, which we refer to as the Non-Ohmic Homogenization Model (NOHM). Boundary conditions at the left end prescribed a pulsatile electrical current, while the right end was prescribed with zero current. To study reverse conduction, the boundary conditions where flipped. The numerical solution of the coupled system of the non-Ohmic cable Eq (11) and kinetic Eq (12) was performed using a standard Galerkin finite-element scheme [20] for the spatial discretization and a Forward Euler scheme for the time discretization implemented in FEiCS [21], see S1 Appendix. Codes are available for download at https://github.com/dehurtado/NonOhmicConduction.
For the sake of comparison, we also consider the case of uniform (Ohmic) gap junction conductivity, i.e.,

$$ \sigma(x, u) = \begin{cases} \sigma_c, & x \in B_{\varepsilon, \delta}^{cyt}, \\ \beta \delta \sigma_g, & x \in B_{\varepsilon, \delta}^{gap}. \end{cases} \qquad (13) $$

Following standard asymptotic-analysis arguments for linear systems [11], one can show that for the case of the piecewise uniform conductivity tensor defined in (13) the effective conductivity tensor takes the form

$$ \hat{\sigma} = \left\{ \frac{1}{\sigma_c} + \frac{1}{\beta \sigma_g} \right\}^{-1}. \qquad (14) $$

We remark that in this case the macroscopic conductivity is not dependent on the voltage gradient. Further, we note that (14) is also obtained as a particular case of the NOHM when the microscopic conductivity (13) is assumed, see S1 Appendix. We refer to the system of Eqs (11) and (12) that considers the uniform effective conductivity tensor (14) as the Linear Homogenization Model (LHM).

## Cellular models of cardiac propagation

To validate the proposed NOHM model we consider the cellular model described in [7, 22], in which a strand of cardiac cells electrically connected by gap junctions are represented using a circuit network, see Fig 2. In the following, we summarize the main aspects of cellular modeling. For the cellular models (CM), a chain of cells is discretized at subcellular level in $n_{div} = 10$ subregions. A generic node $i$ is connected to its neighbor $i + 1$ through a resistor $R^{i,i+1}$. If $i$ and $i + 1$ belongs to the same cell $R^{i,i+1} = R_{myo}$ where $R_{myo}$ is the myoplasmic resistance. Now, if $R^{i,i+1}$ connects nodes from different cells (i.e. a gap junction) its value will be given by $R^{i,i+1} = R_j(V_j)$ where $R_j(V_j)$ is the resistance of the gap junction, which depends non-linearly on the transjunctional voltage $V_j$. For computing the resistance $R^{i,i+1}$ we use

$$ V_j = v^{i+1+(ndiv-1)} - v^{i-(ndiv-1)} $$

where we consider how the experimental relationship $V_j$ versus $g_j$ is obtained through the dual voltage clamp method [23].

Due to the Kirchhoff's law, the current balance yields to the following finite-difference equation

$$ \sigma^{i,i+1} \frac{v^{i+1} - v^i}{\Delta x^2} + \sigma^{i-1,i} \frac{v^{i-1} - v^i}{\Delta x^2} = A_m \left\{ C_m \frac{\partial v^i}{\partial t} + I_{ion}(v^i, w^i) \right\}, \qquad (15) $$

<figure><img src="image_5.png" /><figcaption>Fig 2. Circuit representation for the cellular models.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1007232.g002
where $A_m$ is the surface to volume ratio and $C_m$ the membrane capacitance. The ionic current $I_{ion}$ depends on the transmembrane potential in the node $v^i$ and on the gating variables $w^i$. Local conductivity can be expressed in terms of the resistance as

$$ \sigma^{i,i+1} = \frac{\Delta x}{R^{i,i+1} A_{cell}} $$

with $\Delta x = l_{cell}/n_{div}$. Then

$$ \sigma^{i,i+1} = \begin{cases} \sigma_c & \text{if } i, i+1 \text{ are in the same cell} \\ \delta\sigma_g\beta g_j(V_j) & \text{if } i, i+1 \text{ are in different cells} \end{cases} \qquad (16) $$

and we call the model given by Eqs (15) and (16) CM voltage-gated. If the gap junctions are assumed to be voltage insensitive, then the conductivity reads

$$ \sigma^{i,i+1} = \begin{cases} \sigma_c & \text{if } i, i+1 \text{ are in the same cell,} \\ \beta\delta\sigma_g & \text{if } i, i+1 \text{ are in different cells.} \end{cases} \qquad (17) $$

We refer to the model given by Eqs (15) and (17) as the CM clamped, and we note that it serves as a cellular counterpart to the LHM.

## Conduction experiments in a cardiac strand

In this work, we model the propagation of electrical impulses in a strand with a length $L = 6.4$ mm. Cells were assumed to be cylinders with radius $r_{cell} = 11 \mu m$, length $\varepsilon = 100 \mu m$, and inter-calated-disc length ratio of $\delta = 10^{-4}$. Based on these dimensions, $A_{cell} = 380 \mu m^2$. For the gap-junction conductance model, we set $S = 2$, and assumed $g_{jo} = 2.534 \mu m$ [24], which results in a representative conductivity of $\delta\sigma_g = 6.67 \cdot 10^{-5} Sm^{-1}$. The cytoplasmic conductivity and the membrane capacitance are taken to be $\sigma_c = 0.667 Sm^{-1}$ and $C_m = 1 \mu F/cm^2$, respectively [24]. For the transmembrane ionic current, we considered the Luo-Rudy I model [25]. The surface to volume ratio is given by $A_m = 2RCG/r_{cell}$, where RCG = 2 is the ratio between capacitive and geometrical areas and [24, 26]. Given the time-dependent behavior of the conductance of GJs [5, 27], in our experiments, we consider two limit cases for the temporal state of gap junctions: the instantaneous conductance case, and the steady-state conductance case. For the instantaneous conductance case, we adopt the model of Vogel and Weingart [28] which for the normalized conductance takes the form

$$ g_{j,inst}(V_j) = \begin{cases} \frac{G_j^-}{e^{V_j^-/(1+e^{V_j^-/V_H^-})} + e^{V_j^-/(1+e^{-V_j^-/V_H^-})}} & V_j < 0, \\ \frac{G_j^+}{e^{V_j^+/(1+e^{V_j^+/V_H^+})} + e^{V_j^+/(1+e^{-V_j^+/V_H^+})}} & V_j > 0, \end{cases} \qquad (18) $$

where $V_j$ is the transjunctional voltage, and $G_j^+$, $G_j^-$, $V_H^+$, $V_H^-$ are parameters. For the steady-state case, we assume that the normalized gap-junction conductance follows a Boltzmann
<figure><img src="image_6.png" /><figcaption>Fig 3. Normalized conductance of gap junctions as a function of the transjunctional voltage. (left) Cx43-Cx43 channel, (center) Cx45-Cx45 channel, and (right) Cx43-Cx45 channel. Data extracted from [3]. The (○) and (●) data corresponds to instantaneous and steady state conductance, respectively.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1007232.g003

distribution [29] that reads

$$g_{j,ss}(V_j) = \begin{cases} \frac{1-g_{j,min}^{+}}{p+e^{(A^{+}(V_j-V_{j0}^{+})}} + g_{j,min}^{+}} & V_j < d, \\ \frac{1-g_{j,min}^{-}}{p+e^{(A^{-}(V_j-V_{j0}^{-})}} + g_{j,min}^{-}} & V_j > d, \end{cases} \quad (19)$$

where $p, d, g_{j,min}^{+}, g_{j,min}^{-}, A^{+}, A^{-}, V_{j0}^{+}, V_{j0}^{-}$ are parameters that depend on the type of channel.

To assess the performance of the NOHM model under different conductance distributions, we considered three types of gap-junction channels: the homomeric-homotypic channels Cx43-Cx43 and Cx45-Cx45, and the homomeric-heterotypic channel Cx43-Cx45. The normalized conductance distributions for the instantaneous and steady-state cases of these channels are depicted in Fig 3. The parameters for the instantaneous and steady-state conductance models have been reported in the literature [3], and are summarized in Table 1. The parameter A is computed from z as $A = z/kT$, where $kT = 25.7$ meV is the product of the Boltzmann constant $k$ with the temperature $T$. The effect of gap-junctional coupling on conduction is studied by modulating the maximal gap-junction conductance from $\beta = 100\%$ to $\beta = 0.5\%$. The cardiac strand is excited on one end with an applied transmembrane current using a pacing cycle length of 800 ms and whose amplitudes varied between 10 µA/mm² to 35 µA/mm², which elicits a propagating pulse from left to right. Retrograde-conduction effects are studied by propagating pulses from right to left. Pacing rate dependence was studied by constructing CV restitution curves, using the pacing protocol reported by Gizzi and co-workers [30].

<table><caption>Table 1. Parameters for the conductance distribution of gap junctions, taken from [3]. For $V_{j0}, g_{j,min}, z$ the negative/positive values are presented. The Cx43-Cx45 case considered a modified Boltzmann distribution to improve the fitness to data.</caption><thead><tr><th rowspan="2">Channel</th><th colspan="2">Instantaneous model</th><th colspan="4">Steady-state model</th></tr><tr><th>G<sub>j</sub></th><th>V<sub>H</sub> [mV]</th><th>V<sub>j0</sub> [mV]</th><th>g<sub>j,min</sub></th><th>z</th><th>p</th><th>d</th></tr></thead><tbody><tr><td>Cx43-Cx43</td><td>1.99/2.01</td><td>-175.8/318.4</td><td>-60.8/62.9</td><td>0.26/0.25</td><td>-3.4/2.9</td><td>1</td><td>0</td></tr><tr><td>Cx45-Cx45</td><td>1.99/2.02</td><td>-112.7/135.0</td><td>-38.9/38.5</td><td>0.16/0.17</td><td>-2.5/2.7</td><td>1</td><td>0</td></tr><tr><td>Cx43-Cx45</td><td>1.93/2.0</td><td>-130.0/404.0</td><td>-15.9/149.3</td><td>0.05/0.05</td><td>-2.1/0.7</td><td>0.73</td><td>25</td></tr></tbody></table>

https://doi.org/10.1371/journal.pcbi.1007232.t001
<figure><img src="image_7.png" /><figcaption>**Fig 4.** Impulse conduction features from computational simulations. The propagating wavefront predicted by the CM clamped, CM voltage-gated, LHM and NOHM are compared for three levels of transjunctional coupling: (top row) high coupling $β = 100\%$, (middle row) low coupling $β = 10\%$, and (bottom row) very low coupling $β = 1\%$. In general, the LHM and NOHM drift ahead of their CM counterparts as the GJc is decreased in the case of instantaneous conductance. In contrast, for the case of steady-state conductance the NOHM accurately predicts the CM voltage-gated even for very low coupling levels, whereas the LHM substantially drifts ahead from the CM clamped wavefront.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1007232.g004

## Results

The numerical solution of CM clamped, and CM voltage-gated resulted in dynamical systems with 642 degrees of freedom, using cell subdomains with a length of 10 µm. In contrast, the continuum LHM and the NOHM employed only 65 degrees of freedom, equivalent to a spatial discretization of 100 µm. Fig 4 shows the propagating wavefronts as predicted by the four conduction models studied in this work for the instantaneous conduction and steady-state
conduction cases of the Cx43-Cx43 channel. For high GJc, $\beta$ = 100%, we observe that all four models predicted a very similar wavefront both for the instantaneous and for the steady-state regimes of conductance (Fig 4, top row). When GJc was reduced to low coupling levels, $\beta$ = 10%, continuum models (LHM and NOHM) resulted in propagating waves that drifted ahead of their CM counterparts (CM clamped and CM voltage-gated) in the case of instantaneous conductance. Interestingly, for the case of steady-state conductance, the NOHM accurately predicted the response of the CM voltage-gated, whereas the LHM drifted ahead of the CM clamped (Fig 4, middle row). Remarkably, for very low GJc, $\beta$ = 1%, the NOHM still delivered an accurate wavefront when compared to the CM voltage-gated in the steady-state conductance case, whereas considerable drift occurred in the instantaneous conductance case both for the LHM and the NOHM (Fig 4, bottom row).

Fig 5 shows the CV as a function of the GJc, measured in terms of $\beta$, for the instantaneous and steady-state conduction cases for homomeric-homotypic channels Cx43-Cx43 and Cx45-Cx45. Under instantaneous conductance (Fig 5, left column) all four models converged to CV values between 64–65 cm/s for the high-coupling case $\beta$ = 100%. Cellular models delivered, in general, a very similar behavior. As GJc was reduced, both LHM and NOHM overestimated the CV when compared to CM clamped and CM voltage-gated, respectively. For the case of steady-state conductance (Fig 5, right column), the CM clamped consistently delivered higher values of CV than the CM voltage-gated. Further, the LHM continued to overestimate the CV when compared to the CM clamped as GJc was reduced. Remarkably, the NOHM closely followed the behavior of the CM voltage-clamped, even for very low GJc. Cellular models predicted conduction block for $\beta$ ≤ 0.5% in all cases. Conduction block was captured by the NOHM for the steady-state conductance regime, but not for the instantaneous case. LHM did not predict the conduction block for the range of $\beta$ analyzed. Retrograde wave propagation experiments for Cx43-Cx43 and Cx45-Cx45 channels (not shown in the figure) resulted in CV curves that did not differ from those obtained in normal wave propagation.

Fig 6 shows the dependence of CV on the GJc for the Cx43-Cx45 channel under normal wave propagation (left-to-right direction) and retrograde wave propagation (right-to-left direction). For the instantaneous conductance case (Fig 6, left column), the CM voltage-gated predicted CVs that were slightly higher in retrograde propagation than in normal propagation. Continuum models overestimated the CV when compared to their cellular counterparts. While the LHM predictions were insensitive to the direction of propagation, the NOHM predicted higher CVs for the case of retrograde wave propagation. For the steady-state conductance case (Fig 6, right column), the CM voltage-gated resulted in CVs that were considerably higher in retrograde propagation than in normal propagation, a feature not captured by the CM clamped, which was insensitive to the direction of propagation. The NOHM was able to capture such large differences in CV as well as the conduction block at $\beta$ = 10%, whereas the LHM only delivered reasonable predictions for the retrograde propagation, but considerably overestimated CVs in the normal propagation case, and did not reach conduction block for any of the $\beta$ values analyzed.

The dependence of CV on the pacing rate and propagation orientation in voltage-dependent models of conduction was studied by constructing CV restitution curves for the case of channel Cx43-Cx45, see Fig 7. In all cases, the CM voltage-gate under retrograde propagation resulted in higher CVs than in normal propagation. The NOHM captured this orientation dependence and delivered CV curves with a similar shape but with higher values for the cases of low and very low GJc. For the case of steady-state conductance, the NOHM also captured the conduction block predicted by the CM voltage-gate for low and very low GJc. The cycle length in all simulations was reduced until loss of propagation, which occurred in the range of 310-330 ms.
<figure><img src="image_8.png" /><figcaption>Fig 5. Conduction velocity studies on a cardiac strand and the effect of gap-junction coupling for symmetric conductance distributions: (top row) Cx43-Cx43 channel, (bottom row) Cx45-Cx45 channel, (left column) instantaneous conductance, (right column) steady-state conductance. Black and red colors are used to indicate voltage-independent and voltage-dependent gap-junction conduction, respectively. Voltage-dependent models delivered lower conduction velocities than voltage-independent models of gap-junction conductance, particularly for the steady-state regime.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1007232.g005

The influence of the spatial discretization on the CV is reported in Fig 8. Mesh sizes, interpreted as cell segments in cellular modes, ranging from $\Delta x = 0.01$ mm to $\Delta x = 0.2$ mm were considered for both the instantaneous and steady-state conduction cases of a Cx43-Cx43 channel under high GJc. In both cases, a strong dependence of the CV on the mesh size was observed for the CM clamped and CM voltage-clamped, with higher CVs for larger mesh
<figure><img src="image_9.png" /><figcaption>**Fig 6.** Conduction velocity studies on a cardiac strand and the effect of gap-junction coupling for the Cx43-Cx45 channel with non-symmetric conductance distribution: (left) instantaneous conductance case, (right) steady-state conductance case. Black color denotes voltage-independent models, red and blue colors denote voltage-dependent models. Predictions from gap-junction voltage-independent models CM clamped, and LHM were insensitive to the direction of wave propagation, whereas voltage-dependent models resulted in CVs that strongly depended on the direction of wave propagation for the steady-state conductance case.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1007232.g006

sizes. An attenuated yet considerable mesh dependence was also observed for the LHM and NOHM.

Discussion

In this article, we study the gap-junction-mediated electrical conduction in excitable cardiac tissue through a novel non-ohmic multiscale model. A unique feature of the proposed model is that tissue-level spatial conduction is fully informed by sub-cellular communication mechanisms, specifically by cytoplasmic and gap-junctional conductances. While the upscaling of conduction properties in excitable media has been the subject of some studies in the past using a linear homogenization theory approach [11, 12], our work offers a rigorous mathematical framework that delivers an effective non-linear model of conduction able to represent, at the tissue level, the non-Ohmic conduction that takes place at the sub-cellular level. Although our focus has been on understanding gap-mediated communication between cardiac myocytes, the present model of conduction can be extended to study the electrical propagation phenomenon in other areas of biology, such as the neurosciences, where electrical synapsis occurring in the brain is highly regulated by neural gap junctions [31].

To validate the predictions and understand the unique features of the NOHM, we considered cellular models with and without voltage dependence at the gap junctions (CM clamped and CM voltage-gated) as well as a linear continuum model of conduction (LHM). The
<figure><img src="image_10.png" /><figcaption>Fig 7. Conduction-velocity restitution curves for the homomeric-heterotypic channel Cx43-Cx45 for high, low and very low GJc: (left) instantaneous conductance case, (right) steady-state conductance case. CL = cycle length.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1007232.g007

behavior of all four models of conduction was studied in the propagation of waves in a cardiac
strand [24] with decreasing levels of GJc, see Fig 4. Features that arise in propagating action
potentials under decreasing levels of coupling such as a steeper upstroke and a notch in the
upstroke [5] are predicted both by the LHM and NOHM. Remarkably, this prediction is
achieved at a fraction of the computational complexity involved in cellular models, as the num-
ber of degrees of freedom in continuum models are one order of magnitude smaller. An alter-
native approach is the use of hybrid multiscale models [32], which adaptively partition the
domain to solve macroscopic cable equations in regions with low potential gradients and
impose microscopic equations of conduction in regions of high potential gradients. While
hybrid models can reduce the computational complexity of simulations, they involve a signifi-
cant increase in the number of degrees of freedom when compared to standard homogenized
models. We believe that the NOHM model offers the advantage of delivering accurate predic-
tions while maintaining the computational cost similar to that of standard macroscopic con-
tinuum models. The balance between predictive power and computational cost remains one of
the main hurdles in the development of patient-specific whole-heart simulations [33], which
highlights the importance of developing accurate yet efficient tissue-level models.
<figure><img src="image_11.png" /><figcaption>**Fig 8.** The effect of spatial discretization of the conduction velocity: (left) instantaneous conductance, (right) steady-state conductance. The conduction velocity in cellular models exhibit a stronger dependence on the mesh size than the continuum models of conduction.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1007232.g008

The accuracy of the NOHM and LHM in predicting their cellular counterparts CM clamped and CM voltage-gated was assessed by studying the CV for a wide range of gap-junctional coupling levels for channels with symmetric conductance distributions, see Fig 5. In our work, we considered two limits of the dynamic conductance of gap junctions: the instantaneous and steady-state conductance cases. In the case of instantaneous conductance, the CV predicted by the NOHM was in general higher than the CV predicted by the CM voltage-gated, which was observed to decrease as the GJc was reduced [27]. A similar trend was observed for the LHM when compared with the CM clamped. Further, the NOHM and LHM resulted in similar CV curves. Previous studies have confirmed that the accuracy of LHM in predicting cardiac conduction, as dictated by CM clamped, consistently deteriorates as the GJc is decreased to low levels [13, 16]. Interestingly, for the case of steady-state conduction, substantial differences arise between voltage-dependent (NOHM, CM voltage-gated) and voltage-independent (LHM, CM clamped) models, with the former resulting in considerably lower CVs, see Fig 5 (left column). These results can be explained by noting the shape of the conductance distributions associated to the instantaneous and steady-state cases. In particular, the instantaneous conductance of the Cx43-Cx43 channel displays a flat shape with small variations within the range of transjunctional voltages (Fig 3, left), which resembles an Ohmic electrical response. Thus, in this case, the NOHM is not expected to differ from the LHM, as the conduction in both cases is fairly Ohmic, a behavior observed in our experiments (Fig 5, top left). In contrast, the steady-state conductance distribution for the Cx45-Cx45 channels presents a narrow bell shape (Fig 3, center), which is representative of a marked non-Ohmic electrical behavior. Notably, simulations associated with that conductance distribution are the ones that deliver the most different CV curves when the voltage dependence is included (Fig 5,
bottom right). These results confirm the ability of the NOHM to accurately upscale the volt-
age-dependent behavior of gap junctions, a feature not offered by standard homogenization
models of conduction, such as the LHM. Further, we note here that the difference in CV
between cellular models CM clamped and CM voltage-gated for low GJc has been previously
reported in the literature [5, 27], highlighting the importance of modeling the dynamic con-
ductance of gap junctions for cases of low GJc. Decreased GJc takes particular relevance in the
study of cardiac disease, as the reduction of gap-junctional communication has been correlated
to a marked decreased of CV [34], and slow conduction is considered one of the main mecha-
nisms of sustained reentrant arrhythmias [1, 35].

A unique feature of the NOHM model is its ability to upscale the particular features of volt-
age-dependent conduction mediated by homotypic and heterotypic combinations of homo-
meric connexons. In our simulations for the steady-state regime, action potentials resulting
from channels composed by homotypic Cx43-Cx43 resulted in a considerably higher CV
when compared to simulations considering homotypic Cx45-Cx45 channels (Fig 5 right col-
umn). This result is consistent with observations from dual whole-cell patch clamp experi-
ments, where the lower CV in Cx45-Cx45 channels is explained by the higher sensitivity of
conductance to transjunctional voltage [36]. We note, however, that for the instantaneous
regime, the CV is more sensitive to the number of operational channels and unitary conduc-
tance than to gap junction voltage dependence. Further, here we showed that the asymmetry
of conductance distribution found in heterotypic channels results in propagating action poten-
tials whose CV strongly depends on the direction of propagation (Fig 6). The NOHM success-
fully captured this orientation dependence, as well as it was able to capture the conduction
block predicted by the CM voltage-gated for low and very low GJc in the steady-state conduc-
tion regime. We also studied the effect of pacing rate for the Cx43-Cx45 channel, where both
voltage-sensitive models resulted in similar restitution curves with higher CVs for the case of
retrograde propagation (Fig 7). The orientation-dependent conduction, together with con-
nexin coexpression, may partly explain the differences in CV for normal and retrograde con-
duction that have been observed in the sinoatrial node [37]. It is important to note that
voltage-independent models of conduction (CM clamped, LHM) cannot capture this orienta-
tion dependence, as well as they fail to predict conduction block, resulting in a considerable
overestimation of the CV for the case of normal wave propagation. Future developments
should focus on combining the gap-junction conductance distributions, as several connexin
types are typically co-expressed in cardiac tissue.

We assessed the effect of spatial discretization for the continuum and cellular conduction
models considered in this study. Previous studies have shown that the numerical solution of
continuum electrophysiology models depends on the level of spatial discretization [38, 39]. An
interesting finding of this work is that the mesh dependence found in continuum models is
accentuated for cellular models and that the consideration of voltage sensitivity in the conduc-
tion model does not strongly affect the relation between CV and mesh size (Fig 8). Future
developments of numerical methods for the NOHM may include the consideration of
enhanced spatial interpolation and temporal integration schemes [40, 41] which have shown
to attenuate the mesh dependence of continuum model of cardiac propagation, allowing for
the efficient simulation of larger domains.

Our current work can be extended in several directions. First, the theoretical framework for
the NOHM model should be extended to consider the 3D case of cardiac conduction, includ-
ing the case of anisotropic conduction typically observed in cardiac tissue. We include a heu-
ristic derivation in Remark 2 of the S1 Appendix that points towards this direction. Another
important limitation is the consideration of the transmembrane potential, instead of the intra-
and extra-cellular potentials, in modeling intercellular conduction and ionic ionic currents. In
particular, (3) assumes that transjunctional voltage is the jump in transmembrane potential rather than the jump in intra-cellular potential. We note that such consideration is valid when the extra-cellular potential is constant, an assumption that can be debated in more general contexts of conduction. Future contributions should revisit this assumption by explicitly modeling the extra-cellular potential, i.e., considering bidomain formulations of the continuum electrophysiology problem [9]. Second, we note that the time dependence of gap-junction gating that dynamically modulates the conduction has not been considered in this work. Such a dynamic effect has shown to strongly modulate the conductance response of gap junctions [28, 29]. To date, the time-dependent gating of gap junctions has been incorporated in a few cellular models of cardiac conduction [5, 27], showing the importance of both voltage- and time-dependent dynamics. We note here that it takes several seconds for a gap-junction channel to reach steady-state conductance. Such a time scale can be much longer than the time window where transjunctional voltage is large during normal conduction, i.e., during action potential upstroke. Thus, during normal conduction, the steady-state regime is not expected to occur. In contrast, under cases of poor intercellular coupling, large transjunctional voltage can occur for longer periods, which potentially drive the gap-junction towards a steady-state conduction regime [27]. In this work, we only considered two limiting regimes of the dynamic conductance which yield very different behavior as the steady-state regime typically displays conductance distributions that are more sensitive to transjunctional voltage than those found for the instantaneous-conductance regime. Thus, an interesting avenue of research is the development of multiscale formulations of cardiac tissue conduction that incorporate time-dependent gating dynamics of gap junctions. Third, intercellular communication mechanisms other than gap junctions should be integrated into this theoretical framework. Sodium channels have been reported to co-localize with gap junctions at the intercalated discs, creating an ephestic coupling effect that has been associated to conduction during gap-junction blockage [35]. Further, the spatial distribution of sodium channels around the cellular membrane and on the intercalated discs has been studied using detailed cell-to-cell computational simulations to conclude that channel spatial distribution strongly affect the cardiac conduction [42]. Since the ephestic effect has been considered in homogenization schemes of cardiac conduction in the past by including a cleft-to-ground resistance in the microscopic model of conduction [13, 32], we foresee that future versions of the NOHM could equally incorporate this effect, potentially in 3D formulations with non-uniform distributions of channels. Finally, the applicability of the NOHM model should be tested in the simulation of conduction in the whole heart during diseased conditions [33].

**Supporting information**

**S1 Appendix. Formulation details.** Details and proofs for the asymptotic formulation and the numerical solution for the NOHM model are presented here.
(PDF)

**Author Contributions**

**Conceptualization:** Daniel E. Hurtado, Grigory Panasenko.

**Formal analysis:** Daniel E. Hurtado, Javiera Jilberto, Grigory Panasenko.

**Funding acquisition:** Daniel E. Hurtado.

**Investigation:** Daniel E. Hurtado, Javiera Jilberto.

**Methodology:** Daniel E. Hurtado, Javiera Jilberto.
**Project administration:** Daniel E. Hurtado.

**Resources:** Daniel E. Hurtado.

**Software:** Javiera Jilberto.

**Supervision:** Daniel E. Hurtado.

**Validation:** Daniel E. Hurtado, Javiera Jilberto.

**Visualization:** Daniel E. Hurtado, Javiera Jilberto.

**Writing – original draft:** Daniel E. Hurtado, Javiera Jilberto, Grigory Panasenko.

**Writing – review & editing:** Daniel E. Hurtado.

## References

1. Rohr S, Kucera JP, Kleber aG. Slow Conduction in Cardiac Tissue, I: Effects of a Reduction of Excitability Versus a Reduction of Electrical Coupling on Microconduction. Circulation Research. 1998; 83 (8):781-794. https://doi.org/10.1161/01.res.83.8.781 PMID: 9776725

2. Severs NJ, Coppen SR, Dupont E, Yeh HI, Ko YS, Matsushita T. Gap junction alterations in human cardiac disease. Cardiovascular Research. 2004; 62(2):368-377. https://doi.org/10.1016/j.cardiores.2003.12.007 PMID: 15094356

3. Desplantez T, Halliday D, Dupont E, Weingart R. Cardiac connexins Cx43 and Cx45: formation of diverse gap junction channels with diverse electrical properties. Pflugers Archiv: European journal of physiology. 2004; 448:363-375. https://doi.org/10.1007/s00424-004-1250-0 PMID: 15048573

4. Fast VG, Kléber AG. Microscopic Conduction in Cultured Strands of Neonatal Rat Heart Cells Measured With Voltage-Sensitive Dyes. Circulation Research. 1993; 73(5):914-925. https://doi.org/10.1161/01.res.73.5.914 PMID: 8403261

5. Henriquez AP, Vogel R, Muller-Borer BJ, Henriquez CS, Weingart R, Cascio WE. Influence of dynamic gap junction resistance on impulse propagation in ventricular myocardium: A computer simulation study. Biophysical Journal. 2001; 81(4):2112-2121. https://doi.org/10.1016/S0006-3495(01)75859-6 PMID: 11566782

6. Noble D. Modeling the heart—From genes to cells to the whole organ. Science. 2002; 295(5560):1678-1682. https://doi.org/10.1126/science.1069881 PMID: 11872832

7. Plonsey R, Barr RC. Mathematical modeling of electrical activity of the heart. Journal of Electrocardiology. 1987; 20(3):219-226. https://doi.org/10.1016/s0022-0736(87)80019-5 PMID: 3309111

8. Hurtado DE, Castro S, Gizzi A. Computational modeling of non-linear diffusion in cardiac electrophysiology: A novel porous-medium approach. Computer Methods in Applied Mechanics and Engineering. 2016; 300:70-83. https://doi.org/10.1016/j.cma.2015.11.014

9. Tung LA. Bidomain model for describing ischemic myocardial D-C potential. Massachusetts Institute of Technology; 1978.

10. Colli Franzone P, Pavarino LF, Scacchi S. Mathematical Cardiac Electrophysiology. vol. 13; 2014.

11. Neu JC, Krassowska W. Homogenization of syncytial tissues. Critical Reviews in Biomedical Engineering. 1993; 21(2):137-199. PMID: 8243090

12. Hand PE, Griffith BE, Peskin CS. Deriving macroscopic myocardial conductivities by homogenization of microscopic models. Bulletin of Mathematical Biology. 2009; 71(7):1707-1726. https://doi.org/10.1007/s11538-009-9421-y PMID: 19412638

13. Hand PE, Peskin CS. Homogenization of an electrophysiological model for a strand of cardiac myocytes with gap-junctional and electric-field coupling. Bulletin of Mathematical Biology. 2010; 72(6):1408-1424. https://doi.org/10.1007/s11538-009-9499-2 PMID: 20049544

14. Pennacchio M, Savare G, Colli Franzone P. Multiscale modeling for the electrical activity of the heart. SIAM J Math Anal. 2006; 37(4):1333-1370. https://doi.org/10.1137/040615249

15. Bruce D, Pathmanathan P, Whiteley JP. Modelling the Effect of Gap Junctions on Tissue-Level Cardiac Electrophysiology. Bulletin of Mathematical Biology. 2014; 76(2):431-454. https://doi.org/10.1007/s11538-013-9927-1 PMID: 24338526

16. Costa CM, Silva PAA, Santos RWD. Mind the Gap: A semicontinuum model for discrete electrical propagation in cardiac tissue. IEEE Transactions on Biomedical Engineering. 2016; 63(4):765-774. https://doi.org/10.1109/TBME.2015.2470256 PMID: 26292333
17. Bueno-Orovio A, Kay D, Grau V, Rodriguez B, Burrage K, Interface JRS. Fractional diffusion models of cardiac electrical propagation: role of structural heterogeneity in dispersion of repolarization. *Fractional diffusion models of cardiac electrical propagation: role of structural heterogeneity in dispersion of repolarization.* Journal of the Royal Society, Interface. 2014; 11:20140352. https://doi.org/10.1098/rsif.2014.0352 PMID: 24920109

18. Cusimano N, Gerardo-Giorda L. A space-fractional Monodomain model for cardiac electrophysiology combining anisotropy and heterogeneity on realistic geometries. *Journal of Computational Physics*. 2018; 362(February):409–424. https://doi.org/10.1016/j.jcp.2018.02.034

19. Cherubini C, Filippi S, Gizzi A, Ruiz-Baier R. A note on stress-driven anisotropic diffusion and its role in active deformable media. *Journal of Theoretical Biology*. 2017; 430:221–228. https://doi.org/10.1016/j.jtbi.2017.07.013 PMID: 28755956

20. Hurtado DE, Henao D. Gradient flows and variational principles for cardiac electrophysiology: Toward efficient and robust numerical simulations of the electrical activity of the heart. *Computer Methods in Applied Mechanics and Engineering*. 2014; 273:238–254. https://doi.org/10.1016/j.cma.2014.02.002

21. Logg A, Mardal KA, Wells GN, editors. Automated Solution of Differential Equations by the Finite Element Method. The FEniCS Book. 1st ed. Berlin Heidelberg: Springer-Verlag; 2012.

22. Diaz PJ, Rudy Y, Plonsey R. Intercalated discs as a cause for discontinuous propagation in cardiac muscle: A theoretical simulation. *Annals of Biomedical Engineering*. 1983; 11(3-4):177–189. https://doi.org/10.1007/bf02363285 PMID: 6670783

23. Wilders R, Jongsma HJ. Limitations of the dual voltage clamp method in assaying conductance and kinetics of gap junction channels. *Biophysical Journal*. 1992; 63(4):942–953. https://doi.org/10.1016/S0006-3495(92)81664-8 PMID: 1384745

24. Kucera JP, Rohr S, Rudy Y. Localization of sodium channels in intercalated disks modulates cardiac conduction. *Circulation Research*. 2002; 91(12):1176–1182. https://doi.org/10.1161/01.res.0000046237.54156.0a PMID: 12480819

25. Luo CH, Rudy Y. A Model of the Ventricular Cardiac Action Potential. *Circulation Research*. 1991; 68(6):1501–1526. https://doi.org/10.1161/01.res.68.6.1501 PMID: 1709839

26. Shaw RM, Rudy Y. Ionic mechanisms of propagation in cardiac tissue: Roles of the sodium and L-type calcium currents during reduced excitability and decreased gap junction coupling. *Circulation Research*. 1997; 81(5):727–741. https://doi.org/10.1161/01.res.81.5.727 PMID: 9351447

27. Weinberg SH. Ephaptic coupling rescues conduction failure in weakly coupled cardiac tissue with voltage-gated gap junctions. *Chaos*. 2017; 27(093908). https://doi.org/10.1063/1.4999602 PMID: 28964133

28. Vogel R, Weingart R. Mathematical model of vertebrate gap junctions derived from electrical measurements on homotypic and heterotypic channels. *Journal of Physiology*. 1998; 510(1):177–189. https://doi.org/10.1111/j.1469-7793.1998.177bz.x PMID: 9625876

29. Harris AL, Spray DC, Bennett MV. Kinetic properties of a voltage-dependent junctional conductance. *The Journal of general physiology*. 1981; 77(January):95–117. https://doi.org/10.1085/jgp.77.1.95 PMID: 6259275

30. Gizzi A, Loppini A, Ruiz-Baier R, Ippolito A, Camassa A, Camera AL, et al. Nonlinear diffusion and thermo-electric coupling in a two-variable model of cardiac action potential. *Chaos*. 2017; 27(9). https://doi.org/10.1063/1.4999610 PMID: 28964112

31. Söhl G, Maxeiner S, Willecke K. Expression and functions of neuronal gap junctions. *Nature Reviews Neuroscience*. 2005; 6(3):191–200. https://doi.org/10.1038/nrn1627 PMID: 15738956

32. Hand PE, Griffith BE. Adaptive multiscale model for simulating cardiac conduction. *Proceedings of the National Academy of Sciences*. 2010; 107(33):14603–14608. https://doi.org/10.1073/pnas.1008443107

33. Prakosa A, Arevalo HJ, Deng D, Boyle PM, Nikolov PP, Ashikaga H, et al. Personalized virtual-heart technology for guiding the ablation of infarct-related ventricular tachycardia. *Nature Biomedical Engineering*. 2018; 2(10):732–740. https://doi.org/10.1038/s41551-018-0282-2 PMID: 30847259

34. Dhillon PS, Gray R, Kojodjojo P, Jabr R, Chowdhury R, Fry CH, et al. Relationship between gap-junctional conductance and conduction velocity in mammalian myocardium. *Circulation: Arrhythmia and Electrophysiology*. 2013; 6(6):1208–1214.

35. Tse G, Yeo JM. Conduction abnormalities and ventricular arrhythmogenesis: The roles of sodium channels and gap junctions. *IJC Heart and Vasculature*. 2015; 9:75–82. https://doi.org/10.1016/j.ijcha.2015.10.003 PMID: 26839915

36. Desplantez T, Dupont E, Severs NJ, Weingart R. Gap junction channels and cardiac impulse propagation. *Journal of Membrane Biology*. 2007; 218(1-3):13–28. https://doi.org/10.1007/s00232-007-9046-8 PMID: 17661127
37. Verheijck EE, Van Kempen MJA, Veereschild M, Lurvink J, Jongsma HJ, Bouman LN. Electrophysiological features of the mouse sinoatrial node in relation to connexin distribution. Cardiovascular Research. 2001; 52(1):40–50. https://doi.org/10.1016/s0008-6363(01)00364-9 PMID: 11557232

38. Krishnamoorthi S, Sarkar M, Klug WS. Numerical quadrature and operator splitting in finite element methods for cardiac electrophysiology. International Journal for Numerical Methods in Biomedical Engineering. 2013; 29:1243–1266. https://doi.org/10.1002/cnm.2573 PMID: 23873868

39. Pezzuto S, Hake J, Sundnes J. Space-discretization error analysis and stabilization schemes for conduction velocity in cardiac electrophysiology. International Journal for Numerical Methods in Biomedical Engineering. 2016. https://doi.org/10.1002/cnm.2762 PMID: 26685879

40. Hurtado DE, Rojas G. Non-conforming finite-element formulation for cardiac electrophysiology: an effective approach to reduce the computation time of heart simulations without compromising accuracy. Computational Mechanics. 2018; 61(4):485–497. https://doi.org/10.1007/s00466-017-1473-5

41. Jilberto J, Hurtado DE. Semi-implicit Non-conforming Finite-Element Schemes for Cardiac Electrophysiology: A Framework for Mesh-Coarsening Heart Simulations. Frontiers in Physiology. 2018; 9 (1513):1–12.

42. Horgmo Jæger K, Edwards AG, Mcculloch A, Tveito A. Properties of cardiac conduction in a cell-based computational model. PLoS Computational Biology. 2019; 15(5):e1007042. https://doi.org/10.1371/journal.pcbi.1007042