Conductance heterogeneities induced by multistability in the dynamics of coupled cardiac gap junctions

The description of the electrical activity in the cardiac tissue is complex, and its study represents a very active field of research in applied mathematics and bioengineering.^1–3 The models used to study cardiac dynamics have become more detailed over time, obtaining a better comparison with experimental results, at the expense of a higher computational cost. Besides direct numerical simulations of the models, insight has been obtained through two different, but complementary, venues: one is the use of simplified models that allow faster simulations and a deeper mathematical analysis^4,5 and the other is the use of standard techniques in applied mathematics, such as stability analysis or bifurcation theory, to study the properties of the different cardiac models.

Stability analysis has been successfully applied to solve many problems in disparate areas of biology, such as gene regulatory networks in cells,⁶ coupled neurons,⁷ the study of alternans,⁸ or early afterdepolarizations (EADs)⁹ in a cardiac tissue. In particular, the phenomenon of multistability has been associated with the production of complex patterns, both spatially and temporally.¹⁰ 

In this paper, we will use the aforementioned techniques to study the propagation of the cardiac waves in a system with voltage and time dependent gap junctions (GJs). GJs often appear as clusters of up to many thousands of connexons¹¹ (hemi-channels formed by six connexins). In mammals, most common types of proteins conforming the connexons in cardiac cells are the Cx43 (in atria) and Cx45 (in ventricles). The GJ dynamics has been experimentally characterized by Vogel and Weingart¹² and Desplantez et al.,^13,14 its dependence on temperature,¹⁵ and its voltage regulation¹⁶ has also been recently described. Altered GJ regulation has been associated with various cardiac dysfunctions. For instance, an experimental study by Beardslee et al.¹⁷ has put forward the relation between ischemic hearts and the remodeling of the cardiac electrical conductance through the dephosphorylation of ventricular Cx43. Jalife et al.¹⁸ have shown in mouse hearts the crucial influence of the GJ dynamics with the propensity to trigger arrhythmias. The influence of the GJ through the slowing down of the action potential has also been put forward in a paper by Gutstein et al.¹⁹ Dynamic uncoupling is also a well-known phenomenon in neurons,²⁰ where it has been shown that GJ gating could lead to unidirectional conduction blocks and wave breaks.²¹ In another context, cancer cells usually have downregulated levels of gap junctions. Several lines of evidence suggest that loss of gap junctional intercellular communication is an important step in carcinogenesis.²² 

In the heart, abrupt changes in conductance produce a sink–source mismatch that may result in the dispersion of velocity, often giving rise to a conduction block and re-entry.²³ The spatial changes in conductance are due to remodeling of the cardiac tissue caused by, e.g., fibrosis. Therefore, it is due to a pre-existing heterogeneity of the medium. In a recent publication,²⁴ we observed that the dispersion in velocity might appear in an utterly homogeneous tissue due to the voltage and time dependence of the intercellular conductance. In the present paper, we will further extend the analysis of the spatially dependent conductances and provide a stability analysis of the mathematical model.

The purpose of this paper is threefold: first, we review briefly the main result in Ref. 24 and show how the inclusion of the GJ dynamics in the model for the action potential (AP) propagation may lead under some circumstances to spatially nonuniform electrical conductances; second, we explain why this multistability does occur and study it in a reduced system (considering successively only one, two, and three coupled GJs); third, we propose a simplified spatially extended model that captures the characteristics of the full model but with a computational speed gain of over four orders of magnitude. The results of the latter model compare satisfactorily with the full model. The organization of the paper is as follows. In Sec. II, we provide the full model description for the AP propagation in a cable with GJ dynamics. The simulations of the full model are presented and analyzed in Sec. III. In Sec. IV, we detail the stability analysis of the GJ dynamics and uncover the reason for the multistability observed in Sec. III. In Sec. V, we propose a simplified model, and we compare the results given by the full and simplified models. We found that there is good agreement between them. Finally, discussions of the limitations of this study and future perspectives are given in Sec. VI. In the  Appendix, we provide an explicit calculation of the space constant, an essential parameter for the AP propagation’s characterization.

Let us present the equations that are used to study the long-term effects of the variation of the tissue conductance induced by the GJ dynamics. It was shown in a previous study²⁵ that it is sufficient to couple the gap junction (GJ) dynamics using a monodomain formulation of the cable equation. The same model was used in a recent publication²⁴ by the authors, where they showed numerical evidence of the multistability phenomenon. For completeness, in this section, we present the model used in Ref. 24. Then, in Secs. IV and V, we develop the theory to explain the observed multistability and give a simplified version of the cable equation model that allows studying the long-term effects of gap junctional dynamics much more efficiently.

The model is based on the monodomain formulation of the cable equation,^26,27

where $s$ is the dynamical state vector that collects all the variables involved in the local description of the ionic currents that cross the myocyte membrane. $V$ denotes the transmembrane potential (standard units are in mV), $Im$ is the sum of the ion currents (units are in $μA/cm2$⁠), and $C$ is the cell membrane capacitance per unit area (⁠$≈1μF/cm2$⁠). The term $Iext$ allows for the introduction of an external current as it happens during an external excitation of the cardiac tissue. Here, because we use a monodomain formulation, the external excitation is applied in the intracellular domain. The electrotonic current between adjacent myocardial cells is modeled through the Laplacian term in Eq. (2). Typical values for the conductance when assuming a constant homogeneous tissue lead to a diffusion parameter of $D≈1.5⋅10−3cm2/ms$⁠.²⁸ 

The inclusion of the GJ dynamics has been studied in a modeling paper by Hand and Griffith.²⁹ The latter study combines a multi-scale approach with microstructure details and macroscopic tissue characteristics. In the same vein, Costa et al.³⁰ developed a semi-continuous model to describe the electrical propagation in the cardiac tissue, including the GJ dynamics. Another modeling perspective consists of viewing the heart as a network of different types of cells that are electronically coupled via gap junctional conductance.³¹ 

Due to the introduction of the GJ dynamics into the monodomain model, the diffusion parameter $D=D(x→,t)$ becomes a function of space and time. The relation between the conductance and the GJ dynamical function $g$ is simply written as follows:

where $D¯$ is the fixed nominal value for the inter-cellular diffusion parameter $D¯=1.5⋅10−3cm2/ms$⁠.

