Reverberation of excitation in neuronal networks interconnected through voltage-gated gap junction channels.

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.

INTRODUCTION

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 (g) to transjunctional voltage (V). Single-channel studies have shown that V 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 V 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 V-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 V falls on the gated hemichannel, whereas in fact V falls on both hemichannels and they can gate at the same time. Furthermore, there were several attempts to describe V 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).

Figure 1.

GJ channels contain fast and slow gating mechanisms. (A) Schematic of the GJ channel illustrating fast gates (orange) that close the channel partially to the so-called residual state and the slow gates (blue) that close the channel fully. (B) Schematic of the GJ channel is combined with the principal electrical scheme composed of four variable resistances in series attributed to the fast (γF) and slow (γS) gates. (C) A schematic of possible transitions of the gate, which has open (o) and closed (c) states.

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 V spike of approximately ±100 mV, causing relatively small decay of g, which can accumulate into significant g decays during bursting activity of neurons. Also, we demonstrate that very low g 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 (I) and g 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.

Table 1.

Set of parameters of MC16SM used for simulation of multicellular networks

Gate/Parameter	A	V0	γopen	γclosed	Ropen	Rclosed	Pt	Π
	mV−1	mV	pS	pS	mV	mV
Fast gate	0.15	40	24	3	10,000	10,000	0.00005	−1
Slow gate	0.15	40	24	0	10,000	-	0.00005	−1

The set of parameters of an MC16SM used for simulation of multicellular networks interconnected through electrical synapses formed of GJ channels resembling Cx36 and to assemble Figs. 3 B, 4 (E–G), 5, 6, 7, and 8. In search of minimum g necessary for synchronization and signal transfer of APs in 1-D and 2-D clusters of neurons (g and g; Fig. 5), we varied the number of GJ channels forming electrical synapses.

Table 2.

Set of parameters of S16SM and MC16SM used for simulation of multicellular networks

Gate/Parameter	A	V0	γopen	γclosed	Ropen	Rclosed	Pt	Π
	mV−1	mV	pS	pS	mV	mV
Fast gate	0.15	10	120	10	10,000	10,000	0.00005	−1
Slow gate	0.15	10	120	0	10,000	-	0.00005	−1

The set of parameters of an S16SM and an MC16SM used for simulation of multicellular networks interconnected through electrical synapses formed of GJ channels resembling Cx45 and to assemble Figs. 2, 3 (C and D), and 4 (B–D, F, and G).

Figure 3.

Simulation of cell–cell coupling through soma-somatic GJs. (A) The scheme illustrates two cells, each exhibiting an excitability described by H-H formalism and communicating through GJs operating according with an S16SM or an MC16SM. The left electrode was used to inject current (Iext) under current clamp mode and measure transmembrane voltage in cell 1 (V1). The right electrode was used to record transmembrane voltage in cell 2 (V2). (B) Electrotonic responses evoked in cell 1 (V1) and cell 2 (V2) by an Iext step of −4 pA. Simulations were performed using an MC16SM. (C) Simulation of changes in V1, V2, and g caused by gating of GJs formed of the Cx exhibiting properties similar to Cx45, in response to an Iext step of −18 pA. Simulation was performed using an MC16SM. The left inset in the V time plot shows that before reaching a steady state, V1 increases and V,2 decays, resulting in an increase of V and V gating. Phase differences between APs arising at the end of the Iext step caused a bipolar V spike (inset in V time plot) and, consequently, a small decay of g (inset in g time plot) caused by V gating. (D) Simulation as in A, but performed using an S16SM instead of an MC16SM and assuming that the junction contains 10 GJ channels.

Figure 4.

Changes of GJC in the pair of neurons exhibiting bursting activity. (A–D) Demonstration of changes in V, V, and g over time for GJs formed of hypothetical Cx45 in response to an Iext step of 15 pA. (B) Burst of APs evoked in cell 1 (red) and then transferred to cell 2 (blue). Inset shows that the delay between APs was ∼2 ms. (C) Transjunctional voltage was defined as V = V − V; inset shows bipolar V spikes arising because of a phase difference between APs. (D) Decay of g during the Iext step. (E) Decay of g of GJs formed of connexin exhibiting gating properties resembling those of Cx36. (F) Delay of AP transfer from cell 1 to cell 2 increased by ∼38% in cells expressing Cx45 and at a much smaller level for cells expressing Cx36. (G) Predicted reduction in g during Iext steps of 15 pA in the cell pair expressing various Cx isoforms. All simulations were performed using an MC16SM.

