## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Ephaptic conduction in a cardiac strand model with 3D electrodiffusion

Contributed by Charles S. Peskin, February 6, 2008 (received for review July 29, 2007)

## Abstract

We study cardiac action potential propagation under severe reduction in gap junction conductance. We use a mathematical model of cellular electrical activity that takes into account both three-dimensional geometry and ionic concentration effects. Certain anatomical and biophysical parameters are varied to see their impact on cardiac action potential conduction velocity. This study uncovers quantitative features of ephaptic propagation that differ from previous studies based on one-dimensional models. We also identify a mode of cardiac action potential propagation in which the ephaptic and gap-junction-mediated mechanisms alternate. Our study demonstrates the usefulness of this modeling approach for electrophysiological systems especially when detailed membrane geometry plays an important role.

Cellular electrical activity is central to physiology. It is the basis of sensory perception, communication between neurons, initiation and coordination of skeletal muscle contraction, synchronization of the heart beat, and the secretion of hormones (1). Most mathematical models of cellular electrical activity are based on the cable model, which can be derived from a current continuity relation on a one-dimensional ohmic cable (2⇓–4). As such, its derivation rests on several assumptions (4): ionic concentrations are assumed not to change appreciably over the time of interest, and a one-dimensional picture of cell geometry is assumed to be adequate for purposes of describing cellular electrical activity. These assumptions, however, may not hold in many systems of biological significance, especially in the central nervous system and cardiac tissue, where microhistological features may play an essential role in shaping physiological responses (4, 5).

Previous attempts at generalizing the cable model include refs. 6 and 7. In ref. 6, the authors include electrodiffusion of ions in their framework while retaining the one-dimensional character of the cable model. Their pioneering work seems not to have seen widespread use, possibly because of the numerical difficulty the authors encountered in integrating their system (4). In ref. 7, the authors propose a model that addresses the limitations of the cable model discussed in the previous paragraph. Their attention is, however, restricted mostly to stationary problems (8). Their model is numerically difficult to deal with because of the need to resolve boundary layers associated with space charge layers near the membrane.

In refs. 9 and 10, we proposed a model that addresses the limitations of the cable model as stated above and developed an efficient numerical algorithm that makes dynamic simulations possible. We view biological tissue as 3D space partitioned into intracellular and extracellular spaces by the membrane. In these regions, ionic concentrations and the electrostatic potential obey a system of partial differential equations (PDEs). Transmembrane ionic currents are incorporated in the form of interface conditions to be satisfied on both sides of the membrane.

In this article, we shall apply this modeling methodology to the study of cardiac action potential propagation when gap junction conductance between neighboring cells is severely suppressed. It is widely accepted that gap junctions, which serve as low-resistance connections between cells, are essential to cardiac action potential propagation (11). There has, however, been some controversy as to whether gap junctions are absolutely necessary (12). This controversy has been given renewed attention in recent years in light of the experimental finding that mice engineered not to express gap junctions in the heart (connexin 43 knockout mice) still exhibit cardiac electric conduction, albeit at a reduced velocity (13).

One possible hypothesis that may explain such anomalous conduction is the ephaptic mechanism, first proposed in ref. 14. Simple models of cardiac electrical circuitry have been used to argue that cardiac propagation is possible without gap junctions (15). In ref. 16, the authors address this question using a simple ohmic current model with realistic cardiac ion channel dynamics.

We apply our modeling methodology to perform a more complete study of this issue. As in ref. 16, we consider a linear strand of cardiac cells. Several biophysical parameters are varied to assess their impact on cardiac action potential propagation. Our methodology allows us to explore the effects that membrane geometry, extracellular space, and ionic diffusion have on ephaptic propagation, aspects that cannot be addressed by using conventional models of electrophysiology.

## Model

We give a brief introduction to our mathematical model. For a detailed discussion, we refer the reader to refs. 9 and 10. Let the biological tissue of interest be divided into subregions Ω^{(k)}, indexed by *k*. We denote the membrane separating the regions Ω^{(k)} and Ω^{(l)} by Γ^{(kl)} (Fig. 1).