To simplify the problem, we will work on a one spatial dimensional setting. Therefore, when we discretize space in order to solve the model equations (1) and (2), the index number $i$ refers to the cell number $i$ on the cable. The gap junction between cell number $i$ and cell number $i+1$ has also index number $i$⁠, following our convention. Here, the dynamics of the gap junctions is modeled following the works of Lin et al.³² and Desplantez et al.¹³ Thus, the set of differential equations that govern the GJ dynamics is written as

where $Δϕ=ϕ(i+1)−ϕ(i)$ is the difference in the intra-cellular electrical potential between the two adjacent cells of the $gi$⁠. The local steady-state value in Eq. (4) depends on the local instantaneous $Δϕ$ following this equation:

and the time scale in Eq. (4) is given by $τg=Aτexp⁡[−Bτ|Δϕ|]$⁠. The parameter values entering Eqs. (4) and (5) are taken from the works of Lin and Desplantez.^13,32 In this work, we will consider GJs of type (⁠$Cx43_43$⁠), which are nearly symmetrical with respect with the sign of $Δϕ$ (see Fig. 1). The corresponding parameters are gathered in Table I. We also have the following parameter values: $A=z/26.714(mV)−1$⁠, $Aτ=109,900$ (ms), $Bτ=1/11.8(mV)−1$⁠, and $gi,max(Δϕ=0)=1$⁠. The dependence of $gi,ss$ as a function of $Δϕ$ is displayed in Fig. 1.

We emphasize that because we are using a monodomain formulation, the computation of the term $Δϕ=ϕ(i+1)−ϕ(i)$ reduces to $ΔV=V(i+1)−V(i)$ as a result of assuming that the extra-cellular electrical potential is constant throughout the domain. Another assumption comes from the fact that each cell is electrically isopotential, which is also a simplification.²⁶ 

Recently, several effects that affect the dependence of $gi,ss$ on $Δϕ$ have been reported.^15,16 After an increase in temperature, for instance, it has been observed that the value of $A$ in Eq. (5) increases, $V1/2$ decreases, and gating kinetics accelerate (cf. $τg$ decreases).¹⁵ A variation in the dependence of the conductance with the intercellular potential has also been observed in connexin 45 variants, in which two amino-terminal ends have been altered.¹⁶ 

To study the alteration of the GJ conductance characteristics, we have modified the dependence of $gi,ss$ on the intercellular potential $Δϕ$ in Eq. (5). We rescale $Δϕ$ as $Δϕ′=αΔϕ$⁠. This is equivalent to modifying the parameters of Eq. (5) as follows: $A′=αA$⁠, $V1/2′=V1/2/α$⁠. The larger $α$ is, the narrower and steeper is the curve in Fig. 1. For this reason, we have denoted $α$ as the factor of shrinking, $FS$⁠.

Detailed models of gap junction dynamics usually consider several open and closed configurations of the connexins that compose the gap junctions.^21,33,34 However, the physical meaning of the factor of shrinking, FS, can be better understood assuming a simple two-state system of high and low conductance. Then, we can consider an effective free energy difference between these two states that depends on the intercellular potential as follows:

The conductance will follow a classical Boltzmann sigmoidal equation as written in Eq. (5), with

Identifying term by term between Eqs. (5) and (7), we have that $A=−a/(RT)$⁠, $V1/2=ΔG0/a$⁠. Thus, a change in the shrinking factor corresponds to a change in the factor $a$⁠, i.e., a modification of the voltage dependence of the effective free energy of this two-state system.

For the dynamics of the transmembrane potential, we will use the generic five-variable model by the model of Cantalapiedra et al.,^5,35 with the model parameters fitted to describe human ventricular myocytes.³⁶ A one-dimensional cable of cardiac tissues was simulated using Eqs. (1)–(5). Time and space were discretized using $Δx=0.01$ cm and $Δt=0.01$ ms. We use a simple Euler explicit scheme to solve the set of Eqs. (1)–(5). In particular, Eq. (2) is discretized as follows:

where the superscript $(n)$ refers to variables taken at time step $n$ and index $i$ indicates the spatial localization on the cable. Here, we consider that the size of the myocytes is constant and corresponds to one grid spacing $Δx=100μ$m. Therefore, one GJ is always located between two spatial grid points. The number of cardiac myocytes in the one-dimensional strand of cardiac tissues that we studied was set to $N=300$ and sometimes $N=1000$⁠. We have checked that $Δt$ that we have chosen is sufficiently small to ensure that the numerical results are robust, and the numerical error is bounded and less than 0.5%.

The dynamics of the GJs are studied following the stimulation protocol described next. The strand of tissues is periodically excited at one end of the cable by injecting enough current to elicit an action potential. In particular, the first seven cells of the system are stimulated for 1 ms with an excitation current of $Iext=0.52μA/cm2$⁠. Following this excitation, an action potential (AP) is elicited. Subsequently, the AP propagates toward the opposite end of the cardiac strand of tissues. We repeat the stimulation periodically. For all the simulations considered, we have fixed the stimulation period to $T=480$ ms. Thus, we have chosen a rather short pacing period to speed up simulations. Most of the simulations presented in this paper require many excitations to display non-trivial behaviors. In some instances, we have performed simulations that required as much as 10 000 excitations to reach a steady state for the GJ dynamics. We have further checked that the numerical scheme does not influence the presented results using the Crank–Nicholson³⁷ discretization scheme for some simulations.

This section presents the results of numerical simulations performed with the model developed in Sec. II. We first discuss the question of the parameter selection. In a healthy tissue, the parameters associated with the action potential model and the connexins are such that the conductance is nearly constant in space and time, as it has been shown in the simulations presented in Ref. 38.

Alternatively, under several pathologies, it has been observed that the average conductance can be dramatically reduced^39,40 and also that the time scale associated with the GJ gating may be altered,⁴¹ for instance, due to drug intake,⁴² autoimmune and inflammatory cardiac channelopathies,⁴³ or an increase in fibrosis content.⁴⁴ This paper aims to simulate a diseased tissue and look at long-term effects on the conductance distribution. In particular, we show that non-trivial spatial distributions eventually appear upon stimulating the tissue when modifying the model’s parameters. We have performed simulations with initial values of the conductance set to a maximum of $40%$ of their nominal values; i.e., $gi,max(Δϕ=0)=0.4$⁠. We also have reduced the GJ gating time scale by modifying this parameter to the new value $Bτ=1/5(mV)−1$⁠. A further modification introduces a parameter associated with the steady-state characteristics of the GJ conductance (⁠$gi,ss$⁠). We called this last parameter the factor of shrinking $FS$⁠. It is directly related to the plateau’s width characterizing the GJ dynamics, as shown in Fig. 2. Evidence of varying GJ characteristics can be found in Refs. 15, 16, and 45.