Figure 5.

Synchronization of neuronal clusters. (A and B) Minimum values of g necessary for synchronization of neuronal networks (g). (A) Repeated APs in cell 1 (red) and cell 2 (blue; b and c) evoked by long depolarizing current pulses shown in a. At g = 0, cells exhibited repeated spiking activity at 65 and 95 Hz (b) and g = 0.26 nS (c). (B) Summarized data of g values in 1-D and 2-D clusters of different sizes depending on differences in their firing rates of neurons, ΔFR. (C) Summary of the minimum g necessary for AP transfer between neurons (g) depending on frequency of APs caused by repeated 2-ms-long Iext steps of 30 pA. Values of g were evaluated in the cell pair (1 × 2) and 2-D clusters ranging from 2 × 2 to 7 × 7.

Figure 6.

Unidirectional signal transfer in a cluster of three cells through electrical synapses exhibiting I-V rectification. (A) Schematic of the cells with stimulation of cell 2. (B) g–V plot exhibiting rectification of GJs between cells 1 and 2 and no rectification of GJs between cells 2 and 3. (C) Stimulation of cell 2–initiated AP that spread to cell 2, but not to cell 1, because of reduced g1-2. (D) Stimulation of cell 1–initiated AP that increased g1-2 because of I-V rectification, therefore facilitating AP transfer to cells 2 and 3.

Figure 7.

Formation of reverberation of excitation in a cluster of nine neurons connected through GJs modulatable by an MC16SM. (B) Spread of APs from stimulated cell 1 through all cells shown in A. g values of GJs varied randomly from 0.175 to 0.2 nS and did not exhibit I-V rectification. (C) Initiation and termination of reverberation of excitation (b) by a single stimulus applied to cell 1 at 2 and 146 ms, respectively (a). Cell 5 is located at the core of reverberation and exhibits only electrotonic depolarizations not transforming into APs (c). The first stimulus caused AP, which spreads to cell 2 but not to cell 4 (d). Reverberation of excitation was stopped by the second pulse, which caused depolarization of ∼3 mV in the middle of cell 1 cycle (red arrow in e) and prevented reexcitation of cell 1 by cell 4. (f) Summarized data demonstrating time windows during which positive and negative Iext,cell 1 pulses of minimal amplitude were able to terminate reverberation of excitation.

Figure 8.

Reverberations of excitation in a 15 × 15 cluster of neurons caused by a single stimulus. In A and B, stimulus was applied to cell 1/1 and cell 8/8 (red circles), respectively. I-V–rectifying junctions are shown by light blue segments. Colored bar is calibrated for transmembrane potential.

Figure 2.

Comparison of I and g time plots simulated using an S16SM and an MC16SM. (A) V step of −60 mV applied to one cell of the cell pair. Both cells were voltage clamped. (B) I time plots obtained simulating a junction composed of 500 channels using an S16SM (red) and an MC16SM (blue). (C) g time plots obtained from the ratio I/V. g values evaluated using MC16SM (blue) reflect averaged changes of g values obtained using S16SM (red) and are virtually free of noise (inset).

H-H model

We have modified the original H-H model (Hodgkin and Huxley, 1952) by adding I, described by the following system of ordinary differential equations:

Here, I is the transmembrane current per unit area, I is external stimulus current, V is membrane potential, and C is the membrane capacitance per unit area; V, V, and V are the potassium, sodium, and leak reversal potentials, respectively; g, g, and g 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 α and β are rate constants for the respective ion channels.

The membrane capacity C was chosen as 1 µF/cm2. The corrected values of reversal potentials were chosen as V = −12, V = 115, and V = 10.6 mV, whereas conductances were chosen as g = 36, g = 120, and g = 0.3 mS/cm2 (Hodgkin and Huxley, 1952). The functions α and β 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.

S16SM

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 I and g 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 V. 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 p and p or probabilities to remain in the same state (p and p) are given byandwhere P 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, A characterizes sensitivity to voltage, V is voltage at which K, or K, is equal to 1, and Π is a gating polarity (+1 or −1). V 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 γ, γ, and γ rectify in accordance with the exponential function as proposed in Vogel and Weingart (1998): and where γ0 and γ0 are γ at V = 0, whereas γ0 is γ at V = 0. R, R, and R are rectification coefficients of the open/residual states and fast/slow gates, respectively.

