The intercalated disk (ID) is a specialized subcellular region that provides electrical and mechanical connections between myocytes in the heart. The ID has a clearly defined passive role in cardiac tissue, transmitting mechanical forces and electrical currents between cells. Recent studies have shown that Na+ channels, the primary current responsible for cardiac excitation, are preferentially localized at the ID, particularly within nanodomains such as the gap junction–adjacent perinexus and mechanical junction–associated adhesion-excitability nodes, and that perturbations of ID structure alter cardiac conduction. This suggests that the ID may play an important, active role in regulating conduction. However, the structures of the ID and intercellular cleft are not well characterized and, to date, no models have incorporated the influence of ID structure on conduction in cardiac tissue. In this study, we developed an approach to generate realistic finite element model (FEM) meshes replicating nanoscale of the ID structure, based on experimental measurements from transmission electron microscopy images. We then integrated measurements of the intercellular cleft electrical conductivity, derived from the FEM meshes, into a novel cardiac tissue model formulation. FEM-based calculations predict that the distribution of cleft conductances is sensitive to regional changes in ID structure, specifically the intermembrane separation and gap junction distribution. Tissue-scale simulations predict that ID structural heterogeneity leads to significant spatial variation in electrical polarization within the intercellular cleft. Importantly, we found that this heterogeneous cleft polarization regulates conduction by desynchronizing the activation of postjunctional Na+ currents. Additionally, these heterogeneities lead to a weaker dependence of conduction velocity on gap junctional coupling, compared with prior modeling formulations that neglect or simplify ID structure. Further, we found that disruption of local ID nanodomains can either slow or enhance conduction, depending on gap junctional coupling strength. Our study therefore suggests that ID nanoscale structure can play a significant role in regulating cardiac conduction.
The intercalated disk (ID) is the specialized structure at the site of the cell–cell interface connecting myocytes in the heart, facilitating electrical coupling via gap junctions (GJs) and mechanical coupling via mechanical junctions (MJs; Leo-Macías et al., 2015). Recent evidence has highlighted the structural heterogeneity of the ID (Zhang et al., 1996; Miyamoto et al., 2002; Shimada et al., 2004; Leo-Macías et al., 2015; Pinali et al., 2015; Leo-Macias et al., 2016; Vanslembrouck et al., 2018; Vanslembrouck and Kremer, 2020): the ID consists of highly tortuous plicate regions (which are oriented perpendicular to the myocyte’s axis, comprise membrane “folds,” and contain the majority of MJs [adherens junctions, desmosomes, and composite junctions]) and interplicate regions (which run parallel to the myocyte’s axis and contain most GJs). MJs and GJs have clearly defined roles in the passive conduction of mechanical forces and electrical currents, respectively. However, recent studies show that multiple electrogenic proteins regulating conduction, in particular the voltage-gated Na+ channel, NaV1.5, are localized at the ID (Maier et al., 2004; Lin et al., 2011; Petitprez et al., 2011; Rhett et al., 2012; Agullo-Pascual et al., 2014; Veeraraghavan et al., 2015, 2018; Veeraraghavan and Gourdie, 2016; Vermij et al., 2017), in close proximity to the junctions, specifically GJ-adjacent perinexi and MJ adhesion-excitability nodes, forming specialized nanodomains in the intercellular cleft spaces (Radwański et al., 2018). Further, structural perturbations of the ID have been shown to alter conduction and risk of arrhythmias (Veeraraghavan et al., 2012; Mezache et al., 2020; George et al., 2015; George and Poelzing, 2016; Veeraraghavan et al., 2018), and ID structure is disrupted in patients with atrial fibrillation (Raisch et al., 2018). Collectively, these findings suggest that ID structure plays an important active role in regulating cardiac conduction.
Standard cardiac tissue modeling approaches, e.g., the monodomain model, typically neglect ID structure and nonuniform ion channel distributions. However, prior in silico studies, including our work, that do account for these details predict that NaV1.5 localization at the ID facilitates cell–cell coupling via a mechanism known as ephaptic coupling (Lin and Keener, 2010; Mori et al., 2008; Tsumoto et al., 2020; Kucera et al. 2002; Wei et al., 2016; Wei and Tolkacheva, 2020; Sperelakis and Mann, 1977; Sperelakis and Mann, 2002; Weinberg, 2017; Greer-Short et al., 2017; Nowak et al., 2020; Nowak et al., 2021; Veeraraghavan et al., 2015; Tveito et al., 2017). Ephaptic coupling is mediated by interactions between Na+ currents, INa, on adjacent ID membranes sharing a restricted intercellular cleft space and the associated electrochemical gradients that form in this space. Specifically, during the cardiac action potential, INa at the ID in a depolarizing cell (i.e., prejunctional INa) reduces the electrical potential of the intercellular cleft. This reduction in intercellular cleft potential depolarizes the apposing cell from the extracellular, rather than the intracellular, side of the cell membrane. Additionally, ID-localized INa reduces the local Na+ concentration of the intercellular cleft, and collectively these electrochemical gradients alter both the apposing ID membrane transmembrane potential and postjunctional INa driving force. Prior work has predicted that a narrow intercellular cleft width, which enhances the effects of ephaptic coupling, can either enhance or slow conduction, depending on the timing and magnitude of postjunctional INa (Kucera et al., 2002; Mori et al., 2008; Weinberg, 2017).
While providing important insights into how nonuniform Na+ distribution regulates conduction, most of these modeling studies simply represent the intercellular cleft as a single, uniform compartment. Two recent studies have considered spatial variation within the ID and the cleft: Mori et al. (2008) investigated ephaptic coupling in a radially symmetric strand of cells, and Hichri et al. (2018) modeled the ID Na+ channel distribution as either uniform or in two apposing clusters. However, all of these prior studies make a key simplifying assumption about the ID and intercellular cleft space, specifically that the cleft is a uniform, cylindrical space and the intermembrane separation between apposing ID membranes is uniform throughout the ID. Thus, all prior modeling studies have neglected the structural heterogeneity of the ID and intercellular cleft space (e.g., variations between plicate and interplicate regions). Thus, the mechanisms underlying conduction changes due to ID structural perturbations (Veeraraghavan et al., 2012; Mezache et al., 2020; George et al., 2015; George and Poelzing, 2016; Veeraraghavan et al., 2018) are not well understood.
The goal of this work is to develop an approach to integrate ID structural heterogeneity into a model of cardiac tissue and investigate the impact of this heterogeneity on conduction. For this approach, we first generate realistic 3-D finite element model (FEM) meshes of the ID and intercellular cleft, based on measurements from transmission electron microscopy (TEM) images of mouse ventricles. We then calculate a reduced electrical network, based on the FEM meshes, to represent the cleft, crucially bypassing the computationally costly step of simulating electrostatics on the full 3-D FEM mesh, and incorporate this reduced network into a novel cardiac tissue model formulation that accounts for the heterogeneity at the ID. In our study, simulations predict that ID structural heterogeneity has a significant impact on electrical conduction. Notably, we find that both the monodomain and single-cleft tissue models either over- or underestimate predictions of conduction velocity (CV), depending on the relative strength of GJ coupling. Importantly, the ID structure results in heterogeneous, asynchronous behavior of the electrical potential and Na+ concentration in the intercellular cleft and INa at the ID, and these changes collectively regulate cardiac conduction. These findings improve our understanding of arrhythmia mechanisms associated with pathological ID structural remodeling and suggest the ID structure as a potential new therapeutic target.
Materials and methods
The overall goal of incorporating ID structure into a cardiac tissue model is accomplished in three stages. First, we develop an algorithm to construct a 3-D FEM mesh of the ID and intercellular cleft, reproducing the structure of key ID measurements from TEM images. Importantly, this mesh generation incorporates several orders of magnitude of structural details, from the nanoscale structure of intermembrane separation, up to the microscale structure of the interplicate and plicate regions. Second, we calculate an equivalent electrical network for the conductivity within and out of the intercellular cleft space. We calculate this reduced cleft network by partitioning the full finite element mesh into a tractable number (25–200) of compartments and determining the equivalent electrical conductance between all pairs of adjacent compartments. Finally, we incorporate this cleft network into a cardiac tissue model, in which neighboring myocytes are coupled via GJs and their shared intercellular cleft space. Importantly, by representing ID-localized ion channels that induce heterogeneous electrical polarization and ionic concentration gradients within the intercellular cleft spaces, the tissue model links nanoscale ID structure with macroscale cardiac tissue function.
All animal procedures were approved by Institutional Animal Care and Use Committee at the Ohio State University and were performed in accordance with the Guide for the Care and Use of Laboratory Animals published by the U.S. National Institutes of Health (publication no. 85–23, revised 2011). Male C57/BL6 mice (30 g, 6–18 wk) were anesthetized with 5% isoflurane mixed with 100% oxygen (1 liter/min). After loss of consciousness, anesthesia was maintained with 3–5% isoflurane mixed with 100% oxygen (1 liter/min). Once the animal was stably in a surgical plane of anesthesia, the heart was excised, leading to euthanasia by exsanguination. The isolated hearts were prepared for fixation for TEM. Ventricles were dissected and fixed overnight in 2% glutaraldehyde at 4°C for resin embedding and ultramicrotomy as previously described (Mezache et al., 2020; Veeraraghavan et al., 2018).
TEM and quantification
TEM images of the ID, particularly GJs and MJs, were obtained at 6,000×, 10,000×, and 20,000× magnification on an FEI Technai G2 Spirit transmission electron microscope (Thermo Fisher Scientific). Morphometric quantification was performed using ImageJ (National Institutes of Health) by manually identifying and quantifying the following 11 specific ID measurements. Images at 6,000× magnification were used to quantify total ID cross-sectional length (1). Images at 10,000× magnification were used to quantify the length of plicate and interplicate regions (2 and 3; Fig. 1 B); the amplitude and frequency of plicate folds (4 and 5; Fig. 1 E), and the lengths of GJs and the fraction of membrane comprised of GJs, in both plicate and interplicate regions (6–9; Fig. 1 E). Images at 20,000× magnification were used to quantify the intermembrane distance in the plicate and interplicate regions, identified by regions outside MJs and GJs, respectively (10 and 11; Fig. 1 E). The Wilcoxon rank sum test was used for single comparisons between plicate and interplicate region measurements. A P value of <0.01 was considered statistically significant. Measurements are reported by the mean ± SEM and the first and third quartile range.
FEM mesh generation
Based on measurements from TEM images, we developed an approach to replicate the nanoscale structure of the ID, generating realistic 3-D FEM meshes of the two cell membranes at the ID and the enclosed intercellular cleft space (Fig. 1). The resulting structures are randomly generated but parametrically defined based on the TEM measurements described above (Fig. S1, left). Specifically, these mesh parameters define the key features of the ID structure and fit them to experimental data: ID diameter, plicate and interplicate region length, plicate and interplicate region intermembrane distance, GJ size and distribution in both regions, and plicate fold amplitude and frequency (Table S1).
ID and GJ map generation
The first step in the mesh generation involves developing “maps” that represent the patterns for two key structures, specifically the plicate regions within the ID and location of GJ clusters (Fig. 2). Full details of this map generation approach are described in the Supplemental text at the bottom of the PDF. The ID map represents a cross section of the entire ID, such that the map defines different levels or tiers of plicate regions (defined as representing the (x,y) plane), which are separated by interplicate regions (in the z-direction; see Fig. 1 C). The GJ map represents the locations of individual GJ clusters. The same methodology is used to generate maps of both structures, with different algorithm parameters. In brief, a 2-D map is generated by filtering and thresholding Gaussian random noise (Fig. S2 A), for which the threshold defines two sets of clusters (plicate levels for the ID map, GJ versus membrane for the GJ map). By varying the properties of the Gaussian noise and filter, we vary the properties of the resulting maps (Fig. S2, B and C). Importantly, we adjust mesh generation parameters (Table S1) such that the map properties match the corresponding TEM measurements, specifically measurements of plicate region length, and GJ lengths and distributions in the plicate and interplicate regions. In particular, estimating the map generation parameters required comparing the 2-D map properties with the 1-D measurements obtained from TEM images; we made “slices” of the 2-D maps to obtain a series of 1-D length measurements along each slice to directly compare model and experimental measurements (illustrated in Fig. S3). Additionally, the random nature of the algorithm is a key feature of the mesh generation and overall study, as this enables the generation of multiple IDs with comparable and statistically similar overall properties and allows us to account for experimental variability in ID structure.
Finite element mesh construction
After generating maps of the individual ID structures fit to TEM measurements, we next generate the finite element mesh, which is comprised of plicate regions, interplicate regions, and small connections between these regions (Fig. 2). Each plicate and interplicate mesh is initially defined by two rectangular surfaces, which represent the two apposing cell membranes, and then is populated by GJs, with distribution specific to that ID region and determined independently. As GJs mediate direct electrical communication between cells, at the location of GJs, the two apposing cell membranes connect, which appear as “holes” in the mesh (Fig. 1, C–E). Additionally, the distance between the two membranes is defined, based on TEM measurements, denoted as dP for the plicate region and dIp for the interplicate region. Importantly, we note that these two distances can be independently varied.
We next transform the individual mesh components into the shape of the ID. The plicate meshes are cropped into the shapes of the individual levels, determined by the ID map, and further are also folded using the shape of a 2-D sine wave, to match the amplitude and frequency of the plicate folds (Fig. 1, C–E). The interplicate region meshes are rectangles, with lengths that match the contours of the plicate region levels and heights determined by the TEM measurements of interplicate region length. The interplicate meshes are then curved and rotated to be positioned between the corresponding plicate region meshes. Small mesh pieces then join all adjacent plicate and interplicate regions. This complete mesh represents the two apposing cell membranes. The mesh is then checked and repaired using the meshfix algorithm (Attene, 2010), ensuring that the surfaces completely enclose a volume, representing the intercellular cleft. Finally, we use Gmsh, finite element mesh generation software (Geuzaine and Remacle, 2009), to generate a 3-D tetrahedral mesh of the intercellular cleft space.
Intercellular cleft conductance measurement
Finally, the equivalent conductance between cleft compartments j and k (gc,i,j,k), for the cleft between cells i and i + 1, is calculated by dividing the total current through the mesh, by the voltage difference across the mesh, i.e., gc,i,j,k = Ifem/Vfem. Supplemental text at the bottom of the PDF illustrates this equivalent conductance calculation in a regular rectangular geometry, resulting in the established relationship between electrical conductance, geometry, and electrical resistivity (Fig. S4). We repeat this calculation over all pairs of adjacent compartments in the intracellular cleft to obtain the intracleft conductances gc,i,j,k for all combinations of j, k = 1,…, M and j ≠ k. By definition, gc,i,j,k = 0 for nonadjacent compartments and gc,i,j,j = 0. Additionally, for all compartments on the ID periphery, we calculate the effective conductance between cleft compartment and the bulk extracellular space (gb,i,j) using a similar procedure, replacing the second compartment centroid with the center of the compartment-ID boundary; for all compartments not adjacent to the bulk, gb,i,j = 0.
We investigate the dependence of the cleft conductance measurements on intermembrane separation by repeating the mesh generation and cleft conductance measurements for different values of dP and dIp. Additionally, we perform this process five times to replicate the experimental variability in ID and cleft properties. We also investigate the influence of ID structures by measuring cleft conductances in the absence of either GJs or plicate folds.
The resulting electrical network is formed by nodes representing the average electrical potential of the compartment and edges representing intracleft electrical conductance (Fig. 3 D). We validate this approach by comparing the electrical potential on the reduced cleft electrical network and the corresponding full 3-D cleft finite element mesh (calculated in finite element solver COMSOL). We note that the reduced cleft network represents a significant computational reduction: solving for the 3-D cleft potential requires solving Laplace’s equation on the entire cleft finite element mesh, comprising ∼300,000 finite element tetrahedra (Fig. 3 E), ∼1 min of computational time, while solving for the cleft potential on the reduced cleft network requires solving an M-dimensional linear system, obtained from applying Kirchoff’s current law to the network circuit, requiring <1 s of computational time. Fig. 3 F illustrates the close comparison of the voltage of the cleft network nodes and the corresponding average partition voltage, for all partitions. In the supplemental text, Iteration we quantify the RMSD between the cleft network and cleft finite element solution, and we found that increasing the number of partitions M reduced the error (Fig. S5). Importantly, for ≥50 partitions, we found that the error for the reduced cleft network potential is <10%.
Cleft network tissue model with ID nanodomains
We next developed a novel cardiac tissue model formulation that incorporates the reduced cleft network, thus integrating the effects of ID nanoscale structure into a tissue-scale model. Full details of the cardiac tissue model are provided in the supplemental text, including equations, parameters, and numerical methods. In brief, we simulate a 50-cell cable of ventricular myocytes, with membrane patch dynamics governed by the Luo-Rudy 1 model (Luo and Rudy, 1991). As in our prior work (Weinberg, 2017; Greer-Short et al., 2017; Nowak et al., 2020; Nowak et al., 2021) and work by others (Tsumoto et al., 2020; Kucera et al., 2002; Lin and Keener, 2010; Veeraraghavan et al., 2015; Wei et al., 2016; Wei and Tolkacheva, 2020), we account for nonuniform Na+ subcellular localization by spatially discretizing each cell into axial membrane patches along the lateral membrane and ID membrane patches at site of the cell–cell junctions. In contrast with prior work in which each ID membrane was represented by a single membrane patch, in this study, each ID is discretized into M membrane patches (Fig. 3 G). The extracellular side of each ID membrane patch, representing the intercellular cleft space, is coupled via the reduced cleft electrical network (Fig. 3, D and G, right), as described above. That is, the jth and kth cleft potentials between the ith and i+1th cell (ϕc,i,j and ϕc,i,k, respectively) are coupled with conductance gc,i,j,k. This critical modification incorporates the ID nanoscale structure into the macroscale tissue model.
The intracellular electrical potential (ϕm,i) in each cell is coupled to the pre- and postjunctional intracellular potentials (ϕr,i and ϕl,i, referring to right and left potentials, respectively) with conductance gmyo. In addition to ephaptic coupling via the shared intercellular cleft, pre- and postjunctional intracellular potentials (ϕr,i and ϕl,i+1, respectively) are coupled by a GJ conductance ggap,i. Note that the finite element mesh generates the spatial location of individual GJ clusters, such that the tissue model can represent these distinct electrical connections between the pre- and postjunctional intracellular spaces. However, all of the GJ resistors are connected in parallel, such that the overall GJ conductance between neighboring cells can be represented by a single resistor, with conductance equal to the sum of the individual GJ cluster conductances (i.e., the macroscopic GJ conductance). We also note that while in general the cleft conductances (gc,i,j,k) could differ between different clefts, we assume these values are the same throughout the tissue, i.e., Similarly we assume the same GJ conductance between all adjacent cells, i.e., ggap,i = ggap.
We account for the dynamic [Na+] in each of the intercellular cleft compartments ([Na+]c,j), where [Na+]c,j in the jth cleft compartment is governed by the Na+ current in the jth pre- and postjunctional ID membrane patch, and ionic flux between adjacent compartments (and the bulk, for periphery compartments). The ionic flux between cleft compartments is governed by the electrochemical gradient, i.e., the cleft potential and ionic concentration differences between the adjacent compartments.
We perform simulations with the cleft network tissue model using the cleft conductances obtained from the finite element meshes for different intermembrane separation values. We also vary GJ coupling levels, i.e., values of ggap. For comparison with the cleft network model, we also perform simulations with a single-cleft tissue model, the typical approach of prior studies of ephaptic coupling that neglects ID structural heterogeneity (Kucera et al., 2002; Weinberg, 2017; Lin and Keener, 2010), and also with the standard monodomain tissue model formulation, which does not account for the ID, intercellular cleft space, or nonuniform Na+ channel distribution (Fig. S6). For all the tissue simulations, we pace the leftmost cell (cell 1) at a cycle length of 500 ms. CV was computed by linear regression of the activation times of the intracellular potential of cells 15–35, where the activation time is defined as the time when ϕm crosses above −60 mV.
Online supplemental material
Fig. S1 illustrates mesh generation and cleft conductance calculation algorithms. Fig. S2 shows map generation algorithm and GJ distribution calculations. Fig. S3 shows the fitting of cluster length measurements to 2-D maps. Fig. S4 illustrates Laplace’s equation on a regular geometry and uniform material. Fig. S5 shows the error of the cleft network approximation of the Laplace’s equation solution on the full 3-D cleft FEM mesh. Fig. S6 is a schematic of the electrical circuits for the tissue models compared in the study. Fig. S7 shows interplicate and plicate region conductances for different FEM mesh variants. Fig. S8 shows examples of FEM mesh variants used in the calculations in Fig. 4. Fig. S9 shows mesh incorporating larger GJs at the ID edge. Table S1 lists parameters for cleft FEM mesh generation. Table S2 lists key parameters and values in the cleft network tissue model. Supplemental text at the bottom of the PDF provides full details of the mesh map generation.
ID structural properties
Representative TEM images illustrate the interplicate and plicate regions of the ID (Fig. 1 B) and plicate folds, GJs, and intermembrane distance in the plicate and interplicate regions (Fig. 1 E). Summary measurements of the key ID structures used to generate the finite element meshes are shown in Table 1. In particular, we find several key differences between the plicate and interplicate regions. GJ length is longer and GJs comprise a larger percentage of the cell membrane in the interplicate, compared with the plicate. Intermembrane separation is narrower in the interplicate regions, compared with the plicate, consistent with our recent measurements in atria (Mezache et al., 2020). Our measurements of plicate fold amplitude are comparable with recent studies as well (Vanslembrouck et al., 2018; Vanslembrouck and Kremer, 2020). Additionally, we find that plicate regions tend to be longer than interplicate, although not significantly.
Intercellular cleft conductances
Using the baseline mesh generation parameters obtained from TEM measurements (Table S1), we first demonstrate the generality of the mesh generation process and created five different ID maps and meshes (Fig. 4, A and B). Importantly, each finite element mesh is generated from the same parameters but differs due to randomness inherent in the mesh generation process. Slices through the meshes are visually quite comparable to TEM images but also differ from each other due to differences in the mesh and orientation of the slices, exhibiting differences in plicate and GJ lengths (Fig. 4 C). As described above, we calculate the intracleft conductances between all adjacent cleft compartments for each mesh. We divided each mesh into 200 compartments or partitions, resulting in ∼450 cleft network edges, split between the plicate and interplicate regions. Histograms of the cleft conductances show some small differences between the different meshes (Fig. 4 D); however, overall the cleft conductance distributions are similar, with nearly identical means (black lines).
We next investigated the differences between cleft conductances specifically within the interplicate (gIp) and plicate (gP) regions, and then further to what extent different ID properties altered these conductances. Histograms illustrate that conductances in the plicate are typically larger than in the interplicate, with the mean plicate conductance (gP, orange vertical line) almost twice as large as the mean interplicate conductance (gIp, blue vertical line; Fig. 5, A and E). Note that these conductances represent all of the values from the five meshes in Fig. 4. This overall trend is consistent with higher conductance associated with the wider intermembrane separation in the plicate, compared with the interplicate.
We next investigate to what extent nanoscale ID structures influence cleft conductance, specifically the presence of GJs and plicate folds. We find that GJs decrease conductance in the interplicate (Fig. 5 B), due their high density in this region. Thus, in addition to providing a direct electrical connection between coupled myocytes, GJs also reduce the electrical conductance of the intercellular cleft. In contrast, as GJs are both smaller and rarer in the plicate, the presence of GJs in the plicate minimally influences cleft conductances (Fig. S7 D). We also find that plicate folds greatly decrease plicate conductance (Fig. 5 F), as the folds increase the effective cleft length and thus decrease conductance. As expected, plicate folds do not influence interplicate conductances (Fig. S7 B).
We additionally investigate how changes in intermembrane separation influence the interplicate and plicate conductances. As expected, increasing interplicate intermembrane separation dIp shifts the interplicate conductance histogram to the right (Fig. 5 C), consistent with an overall increase in mean interplicate conductance gIp (Fig. 5 D, blue line). Interestingly, we also observe a small increase in mean plicate conductance gP as dIp increases, due to a few cleft compartments at the interface between interplicate and plicate (Fig. 5 D, orange line). Similarly, we find that increasing plicate intermembrane separation dP shifts the plicate conductance histogram to the right (Fig. 5 G), consistent with an overall increase in mean plicate conductance gP (Fig. 5 H, orange line), with minimal increase in gIp as dP increases (Fig. 5 H, blue line).
Overall, the results show that both intermembrane distance and ID structural features alter local intracleft conductance. Interestingly, calculations also predict that the plicate region is more sensitive to changes in intermembrane separation, compared with the interplicate. In the absence of irregular structure, conductance would be directly proportional to intermembrane separation. However, simulations predict that a 2-fold increase in dP (from 30 to 60 nm) leads to a 3.1-fold increase in gP, while a similar 2-fold increase in dIp (from 15 to 30 nm) leads to a 1.8-fold increase in gIp. Thus, these changes in sensitivity specifically arise due to structural heterogeneities introduced by GJs and plicate folds. Illustrations of the presence/absence of these ID structures and different intermembrane separation are shown in Fig. S8.
ID structure regulation of electrical conduction
We next investigate how nanoscale ID structure and changes in this structure influence tissue level conduction. We first compared our cleft network model with the standard monodomain model and a single-cleft tissue model (Fig. 6). As expected, the monodomain model illustrates the activation of Na+ current during conduction but does not account for changes in cleft potential and [Na+], as these are not accounted for in the formulation (Fig. 6, left column). Consistent with prior work (Lin and Keener, 2010; Weinberg, 2017; Greer-Short et al., 2017; Veeraraghavan et al., 2015), the single-cleft model illustrates cleft hyperpolarization and a transient depletion of cleft [Na+]; additionally, the Na+ current magnitude at the ID is increased, relative to the monodomain model, due to the preferential localization of Na+ current on the ID membrane (Fig. 6, middle column). In contrast with the single-cleft model, the cleft network model illustrates the heterogeneous ID Na+ currents, cleft [Na+], and cleft potentials during conduction (Fig. 6, right column).
Importantly, we note that the cleft network average (thick lines) differs from the single-cleft model, demonstrating that the single-cleft model does not capture the overall dynamics of ID currents, cleft concentrations, and cleft potentials. Specifically, we find that cleft potentials are hyperpolarized to a greater extent and cleft [Na+] refills slower, compared with the single-cleft model, due to the complicated ID geometry and weaker coupling to the bulk. Both of these changes ultimately reduce Na+ current driving force, such that Na+ current at the ID is also reduced, compared with the single-cleft model. In Fig. 6 E, we illustrate the heterogeneity of cleft polarization during conduction. We find that cleft nodes closer to the ID center are more hyperpolarized than those toward the periphery, due to the low conductance path from the center to the bulk. Further, nodes in the interplicate also tend to be more hyperpolarized than in the plicate, due to lower interplicate conductance (compare Fig. 5, A and E).
We next investigate to what extent ID structural differences alter CV (Fig. 7). We measured CV in the cleft network models and compared measurements with the standard monodomain tissue model and a tissue model with a single-cleft compartment. Additionally, in the two cleft models, we varied the intermembrane separation for the single-cleft or specifically in the interplicate regions for the cleft network (i.e., varying dIp), for values within a physiological regimen (15–60 nm) and also considered a pathologically disrupted interplicate (200 nm). We varied interplicate intermembrane separation, as opposed to varying intermembrane separation throughout the ID, to investigate the influence of perturbations that target specific regions of the ID. We note that the same ID and GJ maps were used to generate finite element meshes with varying interplicate intermembrane separation. Experimental measurements of macroscopic GJ conductance vary considerably, typically on the order of 10s to 100s of nanosiemens (Desplantez et al., 2007; Kwak and Jongsma, 1996; McCain et al., 2012; Moreno et al., 1994; Valiunas et al., 2002; Verheule et al., 1997; White et al., 1990; Nielsen et al., 2012; Rüdisüli and Weingart, 1989; Weingart, 1986; Wittenberg et al., 1986). Thus, our simulations consider a range of GJ conductances.
We find that for strong GJ coupling, CV was highest in the monodomain model and lowest in the cleft network model (Fig. 7 A). By neglecting the ID and nonuniform channel distribution, the monodomain model overestimates CV by ∼20–30%. The single-cleft model also overestimates CV by ∼10% by neglecting ID structural heterogeneity. For both cleft models, increased intermembrane separation increases CV, to a slightly greater extent in the cleft network model.
For moderate GJ coupling, CV is similar for all three models, although the monodomain model also predicts a slightly faster CV compared with the cleft models (Fig. 7 B). Interestingly, the cleft models predict different trends for increasing intermembrane separation: CV decreases in the single-cleft model and increases in the cleft network model, although these differences are small. Finally, for weak GJ coupling, we find both quantitative and qualitative differences between model predictions (Fig. 7 C). The single-cleft model predicts CV slower than the monodomain model and CV increasing as intermembrane separation increases. In contrast, CV in the cleft network model is ∼20% faster compared with the monodomain model, and CV decreases as interplicate membrane separation increases. Further, the cleft network model is much more sensitive to changes in intermembrane separation, compared with the single-cleft model, as dIp increasing from 15 to 60 nm results in an ∼20% decrease in CV.
In Fig. 7 D, we plot CV as a function of GJ conductance (ggap) for different values of interplicate intermembrane separation (dIp). This plot in particular highlights that CV is most sensitive to GJ conductance for wide interplicate intermembrane separation, while in contrast, CV is least sensitive to GJ conductance for narrow interplicate intermembrane separation. Additionally, we observe that CV is most sensitive to interplicate intermembrane separation for either low or high GJ conductance, with opposite trends as noted above. Further, the cleft network model CV, normalized relative to the CV in the monodomain model, is most sensitive to interplicate intermembrane separation for low GJ conductance, demonstrating the regimen of the greatest ID structural sensitivity and discrepancy with the monodomain model. Additionally, we note that as dIp increases, for all gap junctional coupling strengths, the cleft network model CV tends to be closer to the monodomain model, i.e., the normalized value approaches 1 (dashed horizontal line), demonstrating that weaker ephaptic coupling in the cleft network model (and specifically in the interplicate region) results in conduction closer to the standard monodomain model.
Collectively, these findings show that CV in the cleft network model exhibits a weaker dependence on GJ coupling, in particular for narrow interplicate membrane separation, compared with standard modeling approaches that neglect the ID or assume a homogeneous ID. Further, we find that ID structural heterogeneity can result in enhanced sensitivity to intermembrane separation, compared with the single-cleft modeling approach assuming a homogeneous ID and intercellular cleft.
We next investigate the mechanism underlying intermembrane separation sensitivity due to ID structural heterogeneity in the cleft network model. To investigate this sensitivity and the differing responses to interplicate expansion, we consider cleft and membrane dynamics during conduction for four cases (Fig. 8): (1) strong GJ coupling and narrow interplicate, (2) weak GJ coupling and narrow interplicate, (3) strong GJ coupling and expanded interplicate, and (4) weak GJ coupling and expanded interplicate. In all cases, we illustrate the intracellular (black) and cleft potentials (blue), cleft [Na+] (red), GJ current (orange), and the Na+ currents at the prejunctional (cyan) and postjunctional (magenta) membranes, in the middle of the tissue (between cells 25 and 26).
For strong GJ coupling and narrow interplicate (Fig. 8 A), the cleft is hyperpolarized in a highly heterogeneous manner (blue; see also Fig. 6 E), and cleft [Na+] (red) is locally depleted, reducing Na+ current driving force heterogeneously. Collectively, these effects result in desynchronized activation of the postjunctional Na+ currents, with variable current magnitude (magenta), which reduces the overall ID Na+ current density (thick magenta). In contrast, for wider interplicate (Fig. 8 C), cleft hyperpolarization is reduced and cleft depletion is attenuated, such that postjunctional Na+ current is more synchronized and less variable, resulting in a larger ID Na+ current density. Additionally, GJ current magnitude (orange) is greater for the wider interplicate. Collectively, the larger Na+ and GJ current result in a faster upstroke and thus faster CV for the wider interplicate.
For weak GJ coupling (Fig. 8, B and D), in both cases, GJ current is greatly reduced and conduction is slower, compared with strong GJ coupling, as expected. Interestingly, expansion of the interplicate exhibits a behavior similar to the previous case, but with opposite effect on conduction: For narrow interplicate (Fig. 8 B), heterogeneous cleft hyperpolarization and cleft [Na+] depletion similarly result in desynchronized postjunctional Na+ current. However, the larger magnitude of cleft polarization results in both an overall earlier activation and longer duration of postjunctional Na+ current, compared with a wider interplicate (Fig. 8 D), which results in an earlier upstroke and faster conduction for the narrower interplicate. Thus, we find that ID structural heterogeneity plays a significant role in regulating conduction: narrow intermembrane separation promotes desynchronized and varying magnitude postjunctional Na+ current, which slows conduction for strong GJ coupling yet enhances conduction for weak GJ coupling.
In this study, we developed an approach to generate realistic FEM meshes of the ID, based on TEM-derived measurements of nano- to microscale structure, which are incorporated into a 1-D tissue model of cardiac conduction. To our knowledge, this is the first study to integrate FEM-based nanoscale structural modeling with a tissue-scale cardiac electrophysiology model. We investigated the effects of ID structure on intracleft conductance on multiple FEM mesh samples. We found that conductance increases with intermembrane distance, with differences between distinct regions (i.e., plicate versus interplicate), and that conductance in interplicate regions is decreased by the nanoscale structure of GJs and associated perinexi (Rhett et al., 2012; Rhett and Gourdie, 2012; Veeraraghavan et al., 2015), while membrane folds have a similar effect in plicate regions.
By incorporating FEM mesh structural data into a tissue model, we demonstrate that changes in regional ID intermembrane distance (i.e., in the interplicate regions) regulate cardiac conduction. In particular, we found that for tissue with strong GJ coupling, interplicate region intermembrane expansion increased CV, while for tissue with weak GJ coupling, the opposite relationship was observed. Importantly, we found that both effects depend on the ID and cleft structural heterogeneity: specifically, heterogeneous cleft hyperpolarization leads to desynchronized activation of the postjunctional Na+ currents. This heterogeneous behavior, for strong GJ coupling, leads to a lower average Na+ current and slower CV. However, in the case of weak GJ coupling, this behavior activates the downstream cell earlier and thus increases CV. More broadly, we find that previous model formulations that neglect ID structural heterogeneity (i.e., the single-cleft model) or the ID entirely (i.e., the monodomain model) can either under- or overestimate predictions of CV. Overall, the cleft network model exhibits a weaker dependence of CV on GJ coupling compared with other modeling approaches. In particular, across a wide range of GJ coupling strengths, we find that CV is less sensitive to GJ coupling for narrower interplicate intermembrane separation, while in contrast, CV sensitivity to GJ coupling is enhanced for wider interplicate intermembrane separation.
Our model predictions demonstrate two critical points with implications for predicting arrhythmia risk in physiological and pathological settings, specifically (a) that incorporating ID structural details into cardiac tissue models can impact conduction predictions (compared with prior approaches), and (b) that perturbations in ID structure (e.g., interplicate intermembrane expansion) can significantly alter conduction. This second point is consistent with prior experimental studies demonstrating that ID disruption can alter conduction and arrhythmia risk. For example, recent work has shown that perinexal structure (interplicate ID) is perturbed in atrial tissue from patients’ atrial fibrillation (Raisch et al., 2018). We also recently showed that vascular leak, induced by inflammatory cytokine vascular endothelial growth factor, disrupts ID nanodomains in both GJ- and MJ-adjacent regions, and ultimately slows conduction and increases susceptibility to arrhythmias (Mezache et al., 2020). We previously showed that disruption of GJ-adjacent perinexi via edema (Veeraraghavan et al., 2012, 2015, 2016; George et al., 2015), inflammation (George et al., 2015), or adhesion-disrupting peptides targeting the Na+ channel subunit (Veeraraghavan et al., 2018) can also slow conduction. In agreement with the majority of GJs residing in the interplicate regions, these conduction-slowing trends are consistent with our findings in the setting of weaker GJ coupling and interplicate expansion (Fig. 7 C). We speculate that differences between experiments and simulations for cases of stronger GJ coupling may arise due to the lack of accounting for intracleft clustering of Na+ channels at GJ-adjacent perinexi (Veeraraghavan et al., 2018; Veeraraghavan et al., 2015; Veeraraghavan and Gourdie, 2016) and MJ-associated nanodomains (Leo-Macias et al., 2016; Mezache et al., 2020) within the different ID regions (Rhett et al., 2012; Veeraraghavan et al., 2015). Incorporating these details is a focus of future work.
We note that the differences between our modeling approach incorporating ID structure and the standard monodomain model over a wide range of GJ coupling strengths counters the argument that gap junctional coupling is solely governing conduction in one regimen, while ephaptic coupling is solely governing conduction in another regimen. Our results suggest that the relationship between these forms of coupling is more nuanced, in the sense that neither is a “switch” that is turned on or off for certain conditions, but rather both mechanisms operate in tandem to maintain robust conduction under different conditions. Specifically, for strong GJ coupling, ephaptic coupling tends to slow conduction, however with relatively weak sensitivity. In contrast, for weaker GJ coupling, ephaptic coupling tends to enhance conduction with greater sensitivity. Similarly, for strong ephaptic coupling, conduction is less sensitive to GJ coupling, while for weak ephaptic coupling, conduction is more sensitive to GJ coupling. Importantly, the relative roles of these forms of coupling are interconnected and dependent.
As noted above, experimental measurements of macroscopic GJ conductance under physiological conditions vary considerably, typically on the order of 10s to 100s of nanosiemens (Desplantez et al., 2007; Kwak and Jongsma, 1996; McCain et al., 2012; Moreno et al., 1994; Valiunas et al., 2002; Verheule et al., 1997; White et al., 1990; Nielsen et al., 2012; Rüdisüli and Weingart, 1989; Weingart, 1986; Wittenberg et al., 1986), such that our simulations for weak GJ coupling correspond with the lower end of this range and those for strong GJ coupling correspond with the upper end of this range. While this variability in measurements may be partially due to differences in species, cell type, and experimental conditions, the individual studies also typically found a wide range in GJ conductance measurements. Additionally, GJ protein expression has been shown to be reduced in the diseased myocardium, in cases of heart failure (Akar et al., 2004, 2007) or myocardial ischemia (Cascio et al., 2005). Collectively, these data suggests that the range of GJ conductance investigated in this study is within a physiologically relevant regimen, with reduced GJ coupling tending to represent more pathological conditions. Our simulation results, particularly for the weak GJ coupling case, may help resolve apparent discrepancies between the time course of GJ remodeling and conduction slowing in the failing heart (Akar et al., 2004, 2007), where edema is known to play a role (Boyle et al., 2007). We therefore expect that CV sensitivity to interplicate intermembrane separation (shown in Fig. 7) is important in modeling these pathologies, especially in a heterogeneously affected tissue. One important direct extension of our current work is to investigate mechanisms of electrical dysfunction in failing hearts. Since we find that interplicate membrane separation has a greater impact on conduction for weaker GJ coupling, our modeling approach should be particularly well suited for modeling pathological cases (such as heart failure or ischemia), for which both GJ coupling and overall ID structure are expected to be perturbed. Importantly, our approach is sufficiently flexible to incorporate ID structural data from failing hearts into FEM and tissue models, and such mechanistic studies are a key focus of future work. The present study is a critical step toward such work, enabling us to develop the necessary imaging, image analysis, and modeling workflows and pipelines.
A number of prior studies quantify ID structure. Our measurements of overall ID size (Bennett et al., 2006; Wilson et al., 2014), intermembrane spacing and GJ cluster size (Severs, 1990), and plicate folds and amplitude (Vanslembrouck and Kremer, 2020) are quite similar to previously published data. In contrast with our approach, other studies have characterized the entire 3-D structure of the ID using 3-D stacking and image segmentation (Leo-Macías et al., 2015; Vanslembrouck et al., 2018; Vanslembrouck and Kremer, 2020). We note that our approach of mesh partitioning and calculating equivalent cleft conductances, which could be then incorporated into the cleft network tissue model, can be applied to a 3-D FEM model based on a 3-D ID reconstruction. However, our approach is more flexible and has a number of specific advantages over a direct 3-D reconstruction. First, we are able to directly modify specific features of the ID mesh geometry, to investigate the role of such features on conduction. A comparable 3-D reconstruction would be required for each experimental case, and as 3-D stacking and segmenting images is quite time-consuming and challenging, this would greatly limit the conditions investigated. Second, our approach can determine ID geometry based on hundreds of images from multiple hearts and IDs, thus accounting for physiological variability, while a 3-D reconstruction provides information only on the specific ID being imaged. Finally, our approach is not limited to TEM images, as we could derive FEM parameters from any imaging modality, provided that there is sufficient resolution to characterize the relevant ID structures.
In a prior modeling study, Hichri et al. (2018) simulated two equally spaced apposing cell membranes, separated by an intercellular cleft space modeled using a uniform FEM mesh. They similarly predicted that cleft potential (ϕc) is variable within the cleft, specifically with larger-magnitude hyperpolarization toward the center (a result also described by Mori et al.  in a 3-D radially symmetric model). This variability in turn drives heterogeneous activation of the postjunctional membrane, similar to our results in Fig. 8, and the authors show that these effects are amplified by both decreased intermembrane distance and Na+ channel clustering. We similarly find that the cleft is more polarized toward the center (Fig. 6). We further predict that the potential distribution is in fact even more variable and heterogeneous due to the irregular structure of the ID meshes and variations in structure that arise due to differences in intermembrane separation in the plicate and interplicate and the presence of GJs and plicate folds. Additionally, our approach to derive an equivalent reduced electrical circuit to represent the cleft conductances eliminates the computationally costly step of simulating electrochemical dynamics on the full 3-D FEM mesh, facilitating simulation of cardiac tissue and predictions of CV.
Interestingly, the mechanism of conduction slowing/enhancing shown in Fig. 8 is similar to the single-cleft model previously predicted by Kucera et al. (2002): for strong GJ coupling, narrow intercellular clefts enhance cleft polarization and thereby reduce Na+ current driving force, which slows conduction, a mechanism described as “self-attenuation.” In contrast, for weak GJ coupling, narrow clefts activate the postjunctional membrane earlier and enhance conduction. Our study builds on this prior work and predicts that cleft heterogeneity also contributes to these mechanisms. Specifically, for strong GJ coupling, self-attenuation and reduced Na+ current also arise due to desynchronized activation of the postjunctional membrane. For weak GJ coupling, narrow clefts enhance cleft polarization heterogeneity that facilitates an overall earlier activation of postjunctional Na+ channels. Additionally, this earlier activation in one region of the postjunctional ID membrane in turn contributes to the activation of additional Na+ channels in other regions, similar to the mechanism of “self-activation” described by Hichri et al. (2018). In addition to illustrating regulation by ID structural heterogeneity, our study illustrates that these mechanisms of conduction regulation can be mediated by intermembrane separation changes in a specific region of the ID, specifically the interplicate. Our ongoing work is investigating how additional ID structural changes impact conduction.
We conclude by acknowledging limitations of our study. Although our model formulation incorporates significant subcellular details in the representation of cardiac tissue, this representation is still a simplification of the complex cardiac tissue structure. Specifically, our tissue model represents a 1-D strand of coupled cells and thus does not account for the 3-D tissue structure of the heart. Importantly, however, our model formulation does represent key subcellular and structural details, specifically nonuniform Na+ channel distribution and ID structure. Further, our approach can modified to account for specific ID structural heterogeneity that is not specifically represented in TEM images. One such example is large GJs that have been identified toward the periphery of the ID (Gourdie et al., 1990, 1991). Although we did not measure these directly in our TEM images, we can investigate their potential role in modulating conduction. To illustrate this approach, we created a FEM mesh with a distribution of larger GJs on the ID periphery (Fig. S9), and interestingly found that these larger GJs reduced connection to the extracellular bulk, and thus amplify the CV sensitivity to changes in intermembrane distance. More broadly, our work can be further extended to incorporate structural details from other imaging modalities, e.g., ion channel clustering within GJ- and MJ-associated nanodomains, as noted above.
Further studies are also needed to investigate the role of additional electrogenic proteins that have been identified to be preferentially localized at the ID, in particular K+ channels (Hong et al., 2012; Veeraraghavan et al., 2015; Vermij et al., 2017); here we focused on the role of Na+ channels and their role in regulating conduction at the ID. Finally, while accounting for the irregular structure of the ID membrane, we assume that the lateral membrane follows an overall cylindrical shape, and prior work has shown that microscale heterogeneity, in particular in cell shape/geometry and GJ protein distribution, can impact conduction delays and lead to heterogeneous conduction (Spach et al., 2000, 2007; Hubbard and Henriquez, 2012). We are interested in how these nano- and microscale heterogeneities can influence each other and impact conduction, in particular in pathological settings and associated structural remodeling. Overall, while these limitations may account for quantitative differences in model predictions, importantly our study demonstrates a novel mechanism in which ID nanoscale structural heterogeneity modulates ID-localized Na+ channels and regulates cardiac conduction.
David A. Eisner served as editor.
This study was supported by funding from National Institutes of Health grant numbers R01HL138003 (to S.H. Weinberg) and R01HL148736 (to R. Veeraraghavan).
The authors declare no competing financial interests.
Author contributions: N. Moise, H.L. Struckman, R. Veeraraghavan, and S.H. Weinberg designed the research studies. N. Moise and S.H. Weinberg performed the numerical simulations. H.L. Struckman conducted the microscopy experiments. H.L. Struckman and C. Dagher performed image analysis. N. Moise, H.L. Struckman, R. Veeraraghavan, and S.H. Weinberg wrote the manuscript. All authors revised, reviewed, and approved of the final version of the manuscript.