Skip to main content

ORIGINAL RESEARCH article

Front. Physiol., 03 January 2024
Sec. Cardiac Electrophysiology
This article is part of the Research Topic Exploring Mechanisms of Cardiac Rhythm Disturbances Using Novel Computational Methods: Prediction, Classification, and Therapy: Volume II View all 6 articles

Nernst-Planck-Gaussian modelling of electrodiffusional recovery from ephaptic excitation between mammalian cardiomyocytes

  • 1Physiological Laboratory, University of Cambridge, Cambridge, United Kingdom
  • 2Department of Veterinary Medicine, University of Cambridge, Cambridge, United Kingdom
  • 3Department of Biochemistry, University of Cambridge, Cambridge, United Kingdom

Introduction: In addition to gap junction conduction, recent reports implicate possible ephaptic coupling contributions to action potential (AP) propagation between successive adjacent cardiomyocytes. Here, AP generation in an active cell, withdraws Na+ from, creating a negative potential within, ephaptic spaces between the participating membranes, activating the initially quiescent neighbouring cardiomyocyte. However, sustainable ephaptic transmission requires subsequent complete recovery of the ephaptic charge difference. We explore physical contributions of passive electrodiffusive ion exchange with the remaining extracellular space to this recovery for the first time.

Materials and Methods: Computational, finite element, analysis examined limiting, temporal and spatial, ephaptic [Na+], [Cl], and the consequent Gaussian charge differences and membrane potential recovery patterns following a ΔV∼130 mV AP upstroke at physiological (37°C) temperatures. This incorporated Nernst-Planck formalisms into equations for the time-dependent spatial concentration gradient profiles.

Results: Mammalian atrial, ventricular and purkinje cardiomyocyte ephaptic junctions were modelled by closely apposed circularly symmetric membranes, specific capacitance 1 μF cm-2, experimentally reported radii a = 8,000, 12,000 and 40,000 nm respectively and ephaptic axial distance w = 20 nm. This enclosed an ephaptic space containing principal ions initially at normal extracellular [Na+] = 153.1 mM and [Cl] = 145.8 mM, respective diffusion coefficients DNa = 1.3 × 109 and DCl = 2 × 109 nm2s-1. Stable, concordant computational solutions were confirmed exploring ≤1,600 nm mesh sizes and Δt≤0.08 ms stepsize intervals. The corresponding membrane voltage profile changes across the initially quiescent membrane were obtainable from computed, graphically represented a and w-dependent ionic concentration differences adapting Gauss’s flux theorem. Further simulations explored biological variations in ephaptic dimensions, membrane anatomy, and diffusion restrictions within the ephaptic space. Atrial, ventricular and Purkinje cardiomyocytes gave 40, 180 and 2000 ms 99.9% recovery times, with 720 or 360 ms high limits from doubling ventricular radius or halving diffusion coefficient. Varying a, and DNa and DCl markedly affected recovery time-courses with logarithmic and double-logarithmic relationships, Varying w exerted minimal effects.

Conclusion: We thereby characterise the properties of, and through comparing atrial, ventricular and purkinje recovery times with interspecies in vivo background cardiac cycle duration data, (blue whale ∼2000, human∼90, Etruscan shrew, ∼40 ms) can determine physical limits to, electrodiffusive contributions to ephaptic recovery.

1 Introduction

1.1 Gap junctional and ephaptic conduction of cardiomyocyte excitation

It has long been established that action potential (AP) conduction between successive cardiomyocytes within myocardial syncytia involves local circuit current flow through low resistance gap junctions (GJs) (Rohr, 2004). The required intercellular conductivities largely arise from connexin Cx40 (unit conductances 120–180 pS) and/or Cx43 channels (unit conductances 40–70 or 90-100 pS depending on phosphorylation conditions) bridging adjoining cells (Moreno et al., 1994; Traub et al., 1994). Ventricular myocyte gap junctions immunostained with anti-Cx43 (and anti-Cx45), but not anti-Cx40 antibodies. Atrial gap junctions immunostained with both anti-Cx40 and anti-Cx43 (and anti-Cx45) antibodies. Double whole cell voltage clamped rabbit atrial and ventricular gap junctions conductance nevertheless showed similar respective, 169 ± 146 and 175 ± 147 nS, average macroscopic junctional conductances. Ventricular myocyte pairs showed a single population of 100 pS, and atrial myocyte pairs, two populations of 100 and 185 pS single gap junction channel conductances (Verheule et al., 1997). Membrane potential change in an active cell thereby electrotonically depolarises the adjacent, initially quiescent, cell within the syncytium, with implications for conduction velocity in turn predisposing to pro-arrhythmic re-entrant circuits. Thus, in both experimental and clinical pro-arrhythmic pathological situations, remodelling gap junction relocation to lateral as opposed to perinexal membrane (Peters et al., 1997; Vetterlein et al., 2006) occurs with the acute ischaemia following coronary occlusion (Severs et al., 2008), advanced chronic ischaemic cardiac disease (Smith et al., 1991), cardiac failure arising from ischaemic, dilated and inflammatory cardiomyopathy (Kitamura et al., 2002; Kostin et al., 2003; 2004), and systemic or pulmonary hypertensive conditions (Kostin et al., 2004). Similarly, albeit likely multifactorial in pathophysiology, atrial fibrillation is associated with altered Cx40 expression and/or its lateral distribution (Polontchouk et al., 2001; Gemel et al., 2014; Lambiase and Tinker, 2015). Finally, though viable, adult Cx40-deficient, Cx40−/−/Cx43+/+ and Cx43-heterozygous, Cx40+/+/Cx43+/−mouse hearts show conduction abnormalities. Double heterozygous Cx40+/−/Cx43+/−showed ECG features suggesting additive Cx40 and Cx43 haploinsufficiency effects on ventricular, but not atrial, conduction (Kirchhoff et al., 2000).

However, propagation of myocardial excitation can persist under circumstances of reduced gap junction coupling. This is exemplified by clinical pathological conditions following ischaemic insult (Rohr, 2004), atrial fibrillation (Chaldoupi et al., 2009) and experimental platforms with loss-of-function, Cx43 (Eloff et al., 2001; Gutstein et al., 2001) or Cx40, genetic modifications (Hagendorff et al., 1999; Bagwe et al., 2005). These findings prompted suggestions of additional or alternative excitation propagation mechanisms including direct cell-to-cell, ephaptic, conduction (Sperelakis and Mann, 1977; Carmeliet, 2019). Widespread such examples occur in peripheral and central nervous systems of numerous species. They can mediate cell-cell, bidirectional or unidirectional excitation (Furshpan and Potter, 1959), either alone (Watanabe and Grundfest, 1961; Bennett, 1977) or combined with other complementary mechanisms (Rash et al., 1996).

1.2 Perinexal regions between cardiomyocytes

In the heart, adjacent cardiomyocytes are joined end-to-end at specialised perinexal regions (Figure 1Aa). These contain structural specialisations compatible with ephaptic function (Rhett and Gourdie, 2012; Rhett et al., 2013; Veeraraghavan et al., 2014). Electron microscopic appearances indicated closely apposed adjacent cell membranes (Veeraraghavan et al., 2015). These showed elevated Cx43 and Na+ channel (Nav1.5) expression. Super-resolution microscopy localised high Nav1.5 densities within 200 nm of the GJ plaques (Veeraraghavan et al., 2015; Hichri et al., 2018; Struckman et al., 2021), likely reflecting recently demonstrated intrinsic Nav1.5 clustering properties (Salvage et al., 2020b), also reported for skeletal muscle Nav1.4 (Sheikh et al., 2001). Smart patch clamping correspondingly indicated greater Na+ current densities (INa) at perinexal than at non-junctional sites (Veeraraghavan et al., 2018). Such high, intercalated disk, Nav1.5 densities might thus support ephaptic conduction.

FIGURE 1
www.frontiersin.org

FIGURE 1. Background geometrical parameters describing the model. (A) The intercalated disc (a) and ephaptic junction (b) with three-dimensional arrangement of the intercalated disc (c) gap-junctions hold the membranes of two adjacent cells in close apposition with Nav 1.5 β-subunits restricting the axial distance in other regions Reproduced from Salvage et al. (2020a), licensed under CC-BY 4.0. (B) Simplified physical scheme describing recovery processes at the ephaptic junction axial distance, w, radius, a containing Na+ and Cl with diffusion coefficients DNa, DCl, following initial Na+ withdrawal consequent upon action potential generation (a). This lists electrodiffusive (b) and other, Na+-K+-ATPase (c) and Kir2.1 (d) mediated membrane transport contributions to ionic composition and consequent membrane potential in the ephaptic space.

Super-resolution imaging additionally demonstrated preferential perinexal β1 (SCN1B) subunit colocation with the Nav1.5 (Veeraraghavan et al., 2018). These could form adhesion trans scaffolds adjacent to the gap junctions across the consequently narrow (∼20 nm) perinexal clefts (Salvage et al., 2020a) (Figure 1Ab). βadp1 peptide selectively inhibited this β1-mediated adhesion and widened guineapig ventricular perinexi. It selectively reduced perinexal but not whole cell INa precipitating arrhythmogenic conduction slowing (Veeraraghavan et al., 2018). These findings together suggested a Nav1.5 mediated ephaptic cell-cell transfer of electrical excitation, potentially offering a novel therapeutic anti-arrhythmic target. This could explain persistent myocardial electrical conduction in mouse models where electrical coupling was decreased by Cx43 knockout (Gutstein et al., 2001).

1.3 Previous modelling of ephaptic activation

Mathematical modelling suggested that ephaptic mechanisms could mediate AP conduction despite compromised gap junctional coupling (Sperelakis and Mann, 1977; Sperelakis, 2002; Mori et al., 2008). Ephaptic activation then provided a slower conduction requiring the nevertheless experimentally reported narrow perinexal axial distance (cleft width) and high Nav1.5 densities. The perinexal space constitutes a restricted extracellular diffusional space permitting formation of microdomains containing reduced [Na+] (Veeraraghavan et al., 2015). These could result from Na+ withdrawal during the AP upstroke of the initially active cell. The consequent negative potential within the resulting Gaussian space could depolarise the membrane of the adjacent, initially quiescent, cell. Should the latter then attain its Nav1.5 activation threshold it initiates an AP (Sperelakis, 2002). Additionally, gated stimulated emission depletion (gSTED) and stochastic optical reconstruction microscopy (STORM) super-resolution microscopy, attribute a significant perinexal K+ conductance to Kir2.1 expression particularly associated with desmosomes (Veeraraghavan et al., 2016; Struckman et al., 2021). Its marked inward rectification property reducing K+ conductance on depolarisation (Difranco et al., 2015; Chen et al., 2018) would enhance ephaptic excitation by minimising electrical shunting of the Na+ transfer. Existence of ephaptic activation between adjacent cardiomyocytes involving capacitive coupling of the closely apposed membranes thus has both experimental and theoretical support (Mori et al., 2008; Lin and Keener, 2010).

1.4 Requirements for ephaptic recovery