Rectification of γ, γ, and γ does not allow for the evaluation of V and V 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 V and V without accounting for rectification:and

Then, we recalculate conductances of gates accounting for rectification and reevaluate a new set of V and V. We repeat the last cycle and at the nth iterative step evaluate conductance of the GJ channel (γ) as well as the difference, Δγ = γ − γ−1 for convergence testing. The value of γ is given by

The iterative procedure is stopped once the tolerance criterion Δγ < (γ−1/1000) is reached. The final values of V and V 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 V, i.e., V = V + V + V + V.

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 γ of the electrical synapse composed of m channels is evaluated, then g is estimated from

MC16SM

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 V 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 g. Thus, we used transitive rather than steady-state probabilities. In doing so, we used a recursive relation, p = p−1 · P, to calculate the transitive probability vector p at a discrete time moment, t, for each of 16 states, where P denotes the transition probability matrix of a discrete-time Markov chain. Each state s can be expressed as a 4-tuple, s = (s; s; s; s), 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 p and p 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 p, can be estimated as follows: p = (1 − p) · p · p · (1 − p).

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 γ. However, because voltage across each gate depends on the state s, i.e., V = V(s) and V = V(s), then γ must also be evaluated at each state, i.e., γ = γ(s). Once the transition probability matrix P is formed and the probability vector p is estimated, then the expected value of junctional conductance g is estimated as

Fig. 2 shows I (B) and g (C) traces generated in response to a V 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. Online supplemental material is available at http://www.jgp.org/cgi/content/full/jgp.201511488/DC1.