To study the effect of the modifications in the GJ dynamics, we have performed simulations in a one-dimensional strand of tissue, stimulating it from one side at a stimulation period of $T=480$ ms during 1000 stimulations. As initial conditions, we consider a uniform value of $g=0.4$⁠. Figure 3(a) displays the spatiotemporal dynamics of the GJ when setting the value of the shrinking factor to $FS=2$⁠. We represent the field $gi(t)$ through a color code [see the colorbar of Fig. 3(a)]. For this value of FS after approximately 300 beats, the conductances show a pattern of alternations that is not regular but rather spatially disordered. Note that all the simulations are performed with a maximum conductance reduced to $40%$ of its nominal value. In Fig. 3(b), we have represented the final values of $gi$ (after 1000 beats) for three values of FS; i.e., FS = 1 (black line), FS = 2 (red line), and FS = 4 (green line). For $FS=1$⁠, the final stage is close to the initial value and a uniform constant $g=0.4$⁠. For FS = 4, we observe again a uniform field but with a value close to $g≈0.1$⁠. For intermediate values, $FS=2$⁠, we observe the spatial alternation of the conductance field. We observe that the spatial correlation length of the disordered pattern [see Fig. 3(a)] is only a few cell’s length.

These results suggest the existence of a transition induced by the increase (or decrease) of the parameter $FS$⁠. This transition occurs between the upper (⁠$g≈0.4$⁠) and the lower (⁠$g≈0.1$⁠) uniform conductances. In between these two limiting cases, we observed a spatially chaotic distribution of the conductance where the two states are mixed. This transition has been studied in a previous publication by the authors.²⁴ In the latter paper, it was shown that the two states competing are, in fact, two limit cycles with a period equals to the forcing period $T=480$ ms. It was also shown that the transition between the two states was sensitive to the initial condition (hysteresis phenomenon) and that the initial noise was also crucial in characterizing this transition. In Sec. IV, we show that the observed phenomenon is a result of the relative stability of the solutions of the GJ dynamics.

Several order parameters can be used to characterize the observed transition in Fig. 3. One possible choice is the spatial average value for the conductance $⟨g⟩$⁠, and another possibility is the proportion of GJs that converge to the upper state, as denoted by $PUP$ in Fig. 4. We characterize the transition as a function of $FS$ and the noise (⁠$Ns$⁠) that is added to the mean initial value of the conductance $gini$⁠. We have set $gi(0)$⁠, with $i$ the spatial index associated with the GJ, as follows:

where $σu$ is a uniformly distributed random variable in the interval $[−1;1]$ and $Ns$ is the parameter characterizing the noise strength in the initial condition. Here, we have used three values for the noise strength $Ns=[10−6,10−4,10−2]$ corresponding to very low, medium, and high initial perturbations. Note that the results for the lowest value of $Ns$ could be obtained equally without adding any initial noise. Indeed, the numerical method and the grid discretization are such that some very small “numerical” noise is always present in the simulations, as explained in more detail in the  Appendix.

In Fig. 4, $gini$ is set to 0.4. For comparison purposes, we have grouped the results of the full model and the simplified model in Fig. 4, even if the simplified model will be explained in detail in Sec. V. Due to the randomness in the initial conditions, we always perform 20 realizations of the initial noise and display the median value over 20 realizations. The striking result from Fig. 4 is that the transition from the upper state to the lower state spans several orders of magnitude of the bifurcation parameter FS when the initial added noise is large enough. We observed a stepped transition, where the intermediate step has an extension of up to four orders of magnitude, as seen in Fig. 4, when $Ns=10−2$⁠. To our knowledge, this is very uncommon in bifurcation theory, and we will provide an explanation of this phenomenon in Sec. IV.

Before turning to the stability problem, let us analyze in more detail the disordered state that is shown in Fig. 3(b). In particular, we have computed the relation between the conductance and the action potential (AP) conduction velocity when we have a non-homogeneous state as it happens for intermediate values of the bifurcation parameter FS.

Figure 5 shows the converged states of two simulations where we observe a dispersed state in the conductance. In Fig. 5(a) where $FS=1.61$⁠, one computes a measured $PUP=0.7275$⁠. We observe alternations of high and low values of the conductance field $g$ with a dominance of UP states. On the contrary, in Fig. 5(b), where we have set FS = 2.19, the lower conductance case is dominant and $PUP=0.18$⁠. Here, we are directly interested in the consequence of this dispersed state onto the AP velocity. Let us recall that according to our notation, the gap junction $gi−1$ lies between cell $i−1$ and cell $i$⁠. The velocity $ci$ corresponds to the cell spacing $dx$ divided by the AP front’s travel time between cell $i−1$ and cell $i$⁠. The correlations between the two fields $gi−1$ and $ci$ are obvious from both Figs. 5(a) and 5(b). We have computed the Pearson correlation coefficient between the conductance and velocity fields, and we get a value of $ρ=0.906$ (a) and $ρ=0.996$ (b). Also, apparent from Figs. 5(a) and 5(b) is that when the local conductance is low, then the local velocity is low and vice versa. However, the relation between the local conductance and the local AP velocity is not strictly local. Indeed, in Fig. 5(a), we observe that the AP velocity corresponding to a cell with high conductance located just before a cell with low conductance has a markedly excess value of speed compared to the others. In the same fashion, in Fig. 5(b), we observe that the AP velocity of a cell with low conductance located before a cell of high conductance has a markedly lower value than its counterparts. These last two observations can be explained physiologically by considering the effect of the electrotonic currents on the AP propagation and the source–sink mismatch. The asymmetry of the electrotonic current for a cell at the edge between low and high conductance explains the phenomenon. It is important to remember that the waves are propagating from left to right in Fig. 5. Figure 5 highlights the non-locality relation between the AP velocity and the tissue conductance. We will return to this non-locality relation in Sec. IV B.