In Ω^{(k)}, the equations satisfied by the ionic concentration *c _{i}* and the electrostatic potential φ are the following:
Here,

**f**

*denotes the flux of the*

_{i}*i*th species of ion, and Eq.

**1**is a statement of ionic conservation.

**f**

*is expressed as a sum of a diffusive and a drift flux.*

_{i}*D*is the diffusion coefficient of the

_{i}*i*th species of ion, and

*qz*is the amount of charge on the

_{i}*i*th species ion, where

*q*is the elementary charge—i.e., the charge on a proton. The coefficient

*qz*/(

_{i}D_{i}*k*

_{B}

*T*) is the mobility of the ionic species (Einstein relation), where

*k*

_{B}is the Boltzmann constant and

*T*is the absolute temperature. Eq.

**2**is a statement of electroneutrality, and ρ

_{0}is the immobile background charge density, which may come from charge contributions from extracellular or cytoskeletal matrix proteins. An alternative would be to replace Eq.

**2**with the Poisson equation, which would then constitute the familiar Poisson–Nernst–Planck system (7). This system, however, is very difficult to simulate numerically because of the presence of Debye layers (to be discussed below).

We now turn to the boundary conditions, satisfied at both the intracellular and extracellular sides of the cell membrane. Across the cell membrane, a jump in electrostatic potential (membrane potential) is maintained, and the cell membrane acts as a capacitor. There is thus a thin layer on both sides of the membrane where electric charge accumulates whose thickness is on the order of the Debye length *r*_{d}, which is approximately 1 nm in physiological systems. In Eq. **2**, we have taken the electroneutrality condition to hold inside and outside of the cell, and this implies that we must treat the presence of Debye layers in the form of boundary conditions.

Ionic current that flows into the membrane either goes across the membrane through ionic channels, transporters, or pumps, or contributes to the change in surface charge. The boundary condition satisfied at the Ω^{(k)} side of the membrane Γ^{(kl)} is
The term *j*_{i}^{(kl)} denotes the transmembrane current from region Ω^{(k)} into Ω^{(l)}. We note by definition that *j*_{i}^{(kl)} = −*j*_{i}^{(lk)}. The variable σ_{i}^{(k)} denotes the contribution of the *i*th species of ion to the surface charge density on the Ω^{(k)} face of the membrane.

