Why strengthening gap junctions may hinder action potential propagation

1 Introduction

Gap junctions are found throughout the body [1–3]. In particular, gap junctions connect astrocytes [4], neurons [5–8], especially during development [9, 10] and heart myocytes [11–13]. More recently, gap junctions were observed in cancer cells [14, 15] and between soma and germline cells [16].

Gap junctions allow ions to pass between cells, and so the voltage in one cell directly influences the voltage in the neighboring cell. Suppose two excitable cells are connected, a “trigger” cell and a downstream neighbor, and both are at rest. If there is a small voltage deflection in the trigger cell, then we will see a proportional deflection in the neighbor. This proportion is often referred to as the coupling coefficient [17]. If the gap junction conductance is low, then an action potential (AP) in the trigger cell may only yield a subthreshold response in the neighbor. This response is called passive propagation since the neighbor's firing currents are never fully activated. On the other hand, if the gap junction conductance is high, then an AP in the trigger cell may yield an AP in the neighbor. We call this active propagation since the neighbor's firing currents amplify the response.

Active propagation of APs across gap junctions is seen in cardiac tissue [11–13] and throughout the nervous system [18–25]. More recently, gap junctions were found to mediate propagating calcium waves through smooth muscle tissues [26], neural progenitor cells [27], possibly astrocytes [28–30], and retinal cells [31].

Several studies have reported that APs propagate most easily across a given network for intermediate gap junction conductances; see Figure 1. Experiments in heart tissue show that reducing the gap junction conductance can allow APs to propagate through expanding tissue [12, 32, 33]. These results are corroborated in computational models of heart tissue, which show that there is an ideal gap junction conductance that maximizes the number of downstream neighbors that APs may propagate across [34] and that maximizes the “safety factor,” a measure of ease of propagation in terms of source-sink mismatch [34–38].

To understand these results, we consider a network of cells connected by gap junctions with possibly many downstream neighbors. If the gap junction conductance between cells is zero, clearly APs cannot propagate through the network. As the gap junction conductance is raised, propagation becomes possible. However, if the conductance is raised too much, current shunted to downstream neighbors may prevent APs from propagating.

Extending the study in [40] and [41], we use a one-dimensional simplification of the FitzHugh-Nagumo equations to analyze when active AP propagation is possible in branching networks. In particular, our framework explains why increasing gap junction conductance may hinder AP propagation beyond pointing to a source-sink mismatch: AP propagation is the easiest when the total gap junction conductance is large enough to overcome the cell's intrinsic conductance below the firing threshold, but small enough to not interfere with the cell's conductance above the firing threshold. We illustrate this phenomenon for AP propagation into a single cell and branching networks. We then verify that our conclusions hold qualitatively for more realistic biophysical models. Our analysis also predicts that if the gap junction conductance reaches a threshold, we may see semi-active propagation where individual cells in the network cannot fire an AP independently from their neighbors, but the network can still propagate APs.

2 Propagation into a single cell

2.1 Model of a single cell

We model a single cell using the equation

where v_T is a parameter with 0 < v_T < 1. We let F(v) stand for the right-hand side of Eq. (1), illustrated in Figure 2. We think of v as a non-dimensionalized membrane potential and will therefore refer to it as voltage. We also think of F(v) as the inward activation current (more accurately, current divided by capacitance) and therefore may interpret F′(v) as the intrinsic conductance of the cell. Equation (1) is one half of the FitzHugh-Nagumo model [42]. We omit the other half, the slow variable that drives v down when it is high. Our focus, for now, is on the initiation, not on the termination of spikes. Single cell simulations in the Supplementary material show that slow gating variables do not appreciably change during AP initiation.

Equation (1) has a stable equilibrium at v = v_R = 0, an unstable one at v = v_T representing the intrinsic threshold, and a stable one at v = v_F = 1 representing the AP firing voltage. If v is perturbed slightly from 0, it returns to 0, but if it is raised above v_T, it rises further to 1 instead. While we technically have a bi-stable differential equation, since we are using this model to focus on AP initiation, we still think of this as a model for excitability. In this study, excitability always means the co-existence of a stable resting state and an unstable threshold, with the property that perturbation away from the stable resting state and past the threshold results in a large excursion to the AP firing voltage. Once the voltage reaches the AP firing voltage, we assume the recovery currents of the full FitzHugh-Nagumo model would take the voltage back to rest.

The graph of F has a local minimum v_min and an inflection point at v_i = (v_T + 1)/3. We assume v_T < v_i, which is equivalent to