The effect of dispersion in the AP velocity has important physiological consequences. It is well known that spatial dispersion of the AP velocity favors re-entries and arrhythmias. With our system, we proceed to a systematic study of the effect of the dispersion on the AP velocity.

Figure 6 condenses our quantitative findings for the relation between the AP velocity and the conductance. In the homogeneous case, we fix the values of the conductance uniformly in the system. On the horizontal scale of Fig. 6, a value of $g=1$ would correspond to a tissue diffusion of $D¯=1.5⋅10−3cm2$/ms. We measure the steady-state AP velocity between grid points 500 and 900 (corresponding to a traveled front distance of 4 cm). We lower the conductance until the conductance is so low (⁠$g≈0.049$⁠) that the AP can no longer propagate through the system. Note that the system size was set to N = 1000 cells to get accurate measurements. The velocity measured in the homogeneous case follows a perfect fitting curve (green line in Fig. 6) given by

The relationship between nerves and cardiac impulses with tissue conductance has been well-studied.^26,46 In particular, in Eq. (10), one should expect that the scaling exponent between the AP velocity and the conductance would be $1/2$⁠. The departure from the $1/2$ exponent can be explained by the fact that we have focused on the very low-speed spectrum of the relationship between the AP velocity and the conductance. For this range, the cardiac tissue’s discreteness from a physiological and numerical point of view matters most.

Let us now turn to the non-homogeneous case. We have seen in Fig. 5 that the local relation between $c$ and $g$ has a strong dispersion. In Fig. 6, we report the velocity computed in the same manner as in the homogeneous case by measuring the AP front’s travel time between grid points 500 and 900 (⁠$⟨c⟩$⁠) after the steady state is reached. Note that $⟨c⟩$ is the harmonic mean of the local velocities between grid points 500 and 900. On the horizontal axis, we compute the average (arithmetic mean) of the conductance between the same points (⁠$⟨g⟩$⁠). To get an extensive range for the dispersion parameter (⁠$pUP$⁠, indicated in the color scale of Fig. 6), we vary the bifurcation parameter $FS$⁠, accordingly.

The overall diffusion parameter $D¯$ is also varied, and the colored circles and diamonds correspond to $D¯=1.5⋅10−3$ and $D¯=0.75⋅10−3cm2/ms$⁠, respectively. Figure 6 shows that the non-homogeneous situation has a dramatic influence on the AP velocity. We observe an overall decrease in the velocities with respect to the homogeneous case. The large fluctuations of the local velocity imply that the velocity computed in the non-homogenous case cannot be easily related to the homogeneous case as evidenced in Fig. 5. In Sec. IV B, we propose an empirical nonlocal relationship between the AP velocity and the conductance values. This relation will be a necessary ingredient for performing the stability analysis in Sec. IV B.

To summarize, one finds that the tissue with non-homogeneous conductance displays a strong dispersion in the local AP velocity. This is especially true when the speed (and conductances) are low. The relationship between the dispersed state and the corresponding average velocity has been measured in our system. An analytical theory (such as, e.g., a mean-field calculation) able to explain our numerical findings of Fig. 6 is beyond the scope of the present paper.

In this section, we aim at explaining, through a series of reduced descriptions, the origin of the multistability phenomenon that was observed in Sec. III (see Fig. 4). The main idea behind the simplified model is to decouple the dynamics of the action potential from the dynamics of the GJ. We want to demonstrate that the onset of the instability leading to the multistability is to be found in the GJ dynamics. We, therefore, computed a typical AP, shown in Fig. 7 and saved it as a data vector with a sample discretization time of $δt=0.01$ ms. In the following developments, whenever we needed the value of the AP in between two sample points, we obtained it using a straightforward linear interpolation between these two points.

Let us start with the simplest model that consists of a single GJ that connects two myocytes (Fig. 8). In this case, we have a non-autonomous first-order differential equation for the GJ dynamics,

where $Δϕ=ΔV=V|2−V|1=V(t−Δt)−V(t)$ is the difference in potential between the two cells. The time delay $Δt$ is the time it takes for the AP to propagate from Cell 1 to Cell 2. Considering that the distance between the two cells is fixed at $δx=100μ$m, one readily computes that $Δt=δx/c(g)$⁠, where $c(g)$ is the function that relates the AP speed with the conductance. This function has been evaluated once for all from the results of simulations of the full model with different constant values of the GJ as indicated by the blue dots and the green fitting curve in Fig. 6. The fitting curve is given by Eq. (10).

Therefore, in this simple model, we have a single differential equation, Eq. (11), that is forced by a periodic function that is perturbed by the value of the gap junction itself (strong perturbation). This makes the analytical solution of Eq. (11) quite intractable, and we rather proceed to solve the equation numerically.

We study the dynamics of the equation by scanning two of its parameters that are the previously defined shrinking factor FS and the initial value of the GJ. In Fig. 9, we collect the results of the converged final state for three different initial conditions [$g(0)=0.1,0.25,0.4$] for the value of FS in the range between 1 and 3. In addition, by launching many simulations for intermediate initial conditions, we were able to compute numerically the unstable manifold (indicated as a dashed blue line in Fig. 9) that separates the two basins of attraction.

Some comments are in order for the analysis of the results presented in Fig. 9. First, the final states plotted in the figure are not fixed points but rather periodic limit cycles (LCs) with a period equal to the period of the forcing (here, $T=480$ ms). To represent the limit cycles, we use the standard Poincaré map section technique that is very handy to study the stability of periodic cycles. Therefore, we have represented in the diagram of Fig. 9 only one representative point of the corresponding limit cycle. Second, Fig. 9 shows two fold bifurcations at approximate values of $FS1=1.46$ and $FS2=1.72$ that defines a bistable region. In this region, we observe the hysteresis phenomenon that is typical of subcritical bifurcations. Indeed, for intermediate values of the parameter FS, we observe the coexistence of two stable limit cycles. The chosen dynamics essentially depends on the initial condition. This bifurcation diagram is reminiscent to what it is observed with the Duffing oscillators.⁴⁷ As an example, let us illustrate in more detail what happens in the bistable region when FS = 1.67.