The term σ_{i}^{(k)} is expressed as
where *D*_{0} = 1.0 μm^{2}/ms is a typical ionic diffusion coefficient and *r*_{d} is the aforementioned Debye length. σ^{(k)} is the total charge per unit area on the membrane surface facing Ω^{(k)} and is the product of *C*_{m}, the capacitance of the membrane, and φ^{(kl)} = φ^{(k)} − φ^{(l)}, the membrane potential. Note that λ_{i}^{(k)} is the fractional contribution of the *i*th species of ion to the surface charge density on face *k* of the membrane (Eq. **4**). The quantity *r*_{d}^{2}/*D*_{0} is the time scale in which λ_{i}^{(k)} relaxes to λ̃_{i}^{(k)}. The form of λ̃_{i}, and also the differential equation that related λ_{i} to λ̃_{i}, as given in Eq. **5** can be obtained by using a matched asymptotic analysis at the membrane as derived in ref. 10. The time scale *r*_{d}^{2}/*D*_{0}, which is on the order of 1 ns, is in fact so fast that λ_{i}^{(k)} ≈ λ̃_{i}^{(k)} for physiological time scales of interest.

The term *j*_{i}^{(kl)} denotes transmembrane currents that flow through ion channels, transporters, or pumps that are located within the cell membrane (17⇓–19). We use the formalism introduced by Hodgkin and Huxley for ion channel currents (2⇓–4), generalized to allow for nonlinear instantaneous current–voltage relations and ion concentration effects. We refer the reader to supporting information (SI) *Text*, Table S1, and Movies S1–S9 for details.

Eqs. **1** and **2** with the interface condition **3** constitute the full system of equations. This is a differential algebraic system in which the electrostatic potential evolves so that the electroneutrality condition **2** is satisfied at each instant. We may derive an elliptic equation that is satisfied by φ, by taking the derivative of **2** in time *t* and substituting **1** for ∂*c*_{i}/∂*t*:
This may be interpreted as a current conservation relation where *a*∇φ is the ohmic contribution and ∇*b* is the diffusive contribution to the electric current. We may identify *a*(**x**, *t*) with the electrolyte ohmic conductivity. We note that the equation satisfied by φ is Eq. **6**, and not the Laplace equation. The electroneutral limit is one in which a small imbalance of large absolute ionic charge densities results in a nonzero Laplacian of the electrostatic potential (see ref. 9 for further discussion of this subtle issue).

In the simulations to follow, we apply the above model to cardiac geometries in which the smallest dimension is the width of an intercellular cleft that we vary from 2 to 50 nm. Since the Debye length is ≈1 nm, this raises the question as to whether electroneutrality is a valid approximation within such a narrow cleft, especially at the lower end of the thickness range that we consider. It turns out, however, that the electroneutral model is still a remarkably good approximation even at this length scale. We have demonstrated this by both asymptotic analysis and high-resolution one-dimensional computations (see refs. 9 and 10).

## Application to Cardiac Physiology

Cardiac muscle consists of a network of cardiac muscle fibers. Each cardiac fiber may be seen as a linear strand of cardiomyocytes, separated from one another by a thin gap known as the intercalated disk. The intercalated disk is believed to be the site of both electrical and mechanical coupling between cells. Gap junctions are primarily located at the intercalated disk, providing a low-resistance pathway of electric current between the cells facing this gap. As discussed in the Introduction, it is not clear as to whether gap junctions are absolutely necessary for cardiac action potential propagation, especially given that connexin 43 knockout mice exhibit cardiac electric conduction (13). In this section, the foregoing model is used to explore the consequences of severely reduced gap junction conductance and the preferential localization of Na^{+} channels to the intercalated disk in a model cardiac muscle fiber.

### Model.

We model a cardiac fiber as a strand of *N*_{c} cardiac cells all of which are assumed to be cylindrical in shape (Fig. 2). We take the radius of the circular cross-section of the cell to be *l*_{R} = 24.7/2 μm and the length of the cell to be *l*_{A} = 157.9 μm (experimentally measured values). Place the *N*_{c} cells so that the cylindrical axes of all of the cells lie along a single line. Take this line to be the *z* axis, and the radial coordinate to be the *r* axis. Label the cells *k* = 1 … *N*_{c} in order of increasing *z* coordinate. Only half of cell 1 and *N*_{c} are included in the computational domain. The intracellular region of cells 1 and *N*_{c} thus meet the outer boundary of the computational domain at the middle of the cell. The width of the gap between cells is *l*_{g}, which we vary in our simulations. The computational domain is taken to be a cylindrical domain of radius (1 + η)*l*_{R} with the *z* axis as the center line. Since the *N*_{c} cells have radius *l*_{R}, the extracellular space corresponds to the region *l*_{R} < *r* < (1 + η)*l*_{R} (which we shall call the extracellular bath) as well as the gaps between cells. We take η = 1 unless indicated otherwise.

We seek radially symmetric solutions, thus obviating the necessity to introduce an angular coordinate. No-flux boundary conditions are imposed at the outer boundary of the computational domain. We consider four ion types in the calculation, Na^{+}, K^{+}, Ca^{2+}, and Cl^{−}. Initial conditions are set to satisfy electroneutrality (see *SI Text* for details). To facilitate comparison with experimental data on mice (13, 20), we shall use the mouse cardiac model of Bondarenko *et al.* (21) as the ion channel model in our simulations with the following modifications. We do not include intracellular calcium handling in the Bondarenko model, since this would require a detailed geometric model of intracellular organelles. In addition, we do not include the background Na^{+} conductance and the background Ca^{2+} conductance, which we have seen induce unwanted spontaneous membrane potential oscillations. The gap junctions are modeled as cytoplasmic pores, the details of which can be found in *SI Text*.

In the following sections, we vary several parameters to examine their effects on cardiac action potential propagation. First, we vary the total gap junction conductance per gap, *G*^{gap}, within the experimentally observed range (20). The next parameter of interest is the gap width *l*_{g}. The exact width of the gap is unknown, and we vary this parameter in the range 2–50 nm. As documented in refs. 16 and 22, recent evidence suggests the preferential localization of Na^{+} channels on the membranes facing the intercellular gap. Denote the proportion of Na^{+} channels expressed in the gap by *p*_{Na}. Following ref. 16, we vary *p*_{Na} [from *p*_{Na} = uniform (no redistribution) to *p*_{Na} = 0.99 in which 99% of Na^{+} channels reside on the membrane facing the gap], while keeping the total number of Na^{+} channels on each cell constant. Finally, we vary η, which controls the size of the extracellular bath.

As documented in refs. 9 and 10, the equations we must simulate are numerically stiff. We use a finite volume discretization in space and an implicit discretization in time to perform the simulations to follow. The details may be found in ref. 10.

### Normal Conduction.

We first tested the model for normal conduction. According to ref. 20, the intercellular conductance along the fiber direction, which we identify as the conductance across gap junctions over the intercalated discs, is *G*^{gap} = *G*_{normal} ≡ 5.58 × 10^{−4} mS. In this and all ensuing simulations, gap junction conductance is distributed uniformly over the membrane facing the intercalated disk so that the gap junction conductance per unit area, *g*^{gap} is given by *g*^{gap} = *G*^{gap}/π*l*_{R}^{2} mS/μm^{2}. We initiate an action potential by transiently increasing the Na^{+} conductance of cell 1.

The conduction velocities (CVs) when *p*_{Na} and *l*_{g} are varied are plotted in Fig. 3*a* (see Movie S1 for a movie of the propagating action potential). When *p*_{Na} = uniform, we see that CV is insensitive to gap width *l*_{g} at ≈36 cm/s. When we redistribute Na^{+} channels so that the density is higher facing the gap, CV decreases as *l*_{g} decreases, and this decrease is greater for higher values of *p*_{Na}.

There can be at least two factors that contribute to the decrease in CV with *l*_{g}. The first is that concentration changes can decrease the driving force for Na^{+} channel currents, which are responsible for membrane depolarization. The narrower the gap, the quicker the depletion of Na^{+} ions in the gaps, thereby reducing the equilibrium potential for Na^{+} ions. The second factor is that the greater resistance of the narrower gap makes it difficult for current to flow through channels facing the gap. This will make the channels facing the gap less effective, thereby decreasing CV. This drop in CV as *l*_{g} decreases is also documented in ref. 16.

According to ref. 13, CV under normal conditions is approximately 40 cm/s in the transverse direction and 60 cm/s in the longitudinal direction. Given that the ion channel model of Bondarenko *et al.* (21) was only calibrated to voltage clamp data, the simulated propagation speed of approximately 30 cm/s may be considered relatively close to the experimentally observed value. We shall henceforth express the simulated CVs in percentages with respect to this value (30 cm/s) as well as in centimeters/second.

The source of the discrepancy between the computed and physiological CVs may be our simplification of taking only a single strand of cardiac cells, thus ignoring the 3D arrangement of cardiomyocytes. In a true cardiac preparation, electric current can go through many pathways to get from one cell to another, thereby reducing the effective resistance between two cells.

### Conduction with Reduced Gap Junction Conductance.

We now take a detailed look at conduction when gap junction conductance is severely reduced. According to ref. 20, the gap junction conductance at the intercalated disk space for connexin 43 (dominant gap junction expressed in cardiac tissue) knockout mice is *G*^{gap} = *G*_{reduced} ≡ 1.10 × 10^{−5} mS. Thus, *G*^{gap} is now reduced to ≈2% of the normal value *G*_{norm}.

In the simulations shown in Fig. 4(see Movies S2–S6), we have taken *l*_{g} = 5 nm and *p*_{Na} = 0.95. We initiate an action potential by transiently increasing the Na^{+} conductance to the membrane of cell 1. We find successful cardiac action potential propagation at ≈10 cm/s (33%). The sequence of events that underlies this propagation is as follows. Consider the propagation of the action potential from cell 1 to cell 2. A depolarizing current activates Na^{+} channels in cell 1 (*t* = 0.8 ms in Fig. 4). The opening of these Na^{+} channels, which are preferentially expressed on the membranes facing the gap, generates a strong current flowing into the gap from the extracellular bath. Since the gap is narrow, a large negative deflection in the extracellular voltage is seen within the gap. We see that there is a voltage gradient within the gap from *r* = *l*_{R} to *r* = 0. Therefore, the voltage at *r* = 0 is most negative with respect to the voltage in the extracellular bath (*t* = 1.6 ms in Fig. 4).

The Na^{+} channels facing the gap on cell 2 sense a depolarized membrane potential because the voltage in the gap is negative with respect to the extracellular bath. This activates the Na^{+} channels on cell 2, leading to its depolarization (*t* = 2.4 ms in Fig. 4). Note that the membrane of cell 2 facing the middle of the gap will experience a greater depolarization due to the voltage gradient that develops in the gap. This cycle of events repeats itself between cells *k* and *k* + 1, resulting in a propagating action potential (*t* = 3.2 ms in Fig. 4). The increase in K^{+} concentration in the gap may also facilitate conduction across the gap by pushing the resting potential of the membrane to a more positive value. Movies of concentration changes of Na^{+}, K^{+}, Ca^{2+}, and Cl^{−} are shown in Movies S3–S6.

### Model Comparison.

In this section, we compare the differences in prediction between different models of cellular electrical activity. We compare CV as a function of gap width, as calculated with the full 3D model, a 3D model in which ionic concentrations are ignored (3D cable model), a 1D model in which ionic concentrations changes are incorporated, and the 1D cable model. The 1D model calculations were performed by taking one mesh for 0 < *r* < *l*_{R} and one mesh for *l*_{R} < *r* < (1 + η)*l*_{R}. In these computations, we let *G*^{gap} = *G*_{reduced} and *p*_{Na} = 0.95. We refer the reader to ref. 10 for a detailed discussion of the above models and their interrelation.

In Fig. 5*a*, CV is plotted against different values of *l*_{g}. There is a significant difference between the profile of the traces between the 1D and the 3D models. The 1D models significantly overestimate CV when the gap width is greater and yield underestimates when the gap width is smaller.

The reason for the discrepancy may be explained as follows. In a 1D model, a single compartment is used to represent the cleft. In this case, the effective resistance for electric current to enter the cleft is overestimated since a single compartment model treats all positions within the cleft as being at the same effective distance away from the lateral membrane. When the cleft is discretized along the radius as is done in the 3D model, the current can easily enter to positions close to the lateral membrane since a voltage gradient can develop. Thus, the effective resistance for the electric current to gain access to the interior of the gap is smaller in a 3D model. In the 1D model, it is more difficult to excite the cleft potential, but once excited, the decay is slower. This is the principal reason why the 3D model shows a peak in CV when *l*_{g} is smaller. This explanation suggests that a larger value of *l*_{R} will accentuate the difference between the 1D and 3D models, which is indeed the case (data not shown).

The peak CV is reached at approximately a 5-nm gap in the 3D model calculations and at approximately a 30-nm gap in the 1D model calculations. It is interesting that in ref. 16, in which a 1D model is used, the authors report peak CVs at approximately *l*_{g} = 30–100 nm, which comes close to the value we obtain using a 1D model. We note, however, that it is not possible to make a straightforward comparison between our results and the results in ref. 16, since many of the parameters, including those of the ionic channel model, are different.

Models with or without ionic concentration effects produce almost identical results. We may conclude that ionic concentration changes play a relatively minor role. This is in line with discussions in ref. 16, where the authors make estimates on the order of magnitude of the effect ionic concentration changes may have on ephaptic coupling. Note, however, that there is a slight deviation between the full 3D and 3D cable models when *l*_{g} is smaller than 5 nm. This can be attributed to the rapid depletion of Na^{+} ions in the gap, which has the effect of reducing the driving force for Na^{+} channels. Thus, we do see a concentration effect, and the full 3D model produces a CV slightly lower than the 3D cable model.

To accentuate the difference with computations that ignore ionic concentration, we performed computations in which we introduce a nonuniform fixed negative charge density within the gap space (see *SI Text* for the exact form of ρ_{0} used). In this case, concentration effects are substantial; see Fig. 5*b*, in which results differ depending on whether the model allows for concentration changes. When there is some fixed negative charge in the gap, a diffusion potential develops in the gap. Without any external stimulus, the electrostatic potential in the gap is negative with respect to the extracellular bath. This tends to accelerate the propagation of the action potential (see Movies S2 and S7).

### Varying Parameters.

First of all, we study the effect of *G*^{gap} on CV as *l*_{g} is varied while keeping *p*_{Na} = 0.95. The most notable point about Fig. 6*a* is that propagation fails completely when *l*_{g} is ≥7 nm if there are no gap junctions connecting two adjacent cells. The propagation velocity is seen to increase with greater gap junction density at the gap, and at *G*^{gap} = 8 × *G*_{reduced}, CV can reach ≈18 cm/s (60%).

We next study the effect of *p*_{Na} and *l*_{g} on CV, while fixing *G*^{gap} = *G*_{reduced}. The plot of CV is given in Fig. 3*b*. We first note that even when *p*_{Na} is low, we do see action potential propagation albeit at a reduced velocity (3∼4 cm/s). This is because cardiac action potentials have a long plateau phase during which the cell is depolarized, and a small gap junctional conductance is enough to inject sufficient current to the next cell to induce depolarization. CV increases with an increase in *p*_{Na}. For each fixed value of *p*_{Na}, CV does not increase (or decrease) monotonically with *l*_{g}. The velocity attains its maximum at approximately *l*_{g} = 5 nm, and this value is at most ≈12 cm/s (40%). The ratio is comparable with that seen in experiments [≈36/62 = 58% in the experimental case (13)]. This nonmonotonic behavior can also be seen in the 1D computations plotted in Fig. 5*a* and is also reported in ref. 16. The explanation for this behavior, given in the following paragraph, is essentially the same as that given in ref. 16.

Suppose that cell A and cell B are separated by a narrow gap and an excited action potential is to propagate from cell A into cell B. As discussed earlier, a narrower gap implies a greater drop in the electrostatic potential in the gap with respect to the extracellular bath. This facilitates action potential propagation by making it easier for ion channels on cell B facing the gap to be activated. However, a narrower gap means that it is more difficult for electric current to flow from the extracellular bath into cell B through the ion channels facing the gap. Cell B will not fully depolarize unless a sufficient amount of current is injected into it. Thus, a gap too narrow will tend to slow down action potential propagation. The propagation velocity slowing effect dominates for smaller (<5 nm) gap widths, and the accelerating effect dominates for larger gap widths. At *l*_{g} = 5 nm, these effects balance and CV reaches a maximum. The length *l*_{g} = 5 nm is approximately equal to the distance between two membranes when bridged by gap junctions (16, 23) and may thus be considered the absolute minimum distance between two opposing membranes (16). It is interesting to note that the maximum CV is achieved at this gap width.

In Fig. 3*b* there is a small kink in the graph when CV is at ≈5–6 cm/sec. In Fig. 7, we plot the time points at which the action potential reaches cell *k* when *l*_{g} = 5, 12, or 30 nm while taking *p*_{Na} = 0.95. At *l*_{g} = 12 nm, we see a biphasic time series. Conduction is rapid between two cells *k* and *k* + 1 when *k* is even but is slow when *k* is odd.

As explained in the foregoing, at *l*_{g} = 5 nm, propagation can be attributed to the ephaptic mechanism. At *l*_{g} = 30 nm, the slow propagation is driven by the plateau phase of the cardiac action potential. At *l*_{g} = 12 nm, these two mechanisms play an equal role. Fast conduction between cells *k* and *k* + 1 when *k* is even is driven primarily by the ephaptic mechanism, and the slow conduction between cells *k* and *k* + 1 when *k* is odd is driven by the plateau phase mechanism (see Movies S8 and S9). The slower modes of propagation at *l*_{g} = 12 and 30 nm depend on L-type calcium current *I*_{Ca}(L), which is in turn influenced by intracellular calcium concentration [Ca^{2+}]_{int}. Given that our model does not incorporate the biophysical machinery of intracellular calcium handling, we tested whether the slower modes of propagation depend strongly on *I*_{Ca}(L) inactivation by [Ca^{2+}]_{int}. In place of [Ca^{2+}]_{int}, we used *f* × [Ca^{2+}]_{int}, 0.01 ≤ *f* ≤ 100 to inactivate *I*_{Ca}(L). Regardless of the value of *f*, we observed these modes of propagation.

In detail, propagation in the intermediate regime (*l*_{g} = 12 nm) proceeds as follows. Consider propagation from cell 1 to cell 2. Cell 1 depolarizes and a strong electric current rushes into the gap, generating a negative electrostatic potential in the gap with respect to the extracellular bath. This drop is not sufficient to depolarize cell 2 and soon dissipates without initiating a depolarization in cell 2. The action potential of cell 1, meanwhile, reaches a plateau phase. Since the electrostatic potential inside of cell 1 is higher than that inside of cell 2, electric current flows into cell 2, slowly depolarizing the cell. Cell 3 also receives current from cell 2, and its membrane potential slowly increases but at a slower rate than that of cell 2. Cell 2 eventually reaches threshold and is fully depolarized. Depolarization of cell 2 induces a negative drop in the electrostatic potential in the gap between cells 2 and 3. This time, the ephaptic mechanism is successful, and cell 3 rapidly depolarizes, unlike what happened between cells 1 and 2. The key difference is that cell 3 was primed during the slow depolarization phase of cell 2. This cycle repeats itself because cell 3 must now activate an unprimed cell 4.

Under reduced gap junction coupling, conduction is primarily due to the ephaptic mechanism when *l*_{g} is small, whereas the gap-junction-mediated mechanism dominates, despite the drastic reduction in gap junction conduction when *l*_{g} is larger. There is a transition zone between the two regimes, in which the ephaptic and gap-junction-mediated mechanisms alternate in propagating the action potential to the neighboring cell.

We finally study the effect of the size of the extracellular bath on CV (Fig. 6*b*). We see that CV tends to decrease with a decrease in η, when the gap width is set to *l*_{g} = 5 nm. This is because a large current flux into the cell is required for rapid conduction. When the extracellular space is limited, as may be the case in cardiac tissue, it becomes more difficult to draw enough current into the gap because of greater resistance.

## Conclusions

In this article, we studied the conduction of cardiac action potentials under severely reduced gap-junctional coupling. The simpler 1D studies analyzing this phenomenon (16) left unexplored the possibility that the details of the geometry and ionic concentration effects in the narrow intercalated disk may significantly alter the characteristics of the proposed ephaptic mechanism (12, 14). We address this issue using a mathematical model of cellular electrical activity that takes into account both 3D geometry and ionic concentration effects. We also use a mouse model to study this problem (21) to make direct comparison possible with experimental data (13, 20).

We demonstrated the possibility of cardiac action potential propagation under severely reduced gap-junction conductance, in agreement with earlier work (12, 16). Our 3D model, however, yields quantitatively different results compared with 1D model studies. We find that in the parameter ranges considered here, a narrower gap width (*l*_{g} = 5 nm) is required for ephaptic propagation to be successful. Nonetheless, the ratio of CVs under normal and reduced gap junctional coupling come reasonably close to the experimentally observed ratio (13). In the course of this study, we also identified a mode of conduction in which ephaptic and gap-junctional propagation alternate. This “mixed” mode of conduction may be searched for experimentally as a signature of ephaptic conduction.

There are several significant limitations in our current study. We took CV as our measure to gauge the correspondence between our computations and experiments performed in refs. 13 and 20. We note, however, that in ref. 13, optical techniques were used to map the propagation of the action potential front on a 2D cardiac surface. To make a better comparison with experimental data, we must go beyond the 1D caricature used here to better represent the 3D connectivity of cardiac cells.

Such a 3D model may help in resolving the discrepancy between the absolute values of the computed and experimentally observed CVs. Since a network of cells has higher connectivity in higher dimensions, it is possible that the ephaptic mechanism is more effective in three dimensions. Another possibility is that although residual gap junction conductance may be reduced on average, there may be a sufficient number of cell–cell links of relatively high gap junction conductance to allow for macroscopic conduction (conductance percolation). It would be important to computationally study such possibilities.

There is much improvement to be made at the level of single-cell geometry. The intercalated disk is not a simple planar plate but is tortuous, and the intermembrane distance may vary significantly from place to place (16). The extracellular matrix proteins that reside in the intercalated disk, the composition of which has yet to be completely characterized, may play a significant role in shaping the electric response, as was suggested in our computations. We have neglected the internal membrane structures present within cardiac cells, which are particularly important given their central role in the calcium dynamics of cardiac cells.

Although such complexities were not addressed in our present article, there is no conceptual difficulty in applying our modeling methodology to this more complex situation. This will, however, require better numerical methods and better anatomical, biochemical, and electrophysiological models that describe the system. This is a direction for future investigation.

## Acknowledgments

We thank the reviewers for helpful comments that improved the manuscript. Y.M. was supported by the McCracken Fellowship and the Dissertation Fellowship of New York University. Y.M. also acknowledges support through a Mathematics of Information Technology and Complex Systems (MITACS, Canada) Team Grant (to Leah Keshet) and a Natural Sciences and Engineering Research Council Discovery Grant (to Leah Keshet). This work was supported in part by National Institutes of Health Grant HL64757 (to G.I.F.).

## Footnotes

- ↵
^{‡}To whom correspondence may be addressed. E-mail: mori{at}math.ubc.ca or peskin{at}cims.nyu.edu

Author contributions: Y.M., G.I.F., and C.S.P. designed research; Y.M. performed research; Y.M., G.I.F., and C.S.P. analyzed data; and Y.M., G.I.F., and C.S.P. wrote the paper.

The authors declare no conflict of interest.

This article contains supporting information online at www.pnas.org/cgi/content/full/0801089105/DCSupplemental.

- Received July 29, 2007.

- © 2008 by The National Academy of Sciences of the USA

## References

- ↵.
- Guyton A,
- Hall J

- ↵
- ↵.
- Keener J,
- Sneyd J

- ↵.
- Koch C

- ↵.
- Shepherd GM

- ↵
- ↵
- ↵.
- Léonetti M,
- Dubois-Violette E,
- Homblé F

- ↵.
- Mori Y,
- Jerome JW,
- Peskin CS

- ↵.
- Mori Y

- ↵
- ↵.
- Sperelakis N

- ↵.
- Gutstein D,
- et al.

- ↵
- ↵.
- Sperelakis N,
- Ramasamy L

- ↵.
- Kucera J,
- Rohr S,
- Rudy Y

- ↵.
- Aidley D

- ↵.
- Hille B

- ↵.
- Kandel E,
- Schwartz J,
- Jessel T

- ↵.
- Yao J,
- Gutstein D,
- Liu F,
- Fishman G,
- Wit A

- ↵.
- Bondarenko V,
- Szigeti G,
- Bett G,
- Kim S,
- Rasmusson R

- ↵.
- Cohen S

- ↵.
- Bers D

## Citation Manager Formats

### More Articles of This Classification

### Physical Sciences

### Related Content

- No related articles found.

### Cited by...

- Microstructure, Cell-to-Cell Coupling, and Ion Currents as Determinants of Electrical Propagation and Arrhythmogenesis
- Ventricular Microanatomy, Arrhythmias, and the Electrochemical Driving Force for Na+: Is There a Need for Flipped Learning?
- Revealing the Concealed Nature of Long-QT Type 3 Syndrome
- Cardiac Response to Low-Energy Field Pacing Challenges the Standard Theory of Defibrillation
- Electrical Coupling and Propagation in Engineered Ventricular Myocardium With Heterogeneous Expression of Connexin43
- Modeling electrical activity of myocardial cells incorporating the effects of ephaptic coupling
- Adaptive multiscale model for simulating cardiac conduction