We combined Hodgkin–Huxley equations and gating models of gap junction (GJ) channels to simulate the spread of excitation in two-dimensional networks composed of neurons interconnected by voltage-gated GJs. Each GJ channel contains two fast and slow gates, each exhibiting current–voltage (I-V) rectification and gating properties that depend on transjunctional voltage (Vj). The data obtained show how junctional conductance (gj), which is necessary for synchronization of the neuronal network, depends on its size and the intrinsic firing rate of neurons. A phase shift between action potentials (APs) of neighboring neurons creates bipolar, short-lasting Vj spikes of approximately ±100 mV that induce Vj gating, leading to a small decay of gj, which can accumulate into larger decays during bursting activity of neurons. We show that I-V rectification of GJs in local regions of the two-dimensional network of neurons can lead to unidirectional AP transfer and consequently to reverberation of excitation. This reverberation can be initiated by a single electrical pulse and terminated by a low-amplitude pulse applied in a specific window of reverberation cycle. Thus, the model accounts for the influence of dynamically modulatable electrical synapses in shaping the function of a neuronal network and the formation of reverberation, which, as proposed earlier, may be important for the development of short-term memory and its consolidation into long-term memory.
In models describing synchronization of neuronal networks or the spread of excitation in excitable structures, cell-to-cell communication through gap junction (GJ) channels composed of the connexin (Cx) protein was typically defined as having constant conductance (see online database of models that contain GJs at https://senselab.med.yale.edu/ModelDB/ModelList.cshtml?id=114107&describe=yes&celldescr=&). It is well established that GJ communication (GJC) exhibits changes caused by gating by voltage and chemical reagents as well as current–voltage (I-V) rectification; however, neither of these properties were accounted for during simulations of networks formed of neurons, cardiomyocytes, or other types of excitable cells. GJC plays an important role in the propagation of excitation in the heart, communication between neurons, metabolic exchange between cells in most organs and tissues, including the lens (which lacks blood circulation), and organ formation during development or regulation of cell proliferation (Bukauskas and Verselis, 2004; González et al., 2007).
Electrical synapses between mammalian neurons are preferentially composed of Cx36, which is commonly expressed by neurons throughout the central nervous system (Connors and Long, 2004; Söhl et al., 2005). Many roles for electrical synapses have been documented, such as metabolic communication, synchronization of neuronal networks, and memory formation (Bennett and Zukin, 2004; Bissiere et al., 2011). Also, electrical synapses rectify and cause asymmetry of signal transfer (Furshpan and Potter, 1959; Auerbach and Bennett, 1969; Phelan et al., 2008; Rash et al., 2013). Modulation of electrical synapses can occur by different factors, including phosphorylation (Ouyang et al., 2005; Cachope et al., 2007; Alev et al., 2008), changes in pH (González-Nieto et al., 2008), intracellular concentration of divalent cations (Palacios-Prado et al., 2013, 2014), and activity-dependent long-term depression (Haas et al., 2011).
A property that appears to be common to all 21 members of the Cx family is the sensitivity of junctional conductance (gj) to transjunctional voltage (Vj). Single-channel studies have shown that Vj causes channels to close to a subconductance (residual) state with fast gating transitions (<1 ms; Bukauskas and Weingart, 1994; Moreno et al., 1994). It was shown that Vj as well as chemical uncouplers can also induce gating transitions to the fully closed state and that these transitions are slow, >10 ms (Bukauskas and Weingart, 1994; Bukauskas and Peracchia, 1997). Gating to different levels via distinct fast and slow gating transitions led to the suggestion that there are two distinct Vj-sensitive gates, termed fast and slow (or loop gating) mechanisms (Bukauskas and Verselis, 2004). Earlier, gating properties of GJ channels were described by using the Boltzmann function (Harris et al., 1981), which assumes that GJ channels have two states: open and fully closed. It was presumed that each hemichannel of the GJ channel gates only when the apposed (or docked) hemichannel (aHC) is open and that all Vj falls on the gated hemichannel, whereas in fact Vj falls on both hemichannels and they can gate at the same time. Furthermore, there were several attempts to describe Vj gating of GJs with a multistate approach, which provided more detailed descriptions of gating processes (Vogel and Weingart, 1998; Ramanan et al., 1999; Chen-Izu et al., 2001). The discovery of fast and slow gates (Fig. 1, A and B) motivated us to develop stochastic 4- and 16-state models of GJ channel gating (Paulauskas et al., 2009, 2012).
However, all previous models were adopted to study GJC between two individual cells. Here, we present a model describing a 2-D network of neurons, each exhibiting excitability according to Hodgkin–Huxley (H-H) formalism and interconnected through a 16-state model of GJ channel gating (16SM). We demonstrate that during the spread of excitation, a delay between action potentials (APs) creates a bipolar Vj spike of approximately ±100 mV, causing relatively small decay of gj, which can accumulate into significant gj decays during bursting activity of neurons. Also, we demonstrate that very low gj allows for synchronization of a network composed of neurons exhibiting different intrinsic firing frequencies. Furthermore, we show that a relatively small heterogeneity in I-V rectification of GJs creates conditions for unidirectional AP propagation, leading to reverberation of excitation, which, as proposed (Lorente de Nó, 1933; Hebb, 1950; Tegnér et al., 2002; Constantinidis and Wang, 2004; Wang et al., 2004, 2011), can serve as a form of short-term memory and can mediate consolidation of short-term memory into long-term memory.
MATERIALS AND METHODS
Model of the neuronal network
We linked up to several hundred cells into a 2-D cluster, each exhibiting excitability described using H-H equations (Hodgkin and Huxley, 1952). Each cell in the cluster communicated with neighboring cells through GJs. Junctional current (Ij) and gj were simulated using stochastic 16SM (S16SM) and Markov chain 16SM (MC16SM). We used the sets of 16SM parameters, which are presented in Tables 1 and 2. All models were implemented in MATLAB.
We have modified the original H-H model (Hodgkin and Huxley, 1952) by adding Ij, described by the following system of ordinary differential equations:
Here, Im is the transmembrane current per unit area, Iext is external stimulus current, Vm is membrane potential, and Cm is the membrane capacitance per unit area; VK, VNa, and Vl are the potassium, sodium, and leak reversal potentials, respectively; gK, gNa, and gl are the maximum values of potassium, sodium, and leak conductances per unit area, respectively; n, m, and h are variables (0 ≤ n, m, h ≤ 1) associated with potassium channel activation, sodium channel activation, and sodium channel inactivation, respectively; and αi and βi are rate constants for the respective ion channels.
The membrane capacity Cm was chosen as 1 µF/cm2. The corrected values of reversal potentials were chosen as VK = −12, VNa = 115, and Vl = 10.6 mV, whereas conductances were chosen as gK = 36, gNa = 120, and gl = 0.3 mS/cm2 (Hodgkin and Huxley, 1952). The functions αi and βi were as follows:
We used Euler’s method for the numerical solution of H-H model equations, which was implemented in MATLAB. A time step of 0.01 ms was used in most numerical experiments.
In an S16SM (Paulauskas et al., 2012), each of four gates is characterized by gating properties and the unitary conductance of the open (o) and closed (c) states (Fig. 1 C). We have used an S16SM approach to describe the operation of each gate, allowing a simulation of Ij and gj changes over time. This model is available online at http://connexons.aecom.yu.edu/Applet.htm, and its application in analyzing gating of GJs under normal and pathologic conditions was reported in Palacios-Prado et al. (2009, 2010) and Skeberdis et al. (2011).
In general, the model defines the probability of an individual gate to remain in the same state or to change its state for any given discrete time interval (Δt) and Vj. Closing one gate changes voltage across the other three gates in series, which will affect their probability of changing the state over the next Δt. The transition probabilities between open and closed states po→c and pc→o or probabilities to remain in the same state (po→o and pc→c) are given by
where Pt adjusts probability of transitions during a discrete time interval, Δt, to equilibrate simulated and experimental times (Paulauskas et al., 2009, 2012). K is the equilibrium constant of slow/fast gate, determined as follows: (Chen-Izu et al., 2001). Here, AS/F characterizes sensitivity to voltage, VS/F,o is voltage at which KS,o↔c or KS,o↔res is equal to 1, and Π is a gating polarity (+1 or −1). VS/F is voltage across fast/slow gate.
There is a solid body of evidence demonstrating I-V rectification of open and residual states of the GJ channel and an unapposed hemichannel (DeVries and Schwartz, 1992; Trexler et al., 1996; Oh et al., 1999; Bukauskas et al., 2002a; Kronengold et al., 2003; Valiunas et al., 2004; Tong and Ebihara, 2006; Rackauskas et al., 2007). We assumed that γS,open, γF,open, and γF,res rectify in accordance with the exponential function as proposed in Vogel and Weingart (1998): and where γF,open,0 and γF,res,0 are γF,open at VF = 0, whereas γs,open,0 is γF,open at VS = 0. RF,open, RF,res, and RS,open are rectification coefficients of the open/residual states and fast/slow gates, respectively.
Rectification of γF,open, γF,res, and γS,open does not allow for the evaluation of VF and VS in one step and necessitates the use of an iterative procedure to estimate the conductance of the GJ channel. In brief, initially we know the states of each gate, which allows for the estimation of VF and VS without accounting for rectification:
Then, we recalculate conductances of gates accounting for rectification and reevaluate a new set of VF and VS. We repeat the last cycle and at the nth iterative step evaluate conductance of the GJ channel (γj,n) as well as the difference, Δγj,n = γj,n − γj,n−1 for convergence testing. The value of γj,n is given by
The iterative procedure is stopped once the tolerance criterion Δγj,n < (γj,n−1/1000) is reached. The final values of VF and VS allow for the evaluation of possible changes in the state of each gate, as well as estimation of voltages across all four gates, the sum of which should be equal to Vj, i.e., Vj = VF,A + VS,A + VF,B + VS,B.
This information is transferred to the next discrete time interval and the process is repeated. In >99% of calculations, the number of iterations was four to seven. Once the effect of I-V rectification on the unitary conductance of each GJ channel γj,i of the electrical synapse composed of m channels is evaluated, then gj is estimated from
To accelerate the simulation and reduce noise, we developed an MC16SM by transforming an S16SM. Instead of evaluating gating of each individual GJ channel, we calculated changes over time of averaged probabilities for each of four gates to be open or closed depending on Vj by using a Markov chain approach. Earlier, we used stationary/steady-state probabilities of Markov chain to model GJ voltage gating (Sakalauskaite et al., 2011; Snipas et al., 2015). However, a steady-state approach is not suitable in this case, because modeling of neuronal networks requires evaluation of dynamic changes in gj. Thus, we used transitive rather than steady-state probabilities. In doing so, we used a recursive relation, pn = pn−1 · P, to calculate the transitive probability vector pn at a discrete time moment, tn, for each of 16 states, where P denotes the transition probability matrix of a discrete-time Markov chain. Each state sj can be expressed as a 4-tuple, sj = (sF,A; sS,A; sS,B; sF,B), where each component denotes the state—open (o) or closed (c)—of the respective gate in the GJ channel. Thus, the model consists of 16 states, which can be numbered in lexicographical order as follows: s1 = (o;o;o;o); s2 = (o;o;o;c); s3 = (o;o;c;o); s4 = (o;o;c;c); s5 = (o;c;o;o); s6 = (o;c;o;c); s7 = (o;c;c;o); s8 = (o;c;c;c); s9 = (c;o;o;o); s10 = (c;o;o;c); s11 = (c;o;c;o); s12 = (c;o;c;c); s13 = (c;c;o;o); s14 = (c;c;o;c); s15 = (c;c;c;o); and s16 = (c;c;c;c).
The transition probabilities among 16 states can be expressed by using the same transition probabilities po→c and pc→o of singular o↔c transitions, as in S16SM. For example, transition probability from s4 = (o,o,c,c) to s6 = (o,c,o,c), which we denote as poocc→ococ, can be estimated as follows: poocc→ococ = (1 − po→c) · po→c · pc→o · (1 − pc→o).
Thus, all basic parameters of the model necessary for building a transition probability matrix P are the same as in the S16SM model, and the same equations are used to estimate K, evaluate I-V rectification, and calculate γj. However, because voltage across each gate depends on the state sj, i.e., VF = VF(sj) and VS = VS(sj), then γj must also be evaluated at each state, i.e., γj = γj(sj). Once the transition probability matrix P is formed and the probability vector pn is estimated, then the expected value of junctional conductance gj is estimated as
Fig. 2 shows Ij (B) and gj (C) traces generated in response to a Vj step of −60 mV (A) using an S16SM (red traces) and an MC16SM (blue traces).
Online supplemental material
Fig. S1 supplements Fig. 4 C. Fig. S2 supplements Fig. 5. Figs. S3–S6 are snapshots of Fig. 8 and Videos 1–4 and provide explanations of conditions under which these videos were recorded. A more detailed description of 16SM is presented in the supplemental materials and methods.
Physiology of fast and slow gates
The demonstration of fast and slow gates (Bukauskas and Verselis, 2004) required consideration of four gates in series, which resulted in the development of an S16SM of gating (Paulauskas et al., 2012). The initial hypothesis of fast and slow gating mechanisms was based on recordings of de novo GJ channel formation observed after providing two patched insect cells (Aedes albopictus, clone C6/36) in contact (Bukauskas and Weingart, 1993, 1994). Based on these data, it was concluded that each hemichannel of the GJ channel possesses a fast gate controlled by Vj that closes the channel from the open state with conductance γj,o to a substate or residual state with conductance γj,res and a slow gate that operates between an open state with conductance γS,open and a fully closed state with conductance γS,closed = 0 (Fig. 1 A; Bukauskas and Weingart, 1994). Furthermore, gating of GJs to the substate were reported in mammalian cells expressing different Cx isoforms (Moreno et al., 1994; Veenstra et al., 1994; Bukauskas et al., 1995; Bukauskas and Peracchia, 1997).
Furthermore, there is a solid body of evidence demonstrating I-V rectification of open and residual states of the single GJ channel (Bukauskas et al., 2002a; Rackauskas et al., 2007) and undocked/unapposed hemichannels (DeVries and Schwartz, 1992; Trexler et al., 1996; Oh et al., 1999; Kronengold et al., 2003; Valiunas et al., 2004; Tong and Ebihara, 2006). I-V rectification of gates is typically reflected in the nonlinear instantaneous gj (gj,inst) dependence on Vj measured at the very beginning of applied Vj steps when Vj gating is not yet activated. Homotypic GJ channels, which exhibit oppositely oriented near-exponential relationships of γF,open or γS,open over Vj in aHCs, demonstrate a relatively small gj,inst dependence on Vj, as shown in Fig. 7 of Paulauskas et al. (2009). However, if cells express different Cx isoforms and form heterotypic GJ channels or are exposed to asymmetric conditions in ionic composition or phosphorylation, which can differently affect I-V rectification of fast and/or slow gates residing in aHCs, then an asymmetry of I-V rectification of GJ channels can be well expressed, causing unidirectional signal transfer and reverberation of excitation in neuronal networks, as we demonstrate later in Figs. 7 and 8.
Simulation of signal transfer through electrical synapses
Studies of electrical cell–cell coupling in brain slices are most often performed by injecting a pulse of constant hyperpolarizing current into one cell and measuring the electrotonic response in a neighboring cell. Initially, we simulated electrotonic signal transfer between two cells connected through soma-somatic junctions (Fig. 3 A). Cell 1 was injected with an external current (Iext) under current clamp conditions, and the changes in transmembrane voltages were estimated in cell 1 (Vm,1) and cell 2 (Vm,2). Fig. 3 B shows that at gj = 0.2 nS, electrotonic responses in cell 1 (V1) and cell 2 (V2) reached steady states after ∼30 ms. The ratio between V2 and V1 is named as a coupling coefficient, K12 = V2/V1. The bottom plot in Fig. 3 A shows the K12–gj relationship simulated for soma-somatic GJs; the gray circle shows averaged K12 (∼0.15) and gj (∼0.2 nS) values measured in rat cerebellar cortex basket cells (Alcami and Marty, 2013). However, for dendro-dendritic GJs located relatively distantly from somas (more than ∼100 µm), an influence of axial resistance of the dendrites (Rd in Fig. S1 B and supplemental materials and methods) should be accounted for. Our preliminary data demonstrate that gj values are increasingly undervalued with an increase in the distance of electrical synapse from somas of communicating neurons (unpublished data).
In testing whether the strength of electrical synapses can cause gating, we initially assumed that cells are connected through the relatively highly Vj-sensitive Cx45, which, like Cx36, is expressed in neurons, but to a lesser degree (von Maltzahn et al., 2004; Chapman et al., 2013). Fig. 3 C shows that an Iext of −18 pA caused changes in Vm,1 and Vm,2, resulting in a Vj of ∼30 mV; simulation was performed using MC16SM. The gj trace demonstrates that gj decayed, and it took ∼600 ms to reach a steady state at ∼50% of the initial gj. With the decrease in gj, Vm,1 increased and Vm,2 decreased (Fig. 3 C, left inset in Vm time plot), thus causing Vj to rise, which accelerated a decay in gj. Furthermore, after an Iext step, the AP developed in cell 1 and was transmitted to cell 2 with a ∼2-ms delay that generated a bipolar Vj pulse (Fig. 3 C, inset on Vj time plot), resulting in a small decrease in gj (arrow in the inset of gj time plot). Assuming that GJs are formed of the less Vj-sensitive Cx36 (VF,o and VS,o were 40 mV instead of 10 mV for Cx45; Bukauskas et al., 2002b; González et al., 2007; Marandykina et al., 2013), then the gj decrease was ∼10-fold lower (unpublished data). Gating parameters of Cx36 and Cx45 are shown in Tables 1 and 2.
Fig. 3 D shows data simulated using S16SM instead of MC16SM; parameters of gates were the same as in Fig. 3 C. It was assumed that the electrical synapse contained 10 GJ channels with a single-channel conductance of 25 pS. This revealed stepwise rather than smooth changes in gj. It is notable that the unitary stepwise changes in gj resulted in slow changes in Vm,1 and Vm,2 (Fig. 3 D, left inset in gj time plot).
To assess the extent to which the Vj values that develop during spiking activity can affect gj, we exposed cell 1 of the cell pair to Iext of +15 pA (Fig. 4 A), which caused a burst of APs at 63 Hz in cell 1 (Fig. 4 B, red trace). Then, excitation was transferred to cell 2 (Fig. 4 B, blue trace) with a ∼2-ms delay (Fig. 4 B, inset, ↔). In this simulation, VF,o = VS,o = 10 mV, which roughly corresponds to Cx45. Fig. 4 C shows repeated Vj values arising mainly during the phase shift between APs in neighboring cells, which we term Vj,AP. The amplitude of the Vj,AP spike can reach approximately ±100 mV, which caused a ∼2% decrease in gj (Fig. 4 D, inset on gj trace). During bursting activity of neurons, small reductions in gj accumulate, giving a ∼28% decrease in gj at the steady state (Fig. 4 D). However, Fig. S1 A (see also supplemental materials and methods) shows that at higher gj values (2 instead of 0.36 nS), the phase shift between APs in neighboring cells decreases, resulting in substantially shorter and smaller-amplitude Vj,AP spikes. Consequently, this caused a ∼30% smaller drop in gj than that shown in Fig. 4 D.
At VF,o = VS,o = 40 mV, which resembles Cx36 properties, a decrease in gj during one Vj,AP spike was ∼0.35%, which accumulated during bursting activity, resulting in a ∼2.5% decrease in gj (Fig. 4 E). Fig. 4 F shows that a greater decrease of gj in cells expressing Cx45 correlates with a larger delay of AP transfer than in Cx36-expressing cells. Fig. 4 G shows how the reduction of gj depends on VF,o and VS,o, which roughly replicate gating of different Cxs at Iext = +15 pA (squares; spiking activity is present). Thus, gj can decrease ∼30% in highly Vj-sensitive Cx45 or Cx57, whereas for less Vj-sensitive Cxs, e.g., Cx36 or Cx26, the decrease in gj was much smaller. Even though electrical synapses are formed preferentially of Cx36 GJs exhibiting low Vj sensitivity under normal conditions, their sensitivity to Vj increases significantly when intracellular Mg2+ ([Mg2+]i) is elevated to 5 mM (Palacios-Prado et al., 2013), which would lead to larger reductions in gj during bursting activity of neurons. Also, a robust increase in sensitivity to Vj was demonstrated for Cx45 and Cx57 during acidification (Palacios-Prado et al., 2009). Hence, the decrease in gj that occurs during bursting activity depends not only on Cx type but also on the surrounding conditions influencing the functional state of the cell.
To examine whether Vj,AP spikes and changes in gj are different for dendro-dendritic junctions, we simulated the cell pair with a principal electrical scheme as shown in Fig. S1 B and the supplemental materials and methods. Parameters of fast and slow gates as well as the value of gj were identical to those of the soma-somatic junction. Fig. S1 C shows that Vj,AP spikes obtained for soma-somatic (red) and dendro-dendritic (green) junctions are similar in amplitude and duration, whereas both cause gj decay similar to that shown in Fig. 4 D.
Role of electrical synapses in synchrony of neuronal networks
Synchronization of neuronal networks can be achieved through chemical (Brunel and Hakim, 1999; Juang and Liang, 2014) as well as electrical (Hormuzdi et al., 2004) synapses. It is known that neurons coupled through electrical synapses can exhibit synchronous spiking activity while being separated by several hundred micrometers (Alcami and Marty, 2013). Otherwise, it is well established that electrical cell–cell coupling between neurons is relatively weak. In neuronal networks connected through dendro-dendritic GJs, gj values are <1 nS (∼0.13 nS in Vervaeke et al. , 0.02–0.2 nS in Venance et al. , and 0.13–1.0 nS in Devor and Yarom ). Higher values of gj were reported for soma-somatic junctions (mean ∼2.1 nS between low-threshold spiking thalamocortical neurons in Gibson et al.  and ∼6 nS in clusters of mesencephalic trigeminal nucleus neurons in Curti et al. ). Relatively small gj values of electrical synapses are most likely caused by a very small single-channel conductance of Cx36 GJ channels (∼6 pS; Moreno et al., 2005) and the small size of junctional plaques (JPs; Rash et al., 2013).
We examined a minimum gj required to achieve phase synchronization/lock in a cluster of two or more cells (gj,min-synchr), each of which exhibits spontaneous activity but differs in firing rates. Fig. 5 A (b) shows that at gj = 0, cell 1 exposed to Iext of 35 pA exhibited APs at 95 Hz, and cell 2 exposed to Iext of 12 pA exhibited APs at 65 Hz. When gj was increased and reached 0.26 nS, cells exhibited phase-locked oscillations (Fig. 5 A, c) at an intermediate frequency of 85 Hz. Fig. 5 A (a) shows gj,min-synchr dependence on the maximum difference in firing rates (ΔFR) among cells in the 2-D clusters; random distributions of firing rates were created by using different Iext values. For each 2-D cluster consisting of n × n cells, five different numerical experiments were performed. Parameters of an MC16SM were selected assuming that GJs are formed of Cx36 (Table 1). The central cell of the cluster was exposed to the highest Iext value, whereas values between minimum and maximum Iext and assignments to other cells were selected randomly. Fig. 5 B shows that gj,min-synchr approaches zero at small ΔFR and increases as ΔFR rises. It also demonstrates that averaged values of gj,min-synchr tend to be higher for bigger clusters but approach steady state at the cluster size of ∼7 × 7. Also, gj,min-synchr values of 1-D or cable-like (1 × n) clusters were at least ∼30% higher than those of 2-D clusters.
At gj = 0, if cell 1 of the cell pair is exposed to 2-ms-long depolarizing pulses of 30 pA at 70 Hz, then it exhibits regular firing activity while other cells remain silent. When gj was increased and reached 0.1625 nS, cell 2 responded with an AP to every AP in cell 1. Summary of Fig. 5 C shows that the minimum gj necessary for one-to-one AP transfer (gj,min-transf) in the cell pair increases with rising frequency of Iext pulses. Also, the values of gj,min-transf increased as the size of the 2-D cluster increased from 2 × 2 to 7 × 7. In all experiments, the cell exposed to repeated stimulation was located maximally close to the center of the cluster. There was a clear convergence to a steady state when the size of the cluster reached 4 × 4. Furthermore, there is well-expressed decay of a cutoff of frequency from ∼95 Hz for the cell pair to ∼70 Hz for 4 × 4 to 7 × 7 clusters with an increase in size of the cluster. Also, at ∼50 Hz, there is a well-expressed minimum of gj,min-synchr for all sizes of clusters, which is most likely related to the enhanced amplitude (∼8%) of APs initiated at the period of excitation of 20 ms.
To evaluate how Vj gating affects synchronization of neuronal networks, we examined a network interconnected through (a) constant resistors, (b) GJs formed of low-Vj-sensitivity Cx36, and (c) GJs interconnected through high-Vj-sensitivity Cx45. Fig. S2 shows that the difference in gmin,synch between Cx36 and the constant resistors is negligible (<2%). However, for Cx45, gmin,synch was more than two times higher. We expect that a similar increase in gmin,synch can occur for Cx36 under intracellular hypermagnesemia (Palacios-Prado et al., 2013, 2014), which increases Vj-sensitive gating of Cx36.
Unidirectional AP transfer between cells connected through GJs exhibiting I-V rectification
Cells, despite well-expressed homeostasis, differ in pH, electrolytic balance, phosphorylation, etc., which may cause an asymmetry in cell-to-cell signaling. Asymmetry in GJC was observed in approximately one third of cases between striatal output neurons (Venance et al., 2004) and between neurons in the inferior olive nucleus (Devor and Yarom, 2002). Recently, we reported that asymmetry in [Mg2+]i leads to rectification of gj,inst (Palacios-Prado et al., 2014). To explore how I-V rectification can affect signal transfer between cells, we simulated AP transfer in three linearly arranged cells connected through an MC16SM (Fig. 6 A). We assumed that the junction between cells 1 and 2 exhibited well-expressed I-V rectification as shown in the red gj,1-2–Vj plot (Fig. 6 B), whereas the junction between cells 2 and 3 was not rectifying (blue gj,2-3 line); I-V rectification coefficients of GJs were equal to 150 and 10,000 mV, respectively. Cell 2 was excited by a 50-pA pulse 1 ms in duration. This caused an AP in cell 2 (Fig. 6 C). Difference in Vm between cells 1 and 2 resulted in short-lasting Vj with an amplitude reaching ∼100 mV. The orientation of I-V rectification was such that at positivity on the cell 2 side, gj,1-2 decreased (red trace), causing a block of AP transfer to cell 1 (green Vm trace), whereas excitations were transferred to cell 3. Stimulation of cell 1 (Fig. 6 D) induced an AP during which cell 1 became more positive than cell 2 and, because of I-V rectification and oppositely oriented Vj, gj,1-2 increased, which facilitated AP transfer to cell 2 and onto cell 3. Thus, relatively weak I-V rectification of GJs can lead to unidirectional signal transfer between neurons.
Reverberation of excitation in a cluster of cells connected through electrical synapses
It is well established that unidirectional signal transfer or an interruption/brake of the wave of excitation can initiate reverberation or reentry, causing heart fibrillation (Davidenko et al., 1995). Here, we present data demonstrating formation of reverberation in a cluster of nine neurons connected through electrical synapses (Fig. 7). Assuming that GJs do not exhibit I-V rectification and gj values of electrical synapses are distributed randomly in a range of 0.175–0.2 nS, then an AP evoked in cell 1 spreads through the cluster and terminates (Fig. 7 B). When I-V rectification was ascribed between cells 1 and 4 and between cells 2 and 5, leaving other GJs lacking rectification, then an Iext pulse of 25 pA (Fig. 7 C, a) initiated an AP in cell 1 that was transmitted to cell 2, but not to cell 4 (Fig. 7 C, d). Consequently, the excitation spread along the yellow circle (Fig. 7 A). After reexcitation of cell 1, the process repeated many times. Interestingly, cell 5 exhibited only low-amplitude electrotonic depolarizations (Fig. 7 C, c), which is typically observed at the core of spiral waves (Davidenko et al., 1992). Reverberation was stopped by a second pulse of 10 pA (Fig. 7 C, a), causing a depolarization of ∼3 mV (Fig. 7 C, e, red arrow) and preventing reexcitation of cell 1. To examine how termination of reverberation depends on the amplitude and timing of Iext, we applied the second pulse at different phases in the reverberation cycle. Gray circles and curves show time windows (Fig. 7 C, f) during which reverberation was stopped using positive or negative Iext pulses. Thus, termination of reverberation requires relatively low levels of excitatory or inhibitory currents locked to a specific phase of reverberation cycle. Once initiated, reverberation does not need inhomogeneity of I-V rectification to continue.
Reverberation of excitation in a 15 × 15 cluster is shown in Fig. 8 (see Fig. S3 and Video 1). A single stimulus applied to cell 1/1 (Fig. 8 A, red circle) caused reverberation of excitation, with its core (Fig. 8 A, dashed circle) identical to that in Fig. 7. During this simulation, it was assumed that all junctions exhibited identical gj, but only two junctions shown by light-blue segments exhibited an asymmetry of I-V rectification equal to that shown in Fig. 6 B. Colored isopotential areas (see calibration bar for transmembrane potential) show the spread of excitation (dashed arrows) across the cluster. A stimulus applied to cell 8/8 (Fig. 8 B, red circle) caused formation of two endless reverberation circles spreading across the cluster. In this experiment, it was assumed that junctions indicated by light-blue segments exhibited I-V rectification (see Fig. S4 and Video 2). Cells located in the center of the reverberation core (Fig. 8, A [cell 2/2] and B [cells 9/9 and 7/9]) exhibited low-amplitude electrotonic responses similar to those shown in Fig. 7 C (c), whereas all other cells in the cluster generated ordinary APs. To examine whether the spread of excitation in a 15 × 15 cluster is affected by boundary conditions, we wrapped the cluster on a virtual torus by connecting all cells in column 1 with cells in column 15 and all cells in row 1 with cells in row 15. Reverberation was initiated by a single stimulus applied to cell 4/8 (red circle); I-V–rectifying junctions are indicated by light-blue segments. Fig. S5 and Video 3 show a returning excitation wave that spreads from right to left and collides with the wave spreading from the reverberation core. These data show that for clusters of ∼15 × 15 and larger, the boundary conditions do not noticeably affect the spread of excitation and formation of reverberation. Our data also show that reverberation can be initiated independently whether the region of I-V rectification is located in the corner of the cluster or more centrally. However, as shown in Fig. S6 and Video 4, the presence of multiple I-V rectification inhomogeneities can cause local deformations of spreading spiral waves when passing these regions. Larger regions of I-V inhomogeneity might break down the front of excitation, leading to formation of a new reverberator (not shown), as was demonstrated in computational model of cardiac arrhythmias (Casaleggio et al., 2014).
In all reported models describing synchronization or phase-locked activity of neuronal networks, as well as the spread of excitation in various excitable structures, GJC is defined as constant. However, it is well established that fast and slow gates of GJ channels are modulated by voltage, chemical uncouplers, and coupling-promoting reagents. These changes occur in a timescale from hundreds of milliseconds to minutes. Furthermore, GJ channels exhibit I-V rectification, which changes the single-channel conductance instantaneously. Hence, the regulation of GJ channels is highly complex, and an influence of Vj gating as well as I-V rectification should be accounted for, especially during high spiking activity of neuronal networks. We acknowledge that 2-D models present a simplified architecture of neuronal networks; however, it was reported that a 2-D lattice describing a network of basket cells is physiologically adequate (Alcami and Marty, 2013). The 2-D lattice is also extensively used to model the tissues of cardiomyocytes (Pertsov et al., 1993; Casaleggio et al., 2014) or in numerical modeling of epileptic seizures (Volman et al., 2011).
In neuronal networks, Vj can change dynamically because of a phase shift between APs in neighboring cells, which we name Vj,AP. One of our major goals was to find whether Vj,APs that last only a few milliseconds can influence GJC through Vj gating and how such influence depends on the Cx isoform. The spikes of Vj,APs, despite their brevity, are high in amplitude (>100 mV) and potentially influence open probability of fast and/or slow gates because of Vj gating and single-channel conductance caused by I-V rectification. Indeed, data shown in Fig. 4 D demonstrate that for Cx45, a single Vj,AP can cause ∼2% decay in gj. During bursting activity of neurons, small gj decays accumulate, leading to significant decays of gj, as is shown in Fig. 4 D.
Although mammalian electrical synapses are preferentially composed of Cx36 GJs with low sensitivity to Vj under normal conditions, their sensitivity to Vj can increase significantly at enhanced [Mg2+]i (Palacios-Prado and Bukauskas, 2012; Palacios-Prado et al., 2014). This suggests that synchronization of neuronal networks may be reduced under ischemic/stroke conditions when an ATP deficit leads to an increase in [Mg2+]i. An even stronger rise in sensitivity to Vj is demonstrated for Cx45 and Cx57 during acidification (Palacios-Prado and Bukauskas, 2009; Palacios-Prado et al., 2010). Because Cx45 is mostly expressed in the inferior olive (Van Der Giessen et al., 2006), whose role in coordination and control of movement is well established, our data predict that gj in this structure should be affected by activity of neurons more than in most other structures that preferentially express Cx36. Because Cx36 is resistant to spike-induced depression (Fig. 4 E) and acidification (González-Nieto et al., 2008), synchronization of neuronal networks during paroxysmal activity or ischemia can be maintained. This may explain why Cx36 predominates over Cx45 and other Cxs in neurons of mammalian brains.
Electrical synapses are expressed in a majority of mammalian brain regions. Along with chemical synapses, they facilitate synchronous or phase-locked activity of neuronal networks. Direct GJC between neurons in hippocampus, neocortex, amygdala–hippocampus–cortical axis, and other areas of the mammalian brain is well established, and knowledge on the regions expressing Cx36 has expanded extensively over the last decade because of application of better antibodies and imaging technologies with enhanced resolution (Bennett and Zukin, 2004; Connors and Long, 2004; Hormuzdi et al., 2004; Saraga et al., 2006; Bissiere et al., 2011). Electrical synapses, by forming dendro-dendritic junctions, participate in synchronizing activity of neuronal networks (Hestrin and Galarreta, 2005) spatially separated by several hundred micrometers (Fukuda and Kosaka, 2000; Saraga et al., 2006). In the somatosensory cortex, synchrony is found among inhibitory interneurons at separations of up to 400 µm (Fukuda, 2007), whereas in the primary visual cortex, this distance was from 500 to 800 µm. In the prefrontal cortex, synchronized firing during working-memory tasks appeared to be within 500 µm (Bissiere et al., 2011).
Synchronization of electrically coupled oscillators was examined theoretically in multicellular networks, and it was shown that relatively low gj can be sufficient to synchronize or phase-lock oscillatory activity of neurons exhibiting different oscillatory activities (Saraga et al., 2006). In different regions of the brain, electrical cell–cell coupling between neurons connected through dendro-dendritic GJs is relatively weak and constitutes ∼0.13–1.0 nS (Devor and Yarom, 2002; Venance et al., 2004; Vervaeke et al., 2010). Such GJC satisfies the conditions necessary for AP transfer between neurons, but it is only a few fold higher than gj,min-transf, which is in the range of ∼0.1–0.3 nS and varies depending on the firing rate of pacemaker cells (Fig. 5 C). However, such GJC is >10-fold stronger than the gj necessary for synchronization neurons if their intrinsic oscillatory activity is contiguous (Fig. 5 B).
It is evident that electrical synapses with a conductance of >10 nS would bypass function of chemical synapses and compromise information processing in neuronal networks. Presumably, this is one of reasons that electrical synapses are formed preferentially of Cx36, which exhibits the smallest single-channel conductance among all members of the Cx family (Moreno et al., 2005); <1% of Cx36 GJ channels clustered in junctional plaques are functional under normal conditions (Curti et al., 2012; Marandykina et al., 2013). An increase in open probability of Cx36 GJs during intracellular hypomagnesemia (Palacios-Prado and Bukauskas, 2012; Palacios-Prado et al., 2014), inhibition of arachidonic acid synthesis, or weak acidification (González-Nieto et al., 2008) can strengthen electrical synapses, thus resulting in hypersynchronization, which is linked to development epileptic seizures (Carlen et al., 2000; Gajda et al., 2003).
Data shown in Fig. 6 demonstrate that relatively weak rectification of GJ channels can change gj instantaneously during short-lasting Vj,AP that can lead to unidirectional signal transfer and consequently to reverberation of excitation (Figs. 7 and 8). Thus, I-V rectification can very quickly reshape the functional state of a neuronal network in response to dynamic changes of bursting activity and initiate self-sustained reverberation. Furthermore, it was proposed that I-V rectification of GJs is defined by an asymmetry of fixed charges, assumed to be charged amino acids, inside the channel pore, and described using the Poisson–Nernst–Plank equation (Oh et al., 1999). Recently, we demonstrated a Mg2+-dependent plasticity of electrical synapses between neurons in the mesencephalic trigeminal nucleus (Palacios-Prado et al., 2013) and in the thalamic reticular nucleus (Palacios-Prado et al., 2013, 2014). Also, we showed a novel form of I-V rectification in Cx36 that depends on Mg2+ binding within the pore (Palacios-Prado et al., 2013, 2014). It was assumed that a domain of extracellular loop 1, containing D47 residue, forms a Mg2+-binding/interaction site, which at higher [Mg2+]i (∼5 mM) creates an exponential I-V rectification, which is oriented in the opposite direction than I-V rectification caused by fixed charges. Superposition of these two forms of I-V rectification results in a hyperbolic Ij,inst-Vj relationship measured uniquely in Cx36 GJ channels (Palacios-Prado et al., 2014). In addition, it was shown that at low [Mg2+]i (∼10 µM) in both cells, the Ij,inst-Vj relationship becomes linear, whereas transjunctional asymmetry in [Mg2+]i can cause an exponential Ij,inst-Vj relationship. This is caused by a presence of I-V rectification in the hemichannel exposed to 5 mM [Mg2+]i and absence of I-V rectification in the hemichannel exposed to 0.01 mM [Mg2+]i. The I-V rectification in electrical synapses has been proposed to affect the dynamic output of neuronal networks, and therefore this novel Mg2+-dependent I-V rectification could explain the switching between firing states and changes in the output of neuronal networks during different metabolic states when [Mg2+]i is strongly affected (Gutierrez and Marder, 2013).
We demonstrate that reverberation of excitation in a multicellular network of neurons interconnected only through electrical synapses and inhomogeneities of I-V rectification can be initiated and terminated by a single excitation of only one cell. Previously, Lorente de Nó (1933) proposed the concept of the continuous reverberation of excitations in the circuit formed from several neurons to explain the vestibulo-ocular reflex. This allows for the generation of a burst of APs in response to a single excitation and was viewed as a mechanism for the active maintenance of working memory. However, this hypothesis is not listed in the latest views on the mechanisms of short- and long-term memories preferentially related to molecular, structural, and functional changes of chemical synapses (Kandel et al., 2014; Bailey et al., 2015). Later, Lorente de Nó’s idea was extended in theoretical studies demonstrating a possibility for reverberation of excitation in neuronal networks interconnected with chemical synapses (Tegnér et al., 2002; Constantinidis and Wang, 2004; Wang et al., 2004, 2011). However, there are no studies demonstrating, either experimentally or theoretically, an involvement of electrical synapses in the realization of reverberation in mammalian brain. On the other hand, reverberation of excitation was demonstrated in crayfish between segmented giant axons connected exclusively through electrical synapses (Watanabe and Grundfest, 1961). The role of reverberation of excitation in memory mechanisms was expanded by the proposal (Hebb, 1950) that memory consolidation in neuronal networks includes two basic mechanisms, post-acquisition reverberation of excitation and structural synaptic plasticity (Ribeiro and Nicolelis, 2004). It was proposed that reverberations in neuronal circuits consistently increased during “slow-wave” sleep phase, which was related to consolidation of short-term memory into long-term memory. Reverberation decreased during wakefulness, which suggests that sensory input might suppress reverberations. All studies relating the role of reverberation in short-term memory and its consolidation into long-term memory are based on circuits in which neurons are connected through chemical synapses (Tegnér et al., 2002; Constantinidis and Wang, 2004; Ribeiro and Nicolelis, 2004; Wang et al., 2004, 2011). Our data show that reverberation of excitation can also form in circuits interconnected through electrical instead of chemical synapses and plays a specific role in a complex mechanisms of memory. In such networks, reverberation can also be terminated by low-amplitude phase locked external stimulation (Fig. 7 C, f), resembling activation of sensory signals during a waking period. We expect that reverberation can also exist in circuits consolidating electrical and chemical synapses, in which the presence of even a single chemical synapse in a chain of neurons connected through electrical synapses can trigger unidirectional signal transfer and consequently reverberation of excitation.
Data presented in Fig. 7 show that the center of reverberation core cells did not exhibit APs but instead displayed multiple low-amplitude depolarizations that appeared to be electrotonic responses from APs in neighboring cells. Such behavior is characteristic only for reverberation and cannot be observed during the high frequency caused by enhanced automaticity. Also, low-amplitude depolarizations are a hallmark of spiral wave dynamics (Davidenko et al., 1992), which have been associated with cardiac arrhythmias (Pertsov et al., 1993). It has been detected that epileptic seizures may also manifest as recurrent spiral waves in cat neocortex (Stacey, 2012). Interestingly, a termination of reverberation requires relatively low levels of excitatory or inhibitory current, whereas active windows are locked to specific phases of the reverberation cycle (Fig. 7 C, f). Again, such a property applies only to reverberation of excitation and not to enhanced automaticity, which cannot be terminated by phase-locked stimulation. Both these criteria can serve to differentiate reverberation from enhanced automaticity in the experiment.
The role of GJs in the structural synaptic plasticity is less known; however, several studies suggest a link between GJ activity and memory consolidation via the process of phosphorylation. Cx phosphorylation can influence GJC in cell cultures (Márquez-Rosado et al., 2012) and hippocampus (Nakai et al., 2014). Indeed, phosphorylation/dephosphorylation could combine both necessary conditions of Hebb’s hypothesis, synaptic plasticity and reverberation of excitation, into a single theoretical mechanism: changes caused by Cx phosphorylation/dephosphorylation can increase I-V rectification, which allows unidirectional signal transfer between cells and, consequentially, enables reverberation of excitation.
We thank Michael V.L. Bennett and Vytautas Verselis for helpful comments and discussions. We thank Angele Bukauskiene, Kevin Fisher, Alis Dicpinigaitis, and Deesha Shah for excellent technical assistance.
This study was funded by a grant (MIP-76/2015) from the Research Council of Lithuania and by National Institutes of Health grant R01NS 072238 to F.F. Bukauskas.
The authors declare no competing financial interests.
Richard W. Aldrich served as editor.