If we set FS = 1.67, i.e., in the bistable region, we observe in Fig. 10 that the final state of the system depends crucially on the initial condition. In this case, there is a saddle point around $g≈0.25445$ that separates the trajectories in the phase space. Insets of Figs. 10(b) and 10(c) show the structure of the two limit cycles. Note that the limit cycle in the inset of (b) is an order of magnitude larger than the other LC. Note also that the time to reach their respective LC is very large. It is of the order of thousands of beats (or periods). This means that the dynamics takes place on a very large time scale compared to one heart beat (⁠$T=480$ ms).

As partial conclusions, in this section, we have shown that a model consisting of a single GJ periodically driven is sufficient to provide bistability. We have determined the range of the parameter FS leading to bistability. The dynamics of the GJ takes place on a very long time scale compared to the basic time scale, i.e., the heart beat period.

Let us now extend the results of Sec. IV A by analyzing a system of two coupled GJs. The system is still relatively simple (see Fig. 11), and we can proceed to a detailed analysis of it. We are especially interested in determining whether the multistability of the two LCs remains present and to what extent.

The equations for the dynamics of the two GJs are as follows:

where $Δϕ1=ΔV1=V|2−V|1=V(t−Δt1)−V(t)$ and $Δϕ2=ΔV2=V|3−V|2=V(t−Δt2−Δt1)−V(t−Δt1)$⁠. Here, again, we use the typical temporal variation of AP as shown in Fig. 7. To proceed to the solution of Eqs. (12) and (13), we need to express the transit times $Δt1$ and $Δt2$ corresponding to the displacement of the AP front from Cell 1 to Cell 2 and from Cell 2 to Cell 3, respectively. To do this, we will use kinematic relations as we did in the one GJ case in Subsection IV A. Here, again, we use the fact that the distance between the cells is held constant to $δx=100μ$m. In the case of several coupled GJs, we have measured in the full simulations that there exist spatial correlations between the AP speeds and the local values of the conductance as shown in Fig. 5. If we want to have a simplified model as close as possible to the full model, we need to introduce the spatial correlations in the AP speed function. Indeed, in the simulations of the full model, we have measured that the AP speed is not a local function of the GJ but rather depends also on the neighboring GJ values. As a first approximation and following what we have observed in Fig. 5, we consider only two GJs in the description of the AP speed. For example (see Fig. 11), for the AP speed between Cell 1 and Cell 2, we have that $c=cf(g1,g2)$⁠, where the subscript $f$ denotes a “forward” approximation. On the contrary, for the AP speed between Cell 2 and Cell 3, we could write the speed as $c=cb(g2,g1)$⁠, where the subscript $b$ denotes a “backward” approximation. We have extracted the data for the AP speed as a function of the GJ $g$-values in the full simulations in order to determine the functions $cf$ and $cb$⁠.

In Fig. 12, we observe a striking similarity between the two functions $cf$ and $cb$ expressing the AP speed as a function of neighboring gap junction values. In this figure, we have measured the speed between the cardiac cell $i−1$ to $i$ across the gap junction GJ$(i−1)$ according to our notation definitions in the computing algorithm. The fact that $cf$ and $cb$ are so similar indicates a strong correlation between $gi−2$ and $gi$⁠, which was already evident from the observation of the full simulation displayed in Fig. 3(b). Indeed, when FS = 2, the alternation of an up–down state in the conductance is the dominant pattern. More importantly, Fig. 12 reveals that a simple local approximation for the AP speed would not be sufficient. Indeed, the AP speed depends on neighboring values of the GJ as it was already apparent from Fig. 5. With this preliminary consideration, we are nearly ready to proceed to the integration of Eqs. (12) and (13). The two kinematic expressions for relating the time delay with the GJ values are as follows:

Because of the large similarity between the functions $cf$ and $cb$⁠, we will use a single function $c$ for the AP speed in the following analysis. Note that AP speed functions are indeed different when varying FS, but in order to keep the model simple, we use the same AP speed function for all the FS values. To continue with the analysis, we can either use a tabular value function defined from the data points or define a fitting function that approximates the data shown in Fig. 12. We will use the latter and define a polynomial approximation for the AP speed as follows:

where the fitting parameters are collected in Table II. Note that for stability reasons of the integration of Eqs. (12) and 13), we cannot just apply a straightforward fit to the data shown in Fig. 12, but we have to impose some constraints on the fitting function. We have found that to ensure stability, we have to add the two conditions $∂gi−1c|{gi−1=0.1,gi=0.1}>0$ and $∂gi−1c|{gi−1=0.1,gi=0.4}<0$⁠. We have used the cftool of Matlab⁴⁸ to perform such a fitting procedure. The resulting fitting function is shown in Fig. 12(a) as a mesh-grid, and the agreement between the data and the function is satisfactory. We computed the coefficient of determination for the fit to be $R2=0.905$⁠.

At this stage, we are ready to perform the study of Eqs. (12) and (13). We have simulated the system of Eqs. (12) and (13) with several fixed values of the parameter $FS$ and many initial conditions for $g1$ and $g2$ in the range [0.05–0.45]. From the simulations, it appears that again, we have multistability, and in this case, we can have up to four different competing attractors coexisting for a fixed value of $FS$⁠. As in the previously discussed single GJ case, the attractors are limit cycles (LCs) with a period corresponding to the period of $V(t)$⁠, which is fixed to $T=480$ ms. Table III gathers the four competing attractors and the numerical values associated with a characteristic point on the LC (Poincaré section).

We study the local stability of the four competing attractors that are identified in Table III in the following manner. We fix the value of the parameter $FS$ and start the simulation with initial values of $g1$ and $g2$ close to the values given in Table III. We proceed to the integration of Eqs. (12) and (13) until we reach a steady state that is defined by a constant Poincaré section in time. We have displayed the results of the local stability analysis in Fig. 13.

From Fig. 13, it is apparent that multistability is present for a much larger range of $FS$ values with respect to the single GJ case considered in Subsection IV A. Indeed, Fig. 13(a) shows that there is multistability for values of $FS$ expanding three orders of magnitude. More precisely, we have indicated in the last column of Table III the range of $FS$ for which we observe local stability of the corresponding solution. It is noteworthy that in the range of $1.24<FS<1.665$⁠, we have coexistence of the four attractors.