However, sustained operation of the resulting ephaptome also requires complete recovery of the altered extracellular ephaptic ionic concentration within each cardiac cycle. On the one hand, this could involve actions of other perinexal membrane molecules. As well as the increased Nav1.5 densities in intercalated disk regions, and Kir2.1 particularly around desmosomes, recent confocal microscopy and super resolution studies reported increased signals reflecting Na+-K+-ATPase occurring around the gap and adherence junctions (Figure 1B) (Noël et al., 1991; Mcdonough et al., 1996; Struckman et al., 2021). These patterns were more marked in atrial relative to ventricular myocytes. Finally, electrodiffusive exchanges particularly involving the principal extracellular ions Na+ and Cl could occur within the ephaptic space, and between the ephaptic space and the remaining extracellular space. Available information bearing on 3-dimensional perinexal geometry (Rhett and Gourdie, 2012; Gourdie, 2019) suggests that gap junction plaques forming relatively large regions of 100s of gap junctions surrounded by perinexal regions containing Nav1.5 clusters (Salvage et al., 2020b) ‘pin’ together the component membranes at discrete sites (Figure 1Ac) (Salvage et al., 2020a). However, they leave an ephaptic space offering pathways for diffusion around them. This permits diffusive interconnections between different ephaptic regions within the intercalated disc, as well as to the bulk extracellular fluid, possibly exerting some restriction to free diffusion. These could contribute to dissipation of the resulting ephaptic membrane potential change. Thus, acute interstitial oedema (AIE)-induced perinexal swelling produced an arrhythmogenic conduction slowing respectively mitigated and exacerbated by Kir2.1 inhibition and Nav1.5 block (Veeraraghavan et al., 2016). In addition, such recovery times have implications for the effectiveness of possible ephaptic conduction contributions at different heart rates (see Discussion).

1.5 Modelling electrodiffusive recovery

The present computational studies assess limiting features of such electrodiffusive recovery processes following AP activation and their variation with both ephaptic junction geometry and diffusional properties, for the first time. The modelling (1) employed experimentally reported values describing ephaptic junction anatomy in mammalian atria, ventricles and purkinje fibres and established initial values of extracellular electrolyte, Na+, Cl, concentrations and their diffusion coefficients (Schoenberg et al., 1975; Severs, 2000; Mori et al., 2008; Richards et al., 2011; Veeraraghavan et al., 2015; 2018). It then (2) reproduced changes in these within the restricted ephaptic space resulting from a ΔV∼130 mV AP upstroke at physiological (37°C) temperatures and (3) followed their subsequent recovery. The latter involved determination of ionic concentrations and their variation with time and through the three dimensional ephaptic geometry (Léonetti, 1998; Koch, 2004) employing a finite element Nernst-Planck modelling of the consequent extracellular Na+ and Cl electrodiffusional processes. (4) The consequent membrane potential changes resulting from the imbalance in ion charges with finite diffusion coefficients within a restricted extracellular space were derived by an adaptation of Gauss’s Flux Theorem. Such a charge difference approach was introduced in previous membrane potential modelling studies from properties of multiple channel components (Fraser and Huang, 2004; 2007). This made it possible (5) to explore the robustness of the solutions and the effect of biological variations in ephaptic dimensions, membrane anatomy, and diffusion restrictions. These features made it possible to consider potential comparative roles of electrodiffusive and membrane transport processes in ephaptic recovery in different, atrial, ventricular and purkinje cardiomyocyte types in relationship to known physiological heart rates in different species.

2 Theory

2.1 Forms of the Nernst-Planck equation describing three-dimensional electrodiffusional fluxes

Solving for the time and spatial dependence of ephaptic ion concentrations, c, electrodiffusionally recovering following initial Na+ withdrawal employed the fundamental time (t) dependent Nernst-Planck equations. This relates flux, J, to concentration c, and electrical potential φ for an ion with diffusion coefficient, D and valence, z assuming Faraday constant, F, gas constant, R, and absolute temperature, T:

J=Dc+zFRTcφ.(1)

Substituting [1] in the conservation relationship:

ct=.J(2)

Gives:

ct=.Dc+zFRTcφ(3)

Applying the product rule and collecting terms to the left hand side of the equation:

ct.DcDzFRTc2φ)DzFRTc.φ=0(4)

Splitting the Nernst Planck Equation into two equations representing positive (+ve) and negative (-ve) ions:

c+vetD+ve2c+veD+vezFRTc+ve2φD+vezFRTφ.c+ve=0(5)
cvetDve2cveDvezFRTcve2φDvezFRTφ.cve=0(6)

2.2 Initial action potential mediated Na+ transfer from the ephaptic space

These equations were solved for an ephaptic space comprising a circularly symmetric gap, radius a and axial separation w, filled with extracellular fluid, free-space permittivity ε0, and water dielectric constant εc, containing principal charge carriers Na+ and Cl. Their initial conditions modelled effects of an AP upstroke withdrawing [Na+] from this space uniformly through the active membrane, perturbing ionic concentration differences within the ephaptic space. Derivation of limiting minimum values for this Na+ transfer employed active membrane surface area:

SAP=πa2.(7)

Assuming specific membrane capacitance, Cm, the Na+ current then discharges capacitance,

πa2Cm.(8)

The charge q transferred across the active membrane, by membrane potential change ΔV is then:

q=ΔVπa2Cm(9)

The charge of a single, univalent (z = +1 or −1) Na+ or Cl ion equals the elementary charge e =1.6021 × 10−19 C. Assuming Avogadro’s number NA = 6.0221367 × 1023 mol-1 (Lide, 1993), the molar Na+ transferred from the ephaptic space is therefore:

MolesofNa+transferred=ΔVCmπa2eNA(10)

For example, for atrial cells with typical ephaptic radius a = 8 × 10−6 m and axial distance w = 20 × 10−9 m, and specific plasma membrane capacitance Cm = 0.01 F m-2, a membrane potential change ΔV = +130 × 10−3 V (Table 1, see also (Draper and Weidmann, 1951; Hilgemann and Noble, 1987)), then transfers 2.7090 × 10−18 mol Na+ from an initially neutral ephaptic space.

TABLE 1
www.frontiersin.org

TABLE 1. Basic parameters for cardiac ephaptic junctions.

2.3 Gaussian determination of the voltage recovery from changes in the ion concentration difference

The Nernst-Planck electrodiffusion analysis then examined spatial and temporal changes in ephaptic space ionic concentrations and their consequences for membrane potential following the initial ephaptic excitation. The ephaptic space was approximated to two closely apposed, active and passive, membranes radius a separated by small ephaptic gap distance w whose margins access the bulk extracellular compartment (Figure 1B),

wa.(11)

The membrane areas (Eq. 12) therefore greatly exceed that of the rim of the ephaptic gap:

2πawπa2.(12)

Following action potential recovery and early rapid axial diffusive movements these might first produce an axially and radially uniform ionic concentration change, and radially uniform membrane potential and electric field changes across the enclosing plasma membranes. The subsequent electrodiffusive recovery primarily involving fluxes of the principal extracellular ions Na+ and Cl then result in the time and position dependent [Na+] and [Cl] of c+ve and c-ve respectively. The consequent time and spatially dependent radial ionic concentration differences {c+ve - c-ve} within the ephaptic space cause correspondingly time and spatially dependent membrane potential alterations, ∆V. Gauss’s Flux Theorem gives the electric flux δΦΕ^ given by the surface integral of the electric field Ê over any closed area S generated by enclosed net charge δq in medium of relative permittivity εr (Lide, 1993):

δΦΕ^=δΕ^.dS=δqε0εr(13)

Consider a surface falling within, and in which most of whose electrical flux traverses, the plasma membranes of relative permittivity εr, enclosing the ephaptic space. As 2πaw << πa2 (Equation (12)), and εr << εc, the position and time-dependent transmembrane membrane potential change ∆V tends to the corresponding φ within the adjoining recovering ephaptic space. The latter geometrically approximates a series of successive coaxial, circularly symmetric, concentric annular volume elements. Each element contains progressively altered {c+ve - c-ve}. These charge membranes of permittivity ε forming the two axial ends of each element. The total membrane cross-sectional area δS in any given annular element of inner radius α and outer radius α + δ α), is:

δS=2πα+δα2πα2=2π2α+δαδα=4παδα+2πδα2(14)

The volume of each annular element in a ephaptic gap of axial width (w) is then:

wδS2(15)

The concentration of charge, [q] within that volume element is the product of the ionic concentration difference, c+vecve, and Faraday’s constant (F), for Na+ and Cl ions of valencies z = +1 and −1 respectively:

q=zFc+vecve(16)

The corresponding quantity of charge present within that element, δq, is then the concentration of the charge q multiplied by its volume. From equations (15) and (16):

δq=wzFδS2c+vecve(17)

In the limit δα0, adjoining annular elements tend towards equal δS, each containing equal charge δq minimizing the proportion of the total electric flux between elements. Most of the flux then traverses the ephaptic membranes at the two ends of each element. For membranes of thickness, ζ, total surface area S, whose electric field is uniform (Goldman, 1943), ΔÊ = φ:

ΔΦÊ=ΔÊ.dS=φζδS(18)

Equating the right sides of Eqs. (13), (18) gives the alteration in transmembrane voltage:

φ=ζδqε0εrδS=δqCmδS,(19)

since:

Cm=ε0εrζ(20)

Finally, substituting the expressions (Eq. 17) for δq:

φ=zFwc+vecve2Cm(21)

This gives:

φ=zFwc+vecve2Cm=zFw2Cmc+vecve(22)

This treatment employing Gauss’s Flux theorem extends previous analyses deriving cell resting potentials from intracellular charge concentration differences (Fraser and Huang, 2004; 2007). The latter studies applied Eq. (19) to a cell of volume ϑ containing homogeneous intracellular electrolyte and protein concentrations Na+i,K+i,Cliand Przi and consequently uniform resting membrane potential φ. The δqδS term becomes the ratio between total charge difference related to the intracellular (i) ion concentrations, ϑFNa+i+K+iCli+zxPrzi, and total membrane surface area S (Fraser and Huang, 2004). The ratio between the resulting total surface membrane capacitance SCm and the ϑ term gives the capacitance of unit cell volume Cϑ yielding a resting potential equation previously derived by alternative means (Fraser and Huang, 2004):

φ=FNa+i+K+iCli+zxPrziCϑ(23)

2.4 Computational implementation of the Nernst-Planck formalisms

Combining Eqs. (5), (6), (22) and the Poisson Equation (Malmberg and Maryott, 1956; Gabriel et al., 1996; Kramlich et al., 2012):

2φ=eε0εccvec+ve(24)

yielded analytic equations for the present specific Nernst-Planck diffusion process over time for numerical solution by finite element analysis (FEA). For the univalent Na+ (z = 1) and Cl (z = −1):

For c+ve (z = +1):

c+vetD+ve2c+veD+veFeRTε0εcc+vecvec+veD+veF2wRT2Cmc+vecve.c+ve=0(25)

For cve (z = −1):

cvetDve2cve+DveFeRTε0εccvecvec+ve+DveF2wRT2Cmc+vecve.cve=0(26)

These can be re-arranged into a form solvable by the matrix-based MATLAB platform:

λ12ct2+λ2ct.λ3c+λ4cλ5=0(27)

Thus, applying the product rule:

.λ3c=λ3.c+λ32c(28)

gives:

λ12ct2+λ2ctλ32cλ3.c+λ4cλ5=0(29)

If λ3 is constant, λ3.c=0, as confirmed in Eqs. (31), (32) below:

λ12ct2+λ2ctλ32c+λ4cλ5=0(30)

The partial differential equations (PDEs) were solved using the Partial Differential Equation (PDE) Toolbox™. The coefficients λ1, λ2…. λ5 can be functions of location (x, y, and, in 3-D, z), and, except for eigenvalue problems, can also be functions of the solution c or its gradient. Thus, comparing terms between equation (25) and (30), for c+ve:

λ1=0λ2=1λ3=D+veλ4=D+veFeRTε0εccvec+veλ5=D+veF2wRT2Cmc+vecve.c+ve(31)

Comparing terms between Eqs. (26), (30), for cve:

λ1=0λ2=1λ3=Dveλ4=DveFeRTε0εccvec+veλ5=DveF2wRT2Cmc+vecve.cve(32)

3 Materials and methods