For neurons and myocytes, this is a natural assumption: The firing threshold is closer to the resting voltage (v = 0) than to the peak voltage (v = 1). So

as illustrated in Figure 2.

2.2 Gap junction connections with neighbors held at fixed voltage

Consider a cell described in the previous section that is connected by gap junctions to one upstream neighbor with voltage v_u and k downstream neighbors with voltage v_d = 0. The upstream and downstream voltages are assumed fixed for now. Only the central cell has dynamics. Simulations (not shown) verify that the assumption v_d = 0 is reasonable as long as the central cell's voltage is below threshold. We assume that all gap junction conductances have the same strength g, and stimulation comes from the upstream neighbor. See Figure 3 for an illustration with k = 3.

Our model now becomes

As is customary, we assume here that gap junctions are governed by Ohm's law—current is proportional to voltage difference. The term g(v_u − v) models the effect of the upstream cell and −gkv models the effect of the downstream cells. With the notation

Eq. (2) becomes

Figure 4A shows F and G in one figure. The graph of G intersects the horizontal axis at (v_u/(k + 1), 0), and the vertical axis at (0, −gv_u).

In general, there may be more upstream cells than just one. There may also be more or fewer downstream cells coupled to the central cell with varying conductances. What matters is that g is the total gap junction conductance of connections between the central cell and upstream cells and kg is the total gap junction conductance of connections between the central cell and downstream cells. Therefore, we will not assume k to be an integer any more, simply thinking of it as the ratio of total downstream conductance over total upstream conductance.

2.3 All neighbors at rest

We first discuss the equilibria of the central cell when the upstream neighbors of the central cell are at rest (v_u = 0) just like the downstream neighbors. The formula for G(v) is now

Fixed points of Eq. (2) are solutions of

If g(k + 1) is not too large, then there are three equilibria, depicted in the left panel of Figure 4B. As g(k + 1) increases, the graph of G becomes steeper. This raises the threshold v_t and lowers the AP firing voltage v_f but does not affect the resting voltage. Eventually, the threshold and firing equilibrium points collide at a voltage which we denote by v_E; see the middle panel of Figure 4B. For larger g(k + 1), there is only one equilibrium at v = 0.

We say that the cell is excitable if there are three equilibria—an unstable threshold surrounded by two stable equilibria. This is consistent with the common notion of an excitable cell, where a small deflection in voltage across a threshold can cause a large excursion to a much higher voltage. Hence, the cell is excitable if

and is not for higher values of g(k + 1). So gap junction connections with neighbors held at rest, if they have high conductance or are numerous enough, can make excitability disappear altogether. We will see below that non-excitable cells can still amplify a response—just not independently from other cells.

2.4 Firing triggered by raising the upstream voltage

As v_u is raised from 0 to some maximum value V_u > 0, the graph of G moves downward. Depending on the values of V_u, k, and g, the threshold and resting voltages may eventually collide and annihilate each other in a saddle-node bifurcation, leaving the AP firing voltage as the only equilibrium. We let v_u,c be the bifurcation point, where v_min < v_u,c < v_i; see Figure 5A.

As the right panel of Figure 5A shows, even when the central cell is not excitable because it is too strongly connected to resting neighbors, raising v_u may restore excitability by making the threshold and AP firing voltage equilibria re-appear via a saddle-node bifurcation. Then, at v_u = v_u,c, the resting and threshold equilibria collide and disappear.

Definition 1. We say that the central cell fires or has a firing response when v_u rises from 0 to V_u if there is a saddle-node bifurcation between the rest and threshold fixed points at some critical value v_u = v_u,c with v_min < v_u,c < v_i.

Definition 2. The critical segment is the segment between v = v_min and v = v_i on the graph of F.

The left panel of Figure 5B shows the critical segment in bold. The collision between the rest and threshold fixed points must occur on the critical segment. Therefore, the following proposition holds.

Proposition 1. Let g > 0 and k ≥ 0. The central cell fires when v_u is raised from 0 to V_u if and only if L_g,k (the graph of G with v_u = V_u) satisfies the following two conditions.

(1) L_g,k lies strictly below the critical segment

(2) L_g,k has slope $<F^{\prime }(v_i)$.

We say L_g,k is admissible if it satisfies the conditions in Proposition 1, as illustrated in the left panel of Figure 5B.

2.5 Region in (g,k)-plane with firing response