Let us illustrate the different basins of attractions of the competing attractors for two selected values of the parameter $FS$ (i.e., $FS=1.4$ and $FS=500$⁠). To create Fig. 14, we simulate Eqs. (12) and (13) starting from a grid of (⁠$101×101$⁠) initial conditions for $g1$ and $g2$ in the range [0.05–0.45]. Each simulation is continued until a steady state is reached (hence, we always have a stable LC). The panels (a) and (b) of Fig. 14 correspond to $FS=1.4$⁠. There, we display the color coded end values of $g1$ and $g2$ after convergence as a function of the initial conditions $g1(0)$ (x axis) and $g2(0)$ (y axis). For the case $FS=1.4$⁠, we observe that the four attractors of Table III coexist as shown in the phase diagram of Fig. 14(c). Depending on the initial conditions, the system will evolve toward one of the four stable LCs. This multistability persists for a large range of FS values. Indeed, panels (d) and (e) of Fig. 14 correspond to $FS=500$⁠. Here, we see that attractors (I)–(III) of Table III coexist as shown in Fig. 14(f). This long-lasting coexistence of multistability while varying the bifurcation parameter FS on such a large range is somewhat uncommon in physical systems and to our knowledge has not been explored yet in detail. Another interesting point to note from Figs. 14(d) and 14(e) is that for large and equal values of $g1(0)$ and $g2(0)$⁠, we are just at the border between two basins of attraction. This means that a very minute noise in the initial conditions will be sufficient to toggle either to attractor (II) or (III).

To summarize this subsection, we have seen that the coupling of two GJs periodically driven also provides multistability and that the range of the parameter FS leading to bistability is greatly enhanced with respect to the single GJ case. Indeed, multistability exists for a range of FS values from approximately 1 to 1000 (three orders of magnitude), which was quite unexpected when compared to the case of the single GJ.

In Secs. IV A and IV B, we have been dealing with the single GJ and two coupled GJ cases. Both cases have shown that the dynamics of the GJ exhibit multistability by modifying the bifurcation parameter FS. However, the range of multistability has been shown to increase by about three orders of magnitude when passing from the single GJ to the two coupled GJ’s cases. This immediately poses the obvious question of what happens when we consider three coupled GJs. Is the multistability range for three GJs comparable to the two GJs case or is it further modified? In this subsection, we address this question. As the setting of the equations for the dynamics of the three GJs case is very similar to the two coupled GJs case, we do not dwell with all the mathematical details in this section. Here, we rather concentrate on the main results. However, it is still good to remind that in the three GJs case, we use the same fitting function [Eq. (16)] for the AP speed as in the two GJs case. One can argue that a more complex function for the AP speed involving three consecutive GJs could have been considered, but to keep the argument simple, we did not examine this possibility here.

When considering three coupled GJs, we have now $23=8$ possibilities for the dynamical attractors. We classify the attractor with a binary notation where digits $1$ and $0$ mean that the final attractor state is with high or low values, respectively. Table IV contains a summary of the different attractors and their local stability. We observe from Table IV that in the range $[1.25<FS<1.58]$⁠, all the eight solutions are locally stable, and therefore, we have coexistence of the eight attractors. The GJ dynamics will select the final state according to the initial conditions.

Figure 15 displays the basins of attraction of the different solutions indexed in Table IV for $FS=1.5$⁠. For this value of the parameter $FS$⁠, we have coexistence of the eight solutions. In order to render the figure easier to analyze, we have regrouped the solutions (II)–(IV) into one group corresponding to all the permutations of the (1, 0, 0) state. In the same fashion, we have regrouped solutions (V)–(VII) into one group corresponding to the permutations of the (1, 1, 0) state. Figure 15 illustrates that the final state of the dynamics depends crucially on the initial conditions of the three GJs (indicated by the three axes of Fig. 15).

As it was the case for the two coupled GJs, we observe that an extensive range of the parameter FS allows for co-existing attractors. The dynamical system consists of several metastable states. Likewise, as the number of cells (and of GJ) increases, the number of metastable states with high and low conductance values increases, resulting in mixed states. Here, we have not pursued this strategy further, and we conclude with the three GJs case. In Sec. V, we have considered a different approach. We take into account the same number of cells as in the simulations of the full model (N = 300) but will consider a simplified description for the action potential.

The simulations of the full model presented in Sec. III were all done with a short one dimensional strand of cardiac tissue $N=300$ cells (corresponding to 3 cm in length). Even with this simple geometry, the simulations are computationally demanding. This is due to the fact that the interesting effect (conductance dispersion) that occurs after a very long transient has elapsed. The reference time scale is the basic forcing period of $T=480$ ms. The typical time scale on which conductance remodeling (“plasticity”) is observed is of the order of thousands or more heartbeats. If one foresees to extend the study to two or three dimensional geometry as it happens for a real heart, we have two options: either use a super-computer to perform the computations or to find a way to speed up the computations by reformulating the model. In this section, we propose a simplified model that retrieves the essential characteristics of the full model but with a gain in computational speed of over four orders of magnitude.

The simplified model has been built upon several assumptions.

The major simplification comes from the decoupling between the AP propagation and the GJ dynamics. This was already assumed in the study of the reduced descriptions in Sec. IV. Again, we will assume that the AP propagates throughout the system with a constant morphology.

The strongest simplifying assumption comes now with the decomposition of the AP morphology into four consecutive straight segments as it is shown in Fig. 16. By doing this, we can exactly integrate the GJ dynamics equations on the four segments and, therefore, avoid the full numerical integration of Eq. (4) with a very small time step $δt=0.01$ (ms). Using this approximation, the time evolution of the membrane potential simply writes as $V=mt$⁠, where $m$ is the corresponding constant slope of the AP taken from Table V and $t$ denotes time. Between two consecutive myocytes in the tissue strand, we have that

where $Δt$ denotes the traveling time for the AP between two consecutive myocytes and $δx$ is the length of the myocyte that coincides with the grid size of the system $δx=0.01$ cm.

Under this last assumption, $Δϕ$ is evaluated and taken constant in each of the four phases of the AP morphology. Therefore, Eq. (4) is reduced to a set of constant coefficient differential equations,

where $g∞$ and $τg$ are time independent, and this equation is readily integrated,

