Role of ephaptic coupling in discordant alternans domain sizes and action potential propagation in the heart

Abstract

Discordant alternans, the spatially out-of-phase alternation of the durations of propagating action potentials in the heart, has been linked to the onset of fibrillation, a major cardiac rhythm disorder. The sizes of the regions, or domains, within which these alternations are synchronized are critical in this link. However, computer models employing standard gap junction-based coupling between cells have been unable to reproduce simultaneously the small domain sizes and rapid action potential propagation speeds seen in experiments. Here we use computational methods to show that rapid wave speeds and small domain sizes are possible when a more detailed model of intercellular coupling that accounts for so-called ephaptic effects is used. We provide evidence that the smaller domain sizes are possible, because different coupling strengths can exist on the wavefronts, for which both ephaptic and gap-junction coupling are involved, in contrast to the wavebacks, where only gap-junction coupling plays an active role. The differences in coupling strength are due to the high density of fast-inward (sodium) channels known to localize on the ends of cardiac cells, which are only active (and thus engage ephaptic coupling) during wavefront propagation. Thus, our results suggest that this distribution of fast-inward channels, as well as other factors responsible for the critical involvement of ephaptic coupling in wave propagation, including intercellular cleft spacing, play important roles in increasing the vulnerability of the heart to life-threatening tachyarrhythmias. Our results, combined with the absence of short-wavelength discordant alternans domains in standard gap-junction-dominated coupling models, also provide evidence that both gap-junction and ephaptic coupling are critical in wavefront propagation and waveback dynamics.

FIG. 1.

Circuit used to model a one-dimensional fiber containing gap-junction and ephaptic intercellular coupling. Black circuit elements are those typically used in standard monodomain models of one-dimensional fibers. Red elements model ephaptic coupling. Blue dashed boxes indicate the locations of cells within this circuit description of the fiber.

FIG. 2.

Wavefront and waveback arrival times (nearly straight green and distinctly curved red traces, respectively) vs arrival location $x$ along the fiber for (a) the GJ system, (b) the EC1 system, and (c) the EC2 system. Waves were launched from the left end of the system at regular intervals. The pacing interval was 207 ms. Data are shown approximately 120 s into the simulations, well after initial transient behavior has dissipated. (d)–(f) $\mathrm{A}\mathrm{P}\mathrm{D}\left(x\right)$ vs $x$ for the last two complete waves shown in (a)–(c) (the last wave is shown by the solid trace and the penultimate wave by the dashed trace). The wavefront and waveback arrival times were defined to be the times the membrane potential $V$ ascended or descended to the level of −65 mV, respectively. The APD was defined to be the difference between these two times.

FIG. 3.

Currents flowing through the gap-junction resistor $R_g$ from left to right in Fig. 1 (blue dashed line) and through the radial resistor $R_r$ to the ground (red solid line) in between cells 194 and 195 during wavefront propagation, in (a) the GJ system and (b) the EC2 system.

FIG. 4.

Action potential wave velocities vs cleft width $w_{\mathrm{c}\mathrm{l}}$ for different strengths of the fast-inward current, as characterized by ${\tau }_{\mathrm{f}\mathrm{i}}^{-1}$ (3.75 ms⁻¹, weakest, and 5.00 ms⁻¹ strongest), for the case of $R_g=395\mathrm{M}\Omega$ and 90% of fast-inward current channels localized on the ends of the cells.

FIG. 5.

(a) Conduction velocity and (b) APD restitution functions for the three systems, GJ (green line), EC1 (blue line), and EC2 (red line), for a cycle length (BCL) of 207 ms. Nonlinear regression was performed to construct the curves shown in both panels, assuming functions of the form $f\left(\mathrm{D}\mathrm{I}\right)={\beta }_1-{\beta }_2\mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{D}\mathrm{I}/{\beta }_3\right)$ in both cases. Additionally, scatterplot data from the three simulations are shown in (b). Data for the GJ and EC1 simulations are hidden underneath the EC2 data. The black dots show the location of the steady-state points in each plot: (a) $(\stackrel{-}{\mathrm{D}\mathrm{I}},\stackrel{-}{\mathrm{C}\mathrm{V}})$ and (b) $(\stackrel{-}{\mathrm{D}\mathrm{I}},\stackrel{-}{\mathrm{A}\mathrm{P}\mathrm{D}})$. The line $\stackrel{-}{\mathrm{A}\mathrm{P}\mathrm{D}}=\mathrm{B}\mathrm{C}\mathrm{L}-\stackrel{-}{\mathrm{D}\mathrm{I}}$ is also shown in (b).

FIG. 6.

Key discordant alternans domain parameters vs pacing cycle length. A description of the calculations of (a) the number of domains $N_{\text{domains}}$, (b) the domain size $l_{\text{domain}}$, (d) the APD coupling length $\xi$, and (e) the node width $l_{\text{node}}$ is provided in Appendix B. The theoretical domain size in (c) is calculated from the formula $(2\pi /\sqrt{3}){\left(2{\xi }^2\Lambda \right)}^{1/3}$ (see the text).

FIG. 7.

Spread of the simulation data about the presumed relationship expressed by Eq. (13), due to the APD’s nonuniformity in $x$, for the EC2 system, with BCL=207 ms. (a) Scatterplot of $\delta \mathrm{A}\mathrm{P}\mathrm{D}\equiv \mathrm{A}\mathrm{P}\mathrm{D}\left(x\right)-a\left(\mathrm{D}\mathrm{I}\right(x\left)\right)$ vs $da\left(\mathrm{D}\mathrm{I}\right(x\left)\right)/dx$ and $d^{2}a\left(\mathrm{D}\mathrm{I}\right(x\left)\right)/d{x}^2$ (in turquoise blue). The plane shown is the least-squares fit to the scatterplot data, obtained using the method described in Appendix B 3. Also shown are the scatterplot data and least-squares fit from (a) plotted vs (b) $da\left(\mathrm{D}\mathrm{I}\right(x\left)\right)/dx$ and (c) $d^{2}a\left(\mathrm{D}\mathrm{I}\right(x\left)\right)/dx$, with dependence on the other derivative subtracted out. Also included are the scatterplots (in blue) of the (d) wave velocity and (e) APD, vs DI. Red curves are the fitted APD and velocity restitution functions.

FIG. 8.

Plots of $V_m$ of a single action potential excited simultaneously at all points on the fiber at time $t=0$, but with slightly different initial conditions in each half of the fiber: $h=0.55$ and $f=0.82$ for $x⩽1.6\mathrm{c}\mathrm{m}$, and $h=0.50$ and $f=0.80$ for $x>1.6\mathrm{c}\mathrm{m}$, for (a) the GJ system, (b) the EC1 system, and (c) the EC2 system. (d)–(f) Close-up of the location of the waveback (in red) for the conditions in (a)–(c) near the center of the fiber.