We use Proposition 1 to derive a detailed description of the set of the parameter pairs (g, k) for which the central cell has a firing response. First, we note that L_g,k cannot be admissible if V_u ≤ v_T. We will therefore assume

Furthermore, since lowering k rotates L_g,k clockwise around the point (0, −gV_u), lowering k makes it easier to fire: For L_g,k to be admissible for a given g, L_g,0 has to be admissible (this makes sense, since lowering k lowers the current shunted to the downstream neighbors). Therefore, we first analyze the case k = 0. Note that L_g,0 passes through (V_u, 0). For L_g,0 to be admissible, g must satisfy

where L_{gmin, 0} is tangent to the critical segment and L_{gmax, 0} has slope $F^{\prime }(v_i)$, as illustrated by the right panel of Figure 5B. Since the slope of L_g,0 is g, then $g_{\text{max}}=F^{\prime }(v_i)$.

Now, we consider k > 0. Suppose g lies between g_min and g_max. Then, L_g,0 is admissible. Other admissible lines, for the same value of g, are obtained by rotating this line around the point (0, −gV_u) in the counter-clockwise direction until it touches the critical segment or reaches slope $F^{\prime }(v_i)$—whichever comes first.

Which of these two does come first depends on g. If g is smaller than some threshold value g_* (discussed further below), then as the line rotates counter-clockwise around (0, −gV_u), it touches the critical segment before it reaches slope $F^{\prime }(v_i)$; see the left panel of Figure 5C. If g is larger than g_*, the line reaches slope $F^{\prime }(v_i)$ before it touches the critical segment; see the right panel of the figure.

From Figure 5C, we see that g_* is determined by the tangent to F at (v_i, F(v_i)), which intersects the vertical axis at −g_*V_u. That is,

Is it possible that $g_{*}\ge g_{\text{max}}=F^{\prime }(v_i)$? According to (6), this means

or

The right-hand side of this inequality is the v-intercept of the tangent to F at (v_i, F(v_i))). So, g_* ≥ g_max only if V_u happens to be very small, just above v_T.

For every g ∈ (g_min, g_max), the slope g(k + 1) of the admissible line L_g,k is bounded above and the upper bound on g(k + 1) translates into an upper bound on k. So, the parameter regime where there is a firing response is of the form

and our next goal will be to understand what the graph of the function k_max looks like.

When g_* ≤ g < g_max, the constraint on the slope of L_g,k is $g(k+1)<F^{\prime }(v_i)$, and therefore

For g_min < g < g_*, k_max(g) is obtained by computing the k where L_g,k is tangent to the critical segment. The two pieces match up continuously at g_*, and we set

Our conclusions are summarized in Figure 6A. When V_u is so small that g_* ≥ g_max, the dashed part of the blue boundary in Figure 6A is simply absent. The dotted and dashed parts of the upper boundary of the blue region in Figure 6A fit together continuously but not differentiably.

In Figure 6A, we also indicate the regions in the (g, k)-plane, where the central cell is excitable when the upstream cell is at rest (v_u = 0). From Eq. 5, the boundary of this region is defined by

We let k_exc(g) stand for the right-hand side of (9) and describe how excitability affects behavior in more detail in section 2.7.

2.6 Ideal gap junction conductance

In Figure 6A, we see that k_max(g) reaches a maximum point at (g_peak, k_peak) at some g ∈ (g_min, g_*). We also see that for all k < k_peak, as g increases AP propagation will first appear and then disappear as illustrated in Figure 1. To understand precisely how (g_peak, k_peak) corresponds to our model, we will examine what happens to L_g,k when it is tangent to the critical segment. In particular, we will follow what happens to the v-intercept (V_u/(k_max(g) + 1)) as g increases. As shown in Figure 6B, initially, the v-intercept falls as g increases until it reaches v_T. As g continues to increase, the v-intercept rises again. Therefore, k_max rises until the v-intercept reaches v_T and then falls. We can solve for (g_peak, k_peak) by setting V_u/(k_peak + 1) = v_T, and therefore

and since $g_{\text{peak}}(k_{\text{peak}}+1)=F^{\prime }(v_T)$,

In other words, AP propagation is the easiest when the total conductance g(k + 1) is the same as the cell's intrinsic conductance when v = v_T. At lower conductances, the gap junction current may not overcome the cell's current below threshold. At higher conductances, the gap junction current may interfere with the cell's firing currents.

2.7 Active, semi-active, and passive propagation