where $B$ is the integration constant. We determine the integration constant by using the initial condition at the beginning of each of the four phases of the AP morphology. We get that $B=g∞−g0$⁠, where $g0$ is the value of $gi$ at the beginning of each time interval. Finally, the expression for the conductance $gi$ as a function of time is given by

The expression in Eq. (20) leads to the time integration of the GJ dynamics in four steps $(td;tp;tr;tf)$ rather than an integration with $Nt=48000$ time steps in the full model. We obtain a piecewise defined function depending on the time interval along the AP morphology. In an algorithmic fashion, we get that if $gi(n)(t)$ denotes the conductance value after the $n$-beat at the GJ located at the position $i$⁠, by using Eq. (20) in the four stages corresponding to the AP morphology, we can compute the value for conductance at the next beat $(n+1)$⁠, $gi(n+1)(t)$ with these four steps,

where the super-indices $0,d,p,r$ indicate the time at the beginning of the corresponding time interval (i.e., $td,tp,tr,tf$⁠). The values for $gi,∞$ and $τi,g$ are also evaluated at the beginning of each of the four time intervals and assumed to be constant in the corresponding interval. The set of Eq. (21) gives a description of the GJ dynamics akin to a discrete map or a cellular automata model,⁴⁹ and this leads to an obvious significative gain in CPU time.

The integration of Eq. (18) requires the knowledge of the AP speed $c$⁠. A first version of the model used the local approximation [Eq. (10)] to relate the local values of the conductance with the AP speed, but this approximation was too crude and the model was unable to reproduce the GJ dispersion. Hence, we considered the non-local relation between the AP speed $c$ and the conductance of neighboring cells given by Eq. (16). This corresponds to an approximation, but it was sufficient for the purpose of the analysis.

In order to use the simplified model as defined in Eq. (21) with the AP speed fitting function defined in Eq. (16), we had to optimize the parameter value of $td$⁠. Indeed, as seen previously, the parameter $td$ defines the depolarization time and is related to the maximum upstroke velocity $V˙M$⁠. This parameter is crucial in the problem, and the GJ dynamics is critically affected by this parameter. There is a large uncertainty on the value of $td$⁠, and we decided to let this parameter to vary freely. We optimize the $td$ parameter in order to bring the results of the simplified model as close as possible to those of the full model.

To quantify the difference between the results of the full and simplified models, we use a combined $L1$ norm,

where $cases$ in Eq. (22) refers to different predictions associated with the two models. In the optimization process, the $L1$ norm will be minimized as a function of $td$⁠. Let us be more explicit about the $cases$ in Eq. (22). They refer to the three values of noise strengths and values for the bifurcation parameter FS up to $104$ as illustrated in Fig. 4.

By setting the parameter to $td=0.791$ ms, we obtain for the conductance transitions as a function of FS a reasonable agreement with the full model. An example for the integration of the simplified model is shown in Fig. 17. The parameters for Fig. 17 are $NS=10−4$ and FS = 2, and one can visually compare the simulation of Fig. 17 with the one displayed in Fig. 3(b). For a quantitative comparison between the simplified and full model, we rather have to look at Fig. 4. One observes that the values of $PUP$ obtained with the simplified model slightly underestimate those of the full model. However, the long intermediate metastable states are well captured by the simplified model as shown in Fig. 4.

We have found that the simplified model reproduces in an acceptable manner the results of the full model. It is rather obvious that a perfect quantitative agreement is not reached, but the simplified model contains approximations of the morphology of the action potential and other approximations that are presumably not accurate enough. However, the gain of four orders of magnitude in computation with the quasi-quantitative agreement is very appealing to pursue this venue further.

This paper has studied the propagation of a cardiac action potential in a one-dimensional fiber with cells electrically coupled through gap junctions. The dynamics of the gap junctions modifies the adjacent cardiac cells’ conductance, which depends on the intercellular voltage difference. Dynamic gap junction conductance has been shown^50,51 to change propagation speeds close to the point of propagation failure. As the AP speed is decreased, the intercellular electrical gradients increase, resulting in reduced conductance and electrical uncoupling. However, one should notice that, due to the slow time scale of the gap junction dynamic, this process is slow, occurring over many heartbeats.

Dynamic electrical uncoupling has been observed experimentally for two coupled cells on which an action potential delay is externally imposed.^52,53 This effect was observed for gap junctions formed by connexins Cx45, which present a stronger dependence of intercellular voltage (a higher value of FS, in the notation we have followed), but not for the connexins Cx40, whose transmembrane voltage dependence curve is broader.¹⁵ 

Here, we have reproduced this dynamical uncoupling at long time scales. Furthermore, we have used a simplified description with a low number of cells to relate the uncoupling to the existence of multistability, thus explaining the previously observed effect²⁴ of the emergence at a long time scale of a very disordered state with inhomogeneous intercellular conductances. The effect on the action potential conduction velocity of the heterogeneities in the conductance distribution has also been addressed.

In addition, we have derived a simplified model for AP propagation that reproduces semiquantitatively the results of the full model at a much lower computational cost. This simplified model will allow producing very long simulations that are needed in the study of the cardiac phenomena occurring at a large time scale.

From a physiological point of view, it is well known that spatial dispersion of repolarization underlies the development of life-threatening ventricular arrhythmias. Spatial dispersion in the cardiac tissue has been associated with ion channelopathies and heterogeneous expression of gap junction proteins.^54,55 These effects are often mitigated by the electrotonic interactions in the cardiac tissue. It would be interesting to incorporate the GJ dynamics, channelopathies, and electrotonic effects together in a mathematical model to test the different scenarios for explaining the emergence of the spatial dynamic heterogeneities.

Dynamic uncoupling is also a well-known phenomenon in neurons, where it is associated with brain plasticity.²⁰ In simulations of a one-dimensional nerve axon and a two-dimensional network of coupled neurons, it has been shown that GJ gating could lead to unidirectional conduction block and wave breaks.²¹ These simulations suggest that a similar effect could be observed in cardiac tissues in situations of decreased electrical coupling, as in ischemia or fibrosis. Our model shows that the GJ dynamics can give rise to heterogeneities in intercellular conductances, even for a perfectly homogeneous tissue. It would be interesting to study how this affects the trigger of arrhythmias.