RESULTS

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 V that closes the channel from the open state with conductance γ to a substate or residual state with conductance γ and a slow gate that operates between an open state with conductance γ and a fully closed state with conductance γ = 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 g (g) dependence on V measured at the very beginning of applied V steps when V gating is not yet activated. Homotypic GJ channels, which exhibit oppositely oriented near-exponential relationships of γ or γ over V in aHCs, demonstrate a relatively small g dependence on V, 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 (V1) and cell 2 (V2). Fig. 3 B shows that at g = 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–g relationship simulated for soma-somatic GJs; the gray circle shows averaged K12 (∼0.15) and g (∼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 (R in Fig. S1 B and supplemental materials and methods) should be accounted for. Our preliminary data demonstrate that g 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 V-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 V1 and V2, resulting in a V of ∼30 mV; simulation was performed using MC16SM. The g trace demonstrates that g decayed, and it took ∼600 ms to reach a steady state at ∼50% of the initial g. With the decrease in g, V1 increased and V2 decreased (Fig. 3 C, left inset in V time plot), thus causing V to rise, which accelerated a decay in g. 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 V pulse (Fig. 3 C, inset on V time plot), resulting in a small decrease in g (arrow in the inset of g time plot). Assuming that GJs are formed of the less V-sensitive Cx36 (V and V were 40 mV instead of 10 mV for Cx45; Bukauskas et al., 2002b; González et al., 2007; Marandykina et al., 2013), then the g 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 g. It is notable that the unitary stepwise changes in g resulted in slow changes in V1 and V2 (Fig. 3 D, left inset in g time plot).

To assess the extent to which the V values that develop during spiking activity can affect g, 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, V = V = 10 mV, which roughly corresponds to Cx45. Fig. 4 C shows repeated V values arising mainly during the phase shift between APs in neighboring cells, which we term V. The amplitude of the V spike can reach approximately ±100 mV, which caused a ∼2% decrease in g (Fig. 4 D, inset on g trace). During bursting activity of neurons, small reductions in g accumulate, giving a ∼28% decrease in g at the steady state (Fig. 4 D). However, Fig. S1 A (see also supplemental materials and methods) shows that at higher g values (2 instead of 0.36 nS), the phase shift between APs in neighboring cells decreases, resulting in substantially shorter and smaller-amplitude V spikes. Consequently, this caused a ∼30% smaller drop in g than that shown in Fig. 4 D.

At V = V = 40 mV, which resembles Cx36 properties, a decrease in g during one V spike was ∼0.35%, which accumulated during bursting activity, resulting in a ∼2.5% decrease in g (Fig. 4 E). Fig. 4 F shows that a greater decrease of g 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 g depends on V and V, which roughly replicate gating of different Cxs at Iext = +15 pA (squares; spiking activity is present). Thus, g can decrease ∼30% in highly V-sensitive Cx45 or Cx57, whereas for less V-sensitive Cxs, e.g., Cx36 or Cx26, the decrease in g was much smaller. Even though electrical synapses are formed preferentially of Cx36 GJs exhibiting low V sensitivity under normal conditions, their sensitivity to V increases significantly when intracellular Mg2+ ([Mg2+]i) is elevated to 5 mM (Palacios-Prado et al., 2013), which would lead to larger reductions in g during bursting activity of neurons. Also, a robust increase in sensitivity to V was demonstrated for Cx45 and Cx57 during acidification (Palacios-Prado et al., 2009). Hence, the decrease in g 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 V spikes and changes in g 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 g were identical to those of the soma-somatic junction. Fig. S1 C shows that V spikes obtained for soma-somatic (red) and dendro-dendritic (green) junctions are similar in amplitude and duration, whereas both cause g 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, g values are <1 nS (∼0.13 nS in Vervaeke et al. [2010], 0.02–0.2 nS in Venance et al. [2004], and 0.13–1.0 nS in Devor and Yarom [2002]). Higher values of g were reported for soma-somatic junctions (mean ∼2.1 nS between low-threshold spiking thalamocortical neurons in Gibson et al. [1999] and ∼6 nS in clusters of mesencephalic trigeminal nucleus neurons in Curti et al. [2012]). Relatively small g 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 g required to achieve phase synchronization/lock in a cluster of two or more cells (g), each of which exhibits spontaneous activity but differs in firing rates. Fig. 5 A (b) shows that at g = 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 g 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 g 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 g approaches zero at small ΔFR and increases as ΔFR rises. It also demonstrates that averaged values of g tend to be higher for bigger clusters but approach steady state at the cluster size of ∼7 × 7. Also, g values of 1-D or cable-like (1 × n) clusters were at least ∼30% higher than those of 2-D clusters.

At g = 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 g 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 g necessary for one-to-one AP transfer (g) in the cell pair increases with rising frequency of Iext pulses. Also, the values of g 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 g 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 V gating affects synchronization of neuronal networks, we examined a network interconnected through (a) constant resistors, (b) GJs formed of low-V-sensitivity Cx36, and (c) GJs interconnected through high-V-sensitivity Cx45. Fig. S2 shows that the difference in g between Cx36 and the constant resistors is negligible (<2%). However, for Cx45, g was more than two times higher. We expect that a similar increase in g can occur for Cx36 under intracellular hypermagnesemia (Palacios-Prado et al., 2013, 2014), which increases V-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 g (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 g1-2–V plot (Fig. 6 B), whereas the junction between cells 2 and 3 was not rectifying (blue g2-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 V between cells 1 and 2 resulted in short-lasting V with an amplitude reaching ∼100 mV. The orientation of I-V rectification was such that at positivity on the cell 2 side, g1-2 decreased (red trace), causing a block of AP transfer to cell 1 (green V 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 V, g1-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 g 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 g, 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).

DISCUSSION

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 V 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, V can change dynamically because of a phase shift between APs in neighboring cells, which we name V. One of our major goals was to find whether Vs that last only a few milliseconds can influence GJC through V gating and how such influence depends on the Cx isoform. The spikes of Vs, despite their brevity, are high in amplitude (>100 mV) and potentially influence open probability of fast and/or slow gates because of V gating and single-channel conductance caused by I-V rectification. Indeed, data shown in Fig. 4 D demonstrate that for Cx45, a single V can cause ∼2% decay in g. During bursting activity of neurons, small g decays accumulate, leading to significant decays of g, as is shown in Fig. 4 D.

Although mammalian electrical synapses are preferentially composed of Cx36 GJs with low sensitivity to V under normal conditions, their sensitivity to V 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 V 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 g 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 g 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 g, 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 g 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 g instantaneously during short-lasting V 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 I-V 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 I-V relationship becomes linear, whereas transjunctional asymmetry in [Mg2+]i can cause an exponential I-V 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.