In Figure 6A, we see several distinct regions defined by k_max(g) and k_exc(g). Notably, k_exc(g) intersects k_max(g) so that a cell may fire without being excitable. Here, we describe these parameter regions.

Definition 3. The propagation region is the region of parameter pairs (g, k) with 0 ≤ k < k_max(g). The active propagation region is the section of the propagation regions where k < k_exc(g), while the semi-active propagation region is the section where k ≥ k_exc(g). The passive propagation region is the region where k ≥ k_max(g).

When (g, k) lies in the active propagation region (illustrated in Figure 6A where the red and blue regions overlap), the central cell is excitable for v_u = 0 and fires when v_u is set to V_u. When (g, k) lies in the semi-active propagation region, the central cell is not excitable for v_u = 0, but as v_u rises from 0 to V_u, first excitability is restored, and then there is a firing response. In the passive propagation region, raising v_u from 0 to V_u still causes the voltage of the central cell to rise but does not trigger firing.

To illustrate the difference between the three parameter regions, Figure 7A shows solutions of (2) with v(0) = 0 and

In the active propagation region, the central cell's voltage rises and converges to the AP firing voltage, indicating an AP would occur in the full FitzHugh-Nagumo model. (Recall that our model includes no explicit mechanism for bringing the voltage back down after firing.) In the semi-active propagation region, the central cell's voltage also rises but drops back down to rest once v_u is set back to 0—no matter how high the central cell's voltage gets. In the passive propagation region, the voltage rises much less and returns to zero soon after v_u is switched back to zero.

2.8 AP height can depend discontinuously on g and k

For g > 0 and k ≥ 0, let v = v(t) be the solution of (2) with v_u = V_u and v(0) = 0. Let

In other words, v_∞ is the voltage obtained by setting the voltage in the upstream cell to V_u, letting the voltage in the central cell equilibrate, starting at v = 0; so, v_∞ is the smallest equilibrium. We view v_∞ as a function of g and k, as illustrated in Figure 7B. Notably, there is a discontinuous jump along the boundary defined by k_max (the dotted blue curve). This boundary is discontinuous in general, as stated in the following theorem.

Theorem 1. v_∞(g, k) has jump discontinuities at parameter pairs (g, k) with g_min ≤ g < g_* and k = k_max(g) and is continuous everywhere else.

Proof. Let g > 0 and k ≥ 0. As before, denote by L_g,k the graph of G with v_u = V_u. Then, v_∞(g, k) is the v-coordinate of the left-most intersection point of L_g,k and the graph of F. As illustrated in Figure 5A, v_∞ depends discontinuously on (g, k) if and only if L_g,k is tangent to the graph of F along the critical segment. The (g, k) for which L_g,k is tangent to the graph of F at the critical segment are precisely the ones that define k_max(g) where g < g_*.

3 Propagation through branching networks

3.1 Model

For motivation, first think about AP propagation through a tree of cells connected by gap junctions, as shown in Figure 8A. Assume that the cells are identical and each cell has exactly one upstream and k downstream neighbors. Propagation is triggered by raising the voltage of the leftmost cell to v = 1.

By symmetry, the voltages in cells that are vertically aligned in Figure 8A are identical. We therefore collapse vertically aligned cells into a single cell and obtain the equivalent chain shown at the bottom of Figure 8A. In this chain, gap junction conductances are not symmetric: g for upstream connections and kg for downstream connections.

As we did earlier, we now drop the assumption that k is an integer, viewing it more abstractly as the ratio of downstream conductance over upstream conductance. We even allow k to be less than 1. This models a situation in which a signal originates from a large array of cells, propagating into a narrowing network: Think of a signal propagating through a tree from the tips of the outermost branches to the root, as shown in Figure 8B. We therefore consider the collapsed chain in Figure 8A as a mathematical idealization of many different kinds of branching networks.

We make the additional simplification that each cell responds only to the voltage upstream of it and behaves as though its downstream neighbors were held at zero voltage; its voltage therefore rises to v_∞, as defined in Eq. 12. Propagation through the chain amounts to iterating the map

which maps [0, 1] into [0, 1].

3.2 AP propagation determined by sequence of max voltages

A complication in our analysis is that the map φ need not be continuous. To see this, first note that for any equilibrium v_∞ of eq. (2), that is, for any solution v_∞ of

we can solve for v_u in terms of v_∞:

We denote the right-hand side of this equation by ψ(v_∞); Figure 9A shows an example. The graph of φ is obtained by swapping the v_∞- and v_u-axes, and using the minimum value of v_∞ for a given v_u when there is ambiguity; see Figure 9B.

Because of the discontinuity, we will not use the standard theory of iterated one-dimensional maps to understand the effect of iterating φ, but explicitly construct the iterates. Starting with v₀ = 1, we define v₁, v₂, … by v_j+1 = φ(v_j) and j = 0, 1, 2, …. Given v_j, v_j+1 is the smallest solution of F(v) + g(v_j − v) − gkv = 0 or

The graph of the right-hand side of this equation is the line L_g,k with V_u = v_j, which passes through the point (v_j, gkv_j). To find v_j+1 from v_j, we first draw the line P(v) = gkv as a reference and then draw L_g,k through the point (v_j, gkv_j). The v-coordinate of the smallest intersection of L_g,k with F is v_j+1. The construction of v_j+1 from v_j is illustrated by Figures 10A–C.

From Figure 10, we read off the following theorem on the convergence of the v_j. In this theorem, v_E is defined as in Section 2.3 (Figure 4).

Theorem 2. Suppose $gk\le F^{\prime }(v_E)$ and v₊ is the maximum solution to F(v) = gkv. Then $\underset{j\to \infty }{lim}{v}_j=v_{+}$ if the line L_g,k through (v₊, F(v₊)) is admissible or is tangent to F(v) along the critical segment. In all other cases, $\underset{j\to \infty }{lim}{v}_j=0$.

If $\underset{j\to \infty }{lim}{v}_j=v_{+}$, we say that there is persistent propagation. If the line L_g,k through (v₊, F(v₊)) is tangent to F(v), we say it satisfies the tangency condition.

3.3 The region of persistent propagation

Our aim now is to understand for which values of g and k—that is, for which total upstream gap junction conductances and degree of expansion in our branching network—there can be persistent propagation. The result will be that the region of persistent propagation in the (g, k)-plane is of the form g ≥ g_min and 0 ≤ k ≤ k_prop(g), as shown in Figure 11. In particular, there is an “ideal” conductance which maximizes the k that can support persistent propagation.

The geometric construction described by Figures 10A–C and theorem 2 mean that there is persistent propagation precisely if the line P(v) = gkv (depicted in black in Figure 10E) intersects or at least touches F(v) at some positive value v₊ ≥ v_E and with v_u = v₊ (depicted in red) lies below the critical segment or at most touches it in a single tangency point.

To better understand the shape of the region of persistent propagation (Figure 11), we fix g and ask how Figure 10E changes as k is increased (illustrated in Figure 12). The slopes of both P(v) (black) and L_g,k (red) increase with k. Because the slope of P(v) increases, v₊ decreases. Therefore, the vertical intercept of L_g,k (−gv₊) also increases. This implies that L_g,k rises and is closer to the critical segment if it was admissible before.

We conclude that for a fixed g, the range of values of k that allow persistent propagation is either empty or an interval of the form [0, k_prop(g)] with k_prop(g) ≥ 0. In particular, if there is persistent propagation for any k at all, there is persistent propagation for k = 0. Since the vertical intercept of L_g,0 (red) is −g, the minimum possible g is g_min, where L_{gmin, 0} is tangent to the critical segment (Figure 10F).

For a given g > g_min, we obtain the range of values of k for which there is propagation by watching how P(v) = gkv and L_g,k (black and red lines) rotate as k increases, starting at k = 0. Figure 12 shows two examples. This leads to the following result.

Theorem 3. The region of persistent propagation in the (g, k)-plane is defined by

where L_{gmin, 0} is tangent to the critical segment and k_prop is continuous with k_prop(g_min) = 0. For small g, k_prop is defined by the tangency condition (described in Theorem 2 and illustrated in Figure 10E). For large g, $k_{\text{prop}}(g)=F^{\prime }(v_E)/g$.

Since k_prop is 0 at g_min and tends to zero as g → ∞, there is always a value of g that is “optimal” for propagation in the sense that k_prop is maximal. Figure 11 shows the region of propagation for v_T = 0.2. Numerical simulations suggest that the peak still occurs approximately where $g(k+1)=F^{\prime }(v_T)$. This reflects our results for a single cell, that AP propagation is the easiest when the total gap junction conductance and cell's intrinsic conductance at the threshold balance.

4 Comparison with more realistic models