We have considered a simple Hodgkin–Huxley (HH) type gate model for the description of GJ dynamics. However, more detailed models exist that consider transitions between low and high conductance conformations of the connexins that form the GJ. These more complex models are based on a Markov formalism from 4 up to 16 states.^21,33,34 In those models, the transition probabilities depend on the intercellular potential. The conductance is then obtained as the average of a stochastic ensemble of GJs. In general, these Markov models give similar results to a HH model. However, they naturally incorporate the fluctuations in the GJ conductance. This is an important point because, as we have observed in this paper, the amount of noise in the conductances’ initial states is paramount in determining the final GJ states. Thus, a straightforward extension of our work would be to consider the effects of dynamic noise. Two complementary ways are possible: either by considering a set of stochastic GJs or adding an ad hoc noise term to the HH equations.

We are aware of the limitations of the simple mathematical model presented here. In particular, we have not considered the ephaptic coupling,^51,56 which significantly affects low electrical coupling. There are also two important points that we leave for the future. One is a better understanding of heterogeneous conductance’s effect on AP propagation speed, probably through mean-field homogenization.²⁶ The other is the generalization of the model to two (or three) dimensions, where the possible proarrhythmic effect of the conductance heterogeneities can be assessed. This is, however, not a trivial point due to coupling anisotropy among cells that would require to consider a more detailed model with intracellular spatial discretization.

J.B., I.R.C., A.P., and B.E. acknowledge the financial support by MICINN (Spain) under Grant No. SAF-2017-88019-C3-2-R. J.B. wishes to thank Professor F. B. Sachse and Professor J. Keener (University of Utah) for valuable suggestions and comments on a previous version of the manuscript. D.L. acknowledges partial financial support from FONDECYT (No. 1180905) and from Centers of Excellence with BASAL/ANID financing (No. AFB180001), CEDENNA.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The space constant $λ$ is an important concept in the description of the cardiac action potential propagation. It corresponds to the characteristics of the spatial length scale associated with the action potential wave. In the cardiac syncytium description, one often observes a phenomenon called “current redistribution” between several portions that constitute the cardiac tissue. This redistribution of currents from the intracellular to extracellular space takes place over a spatial extent of a few $λ$⁠.⁵⁷ The space constant $λ$ is often assumed to be approximately 20 times longer than the individual myocyte length,⁵⁸ although the space constant varies depending on the type of cardiac myocyte. Here, we will provide a theoretical estimate of the space constant $λ$ for our model and see how $λ$ is involved in the description of the wave propagation in a simple monodomain formulation.

To start with, we need to define the resistance times unit area $Rm$ (units are $kΩ.cm2$⁠) associated with our specific membrane model,

where $Iion$ is the ionic current through the membrane and $V$ is the transmembrane potential; units are $μA/cm2$ and mV, respectively. The membrane resistance times unit area $Rm$ is dependent on the voltage, but we evaluate its value close to the resting potential $Vrest$⁠. Specifically, in the four current model, we obtain straightforwardly that only one current $Iso$ is non-vanishing close to the resting potential and that

the numerical values used here give that $Rm=7.028kΩ.cm2$⁠. From there, we can readily compute the time constant of the cable equation using the RC time constant $τ=RmCm$⁠, where $Cm$ is the membrane capacitance per unit area and its numerical value is $Cm=1μF/cm2$⁠. Here, we compute that $τ=7.028$ ms, and we have checked that this value is indeed well reproduced in the numerical simulations. To compute the space constant, we can use the definition of the diffusion parameter in the monodomain equation that we have also set as a parameter of the problem,

recalling that we have set $D=1.54⋅10−3cm2$/ms and if we choose $g=0.4$ (dimensionless) as a starting value for the GJ values, we readily obtain for the space constant that $λ≈6.5810−2$ cm. This value is between six to seven myocyte lengths, and it does coincide with what is observed in some test numerical simulations.

Let us illustrate how the space constant shows up in the numerical simulations. We will simulate Eqs. (1) and (2) with a constant value of $g=0.4$ (the dynamics of $g$ is frozen here). We are interested in showing that even with a homogeneous system, the space constant will emerge from the simulations. We compute the local velocity of the propagating AP as follows: at each grid point (separated by a distance $δx=0.01$ cm), we register the time in the depolarization phase when the voltage reaches 60% of its maximum value. Then, the local velocity is computed by $δx/Δt$⁠, where $Δt$ is the transit time of the AP front between successive grid points.

Figure 18 displays the AP propagation results when the GJ is held constant to $g=0.4$ throughout the simulation. It is interesting to note that the AP speed is not constant but exhibits small periodic variations. This phenomenon is known as saltatory conduction. This effect becomes stronger when the conductance is further decreased.⁵⁹ Figure 18(b) shows the auto-correlation of the AP speed, and we observe a marked secondary peak for a spatial lag of eight grid units in agreement with the space constant $λ$ of the problem. The small difference between $λ$ and the lag at which the secondary maximum appears is presumably due to the numerical scheme that we use to integrate the equation. A different numerical scheme will give slightly different results.

Figure 19 displays the maximum upstroke velocity $V˙M$ computed at each grid point in the simulation. Recall that the AP upstroke is a significant factor in the GJ dynamics. Figure 19(a) shows that $V˙M$⁠, in the same manner as the AP speed, oscillates in space. This confirms that even with a homogeneous setting, the system will generate fluctuations again at the scale of the space constant $λ$⁠. Note that the fluctuations in the case of $V˙M$ are one order of magnitude smaller than the fluctuations of the AP speed. Indeed, Fig. 19(a) shows that only the fourth precision digit is affected in this case. Also, important to comment is that the values of $V˙M$ calculated by the present model are somewhat smaller than the known values of $V˙M$ for human cardiac myocytes.⁶⁰ Figure 19(b) confirms that the spatial periodicity of $V˙M$ is approximately equal to eight grid units, as was the case for the periodicity of the AP speed shown in Fig. 18.

As a partial conclusion, we have shown that in the case of simulation where the GJ dynamics was frozen, one still observes spatiotemporal fluctuations in the AP speed and shape. The velocity fluctuations are one order of magnitude larger compared with the distortion of the AP shape. The fact that we observe spontaneous fluctuations due to the AP propagation is important. These fluctuations are the seeds for the separation of the GJ dynamics that we observe in the simulations as shown in Fig. 3(b).