The Nernst-Planck diffusion of Na+ and Cl, and their subsequent patterns of temporal and spatial ion concentration and membrane potential differences were examined in a computational model of cardiac ephaptic junctions. Their anatomy was replicated using published values reported for mammalian atrial, ventricular, and purkinje cardiomyocytes. Further computations were also applied following proportional changes in the axial distances and radii of the junctions examined, and diffusion coefficients defining the underlying flux processes. The analysis modified a previous approach examining diffusional fluxes in Ca2+ microdomain formation at skeletal muscle triad junctions (Bardsley et al., 2021). This used the matrix-based MATLAB language, implemented on the MATLAB platform, in a finite element analysis (FEA) of ion diffusion over time. The partial differential equations (PDEs) involved were solved using the Partial Differential Equation (PDE) Toolbox™, modifying the relevant code to model Nernst-Planck diffusion.

The studies first involved meshing the geometry and selection of physical conditions (initial and boundary conditions). This permitted generation of a set of solution matrices for subsequent processing and figure presentation. The studies used the MATLAB matrix-based programming platform and language (version R2021b win64 9.11.0.1769968, version released 22 September 2021, MathWorks, Cambridge, UK) (https://www.mathworks.com/discovery/what-is-matlab.html). The programmes performed the required data array manipulations and generated all the graphics shown in the results section excluding Figure 8 which was produced using SciDAVis (version 2.4, version released fifth May 2021, developed by Miquel Garriga, Arun Narayanankutty, Dmitriy Pozitron, Russel Standish). This modelling was performed using an IBM compatible computer (CPU: Intel® Core™ i5-8265U CPU@1.60GHz; installed RAM: 8 GB, running Windows (Microsoft, Washington, USA) 11, Home 64-bit version 21H2).

The ephaptic junction geometry was generated within MATLAB and then divided into a set of tetrahedral sub-domains (elements) joined along their edges and at vertices (nodes), to be used in a FEA solving a set of PDEs. This culminated with a representation of the global solution whilst retaining local effects. We made use of both Neumann-Type and Dirichlet-Type boundary conditions for the modelling of ephaptic junction recovery which provide values of the field variables at the boundaries of our geometry. FEA equations are used to calculate the value of the field variable at nodes, interpolation of these values is then used to calculate values for non-nodal points (Rapp, 2017). Where two tetrahedra connect (vertices and edges) the value of the field variable is the same between the two elements. This ensures continuity between them, preventing gaps in the solution. Each element represents our set of PDEs, which are time-dependent, with a set of ordinary differential equations (ODEs) that approximate the solution to the original equations (Ramirez, 1997). This set of ODEs is then numerically integrated by MATLAB solvers to provide a solution, for which a source listing is provided in the Supplementary file.

Meshing and the FEA performed used the PDE Toolbox™ (version 3.7 installed on 27 September 2021 by MathWorks) within MATLAB. As indicated above, the PDE Toolbox™ requires PDEs to be in the form:

λ12ut2+λ2ut.λ3u+λ4u=λ5(27)

Where the solution u is the concentration of a given ion, and λ1, λ2…. λ5 are inputted coefficients. This automatically meshes the ephaptic junction geometry given a maximum mesh size (maximum edge length), providing a platform for implementing our set of PDEs. It then outputs and stores the solutions as matrices. These could then be accessed for processing of the raw-data and subsequent presentation in MATLAB.

4 Results

4.1 Na+ transfer across the active membrane: Initial conditions for computational modelling

The finite element analysis (FEA) employing the PDE Toolbox™ modelled electrodiffusive ion fluxes through stepsizes representing both its spatial and temporal finite elements. The ephaptic geometry was meshed into tetrahedral elements joined along their edges and at their vertices (nodes) where the side length, defining the mesh size, determined the spatial resolution of the computational solutions. In turn, the time step adopted for the solver gives the stepsize interval between successive iterations, therefore determining the temporal resolution. In general, smaller stepsize generate more accurate and stable solutions whilst requiring greater computer resources. Computational solutions were also obtained using combinations of varied spatial mesh and temporal stepsize to ensure appropriate optimisations between mathematical conditions.

Table 1 summarises the adopted, previously reported, physical and geometrical properties of the modelled atrial, ventricular and purkinje fibre ephaptic spaces and relevant physical constants. Values of atrial, ventricular and purkinje cardiomyocyte fibre radius, a, were obtained from previous experimental reports on hearts from large mammals (sheep). The selected values additionally permitted parallel investigations of the effect of variations in a on the computed solutions (Schoenberg et al., 1975; Severs, 2000; Richards et al., 2011) and comparisons with previous modelling results on ephaptic activation. Thus the value for ventricular radius was close to that used in recent reports (Mori et al., 2008). Values of axial ephaptic distance w employed recent, 15–30 nm, cell-cell distance measurements (Veeraraghavan et al., 2015; 2018). In any case, the investigations explored the effects of twofold alterations in w on the solutions. Similarly, details of membrane geometry that could affect Na+ and Cl diffusion were reflected in calculations varying their respective diffusion coefficients DNa and DCl. Thus, the wide range of a, w, DNa and DCl investigated permitted a broad clarification of the effects of these parameters on the electrodiffusive recovery process. These parameters defined the geometrical boundaries of the space under analysis, specifically: the active membrane performing the initiating AP, the passive membrane responding to the resulting ephaptic conduction and the rim through which the ephaptic space, in which the consequent electrodiffusive processes under analysis take place, communicates with the surrounding extracellular space.

The computational solutions involved two stages. First, they were initiated by modelling an AP mediated Na+ transfer adopting a limiting minimum value for the Na+ transfer from the ephaptic space of radius a and axial width w. The magnitude of this Na+ transfer was determined by the Na+ flux produced by activation of evenly distributed, voltage-gated Nav1.5 channels responsible for driving the AP upstroke in the active membrane. This formal analysis assumed a membrane potential change ΔV = +130 mV (Draper and Weidmann, 1951; Hilgemann and Noble, 1987) across a plasma membrane, without infoldings, of specific capacitance Cm = 1 μF/cm2. Equations (7)–(10) permitted calculations of the charge q and moles of Na+ transferred across the active membrane of the ephaptic space. The biophysically derived value of Na+ flux was converted into a flux density referred to unit area of active ephaptic membrane surface. This provided the first, spatial, Neumann, boundary condition of the active ephaptic membrane. Boundary conditions at the ephaptic rim and the passive membrane were here set to zero indicating a zero ionic flux. The computational runs then determined the resulting axial and radial gradients of [Na+] and [Cl] following this initial Na+ transfer.

Secondly, following the above 0.7 ms duration step, 2.7090 × 10−18 mol, Na+ transfer from the ephaptic space representing the AP upstroke, a Nernst-Planck electrodiffusion analysis examined limiting properties of ephaptic space [Na+] and [Cl] recovery, and consequent membrane potential changes in the passive membrane. The Neumann boundary conditions of the active and passive membranes were now both set to zero representing zero flux during this recovery. The ephaptic rim was assigned Dirichlet boundary conditions holding [Na+] and [Cl] at 153.1 mM and 154.8 mM respectively representing a well stirred bulk extracellular fluid of constant concentrations. Figure 2 illustrates early changes in atrial cells, ephaptic radius a = 8 × 10−6 m (Richards et al., 2011), axial distance w = 20 × 10−9 m (Veeraraghavan et al., 2018), and specific membrane capacitance Cm = 0.01 F m-2 and a ΔV = +130 × 10−3 V membrane potential change (Table 1). The computation employed a 400 nm mesh size and 0.01 ms stepsize interval. The initial 400 nm mesh size corresponded to 5% of the ephaptic radius, and the 0.01 ms stepsize interval entailed 100 iterations per ms time modeled. Figures 2A–C provides heat maps of axial (vertical axis) and radial (horizontal axis) [Cl] (A) and [Na+] at mM (B) and the resulting charged ion concentration differences at higher μM, resolution (C) immediately following the Na+ transfer. Radial [Na+] gradients here were expectedly absent. Axial [Na+] and charged ion concentration differences, [q], were also extremely small and only detectable at high colour magnifications. The very small size of these gradients and the time-period of recovery to allow dissipation of these gradients permitted adoption of uniform axial concentration profiles in the initial conditions.

FIGURE 2
www.frontiersin.org

FIGURE 2. Midline slice heat maps of concentration and voltage profiles close to the end of the action potential upstroke. (A–C) Heat maps of [Cl] (A) and [Na+] (B) and resulting ion concentration and membrane potential differences (C) at different axial (vertical axis) and radial positions (horizontal axis) in the atrial ephaptic space. Slices obtained at end of a 0.7 ms period of sodium extrusion representing the cardiac action potential upstroke. Note the very small axial gradients generated. Given the size of these gradients and the time-period of recovery to allow dissipation of these gradients we have assumed a uniform concentration profile for the initial conditions. Computational mesh size 400 nm. Stepsize interval 0.01 ms. (D, E) Series of midline slice heat maps at successive time points of recovery (left axis) over the first millisecond of recovery beginning from the initial conditions derived from the Na+ transfer established in (A–C). (D) [Na+] and (E) ionic concentration differences. Scale bars on right: concentration differences (D, E) and membrane potential changes (E) at the passive membrane. Mesh size 400 nm; stepsize interval 0.001 ms.

The latter conditions were confirmed in limiting properties of ephaptic space [Na+] and [Cl] recovery and consequent passive membrane potential changes obtained over the first millisecond of the recovery with parameters in Table 1 at mesh size 400 nm and stepsize interval 0.001 ms. Figure 2D illustrates results of the Nernst-Planck electrodiffusion analysis for the limiting properties of the recovery of ephaptic space [Na+] and [Cl]. Figure 2E portrays the consequent membrane potential changes in the passive membrane following the Na+ transfers produced by the initial ephaptic excitation. These confirm undetectable axial concentration gradients in contrast to time-evolving radial changes. Thus, additional to portraying changes within the first ms the following Na+ transfer, these findings suggest that the recovery process could be simulated assuming uniform concentration profiles for the initial conditions. Such additional simulations indeed yielded demonstrated concordant results (Supplementary Figure S1). This permitted a uniform concentration profile approximation in subsequent computations over longer time intervals to reduce computational time and resource requirements.

4.2 Electrodiffusive recovery in the atrial ephaptic space

Modelling subsequent recovery timecourses accordingly adopted axially uniform initial ephaptic ion concentrations. Neumann boundary conditions of the active and passive membrane were now set to zero reflecting zero flux across these during recovery. Dirichlet boundary conditions set the remaining well stirred extracellular space [Na+] and [Cl] at 153.1 mM and 154.8 mM respectively (Trautwein et al., 1962; Dijkman et al., 1997; Lai et al., 2007). The computation, exemplified in section D in the supplementary file, applied the atrial ephaptic parameters a = 8,000 nm, w = 20 nm (Table 1). It extended over a period covering 99.99975% recovery to resting [Na+] at the ephaptic centre in turn corresponding to a 99.94% ion concentration recovery. Cross-sectional (Figure 3A) heatmaps of [Na+] (Figures 3Aa,b) and ion concentration differences (Figure 3C) demonstrated uniform [Na+] through the radial plane at the outset (Fig. 3Aa). Maps shown following recovery at 80 ms (Figures 3Ab,c) demonstrated extremely small residual differences symmetrical along the radial direction. Midline sectional heatmaps (Figures 3B,C) of [Na+] (Figure 3B) and ion concentration differences (Figure 3C) at critical, 1.25 ms, 10 ms, and 80 ms, time-points during the recovery (a, b, c) demonstrated uniform values along the axial direction, and variations in the radial gradients during recovery.

FIGURE 3
www.frontiersin.org

FIGURE 3. Heat maps of atrial ephaptic recovery. (A) Radial distributions of [Na+] (a, b) and of charged ion concentration differences (c) at the start (a) and end (b, c) of the recovery period. (B, C) Midline slice heat maps of [Na+] (B) and of charged ion concentration differences (C) at critical time points (1.25 ms, 10 ms, and 80 ms) in the atrial ephaptic space during recovery.

Graphical displays at successive logarithmically incremented time points during recovery provided more detailed representations of midline slice [Na+] (Figure 4Aa) and ionic charge concentration differences (Figure 4Ab), the corresponding membrane potential changes (Figure 4Ab) and their spatial profiles. The corresponding quantified graphs showed marked changes both in absolute values and spatial profiles of [Na+] and the corresponding membrane potential profiles with time (Figure 4Ba). The quantification also yielded detailed monotonic recovery timecourses of the [Na+] and membrane potential changes at the ephaptic rim, ephaptic centre and half-way between the two (Figure 4Bb). They also compared the ionic concentration differences with the membrane potential change (Figure 4Bc). Note that the Dirichlet boundary conditions would result in a constant ephaptic rim [Na+] = 153.1 mM and zero membrane potential changes. Within the ephaptic space the recovery followed an approximately exponential trend and was effectively complete within the 80 ms modelling interval.

FIGURE 4
www.frontiersin.org

FIGURE 4. Time and spatial dependences of the atrial ephaptic recovery process. (A) Midline slice heat maps at successive time points of recovery representing (a) [Na+] and (b) ionic concentration differences and membrane potential changes. (B) Quantification of the temporal and spatial recovery. (a) [Na+] spatial profiles with time; Recovery timecourses of (b) [Na+] and membrane potential and (c) ionic concentration differences at the ephaptic rim, ephaptic centre and half-way between the two.

Significant [Cl] changes were contrastingly absent whether at the beginning (Figures 5A,C), in the course of (Figure 5D), or at the end of the recovery intervals (Figures 5B,E), whether in the radial direction (Figures 5A, B) or along axial midline sections (Figures 5C–E). Successive heatmaps systematically representing radial and axial [Cl] distributions through the timecourse of the recovery confirmed these findings (Figure 5F). Use of expanded colour scales focussing on changes in rather than absolute concentration indicated [Cl] changes <1.5 × 10−7 mM, of <0.0000001% (Figure 5G), around otherwise uniformly stable ∼154.8 mM concentrations, through the present and all subsequent computations.

FIGURE 5
www.frontiersin.org

FIGURE 5. Temporal and spatial [Cl] gradients through atrial ephaptic recovery. (A, B) Radial gradients at the outset (A) and 80 ms into (B) the recovery period. (C–E) [Cl] through midline slices taken at successive timepoints (left axis) at the outset (C), 10 ms (D) and 80 ms (E) into the recovery period. (F, G) Sequence of midline slices at successive time points through the recovery period (left axis): (F) uses the same colour scale as (A–E). (G) uses an expanded colour range to permit demonstration of small changes in [Cl] distribution.

Evaluations of the resolution and stability of the computational solutions and their underlying mathematical functions under the chosen simulation conditions examined the effects on these of altered stepsize involving both mesh sizes and stepsize intervals. These were varied from the initially selected 0.01 ms stepsize interval and 400 nm mesh size used above (Supplementary Table S1). First, mesh size was successively varied at a fixed 0.01 ms stepsize interval, extending to a fourfold increased, 1,600 nm, mesh size. Secondly, stepsize interval was varied with fixed 400 nm and then with a substantially increased 1,200 nm mesh size. In the latter case, the range of stepsize intervals was extended to the range 0.005–0.08 ms. All these manoeuvres yielded concordant results, exemplified in results from employing the widest mesh size = 1,200 nm and stepsize interval = 0.08 ms (Supplementary Figure S2).

4.3 Electrodiffusive recovery in the ventricular ephaptic space

Figure 6 displays corresponding midline slice [Na+] (Figure 6Aa), ionic charge concentration differences (Figure 6Ab), corresponding membrane potential changes (Figure 6Ab), and their spatial profiles at successive logarithmically incremented time points during recovery in a ventricular ephaptic space of 12,000 nm radius (Severs, 2000). These similarly yielded absolute values and spatial profiles of spatial [Na+] and membrane potential profiles (Figure 6Ba). The quantification also yielded detailed monotonic recovery timecourses of the [Na+] and membrane potential changes at the ephaptic rim, ephaptic centre and half-way between the two (Figure 6Bb). They also compared the ionic concentration differences with the membrane potential change with time (Figure 6Bc). The recovery trends were similar to that of the atrial ephaptic space but extending over increased 180 ms recovery periods.

FIGURE 6
www.frontiersin.org

FIGURE 6. Time and spatial dependences of the ventricular ephaptic recovery process. (A) Midline slice heat maps at successive time points of recovery representing (a) [Na+] and (b) ionic concentration differences and membrane potential changes. (B) Quantification of the temporal and spatial recovery: (a) [Na+] spatial profiles with time and (b, c) recovery timecourses of (b) [Na+] and membrane potential, and (c) ionic concentration differences at the ephaptic rim, ephaptic centre and half-way between the two, respectively.

4.4 Geometrical and diffusional effects on ephaptic recovery

We next explored the effects of further variations involving the important biological parameters, the ephaptic radius, a, the axial distance, w, and the relevant Na+ and Cl diffusion coefficients, DNa and DCl (Table 2) The computations adopted stepsize all falling within the range validated above. A 0.01 ms stepsize interval and 400 and 800 nm mesh sizes were used when modelling the respective atrial and ventricular parameters, with larger stepsizes when modelling purkinje cell properties. First, varying a above and below reported atrial and ventricular radii modelled above explored the sensitivities of the solutions to decreasing or increasing ephaptic radii representing either varying actual ephaptic radii or the effects of membrane folding in the intercalated disk. The latter extended to ephaptic junctions between purkinje cells completing the data set for different cardiac cell types. The computations completed a data set of results of halving the atrial, doubling the ventricular ephaptic space and finally including the purkinje ephaptic space. This completed a range of radii, a, of 4,000 nm, 8,000 nm, 12,000 nm, 24,000 nm, and 40,000 nm respectively. In modelling the purkinje ephaptic space 0.4 and 0.8 ms stepsize intervals yielded indistinguishable results. Figure 7 demonstrates times required to return to 99.99975% of resting [Na+] and 99.94% and 99.95% at the ephaptic centre respectively. These were, in order of increasing ephaptic radius, 20, 80, 180, 720 and 2000 ms. They fell along a function suggesting that recovery time was directly proportional to the square of the ephaptic radius (a2).

TABLE 2
www.frontiersin.org

TABLE 2. Effects of ephaptic geometry and diffusion properties on ephaptic recovery times.

FIGURE 7
www.frontiersin.org

FIGURE 7. Quantification of temporal and spatial ephaptic recovery at varying a and constant w, DNa and DCl. (A) [Na+] spatial profiles with time; (B) Recovery timecourses of [Na+] and membrane potential and ionic concentration differences at the ephaptic rim, ephaptic centre and half-way between the two. Values of a varied through a = 4000 (a), 8000 (b), 12000 (c), 24000 (d) and 40000 nm (e) respectively.

Secondly, previous reports had suggested possible functional roles of a fixed w term resulting from mechanical couplings between opposing Nav1.5β subunits in the respective active and passive membranes (Veeraraghavan et al., 2018; Salvage et al., 2020a). We accordingly modelled atrial (Supplementary Figure S3) and ventricular (Supplementary Figure S4) junctions with axial distances of w = 5, 10 nm and 40 nm respectively, comparing results with those from the control axial distance w = 20 nm. Supplementary Figure S3 demonstrates absolute [Na+] varying with w, reflecting the constant initial molar Na+ transfer from the ephaptic volume. For example, doubling the axial distance halved the ionic concentration differences at any defined spatio-temporal point. Nevertheless, trends in the spatial and temporal [Na+] recovery profiles were proportionally similar, and the corresponding absolute voltage recovery profiles were indistinguishable between different values of w. Hence ephaptic recovery is not greatly affected by the changes in axial distance explored.

Thirdly, varying the diffusion coefficients could assess for the effects on any restrictions to diffusion attributable to the presence of proteins in the ephaptic space, and non-uniformities in w that might arise from regions where the component ephaptic membranes are more closely associated as at gap junction plaques. We explored the effects of halving and doubling DNa and DCl in the presence of the adopted atrial (Supplementary Figure S5) and ventricular (Supplementary Figure S6) ephaptic space parameters of a and w. Supplementary Figure S5,6 demonstrate similar initial conditions; these were followed by [Na+] and voltage recoveries. However, their time periods varied through values of 160, 80 and 40 ms in the atria and 360, 180 and 90 ms in the ventricle for the runs with halved, normal and doubled DNa and DCl. Hence the time-period for recovery is affected by the relative value of the diffusion coefficients with a doubling of these halving the recovery time-period for recovery. Furthermore, equivalent time points showed similar patterns of radial recovery. For example, 1.25 ms within an 80 ms recovery period gave a similar recovery as at 2.5 ms out of an 160 ms recovery period.

Table 2 summarises the results from the simulations completed here, providing recovery time data for comparison with reported cardiac cycle lengths made in the Discussion. They predicted ephaptic recovery times of order 40 and 180 ms in atrial and ventricular cardiomyocytes. High limits from doubling ventricular radius or halving diffusion coefficient gave respective 720 and 360 ms recovery times. Finally, the larger diameter purkinje fibres gave recovery times of order 2000 ms.

4.5 Dependences of recovery half times upon a, w and DNa and DCl

Figure 8 summarises these effects in terms of the dependences of recovery half times, t1/2, upon a, w and DNa and DCl respectively. The computations were consistently performed with a 800 nm mesh size and stepsize intervals of 0.01 and 0.02 ms for the atrial and ventricular ephaptic spaces respectively. Results were in agreement with those using 1,200 nm and 0.08 ms mesh sizes used to study the purkinje cell ephaptic space. This first demonstrated relationships between t1/2, upon a (Figure 8A) that yielded a linear t1/2, - a2 relationship (Figure 8B). Secondly, t1/2 showed little variation with values of w on the millisecond scale (Figure 8C), although there were sub-millisecond effects of increasing w that appeared to plateau around w > 20 nm (Figures 8E,F). Finally, the linear double logarithmic plots of t1/2 against the relative diffusion coefficients (Figure 8D) yielded lines of best fit of the form:

logt1/2=mlogrelativediffusioncoefficient+c(33)

FIGURE 8
www.frontiersin.org

FIGURE 8. Dependences of recovery half times upon a, w and DNa and DCl. Relationships between recovery half times and: (A, B) ephaptic radius, a (A), and a2 (B), (C) atrial (black lines) and ventricular (red lines) axial distance, w, and (D) in double logarithmic plots, relative diffusion coefficients. (E, F) Small trends in atrial (E) and ventricular recovery half times (F) with axial distance, w, following magnification of the ordinate.

For the atrial data, determination of the constants m and c, gave:

logt1/2=1.15461.1013×logrelativediffusioncoefficient(34)

For which,

t1/2=14.2758×relativediffusioncoefficient1.1013(35)

For the ventricular data,

logt1/2=1.56651.1220×logrelativediffusioncoefficient(36)

For which,

t1/2=36.8553×relativediffusioncoefficient1.1220(37)

5 Discussion

5.1 Gap junction and ephaptic conduction of the cardiac action potential

Conduction of cardiac myocardial excitation involves action potential (AP) generation in active cardiomyocytes electrotonically depolarising membranes of adjacent, initially quiescent, cells within the syncytium, thereby activating their Nav1.5 Na+ channels. Additional to gap junction mediated electrotonic current flow (Rohr, 2004), increasing evidence implicates ephaptic excitation in activation of such conduction (Sperelakis and Mann, 1977; Carmeliet, 2019). Ephaptic conduction has a widespread occurrence; specific examples transmit either bidirectional or unidirectional cell-cell excitation (Furshpan and Potter, 1959) alone (Watanabe and Grundfest, 1961; Bennett, 1977) or combined with other electrical or chemical mechanisms (Rash et al., 1996). In the heart, it could occur at specialised perinexal regions separating adjacent cardiomyocytes (Rhett and Gourdie, 2012; Rhett et al., 2013; Veeraraghavan et al., 2014). It could either complement, or in certain pathological circumstances replace, conduction through well-established local circuit, gap junction (GJ) current pathways (Rohr, 2004) involving ventricular connexin, Cx43 and atrial, Cx43 and Cx40 (Verheule et al., 1997). Syncytial myocardial conduction may then persist despite compromised gap junction coupling following clinical ischaemic insult (Rohr, 2004), atrial fibrillation (Chaldoupi et al., 2009) and in experimental loss-of-function, Cx43 (Eloff et al., 2001; Gutstein et al., 2001) or Cx40, genetic platforms (Hagendorff et al., 1999; Bagwe et al., 2005).

5.2 Ephaptic activation of conduction between adjacent cardiomyocytes

Existence of ephaptic activation of conduction between cardiomyocytes has both experimental and theoretical support (Léonetti, 1998; Koch, 2004; Moise et al., 2021). Electron microscopy of perinexal region specialisations demonstrated closely apposed adjacent component cell membranes (Veeraraghavan et al., 2015) with increased Nav1.5 in addition to Cx43 expression. Super-resolution microscopy demonstrated dense Nav1.5 clusters within 200 nm of GJ plaques (Veeraraghavan et al., 2015; Hichri et al., 2018; Salvage et al., 2020b). Smart patch clamping indicated elevated junctional over nonjunctional INa (Veeraraghavan et al., 2018). Elevations in perinexal β1 (SCN1B) expression, co-localised with that of the Nav1.5 potentially providing adhesion trans scaffolds ensuring narrowed (∼20 nm) perinexal clefts adjacent to the gap junctions (Veeraraghavan et al., 2018; Salvage et al., 2020a). βadp1 peptide inhibited this β1-mediated adhesion, widened guineapig ventricular perinexi, reduced perinexal but not whole cell INa and pro-arrhythmically slowed AP conduction (Veeraraghavan et al., 2018).

Mathematical modelling of linear strands of cardiac cells under different membrane geometrical, extracellular space, ionic concentration and diffusion conditions suggested that AP propagation can persist even with compromised GJ coupling (Sperelakis and Mann, 1977; Sperelakis, 2002; Mori et al., 2008). They also raised the possibility of alternating ephaptic and gap-junction-mediated activation (Mori et al., 2008). Such features extended to 3-dimensional myocyte sheets incorporating functioning GJs and more realistic inhomogeneities in extracellular spacing and ionic channel distribution (Lin and Keener, 2010). Recent finite element replications integrated experimentally determined nanoscale transmission electron microscopy intercalated disk structure with measurements of intercellular cleft electrical conductivity. Intermembrane separation and gap junction distribution and their heterogeneity then significantly affected spatial characteristics of electrical polarization within the intercellular cleft modifying the postjunctional Nav1.5 activation (Moise et al., 2021).

5.3 Ephaptic recovery following conduction

Fewer studies have examined recovery processes following ephaptic excitation, essential for any subsequent excitation as in repetitive excitation. The present formalised computational analysis explored limiting electrodiffusive contributions within the ephaptic space, and between ephaptic space and remaining extracellular space to ephaptic membrane recovery for the first time. It then explores effects upon these of ephaptic junction geometrical and diffusional properties. It extends a charge difference approach utilising Gauss’s flux theorem previously introduced to determine membrane voltages across and within arbitrary numbers of compartments separated by defined membrane capacitances (Fraser and Huang, 2004; 2007; Fraser et al., 2011; Pedersen et al., 2011). These are here exemplified by the intracellular compartments of the participating cardiomyocytes, the restricted ephaptic, and the remaining free, extracellular space. This approach determined membrane voltage from concentration differences resulting from recovery fluxes of finitely mobile charge within a restricted ephaptic space. It is particularly useful analysing an ephaptic recovery involving electrochemically driven ion fluxes producing ion concentration and consequent potential changes dependent upon both compartment volumes and their intervening membrane capacitances. Previous computational modelling contrastingly involved summing of ionic current as opposed to ion concentration components (Hodgkin and Huxley, 1952; McAllister et al., 1975; Huang and Peachey, 1992; Luo and Rudy, 1994; Puglisi and Bers, 2001). It provides changes in activation transmembrane potential differences between compartments (Sperelakis, 2002; Mori et al., 2008; Lin and Keener, 2010; Veeraraghavan et al., 2015) but not absolute potentials within compartments. This requires at least one compartment to provide a reference potential, and frequently cannot explicitly simulate concentration changes.

5.4 Modelling electrodiffusive processes underlying ephaptic recovery

5.4.1 Initiation of electrodiffusive modelling

The present analysis first adopted experimentally established experimental values describing mammalian atrial, ventricular and purkinje cardiomyocyte ephaptic junction anatomy, extracellular, Na+ and Cl, electrolyte concentrations and their free diffusion coefficients, DNa and DCl (Passiniemi, 1983) and AP upstroke magnitudes (Schoenberg et al., 1975; Severs, 2000; Mori et al., 2008; Richards et al., 2011; Veeraraghavan et al., 2015; 2018). Secondly, computational analysis was initialised by determining a lower limit of transmembrane Na+ withdrawal by the AP of the initially active cell by a ΔV∼130 mV AP upstroke at physiological (37°C) temperatures. This would generate domains of reduced [Na+] within the restricted diffusional perinexal ephaptic space. This calculation likely provides a lower limit value as it assumes a simple circular geometry of unfolded contiguous membranes of equal radial diameter a. Available quantitative information bearing on possible transmembrane fluxes within the ephaptic space does not provide precise experimental values for Nav1.5 channel numbers within the ephaptic space itself. Combining reported peak cardiomyocyte INa = 20–140 pA pF-1 ranges (Casini et al., 2009; Rougier et al., 2019) and a ∼120 pF total (guinea-pig) myocyte membrane capacitance (Wan et al., 2003), yields a ∼2,400-16,800 pA range per myocyte. This suggests a typical total 1,200-8,400 Nav channel number per myocyte assuming reported ∼2 pA unit currents reported in neuroblastoma and SCN5A-transfected HEK293 cells (Aldrich and Stevens, 1987; Clatot et al., 2017). Some reports localise ∼50% of cardiomyocyte NaV channels to intercalated disc plasma membranes (Petitprez et al., 2011). This suggests between 600-4,200 active Nav1.5 in each membrane flanking the ephaptic space. A value of 2,400 channels in the middle of this range for atrial cardiomyocytes gives a 4,800 pA = 4,800 × 10−12 C s-1 current at the intercalated disk. The resulting 4.9748×1014mol s1 molar Na+ transfer rate predicts a 3.4824 × 10−17 mol Na+ removal from the ephaptic space during the action potential upstroke assuming ∼0.7 ms mean Na+ channel open times (Sigworth and Neher, 1980). This exceeds the limiting values adopted here.

5.4.2 Modelling the electrodiffusive recovery

Thirdly, finite element computation utilising three-dimensional Nernst-Planck electrostatic equations determined the extracellular Na+ and Cl electrodiffusional fluxes. This yielded the time-dependence of axial and radial ionic concentrations within the three-dimensional ephaptic space both immediately following the AP and its subsequent recovery (Léonetti, 1998; Koch, 2004). Computations of the timecourse and spatial distribution of the resulting [Na+] changes adopted a bilaterally and radially symmetrical atrial ephaptic space flanked by uniformly contiguous directly opposing membranes of equal area and specific capacitance, corresponding to equal total capacitances (Schoenberg et al., 1975; Severs, 2000; Mori et al., 2008; Richards et al., 2011). Immediately following Na+ transfer, radial [Na+] gradients were absent in [Na+] heat maps. Axial gradients in [Na+] and charged ion concentration differences were detectable only at high colour magnifications. Such conditions were confirmed in computations obtained over the first millisecond of recovery.

Fourthly, the consequent membrane potential changes induced by the Na+ withdrawal were then approximated applying Gauss’s Flux theorem. The majority of the electric flux, ΦÊ resulting from the resulting charging of the ephaptic space would traverse its bounding membranes. Thus, these two closely apposed, active and passive, membranes were separated by a small gap distance w << a. Hence membrane areas accounted for most of the gaussian surface surrounding the ephaptic space in comparison to its rim accessing the extracellular compartment: πa2>>2πad. The simulations then sought to determine voltage change V in the passive membrane and its subsequent recovery, produced by AP activation in the active membrane of the ephaptic junction.

For an ephaptic space of uniform radial extent comprising closely apposed participating membranes each of equal area or total capacitance, the latter membrane would fire an AP were V to exceed the threshold amplitude, V>ΔVth,. More general treatments for the latter firing condition could consider passive and active opposing membranes of differing surface areas, Sth and SAP respectively, and similar Cm. An action potential, amplitude ΔVAP, would now transfer from the ephaptic space, total charge,

q=SAPCmVAP(38)

In the limit corresponding to small gap distance w→ 0, the two membranes will show voltage change following the excitation condition, illustrating the importance of the relative areas of the active and passive membranes, and their implications for a capacity for uni- or bi-directional conduction:

ΔV=SAPVAPSAP+Sth>Vth(39)

Characterisation of the subsequent recovery phase simulated Na+ diffusion from the remaining extracellular space into the radially symmetrical atrial cardiomyocyte ephaptic space. Communication of ephaptic space with a well stirred extracellular fluid was reflected in Dirichlet boundary conditions holding constant the ephaptic rim to established, 153.1 mM [Na+] and 154.8 mM [Cl], extracellular values (Trautwein et al., 1962; Dijkman et al., 1997; Lai et al., 2007). The resulting cross sectional and midline sectional heatmaps confirmed uniform, and close to uniform radial [Na+] and ion concentration differences at the outset and following recovery. There were uniform axial values with variations in the radial gradients at critical time-points during recovery. Systematic displays and quantified graphs of the recovery timecourses demonstrated marked changes in absolute values and spatial profiles of midline slice [Na+], ionic concentration differences and the corresponding membrane potential profiles with time. Their recovery timecourses at the ephaptic rim, ephaptic centre and half-way between the two assumed approximately exponential trends to their final recovered values. In contrast, there were no significant changes in the counterion [Cl] levels through the simulations. Further computations successively varying the parameters of mesh size and stepsize interval holding the remaining parameter constant corroborated these analyses. They demonstrated that all the computations fell within a stepsize parameter space producing consistent and convergent results. Finally, calculations extended to the ventricular ephaptic space demonstrated similar recovery trends that took place over longer recovery periods.

5.4.3 Effects of ephaptic junction anatomy and diffusive properties

These estimations began with considering a simplified and formalised system. However, in vivo ephaptic junctions are complex 3-dimensional structures comprising plicate and inter-plicate regions resulting in varying rather than constant ephaptic axial distances (Veeraraghavan et al., 2015; 2018; Vermij et al., 2017). These as well as diffusion restrictions reflecting contained proteins and nonuniformities in ephaptic thickness at gap junction plaques could modify recovery. These complexities would all tend to predict higher initial overall Na+ fluxes attributed to AP activation by the active membrane, higher effective radial diameters and lower effective axial distances and Na+ and Cl diffusion coefficients within the ephaptic gaps. All these individually and together would result in more prolonged recovery times than predicted by a formal model, which therefore provide lower limits for such recovery times. Nevertheless, the extent to which these variations might influence the findings was assessed by exploring the effects of systematically increasing and decreasing values of cardiomyocyte radii, axial distances, and adopted Na+ and Cl diffusion coefficient values in turn, whilst the remaining parameters were held constant. These additionally provided useful quantitative insights into the physical consequences of specific properties in the biological systems. These could prompt and guide further detailed experimental study. These assessments first explored the effects of varying ephaptic radius a. They thereby assessed the sensitivities of the solutions to ephaptic radius, provided a formal representation of membrane folding in the intercalated disk, and extended the computations to the significantly larger diameter purkinje cells completing the data set for different cardiac cell types. The findings suggested that recovery half-lives were directly proportional to the square of the ephaptic radius. Contrasting results resulted from varying the axial distances w. Despite significant effects on absolute [Na+] reflecting the constant initial molar Na+ transfer from the ephaptic volume, this exerted little significant effect upon voltage recovery profiles. Finally, halving and doubling DNa and DCl markedly affected recovery half-lives through double logarithmic relationships.

5.5 Physiological significance of the computed properties

Finally, the resulting spatial and temporal electrodiffusion properties of the modelled recovery could be compared with known physiological heart rates, in comparing its contributions to those of the other biomolecules associated with such clusters. Collation of recovery times from these investigations provided limiting indications of the extent to which electrodiffusion processes might contribute to ephaptic recovery in different cardiomyocyte types, pacing conditions, and animal species. The initial modelling conditions predicted ephaptic recovery times of order 40 and 180 ms in atrial and ventricular cardiomyocytes. High limits from doubling ventricular radius or halving diffusion coefficient gave respective 720 and 360 ms recovery times. The larger diameter purkinje fibres gave recovery times of order 2000 ms. All these values could be compared with resting and exercising heart rates amongst mammalian species. For example, resting heart rates of blue whale ∼30; elephant, ∼30; human, ∼70; mouse, ∼650; shrew, ∼835, Etruscan shrew, ∼1,511 bpm, correspond to cardiac cycle durations of ∼2000, ∼2000, ∼850, ∼90, ∼70 and ∼40 ms respectively (Benedict and Lee, 1936; Fons et al., 1997; Ho et al., 2011; Goldbogen et al., 2019; Ponganis, 2021). Electrodiffusion could thus potentially either entirely, or partially but significantly, account for such recovery in atrial and ventricular cardiomyocytes respectively, to extents dependent on the species heart rates.

In contrast, purkinje cardiomyocyte recovery times were clearly incompatible with ephaptic transmission in an absence of mechanisms additional to electrodiffusional processes to aid [Na+] recovery within the ephaptic space. Thus, the increased purkinje cell diameter whilst enhancing conduction velocity by increased local circuit gap junction current could compromise ephaptic conduction by reducing its electrodiffusive recovery. Thus, cable theory predicts that conduction velocity attributable to gap junction mediated local circuit currents (King et al., 2013) as opposed to ephaptic conductance contributions depends on membrane capacitance, and cell-cell intracellular resistance, of unit fibre length (see Figure 5 of (Edling et al., 2019)) and that the latter varies with square root of the fibre radius (Huang, 2021). Hence in larger diameter cardiomyocytes, in contrast to their prolonged recovery times compromising ephaptic conduction contributions, the resulting increased gap junction conduction positively contributes to action potential propagation velocity. The latter is reflected in the higher conduction velocities shown by larger diameter Purkinje cardiomyocytes. This would contrast with the larger expected relative ephaptic contributions to conduction velocity in smaller diameter atrial fibres.

6 Conclusion

In conclusion, the present analysis sought fundamental physical insights into ephaptic conduction and its recovery. It combined electrochemical and Gaussian analyses of diffusive and transmembrane voltage events, in contrast to previous approaches involving circuit modelling solely of electric currents, for a formalised ephaptic space. It initially adopted limiting ephaptic radial diameters, axial distances and Na+ and Cl diffusion coefficients that likely provide lower limiting recovery times. However, it then explored the effects of more complex ephaptic junction organization through systematically varying each of these parameters, to provide more broadly applicable physical insights. These insights could prompt further experimental and theoretical studies scrutinizing the applicability of electrodiffusive recovery to cardiac ephaptic transmission in vivo (Struckman et al., 2021). The latter could first determine effects on effective ephaptic radii and diffusion coefficients of membrane plication, nonuniformities in axial distance, diffusive obstructions and local factors affecting diffusion coefficients. One might also speculate that contractile activity could exert mechanical effects upon ephaptic geometry contributing recovery through hydrostatic effects additional to the electrodiffusion processes replacing Na+ from the bulk extracellular space analysed here (Salvage et al., 2020a).

Second, the relative contributions of electrodiffusive and actions of other membrane molecules potentially making up a functional ephaptome merit determination. The Nav1.5 occur in clusters within perinexal regions and surround gap junction plaques with large regions of gap junction clusters (Rhett and Gourdie, 2012; Rhett et al., 2013; Gourdie, 2019) ‘pinning’ together the component membranes at discrete sites (Salvage et al., 2020a). These are separated by an ephaptic space permitting diffusion around them. Besides electrodiffusion, ephaptic [Na+] or positive charge could also be restored by transport pathways extruding intracellular Na+ or other positive charge into the ephaptic space. Confocal microscopy and super resolution studies indicated increased Na+-K+-ATPase densities (Noël et al., 1991; Mcdonough et al., 1996) particularly around gap and adherence junctions (Struckman et al., 2021), particularly in atrial relative to ventricular myocytes. Cycles of electrogenic intercalated disk Na+-K+-ATPase activity each would add 3Na+ and withdraw 2K+ from the ephaptic space (Vermij et al., 2017).

In addition, gated stimulated emission depletion (gSTED) and stochastic optical reconstruction microscopy (STORM) super-resolution microscopy attribute much perinexal K+ conductance to Kir2.1 (Difranco et al., 2015; Chen et al., 2018) clustered around desmosomal regions (Veeraraghavan et al., 2016). Their strongly inwardly rectifying property would close depolarised channels, reducing outward current during AP depolarisation, minimising electrical shunting of the Na+ transfer. Contrastingly, they would permit positive K+ flux into the recovering ephaptic space across repolarised membranes. Thus, acute interstitial oedema-induced perinexal swelling causes a pro-arrhythmic slowed conduction respectively relieved and accentuated by Kir2.1 inhibition and Nav1.5 block (Veeraraghavan et al., 2016).

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

JAM: Data curation, Investigation, Validation, Writing–original draft, Methodology, Software. OJB: Data curation, Investigation, Methodology, Software, Writing–original draft, Visualization. SS: Conceptualization, Formal Analysis, Funding acquisition, Supervision, Writing–review and editing. APJ: Conceptualization, Formal Analysis, Funding acquisition, Writing–review and editing, Resources, Visualization. HRM: Conceptualization, Investigation, Project administration, Software, Supervision, Validation, Writing–original draft. CL-HH: Conceptualization, Investigation, Project administration, Supervision, Validation, Writing–original draft, Data curation, Formal Analysis, Funding acquisition, Writing–review and editing, Visualization.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. We thank the British Heart Foundation (PG/14/79/31102 and PG/19/59/34582 and Cambridge Centre for Research Excellence) for its generous support.

Conflict of interest

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

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

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

References

Aldrich R. W., Stevens C. F. (1987). Voltage-dependent gating of single sodium channels from mammalian neuroblastoma cells. J. Neurosci. 7, 418–431. doi:10.1523/JNEUROSCI.07-02-00418.1987

PubMed Abstract | CrossRef Full Text | Google Scholar

Bagwe S., Berenfeld O., Vaidya D., Morley G. E., Jalife J. (2005). Altered right atrial excitation and propagation in connexin40 knockout mice. Circulation 112, 2245–2253. doi:10.1161/CIRCULATIONAHA.104.527325

PubMed Abstract | CrossRef Full Text | Google Scholar

Bardsley O. J., Matthews H. R., Huang C. L. H. (2021). Finite element analysis predicts Ca2+ microdomains within tubular-sarcoplasmic reticular junctions of amphibian skeletal muscle. Sci. Rep. 11, 14376. doi:10.1038/s41598-021-93083-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Benedict F. G., Lee R. C. (1936). The heart rate of the elephant. Proc. Am. Philos. Soc. 76, 335–341. Available at: http://www.jstor.org/stable/984548.

Google Scholar

Bennett M. (1977). “Electrical transmission: a functional analysis and comparison to chemical transmission,” in Handbook of Physiology. Editor E. Kandel (Bethesda, MD: American Physiological Society), 357–416.

CrossRef Full Text | Google Scholar

Carmeliet E. (2019). Conduction in cardiac tissue. Historical reflections. Physiol. Rep. 7, e13860. doi:10.14814/phy2.13860

PubMed Abstract | CrossRef Full Text | Google Scholar

Casini S., Verkerk A. O., van Borren M. M., van Ginneken A. C., Veldkamp M. W., de Bakker J. M., et al. (2009). Intracellular calcium modulation of voltage-gated sodium channels in ventricular myocytes. Cardiovasc Res. 81, 72–81. doi:10.1093/cvr/cvn274

PubMed Abstract | CrossRef Full Text | Google Scholar

Chaldoupi S. M., Loh P., Hauer R. N. W., De Bakker J. M. T., Van Rijen H. V. M. (2009). The role of connexin40 in atrial fibrillation. Cardiovasc. Res. 84, 15–23. doi:10.1093/cvr/cvp203

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen K., Zuo D., Liu Z., Chen H. (2018). Kir2.1 channels set two levels of resting membrane potential with inward rectification. Pflugers Arch. Eur. J. Physiol. 470, 599–611. doi:10.1007/s00424-017-2099-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Clatot J., Hoshi M., Wan X., Liu H., Jain A., Shinlapawittayatorn K., et al. (2017). Voltage-gated sodium channels assemble and gate as dimers. Nat. Commun. 8, 2077. doi:10.1038/s41467-017-02262-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Difranco M., Yu C., Quiñonez M., Vergara J. L. (2015). Inward rectifier potassium currents in mammalian skeletal muscle fibres. J. Physiol. 593, 1213–1238. doi:10.1113/jphysiol.2014.283648

PubMed Abstract | CrossRef Full Text | Google Scholar

Dijkman M. A., Heslinga J. W., Allaart C. P., Sipkema P., Westerhof N. (1997). Reoxygenated effluent of Tyrode-perfused heart affects papillary muscle contraction independent of cardiac perfusion. Cardiovasc. Res. 33, 45–53. doi:10.1016/S0008-6363(96)00173-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Draper M. H., Weidmann S. (1951). Cardiac resting and action potentials recorded with an intracellular electrode. J. Physiol. 115, 74–94. doi:10.1113/jphysiol.1951.sp004653

PubMed Abstract | CrossRef Full Text | Google Scholar

Edling C. E., Fazmin I. T., Saadeh K., Chadda K. R., Ahmad S., Valli H., et al. (2019). Molecular basis of arrhythmic substrate in ageing murine peroxisome proliferator-activated receptor γ co-activator deficient hearts modelling mitochondrial dysfunction. Biosci. Rep. 39, BSR20190403. doi:10.1042/BSR20190403

PubMed Abstract | CrossRef Full Text | Google Scholar

Eloff B. C., Lerner D. L., Yamada K. a., Schuessler R. B., Saffitz J. E., Rosenbaum D. S. (2001). High resolution optical mapping reveals conduction slowing in connexin43 deficient mice. Cardiovasc. Res. 51, 681–690. doi:10.1016/s0008-6363(01)00341-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Fons R., Sender S., Peters T., Jürgens K. D. (1997). Rates of rewarming, heart and respiratory rates and their significance for oxygen transport during arousal from torpor in the smallest mammal, the Etruscan shrew Suncus etruscus. J. Exp. Biol. 200, 1451–1458. doi:10.1242/jeb.200.10.1451

PubMed Abstract | CrossRef Full Text | Google Scholar

Fraser J. A., Huang C. L.-H. (2004). A quantitative analysis of cell volume and resting potential determination and regulation in excitable cells. J. Physiol. 559, 459–478. doi:10.1113/jphysiol.2004.065706

PubMed Abstract | CrossRef Full Text | Google Scholar

Fraser J. A., Huang C. L.-H. (2007). Quantitative techniques for steady-state calculation and dynamic integrated modelling of membrane potential and intracellular ion concentrations. Prog. Biophys. Mol. Biol. 94, 336–372. doi:10.1016/j.pbiomolbio.2006.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Fraser J. A., Huang C. L.-H., Pedersen T. H. (2011). Relationships between resting conductances, excitability, and t-system ionic homeostasis in skeletal muscle. J. Gen. Physiol. 138, 95–116. doi:10.1085/jgp.201110617

PubMed Abstract | CrossRef Full Text | Google Scholar

Furshpan E. J., Potter D. D. (1959). Transmission at the giant motor synapses of the crayfish. J. Physiol. 145, 289–325. doi:10.1113/jphysiol.1959.sp006143

PubMed Abstract | CrossRef Full Text | Google Scholar

Gabriel S., Lau R. W., Gabriel C. (1996). The dielectric properties of biological tissues: III. Parametric models for the dielectric spectrum of tissues. Phys. Med. Biol. 41, 2271–2293. doi:10.1088/0031-9155/41/11/003

PubMed Abstract | CrossRef Full Text | Google Scholar

Gemel J., Levy A. E., Simon A. R., Bennett K. B., Ai X., Akhter S., et al. (2014). Connexin40 abnormalities and atrial fibrillation in the human heart. Curr. Ther. Res. - Clin. Exp. 76, 159–168. doi:10.1016/j.yjmcc.2014.08.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldbogen J. A., Cade D. E., Calambokidis J., Czapanskiy M. F., Fahlbusch J., Friedlaender A. S., et al. (2019). Extreme bradycardia and tachycardia in the world’s largest animal. Proc. Natl. Acad. Sci. U. S. A. 116, 25329–25332. doi:10.1073/pnas.1914273116

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldman D. E. (1943). Potential, impedance and rectification in membranes. J. Gen. Physiol. 27, 37–60. doi:10.1085/jgp.27.1.37

PubMed Abstract | CrossRef Full Text | Google Scholar

Gourdie R. G. (2019). The cardiac gap junction has discrete functions in electrotonic and ephaptic coupling. Anat. Rec. Hob. 302, 93–100. doi:10.1002/ar.24036

CrossRef Full Text | Google Scholar

Gutstein D. E., Morley G. E., Tamaddon H., Vaidya D., Schneider M. D., Chen J., et al. (2001). Conduction slowing and sudden arrhythmic death in mice with cardiac-restricted inactivation of connexin43. Circ. Res. 88, 333–339. doi:10.1161/01.res.88.3.333

PubMed Abstract | CrossRef Full Text | Google Scholar

Hagendorff A., Schumacher B., Kirchhoff S., Lüderitz B., Willecke K. (1999). Conduction disturbances and increased atrial vulnerability in connexin40-deficient mice analyzed by transesophageal stimulation. Circulation 99, 1508–1515. doi:10.1161/01.CIR.99.11.1508

PubMed Abstract | CrossRef Full Text | Google Scholar

Hichri E., Abriel H., Kucera J. P. (2018). Distribution of cardiac sodium channels in clusters potentiates ephaptic interactions in the intercalated disc. J. Physiol. 596, 563–589. doi:10.1113/JP275351

PubMed Abstract | CrossRef Full Text | Google Scholar

Hilgemann D. W., Noble D. (1987). Excitation-contraction coupling and extracellular calcium transients in rabbit atrium: reconstruction of basic cellular mechanisms. Proc. R. Soc. Lond. Ser. B. Biol. Sci. 230, 163–205. doi:10.1098/rspb.1987.0015

PubMed Abstract | CrossRef Full Text | Google Scholar

Ho D., Zhao X., Gao S., Hong C., Vatner D. E., Vatner S. F. (2011). Heart rate and electrocardiography monitoring in mice. Curr. Protoc. Mouse Biol. 1, 123–139. doi:10.1002/9780470942390.mo100159

PubMed Abstract | CrossRef Full Text | Google Scholar

Hodgkin A. L., Huxley A. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117, 500–544. doi:10.1113/jphysiol.1952.sp004764

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang C. (2021). Keynes and aidley’s nerve and muscle. 5th ed. Cambridge: Cambridge University Press.

Google Scholar

Huang C. L., Peachey L. D. (1992). A reconstruction of charge movement during the action potential in frog skeletal muscle. Biophys. J. 61, 1133–1146. doi:10.1016/S0006-3495(92)81923-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Hutchison J. S., Ward R. E., Lacroix J., Hébert P. C., Barnes M. A., Bohn D. J., et al. (2008). Hypothermia therapy after traumatic brain injury in children. N. Engl. J. Med. 358, 2447–2456. doi:10.1056/NEJMoa0706930

PubMed Abstract | CrossRef Full Text | Google Scholar

King J. H., Huang C. L.-H., Fraser J. A. (2013). Determinants of myocardial conduction velocity: implications for arrhythmogenesis. Front. Physiol. 4, 154. doi:10.3389/fphys.2013.00154

PubMed Abstract | CrossRef Full Text | Google Scholar

Kirchhoff S., Kim J. S., Hagendorff A., Thönnissen E., Krüger O., Lamers W. H., et al. (2000). Abnormal cardiac conduction and morphogenesis in connexin40 and connexin43 double-deficient mice. Circ. Res. 87, 399–405. doi:10.1161/01.RES.87.5.399

PubMed Abstract | CrossRef Full Text | Google Scholar

Kitamura H., Ohnishi Y., Yoshida A., Okajima K., Azumi H., Ishida A., et al. (2002). Heterogeneous loss of connexin43 protein in nonischemic dilated cardiomyopathy with ventricular tachycardia. J. Cardiovasc. Electrophysiol. 13, 865–870. doi:10.1046/j.1540-8167.2002.00865.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Koch C. (2004). Biophysics of Computation: information processing in single neurons. London: Oxford University Press.

Google Scholar

Kostin S., Dammer S., Hein S., Klovekorn W. P., Bauer E. P., Schaper J. (2004). Connexin 43 expression and distribution in compensated and decompensated cardiac hypertrophy in patients with aortic stenosis. Cardiovasc. Res. 62, 426–436. doi:10.1016/j.cardiores.2003.12.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Kostin S., Rieger M., Dammer S., Hein S., Richter M., Klövekorn W. P., et al. (2003). Gap junction remodeling and altered connexin43 expression in the failing human heart. Mol. Cell. Biochem. 242, 135–144. doi:10.1023/A:1021154115673

PubMed Abstract | CrossRef Full Text | Google Scholar

Kramlich A., Bohnert J., Dössel O. (2012). “Transmembrane voltages caused by magnetic fields – numerical study of schematic cell models,” in, 337–342.

CrossRef Full Text | Google Scholar

Lai Y.-J., Chen Y.-Y., Cheng C.-P., Lin J. J.-C., Chudorodova S. L., Roshchevskaya I. M., et al. (2007). Changes in ionic currents and reduced conduction velocity in hypertrophied ventricular myocardium of Xin alpha-deficient mice. Anadolu Kardiyol. Derg. 7, 90–92.

PubMed Abstract | Google Scholar

Lambiase P. D., Tinker A. (2015). Connexins in the heart. Cell Tissue Res. 360, 675–684. doi:10.1007/s00441-014-2020-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Léonetti M. (1998). On biomembrane electrodiffusive models. Eur. Phys. J. B 2, 325–340. doi:10.1007/s100510050256

CrossRef Full Text | Google Scholar

Lide D. R. (1993). CRC handbook of chemistry and physics: a ready-reference book of chemical and phyical data/editor-in-chief, David R. Lide; associate editor, H.P.R. Frederikse. 76th ed. Boca Raton, Fla, London: CRC Press.

Google Scholar

Lin J., Keener J. P. (2010). Modeling electrical activity of myocardial cells incorporating the effects of ephaptic coupling. Proc. Natl. Acad. Sci. U. S. A. 107, 20935–20940. doi:10.1073/pnas.1010154107

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo C. H., Rudy Y. (1994). A dynamic model of the cardiac ventricular action potential: II. Afterdepolarizations, triggered activity, and potentiation. Circ. Res. 74, 1097–1113. doi:10.1161/01.RES.74.6.1097

PubMed Abstract | CrossRef Full Text | Google Scholar

Malmberg C. G., Maryott A. A. (1956). Dielectric constant of water from 0 to 100 C. J. Res. Natl. Bur. Stand. 56, 1. doi:10.6028/jres.056.001

CrossRef Full Text | Google Scholar

McAllister R. E., Noble D., Tsien R. W. (1975). Reconstruction of the electrical activity of cardiac Purkinje fibres. J. Physiol. 251, 1–59. doi:10.1113/jphysiol.1975.sp011080

PubMed Abstract | CrossRef Full Text | Google Scholar

Mcdonough A. A., Zhang Y., Shin V., Frank J. S. (1996). Subcellular distribution of sodium pump isoform subunits in mammalian cardiac myocytes. Am. J. Physiol. - Cell Physiol. 270, C1221–C1227. doi:10.1152/ajpcell.1996.270.4.c1221

PubMed Abstract | CrossRef Full Text | Google Scholar

Moise N., Struckman H. L., Dagher C., Veeraraghavan R., Weinberg S. H. (2021). Intercalated disk nanoscale structure regulates cardiac conduction. J. Gen. Physiol. 153, e202112897. doi:10.1085/jgp.202112897

PubMed Abstract | CrossRef Full Text | Google Scholar

Moreno A. P., Sáez J. C., Fishman G. I., Spray D. C. (1994). Human connexin43 gap junction channels: regulation of unitary conductances by phosphorylation. Circ. Res. 74, 1050–1057. doi:10.1161/01.RES.74.6.1050

PubMed Abstract | CrossRef Full Text | Google Scholar

Mori Y., Fishman G. I., Peskin C. S. (2008). Ephaptic conduction in a cardiac strand model with 3D electrodiffusion. Proc. Natl. Acad. Sci. U. S. A. 105, 6463–6468. doi:10.1073/pnas.0801089105

PubMed Abstract | CrossRef Full Text | Google Scholar

Noël F., Wibo M., Godfraind T. (1991). Distribution of alpha 1 and alpha 2 (Na+,K+)-ATPase isoforms between the junctional (t-tubular) and non-junctional sarcolemmal domains of rat ventricle. Biochem. Pharmacol. 41, 313–315. doi:10.1016/0006-2952(91)90494-P

PubMed Abstract | CrossRef Full Text | Google Scholar

Passiniemi P. (1983). Accurate tracer diffusion coefficients of Na+ and Cl− ions in dilute aqueous sodium chloride solutions measured with the closed capillary method. J. Solut. Chem. 12, 801–813. doi:10.1007/BF00653183

CrossRef Full Text | Google Scholar

Pedersen T. H., Huang C. L.-H., Fraser J. A. (2011). An analysis of the relationships between subthreshold electrical properties and excitability in skeletal muscle. J. Gen. Physiol. 138, 73–93. doi:10.1085/jgp.201010510

PubMed Abstract | CrossRef Full Text | Google Scholar

Peters N. S., Coromilas J., Severs N. J., Wit A. L. (1997). Disturbed connexin43 gap junction distribution correlates with the location of reentrant circuits in the epicardial border zone of healing canine infarcts that cause ventricular tachycardia. Circulation 95, 988–996. doi:10.1161/01.CIR.95.4.988

PubMed Abstract | CrossRef Full Text | Google Scholar

Petitprez S., Zmoos A. F., Ogrodnik J., Balse E., Raad N., El-Haou S., et al. (2011). SAP97 and dystrophin macromolecular complexes determine two pools of cardiac sodium channels Nav1.5 in cardiomyocytes. Circ. Res. 108, 294–304. doi:10.1161/CIRCRESAHA.110.228312

PubMed Abstract | CrossRef Full Text | Google Scholar

Polontchouk L., Haefliger J. A., Ebelt B., Schaefer T., Stuhlmann D., Mehlhorn U., et al. (2001). Effects of chronic atrial fibrillation on gap junction distribution in human and rat atria. J. Am. Coll. Cardiol. 38, 883–891. doi:10.1016/S0735-1097(01)01443-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Ponganis P. J. (2021). A physio-logging journey: heart rates of the emperor penguin and blue whale. Front. Physiol. 12, 721381. doi:10.3389/fphys.2021.721381

PubMed Abstract | CrossRef Full Text | Google Scholar

Puglisi J. L., Bers D. M. (2001). LabHEART: an interactive computer model of rabbit ventricular myocyte ion channels and Ca transport. Am. J. Physiol. - Cell Physiol. 281, C2049–C2060. doi:10.1152/ajpcell.2001.281.6.c2049

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramirez W. F. (1997). “Chapter 8 - solution of partial differential equations,” in Computational methods in process simulation. Editor W. F Ramirez. Second Edition (Oxford: Butterworth-Heinemann), 353–430.

CrossRef Full Text | Google Scholar

Rapp B. E. (2017). “Chapter 32 - finite element method,” in Microfluidics: modelling, Mechanics and mathematics micro and nano technologies. Editor B. E. Rapp (Oxford: Elsevier), 655–678.

CrossRef Full Text | Google Scholar

Rash J. E., Dillman R. K., Bilhartz B. L., Duffy H. S., Whalen L. R., Yasumura T. (1996). Mixed synapses discovered and mapped throughout mammalian spinal cord. Proc. Natl. Acad. Sci. U. S. A. 93, 4235–4239. doi:10.1073/pnas.93.9.4235

PubMed Abstract | CrossRef Full Text | Google Scholar

Rhett J. M., Gourdie R. G. (2012). The perinexus: a new feature of Cx43 gap junction organization. Hear. Rhythm 9, 619–623. doi:10.1016/j.hrthm.2011.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Rhett J. M., Veeraraghavan R., Poelzing S., Gourdie R. G. (2013). The perinexus: sign-post on the path to a new model of cardiac conduction? Trends cardiovasc. Med. 23, 222–228. doi:10.1016/j.tcm.2012.12.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Richards M. A., Clarke J. D., Saravanan P., Voigt N., Dobrev D., Eisner D. A., et al. (2011). Transverse tubules are a common feature in large mammalian atrial myocytes including human. Am. J. Physiol. - Hear. Circ. Physiol. 301, H1996–H2005. doi:10.1152/ajpheart.00284.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Rohr S. (2004). Role of gap junctions in the propagation of the cardiac action potential. Cardiovasc. Res. 62, 309–322. doi:10.1016/j.cardiores.2003.11.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Rougier J.-S., Essers M. C., Gillet L., Guichard S., Sonntag S., Shmerling D., et al. (2019). A distinct pool of Nav1.5 channels at the lateral membrane of murine ventricular cardiomyocytes. Front. Physiol. 10, 834. doi:10.3389/fphys.2019.00834

PubMed Abstract | CrossRef Full Text | Google Scholar

Salvage S. C., Huang C. L.-H., Jackson A. P. (2020a). Cell-adhesion properties of β-subunits in the regulation of cardiomyocyte sodium channels. Biomolecules 10, 989. doi:10.3390/biom10070989

PubMed Abstract | CrossRef Full Text | Google Scholar

Salvage S. C., Rees J. S., McStea A., Hirsch M., Wang L., Tynan C. J., et al. (2020b). Supramolecular clustering of the cardiac sodium channel Nav1.5 in HEK293F cells, with and without the auxiliary β3-subunit. FASEB J. 34, 3537–3553. doi:10.1096/fj.201701473RR

PubMed Abstract | CrossRef Full Text | Google Scholar

Satoh H., Delbridge L. M., Blatter L. A., Bers D. M. (1996). Surface:volume relationship in cardiac myocytes studied with confocal microscopy and membrane capacitance measurements: species-dependence and developmental effects. Biophys. J. 70, 1494–1504. doi:10.1016/S0006-3495(96)79711-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Schoenberg M., Dominguez G., Fozzard H. A. (1975). Effect of diameter on membrane capacity and conductance of sheep cardiac Purkinje fibers. J. Gen. Physiol. 65, 441–458. doi:10.1085/jgp.65.4.441

PubMed Abstract | CrossRef Full Text | Google Scholar

Severs N. J. (2000). The cardiac muscle cell. BioEssays News Rev. Mol. Cell. Dev. Biol. 22, 188–199. doi:10.1002/(SICI)1521-1878(200002)22:2<188::AID-BIES10>3.0.CO;2-T

CrossRef Full Text | Google Scholar

Severs N. J., Bruce A. F., Dupont E., Rothery S. (2008). Remodelling of gap junctions and connexin expression in diseased myocardium. Cardiovasc. Res. 80, 9–19. doi:10.1093/cvr/cvn133

PubMed Abstract | CrossRef Full Text | Google Scholar

Sheikh S. M., Skepper J. N., Chawla S., Vandenberg J. I., Elneil S., Huang C. L. H. (2001). Normal conduction of surface action potentials in detubulated amphibian skeletal muscle fibres. J. Physiol. 535, 579–590. doi:10.1111/j.1469-7793.2001.t01-1-00579.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Sigworth F. J., Neher E. (1980). Single Na+ channel currents observed in cultured rat muscle cells. Nature 287, 447–449. doi:10.1038/287447a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith J. H., Green C. R., Peters N. S., Rothery S., Severs N. J. (1991). Altered patterns of gap junction distribution in ischemic heart disease: an immunohistochemical study of human myocardium using laser scanning confocal microscopy. Am. J. Pathol. 139, 801–821.

PubMed Abstract | Google Scholar

Sperelakis N. (2002). An electric field mechanism for transmission of excitation between myocardial cells. Circ. Res. 91, 985–987. doi:10.1161/01.RES.0000045656.34731.6D

PubMed Abstract | CrossRef Full Text | Google Scholar

Sperelakis N., Mann J. E. (1977). Evaluation of electric field changes in the cleft between excitable cells. J. Theor. Biol. 64, 71–96. doi:10.1016/0022-5193(77)90114-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Struckman H., Moise N., Dunlap I., Chen Z., Weinberg S., Veeraraghavan R. (2021). Abstract 12293: intercalated disk nanostructure and molecular organization underlie conduction differences between atria and ventricles. Circulation 144, A12293. doi:10.1161/circ.144.suppl_1.12293

CrossRef Full Text | Google Scholar

Traub O., Eckert R., Lichtenberg-Frate H., Elfgang C., Bastide B., Scheidtmann K. H., et al. (1994). Immunochemical and electrophysiological characterization of murine connexin40 and -43 in mouse tissues and transfected human cells. Eur. J. Cell Biol. 64, 101–112.

PubMed Abstract | Google Scholar

Trautwein W., Kassebaum D. G., Nelson R. M., Hecht H. H. (1962). Electrophysiological study of human heart muscle. Circ. Res. 10, 306–312. doi:10.1161/01.RES.10.3.306

PubMed Abstract | CrossRef Full Text | Google Scholar

Veeraraghavan R., Hoeker G. S., Laviada A. A., Hoagland D., Wan X., King D. R., et al. (2018). The adhesion function of the sodium channel beta subunit (β1) contributes to cardiac action potential propagation. Elife 7, e37610. doi:10.7554/eLife.37610

PubMed Abstract | CrossRef Full Text | Google Scholar

Veeraraghavan R., Lin J., Hoeker G. S., Keener J. P., Gourdie R. G., Poelzing S. (2015). Sodium channels in the Cx43 gap junction perinexus may constitute a cardiac ephapse: an experimental and modeling study. Pflügers Arch. - Eur. J. Physiol. 467, 2093–2105. doi:10.1007/s00424-014-1675-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Veeraraghavan R., Lin J., Keener J. P., Gourdie R., Poelzing S. (2016). Potassium channels in the Cx43 gap junction perinexus modulate ephaptic coupling: an experimental and modeling study. Pflugers Arch. Eur. J. Physiol. 468, 1651–1661. doi:10.1007/s00424-016-1861-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Veeraraghavan R., Poelzing S., Gourdie R. G. (2014). Intercellular electrical communication in the heart: a new, active role for the intercalated disk. Cell Commun. Adhes. 21, 161–167. doi:10.3109/15419061.2014.905932

PubMed Abstract | CrossRef Full Text | Google Scholar

Verheule S., van Kempen M. J., te Welscher P. H., Kwak B. R., Jongsma H. J. (1997). Characterization of gap junction channels in adult rabbit atrial and ventricular myocardium. Circ. Res. 80, 673–681. doi:10.1161/01.res.80.5.673

PubMed Abstract | CrossRef Full Text | Google Scholar

Vermij S. H., Abriel H., van Veen T. A. B. (2017). Refining the molecular organization of the cardiac intercalated disc. Cardiovasc. Res. 113, 259–275. doi:10.1093/cvr/cvw259

PubMed Abstract | CrossRef Full Text | Google Scholar

Vetterlein F., Mühlfeld C., Cetegen C., Volkmann R., Schrader C., Hellige G. (2006). Redistribution of connexin43 in regional acute ischemic myocardium: influence of ischemic preconditioning. Am. J. Physiol. - Hear. Circ. Physiol. 291, H813–H819. doi:10.1152/ajpheart.01177.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Wan X., Bryant S., Hart G. (2003). A topographical study of mechanical and electrical properties of single myocytes isolated from normal Guinea-pig ventricular muscle. J. Anat. 202, 525–536. doi:10.1046/j.1469-7580.2003.00187.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Watanabe A., Grundfest H. (1961). Impulse propagation at the septal and commissural junctions of crayfish lateral giant axons. J. Gen. Physiol. 45, 267–308. doi:10.1085/jgp.45.2.267

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: ephaptic conduction, cardiomyocytes, action potential propagation, sodium channels, electrodiffusion

Citation: Morris JA, Bardsley OJ, Salvage SC, Jackson AP, Matthews HR and Huang CL-H (2024) Nernst-Planck-Gaussian modelling of electrodiffusional recovery from ephaptic excitation between mammalian cardiomyocytes. Front. Physiol. 14:1280151. doi: 10.3389/fphys.2023.1280151

Received: 19 August 2023; Accepted: 01 December 2023;
Published: 03 January 2024.

Edited by:

Fernando Soares Schlindwein, University of Leicester, United Kingdom

Reviewed by:

Seth H. Weinberg, The Ohio State University, United States
Alvaro Macias, Spanish National Centre for Cardiovascular Research, Spain

Copyright © 2024 Morris, Bardsley, Salvage, Jackson, Matthews and Huang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Christopher L-H. Huang, clh11@cam.ac.uk

These authors share first authorship

These authors share last authorship

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.