We apply our analysis to see if we can predict AP propagation in simulations using more realistic models of a neuron and a cardiac myocyte. We use a four-dimensional Hodgkin-Huxley-type model of the neuronal axon [39], which serves as a base model for hippocampal axons connected by gap junctions in several studies [43, 44]. For the cardiac myocyte, we use the eight-dimensional Luo-Rudy model [45], which captures many of the complexities in the cardiac action potential [46]. We used Anaconda Jupyter notebooks [47] with Python 3 to run simulations and apply our analysis. Numerical simulations were conducted using solve_ivp in the Scipy Integrate library (version 1.10.1, [48]) using implicit differentiation specified by the ‘BDH' option. The code for all simulations and analysis are available in the Supplementary material.

4.1 A neuronal model

To start, we compare the single cell analysis of Section 2 with simulation results for a model of the neuronal axon. We consider a single cell with downstream neighbors held at rest (v_R = 0 in [39]). To apply our analysis, we construct the activation curve for this model, F_HH, by taking the right-hand-side of the voltage differential equation and replacing m with m_∞(v), and keeping n and h (which are slower than v near v_R) fixed at resting values. Using F_HH, we calculate v_F, the AP firing voltage of the disconnected cell, and v_E, the maximum possible threshold voltage (see Figure 4 and the Supplementary material).

Using the full four-dimensional model, we run simulations of a single cell connected to fixed up- and downstream neighbors (Figure 3) for many (g, k) pairs under two different conditions:

1. The upstream cell is held at v = v_F for the entire simulation.

2. The upstream cell is held at v = v_F until the cell's voltage reaches v_E.

For each simulation, we calculate the maximum voltage of the cell over the entire simulation for conditions 1 and 2. If the cell's voltage rises beyond v_E under both conditions, then we say that active propagation is possible. However, if the cell's maximum voltage goes beyond v_E only under the first condition, then only semi-active propagation is possible since the cell cannot maintain a full AP on its own. These results are shown in Figure 13A.

We then calculate k_max(g) and k_exc(g) using F_HH, and compare them directly to the simulation results. Figure 13A shows that k_max(g) for a single cell matches the simulated active propagation region boundary quite well for small g, with a discontinuous jump in the voltage height up to g_*. While the simulated boundary for the semi-active propagation is sharp as predicted in our framework, the simulated boundary lies to the left of the predicted boundary. The boundary matches the simulated active propagation regions, especially for g < g^*, which confirms that fast activation currents are primarily responsible for AP initiation in a single cell.

For the branching network, we make modifications in both the analysis and the simulations. In the simulations, we model 20 layers of cells and record the maximum voltage of the penultimate cell to indicate whether APs propagated through the network. We also only hold downstream neighbors for the last cell fixed at v = 0. In our analysis, we therefore adjust the assumption that downstream neighbors are held at v = 0. Instead, we make a cell's downstream neighbor voltages proportional to the cell's voltage [49], which follows the steady-state condition when there are only passive currents. The details of this modification are explained in the Appendix in the Supplementary material.

In Figure 13B, we see that the simulated active propagation region for the branching network lies entirely within the predicted region and has a similar shape. That k_prop over-estimates the boundary for AP propagation is to be expected since our analysis does not take the duration of an AP into account. (Though higher-order factors may shape the region as well, see section 5.) Regardless, we may still use the curve $g(k+1)=F^{\prime }(v_T)$ to estimate the peak, where the total gap junction conductance matches the cell's conductance at the intrinsic threshold. In Figure 13B, we see that the peak of the simulated active propagation region lies to the right of the peak curve. As predicted, there is a sharp boundary between the simulated semi-active propagation and active propagation regions. However, k_exc(g) lies to the right of the simulated boundary. Figure 14 illustrates the behavior seen in each region of Figure 13B for k = 2, and Figure 15 gives an example of AP propagation on both sides of the simulated semi-active propagation boundary. In particular, in Figure 15, we see APs propagate in the branching network only if cells in multiple layers are depolarized enough to fire.

4.2 A myocyte model

We compare simulations using the full eight-dimensional Luo-Rudy myocyte model [45] to our analysis in a similar manner to the neuron model above. We construct the activation curve F_LR by replacing m with m_∞(v) and holding all other gating variables (which are slower than v near v_R) fixed at their resting values. For ease of calculation, we also re-center F_LR(v) so that the resting fixed point is v_R = 0. Using the same simulation setup as in section 4.1, in Figures 13C, D, we again see that simulated AP propagation into a single cell matches the predicted boundaries very well, especially for small g. Similarly, for the branching network, the simulated active propagation region lies entirely within the predicted boundary but does not quite have a similar shape. This may be due to the fact that slower calcium currents in the Luo-Rudy myocyte model allow APs to propagate for low g [38]. Similar to the neural model, the peak of the simulated active propagation region lies to the right of the peak curve. There is a sharp boundary between the semi-active and active propagation regions, which lies to the left of the predicted boundary.

5 Discussion

We explained why AP propagation through a branching network of excitable cells connected by gap junctions depends non-monotonically on gap junction conductance. We did this by reducing models to one-dimensional firing currents, which allows us to visualize how the firing currents and gap junction currents interact. From this visualization, we estimate that AP propagation is easiest when the total gap junction conductance matches the cell's conductance at its intrinsic threshold. We also found a discrete behavioral region where cells may propagate APs but are so strongly connected that they are no longer individually excitable.

We believe this framework can give a simple and efficient method for understanding where AP propagation may occur and may add guidelines on whether to raise or lower gap junction conductance to enhance propagation. For example, cells in realistic networks have a heterogenous number of connections. Our framework may predict when AP propagation failure and subsequent re-entry occurs at cells with higher connectivity [12, 43, 50]. Currently, simulations are used to determine whether APs will propagate through any given network. For instance, whole heart simulations are a promising diagnostic tool in heart patients but can be computationally costly [51]. While cellular automata may alleviate the computational cost [52], our framework can inform phenomenological rules and greatly reduce the parameter space that needs to be searched. Finally, our framework outlines how the AP height, threshold, and subthreshold response changes with gap junction conductance [53]. These measures may be useful in determining the state of a network and predicting network-wide phenomena such as synchronization.

Moreover, this framework may be used to understand how network behavior may change with changes in connectivity, gap junction conductance, and firing currents. Cardiac arrhythmia can be brought about by pathogenic stressors [54, 55], adrenergic stimulation [56], or myocardial ischemia (buildup of plaques) associated with decreased gap junction coupling [57, 58]. Abnormal gap junction ubiquitination can affect both the heart and nervous system [59]. Gap junctions may be sensitive to toxins [60]. Epilepsy has been tied to gap junctions connecting networks of pyramidal cell axons, interneurons, and astrocytes, each of which play a different role in network excitability [61–64]. Furthermore, gap junctions are indicated in Parkinson's disease [65], Alzheimer's disease [66], nerve injury [67–69], and degenerative diseases in the retina [70]. Modulating gap junction conductance may help with several of these conditions; however, questions remain on how much gap junctions can be modulated without producing adverse side effects [71, 72]. An alternative may be to modulate firing currents through gene therapy [73]. There are also many conditions where external stimulation is applied as a part of therapeutic treatment [51, 74–77]. We hope that our framework can give insight on appropriate ranges to maintain proper function in a variety of networks.

Our study is closely related to several other studies. We directly extend the theory presented by Keener [40], which addresses the conditions where AP propagation may occur in a chain of cells with no branching. Our study also extends the study of Kouvaris et al. [41], who analyzed a bistable system where both the high and low resting voltages are fixed. This matches our analysis for AP propagation into a single cell, so their bifurcation diagram matches our bifurcation diagram for a single cell (Figure 6). We take a more geometric approach, which clearly shows that the analysis depends only on the qualitative shape of the activation curve. This allows us to analyze the peak value of k for AP propagation. We also show how connectivity affects the AP firing voltage, which then determines whether persistent propagation is possible. Wang and Rudy first showed a biphasic relationship between gap junction conductance and an expansion ratio using simulations of cardiac tissue [34]. Several experimental and computational studies show a biphasic relationship between gap junction conductance and the safety factor [34, 38, 78, 79]. While our framework can predict where AP propagation is possible, the safety factor measures the source-sink mismatch from a recording or simulation based on the incoming current vs. outgoing current in a single cell. Finally, there are other studies that show an ideal gap junction conductance for propagation using a mean-field model [80], in a network using small-world topology [81], and an ideal cable diameter for AP propagation [82, 83].

We note that our framework only provides simple guidelines outlining where AP propagation may occur. While this theory can give a first estimate on when to expect AP propagation, there are many ways we could increase the accuracy of this prediction. For instance, we could link different versions of φ(v) to model AP propagation through a network with heterogenous numbers of downstream neighbors. We may also improve the prediction by taking the connectivity of downstream neighbors more fully into account. In fact, AP propagation may not only depend on the number of immediate downstream neighbors, but on the number of neighbors several steps away [43]. Future studies may incorporate multiple layers of downstream neighbors. The predicted AP propagation regions in our framework overestimate the regions found in the original models. Preliminary results show that for the Fitzhugh-Nagumo model, the simulated boundary may be better predicted by taking into account the time to reach the threshold. However, a correction based on the timing still does not quantitatively predict the region in higher-order models. This may be due to higher-dimensional interactions within the cell, such as calcium currents in the Luo-Rudy myocyte model which may allow APs to propagate for low g [38].

Our framework studies propagation assuming the first cell is firing. It does not address getting the first cell to fire in the first place but rather reproduces the situation where the cell is voltage clamped. As gap junction conductance increases, it may take an increasing amount of current to bring the first cell to the desired voltage. In some cases, we may need to depolarize multiple downstream neighbors in the network to get AP propagation, similar to what is seen in semi-active propagation.

Overall, reduction of cell models to one-dimensional firing currents gives insight into how gap junction currents allow active propagation, semi-active propagation, or passive propagation. The shape of the firing currents explains how the threshold and peak firing voltage change with gap junction conductance and allows us to approximate an ideal gap junction conductance for AP propagation. This framework gives further insight into the mechanisms of these networks - which play a key role in both normal biophysiological function and disease.

Data availability statement

The original contributions presented in the study are included in the article and Supplementary material. Further inquiries can be directed to the corresponding author.

Author contributions

EM developed the original concept and analysis. EM and CB revised the analysis and visualization and wrote and edited the manuscript.

Funding

We thank Tufts University for paying for the open access page charges through a Faculty Research Awards Committee (FRAC) grant. EM was partially funded by a National Science Foundation (NSF) Mathematical Sciences Postdoctoral Research Fellowship, NSF Award Number 0802889.

Acknowledgments

The authors thank Nancy Kopell for many insightful comments on this study, and students who participated in this research, especially Abel Gonzalez and Mitra Kermani.

Conflict of interest

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

Publisher's note

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

Supplementary material

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

Presentation 1. Appendix.

Data Sheet 1. Code for all simulations and analysis.

References

Glossary

active response: a cell is excitable and fires as v_u increases from 0 to V_u

admissible: L_g,k lies below the critical segment

AP: action potential

critical segment: segment on the graph of F between v_min and v_i

excitable cell: a cell with a stable resting voltage and unstable threshold when v_u = 0

fire: in a single cell, a saddle-node bifurcation occurs on the critical segment as v_u increases from 0 to V_u

F(v): inward activation current

g: total upstream gap junction conductance

g_min: minimum g where firing is possible

g_max: maximum g where firing is possible in a single cell

g_peak: g where k_max(g) attains its maximum

g_*: value of g where the condition determining k_max(g) changes

G(v): outward gap junction current (see Eq. 3)

k: ratio of downstream gap junction conductance over upstream conductance

k_exc(g): maximum k where a single cell is excitable

k_max(g): maximum k where a single cell can fire

k_peak: maximum value of k_max(g)

k_prop(g): maximum k where $\underset{j\to \infty }{lim}{v}_j=v_{+}$ as a function of g

k_*: k_max(g_*)

L_g,k: graph of G when v_u = V_u

passive response: single cell doesn't fire as v_u increases from 0 to V_u

semi-active response: single cell isn't excitable but fires as v_u increases from 0 to V_u

v_d: voltage of a downstream cell

v_E: value v where a saddle-node bifurcation occurs as g(k + 1) increases and v_u is fixed at 0, rendering the cell unexcitable

v_F: intrinsic AP height of single cell without gap junction connections

v_f: AP height of a single cell with gap junction connections

v_i: inflection point of F(v)

v_j: voltage of cells in the jth layer of a branching network

v_min: v where F(v) attains its local minimum

v_R: the resting equilibrium point of a single cell

v_T: intrinsic threshold of single cell without gap junction connections

v_t: threshold of a single cell with gap junction connections

v_u: voltage of an upstream cell

V_u: maximum voltage or AP height of an upstream cell

v_u,c: bifurcation value of v_u as a single cell fires

v_∞: lowest fixed point in a single cell when v_u = V_u

v₊: maximum intersection of F and P(v) = gkv if $gk\le F^{\prime }(v_E)$

φ(v_u): v_∞ as a function of v_u

ψ(v): v_u that will produce a fixed point at v