Network synchrony depends on gap junctions strength

To study the effect of gap junction plasticity, we developed a network of coupled inhibitory and excitatory neurons in the fluctuation-driven state (Fig 1A). The Izhikevich model was used for the inhibitory neuron population to fit the fast-spiking inhibitory neuron firing pattern [63]. Excitatory neurons are modelled by leaky integrate-and-fire models. As in [24], the excitatory neurons act as low pass-filters for their inputs while the FS neurons have a sub-threshold resonance in the gamma range [42–46]. To demonstrate this, we injected an oscillatory current of small amplitude in a single cell and recorded the amplitude response for different oscillatory frequencies. Excitatory neurons better respond to low frequency inputs, while FS neurons respond maximally for gamma inputs (Fig 1B). This is in line with the experimental evidence of Cardin et al. showing that FS-specific light stimulation amplifies gamma-frequencies [33].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Network synchrony depends on gap junction strength.

(A) The network consists of excitatory (E) and inhibitory (I) neurons. The neurons are coupled in an all-to-all fashion with chemical synapses. The inhibitory neurons are also connected by gap junctions (jagged green line). (B) Voltage response of one single excitatory (red line) / inhibitory (blue line) neuron to a sub-threshold oscillatory input current (see Methods). Excitatory neurons act as low-pass filters, whereas the inhibitory neurons show a resonance frequency in the gamma range. This resonance is in agreement with the network wide response observed by Cardin et al. 2009, when FS neurons are stimulated in the gamma range (black line, figure redrawn from [32] figure 3d). (C) Simulation of a pair of electrically coupled neurons N1 and N2, where N1 is voltage-clamped (red) such that it is hyperpolarized (light blue) and the potential of N2 is measured for different value of gap junction strength (γ = 3 and γ = 5). (D) Power of the main frequency component in the Fourier domain of the population activity (PA) of inhibitory neurons. The blue area denotes the lack of oscillations AI whereas the red area SR shows periodic oscillations in the spiking activity of inhibitory neurons. (E) Oscillation frequency of the network activity. The white area represents a region where the network is not oscillating and has no oscillation frequency. (F) Histogram of the oscillation frequency of population spiking activity. The values are contained in the γ range, from 30 to 60 Hz. (G) Ratio of bursting A_bursting over spiking A_spiking activity, averaged over 2 seconds. Bursting activity prevails in the light region and sparse firing dominates in the dark region. For the following Figures 1H and 1I, 100 ms of data is represented. (H) Raster plots of 100 FS neurons (blue) and 100 pyramidal neurons (red) for two values of the gap junction coupling, where dots represents spiking times and each line represents a neuron (note that the network E/I proportion is actually 80%/20%). Top raster plot shows asynchronous activity for low gap junction coupling and bottom raster plot shows synchronous activity in inhibitory and excitatory neuron populations, for strong gap junction coupling. (I) Membrane voltage traces of individual inhibitory neurons (dark blue) and population average (light blue, down-shifted) for different values of the gap junction coupling. Bursts appear for strong gap junction coupling on the peaks of the membrane voltage oscillations.

https://doi.org/10.1371/journal.pcbi.1006025.g001

All neurons have chemical synapses but only inhibitory neurons are also coupled via gap junctions (Fig 1A). The gap junctions are modelled such that a voltage hyperpolarization (depolarization) in one neuron induces a voltage hyperpolarization (depolarization) in the connected neuron. The current contribution of gap junction coupling is proportional to the difference of voltages between the coupled neurons, multiplied by the gap junction strength γ (Fig 1C). Moreover, when one neuron spikes, it emits a spikelet in the coupled neuron. We model this by a positive inhibitory to inhibitory electrical coupling, which we add on top of the negative inhibitory to inhibitory chemical coupling (see Materials and methods).

In order to understand the effects of gap junction plasticity, we initially considered the network without plasticity. We first explored the network behavior for different values of the mean gap junction strengths γ and mean external drive to the inhibitory neurons ν_I. As demonstrated by [24], our network exhibits two regimes (Fig 1D): an asynchronous irregular (AI) regime and a synchronous regular regime (SR). The AI regime occurs for networks with weak external drive and weak gap junctions. In this regime the network is in the fluctuation driven regime so that the neurons spike due to variations in their input. The SR regime occurs for strong external drive and strong gap junctions. This regime leads to the emergence of gamma oscillations. Mathematically, the network undergoes a Hopf bifurcation [24, 39]. The oscillations arise as the network directly inherits the resonance properties of the individual neurons. This is mediated through the gap junction coupling which effectively allows positive coupling through their spikelets. Moreover, the gap junctions reduce sub-threshold voltage differences between neurons which promotes synchrony. The excitatory neurons are not necessary for the oscillations but they amplify the dynamics (see [24] for mathematical derivations). When placed in the SR regime, the network oscillates in the gamma-range at a frequency near the single neuron resonance frequency (Fig 1E and 1F). In addition, we observe that the spiking activity is characteristic to the network regime, with bursting activity in the synchronous regime and spikes in the asynchronous regime (Fig 1G–1I).

To summarize, increased gap junction coupling and input drive into the network promotes gamma oscillations. To explain the relationship between network activity and gap junction plasticity, we first model the simplest case of plasticity between a pair of electrically coupled neurons. We then apply the plasticity rule to a population of neurons and investigate the effects on the network dynamics.

Model of gap junction plasticity: Bursting induces gLTD, spiking gLTP

To determine how gap junction plasticity can alter network dynamics, we developed a model of the plasticity based on experimental observations. [55] have shown that bursts in one or both neurons in an electrically coupled pair lead to long-term depression (gLTD). Therefore, we modeled gLTD as a decrease in the gap junction strength that is proportional to the amount of bursting. The constant of proportionality, α_gLTD serves as the learning rate. To infer α_gLTD, we reproduced the bursting protocol in Haas et al., where a neuron bursting for a few milliseconds, 600 times for 5 minutes, leads to 13% decrease (Fig 2A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Model of gap junction plasticity: Bursting induces gLTD, spiking gLTP.

(A) Bursting protocol replicated from Haas et al. [16]. A current (red line, top panel) of 300 pA for 50 ms at 2 Hz and of -80 pA otherwise is injected into a pair of coupled neurons induces repeated bursting (blue line, middle panel, voltage trace). To quantify the amount of bursting, we low-pass filtered (b_i) the voltage trace, threshold it at θ_burst = 1.3 (discontinued dark line), and integrate. Light blue areas represent the periods during which bursts are detected and therefore gap junctions are depressed. (B) When neurons N1 and N2 spike sparsely (top panel, dark blue, first part of the stimulus), gap junctions are potentiated (bottom panel, green line, first part of the simulation), whereas when they are bursting, gap junctions are depressed (second part of the simulation). (C) Green dots show steady-state values of the mean gap junction coupling for the gLTP with soft bounds, for different values of the network drive along the y-axis. For slow gLTP, the steady-state can be found in the AI regime, where the power of the oscillations of the population spiking activity is low (blue area). (D) Network architecture: A step excitatory drive is fed to the network of E and I neurons (same network detailed on Fig 1, with plastic gap junctions) inducing gamma oscillations. The activity of the network is read out by a downstream population of 200 regular spiking cells. (E) Top panel, step excitatory drive fed to the networks. Second panel, evolution of the mean gap junction coupling. As the excitatory drive is delivered, a gamma oscillation appears, leading to an increase in bursting activity which is followed by a depression of the gap junctions, until the new fixed point is reached. Bottom panels, raster plots of the inhibitory neurons (blue, I1), excitatory neurons (red, E1) and read-out neurons (red, RON). 6 s of data is represented. (F) Top panel, step excitatory drive. Other panels, population activity of the read-out neurons in red, evolution of the mean gap junction coupling in light blue. Second panel, simulation with plastic gap junctions. The read-out neurons are the most active during the transient oscillations. Third panel, static gap junction coupling. The read-out neurons are active as long as the excitatory drive is high. Bottom panel, no gap junction coupling. The read-out neurons are not active. 10 s of data is represented.

https://doi.org/10.1371/journal.pcbi.1006025.g002

Activity-dependent gap junction long-term potentiation (gLTP) has not been reported experimentally yet in the mammalian brain. There is evidence for activity dependent short-term potentiation in vertebrates [53, 64]. However, without potentiation, all gap junctions would likely become zero with time. To address this concern, we hypothesize that gap junctions can undergo gLTP and we modeled it such that single spikes induce gLTP by a constant amount given by the potentiation learning rate α_gLTP (Fig 2B, first half). Furthermore, we considered activity-independent gLTP rules in the supplementary materials (S1 Fig).

Gap junction plasticity regulates network-wide oscillations

Our plasticity model therefore potentiates gap junctions under spiking activity and depresses under bursting activity. Therefore, we wondered how gap junction plasticity can alter network dynamics. We previously quantified the amount of spiking versus bursting in our network for different levels of fixed gap junction strength and mean drive. For low levels of both, the network is spiking whereas for high levels of both the network is bursting. The spiking to bursting transition (Fig 1G) corresponds to the bifurcation (Fig 1D) from asynchronous irregular to synchronous oscillations at gamma frequency. When inhibitory neurons are oscillating, they fire a burst of spikes at the peak of the oscillations (Fig 1I, γ = 5). Therefore, when gap junctions are plastic, the network steady state can be found on the side of the bifurcation that balances the amount of potentiation due to spiking activity with the amount of depression due to bursting activity. The depression learning rate is inferred from Haas et al., while the potentiation learning rate is left as a free parameter.

We found that a strong relationship exists between gap junction plasticity and network synchrony. When the network is in the AI regime, characterized by low prevalence of bursting activity, gap junction potentiation dominates. However, for a strong mean coupling strength, the emergence of oscillations is associated by high bursting activity which leads to depression of the gap junctions. Therefore gap junction plasticity in our network maintains a tight balance between asynchronous and synchronous activity. Depending on the value of α_gLTP, the position of the plasticity fixed point lies either in the asynchronous regime (low α_gLTP, Fig 2C) or in the synchronous regime (high α_gLTP). For high values of α_gLTP, potentiation is fast while for low values, the potentiation is slow.

Gap junction plasticity allows for sparse but salient information transfer

We wondered how gap junction plasticity would interact with time-varying inputs. For the following experiment we consider slow gLTP. First, we let the network reach its steady state with a low level of drive (Fig 2E, beginning). As previously observed, the mean gap junction strength reaches a value which sets the network near the AI/SR transition. Then, we proceeded by injecting an additional constant current to the network. This new current baseline induces network level oscillations (Fig 2E, transition). However, over time the mean gap junction strength decays due to the gap junction plasticity mechanism. This gap junction depression is followed by a loss of synchrony and the network reaches its new steady state (Fig 2E, end), again near the border of asynchronous and synchronous regimes.

We measured the response of read-out neurons which receive projections from the excitatory and inhibitory neurons in our network (Fig 2D). At the onset of the current step, the network undergoes transient oscillations. When the gap junctions are plastic, the downstream neurons increase their spiking activity only for a few hundred milliseconds during the transient oscillations and then became almost quiescent again (Fig 2F, second panel). This contrasts with the simulation of a static network where the downstream keep a high firing rate (Fig 2F, third panel).

These results suggest that synchronous activity is a powerful signal to provoke spiking in downstream neurons. But oscillations and high firing rates of downstream neurons are also metabolically costly [65]. With transient oscillations however, the downstream neurons only sparsely fire when the stimulus changes but not when it is predictable. Thus, the regulation of oscillations mediated by gap junction plasticity allows for sparse but salient information transfer.

Gap junction plasticity enhances the ability of sub-populations of neurons to synchronize

We now sought to study the functional implications of fast gLTP. As stated before, this synchronizes the network into gamma oscillations. Synchronization between networks is considered to be one possible mechanism of information transfer [66–69]. We wondered whether gap junction coupling could mediate cross-network synchronization, and how gap junction plasticity would regulate this synchronization. To test this hypothesis, we considered two subnetworks having different oscillation frequencies and coupled by gap junctions (Fig 3A). A fast network oscillates at a gamma frequency and therefore is called the gamma-network. Then, a slow-network oscillates at a slower frequency as the membrane time constant of its inhibitory neurons is chosen to have a larger value. Indeed, previous analyses show that the network frequency in our model is inherited from the single neuron resonance frequency of inhibitory neurons [24, 70]. As a result, increasing the membrane time constant of the inhibitory neurons results in a decrease of the network oscillation frequency (Fig 3B–3D). Cross-network gap junctions reduce the frequency and phase difference between the gamma- and slow-network (Fig 3E and 3F) and larger differences of subnetwork resonant frequencies require a larger number of cross-network gap junctions for the networks to oscillate in harmony (Fig 3E and 3G). Their common frequency lies between the resonant frequencies of the decoupled networks. Importantly, cross-network synchronization requires the subnetworks to be in phase. If the gamma- and slow-network do not share enough gap junctions, there is little mutual information and no correlation in their population activities (Fig 3H and 3I), despite having a common oscillation frequency in some cases ([Δf_res = 0; number of shared GJs = 0] on Fig 3I). However, for small differences in the subnetworks resonant frequency Δf_res, increasing the number of shared gap junctions induces the oscillations to lock together. The networks oscillate in phase (Fig 3F, end of first row) as reflected in their mutual information (Fig 3H, dark blue area) and their correlation (Fig 3I, dark red area). In summary, two networks in the SR regime with different resonance frequencies and/or out-of-phase can synchronize if they are coupled by gap junctions. Furthermore, a large number of shared gap junctions is required for large differences of resonant frequency.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Subnetworks having different frequency preferences can synchronize their activity if they share gap junctions.

(A) Both subnetworks have the same topology with all-to-all connected inhibitory and excitatory neurons. Inhibitory neurons have static gap junctions (GJs). The Gamma Network (GN) is connected to the Slow Network (SN) with a varying number of gap junctions. The time constant of the SN inhibitory neuron membranes is varied. (B) Frequency-transfer characteristics of one single inhibitory neuron to a sub-threshold oscillatory input current (see Methods) for different values of its membrane time constant τ_v. The sub-threshold resonance frequency decreases as τ_v increases. Data of Cardin et al. 2009 is also represented (black line, figure redrawn from [32] figure 3d). (C) Changing the single neuron sub-threshold resonance modifies the network oscillation frequency. Mean inhibitory membrane potential for τ_v = 17 ms (continuous line) and τ_v = 55 ms (dashed line). 100 ms of data is represented. (D) Relationship between single neuron resonance (black line) and network oscillation frequency (gray line). For the following figures E and F, for the tuples (Δf_res; Number of shared GJs), the upper (lower) triangle represents the value in the SN (GN). For panels E, F, H, I, the x-axis represents the number of cross-network gap junctions between the GN and SN. The y-axis represents the difference of resonance frequency between the GN and the SN. (E) Oscillation frequencies. We observe that the GN and the SN adopt the same oscillation frequency for low Δf_res and high number of shared gap junctions. (F) Phase differences between population activities of the GN and the SN, when they share the same frequency. Lighter squares denote parameters for which the phase difference is lower. The GN and the SN are considered in phase when the phase difference is zero. Dark blue squares describe a region that is excluded because the GN and the SN do not oscillate at the same frequency, therefore cannot be in phase. (G) Raster plots, where dots represent spiking times and each line represent a neuron, for small (first column) and large (second column) differences in Δf_res. For all raster plots, from top to bottom are represented excitatory and inhibitory neurons from the SN, then inhibitory and excitatory neurons from the GN. 100 neurons are shown for each population. When no gap junctions are shared (bottom row), both networks do not synchronize and are out-of-phase. With 40 shared gap junctions (top row), the networks synchronize and are in phase for small values of Δf_res. 100 ms of data is represented. (H) Mutual information between the PAs of the GN and SN. The increase in mutual information for the top row, where Δf = 21Hz, can be due to the fact the SN oscillates at half the frequency of the GN (which oscillates around 40Hz). (I) Pearson’s correlation of the PAs of the GN and SN. Comparing with panel H, there is high correlation when the GN and the SN are in phase.

https://doi.org/10.1371/journal.pcbi.1006025.g003

As gap junctions can synchronize two oscillating populations of neurons, we wondered whether the same synchronization would occur with one population in the AI regime. First, we initialized the gamma-network in the AI regime while the slow-network was initialized in the SR regime (Fig 4A). After coupling the gamma- and slow-network together, we found that, while the oscillation frequency of the gamma- and slow-network matched (Fig 4B), the two networks could not synchronize. The networks were always out-of-phase with very weak correlation between the population activities (Fig 4C and 4D). The results were similar if the gamma- and the slow-network were initialized in the reverse synchronous and asynchronous parameter regimes, respectively (not shown). Cross-network synchronization is not robust when one network is not oscillatory.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Gap junction plasticity lets networks recover synchronization.

For all panels, the x-axis represents the number of cross-network gap junctions between GN and SN. The y-axis represents the difference of resonance frequency between the GN and the SN. The gap junctions are static from panels A to D and plastic from panels E to H. Values for the Gamma Network (resp. Slow Network) are represented by the lower (upper) triangles. The GN (SN) has weak (strong) initial mean GJ coupling. Shared GJs are initialized with mean coupling strength in the middle between those of the GN and the SN. (A) Oscillation power. The GN, with weak GJ coupling, shows weak or no oscillations. (B) Oscillation frequency. We observe that the GN and the SN oscillate at the same frequency only for high number of shared GJs. (C) Phase differences between PAs of the GN and the SN (as for Fig 3H). The GN and the SN stay mostly out-of-phase. (D) Correlation of the PAs of the GN and the SN. Except for the particular case where Δf_res = 0 and the number of shared GJs is high, the PAs of the GN and the SN show no correlation. (E) Oscillation power. Comparing with panel A, we observe that the oscillation power seems to match in both networks, with mostly the oscillation power of the GN (initially weak) increasing to the SN’s levels (initially strong). (F) Oscillation frequency. Comparing with panel B, we observe an extension of the region where the GN and the SN oscillate at the same frequency. (G) Phase differences between PAs of the GN and the SN. We observe here a large region where the GN and SN are in-phase. (H) Correlation of the PAs of the GN and the SN. Comparing with panel D, we observe a large extension of the region where both networks are synchronized.

https://doi.org/10.1371/journal.pcbi.1006025.g004

Given these constraints on cross-network synchronization, we wondered if gap junction plasticity could remedy the situation and allow for robust cross-network synchronization. To test this hypothesis, we repeated the simulation protocols with the gamma- and slow-network initialized in the asynchronous and synchronous regimes (respectively) and with plastic gap junctions. Here we considered the case where the gLTP rates were slow. As shown previously, gap junction plasticity regulates oscillations such that the network in the asynchronous irregular regime transitions to the oscillatory regime (Fig 4E). The oscillation frequencies of these two networks match (Fig 4F). Strikingly, even with a large resonant frequency difference, the gamma- and slow-network now synchronize through a small number of shared gap junctions (Fig 4G and 4H). This indicates that gap junction plasticity allows for cross-network synchronization that is robust to the underlying neuronal parameters for small numbers of shared gap junctions.

Gap junction plasticity allows for robust information transfer

We hypothesized that cross-network synchronization mediated by plasticity allows information transfer. To investigate this, we considered a similar network architecture as previously studied, with two networks, an input-network and an output-network. The input-network receives an input projected by random weights to its neurons. The output-network is connected to the input-network with a small number of gap junctions and inhibitory chemical synapses.

First, to demonstrate the information transfer capability of the network, we consider static gap junctions with oscillatory inputs to the input-network. The stimulus information is transmitted to the output-network via the frequency modulation of the synchronized oscillations and not by spike transmission nor amplitude modulation (Fig 5A–5D). When sharing gap junctions, the input- and output-network synchronize together (Fig 5A) and their spiking activity is locked (Fig 5B). As the amplitude of the input signal increases, the spiking activity increases in the input-network but not in the output-network (Fig 5C). For a network in the SR, there is a positive correlation between the signal amplitude and the network oscillation frequency (Figs 1E and 5D). This frequency modulation is transferred from the input- to the output-network. Thus, the input amplitude can be estimated from the oscillation frequency of the output-network, despite the absence of chemical synapses between the input-network and the output-network (Fig 5E). However, this synchrony code is only possible for signals below a certain frequency (Fig 5F and 5G). Indeed, the instantaneous oscillation frequency is estimated by measuring the period between consecutive peaks of the population activity. For example, oscillations at 50 Hz have a period of 20 ms. Variations happening within those 20 ms are compressed to a single period value and thus are not transferred via frequency modulation. Mechanisms for estimating the input value from the oscillation frequency of the output-network are discussed further in the methods section. Finally, we tested if this synchrony code was valid for non-oscillatory signals (Fig 5H). We found that non-oscillatory, slowly varying random signals could also be robustly transmitted from the input- to the output-network with gap junction coupling (Fig 5I).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Gap junction coupling allows networks to transmit information and gap junction plasticity improves robustness of the transfer.

(A Voltages traces of inhibitory neurons in the input-network (IN) in light blue and in the output-network (ON) in purple, when networks share no GJs (first rows) or 40 GJs (bottom rows). Despite not directly receiving the input signal, the ON synchronizes its activity with the IN. For panels B to I, the networks share 40 GJs. 50 ms of data is represented. For the following figures 5B, 5C and 5H, 1 s of data is represented. (B) Input signal in red, number of spiking events of inhibitory neurons of the IN in light blue and of the ON in purple, for time bins of 0.1 ms. (C) Input signal in red, number of spiking events of inhibitory neurons of the IN in light blue and the ON in purple, for time bins of 25 ms. (D) Input signal amplitude A_i as function of the corresponding PA peak interval T_i for input signals oscillation at 4 Hz with mean varying from 0 to 1000 (See Methods). (E) Input signal in red and decoded input signal in purple. The PA peak interval T_i is used to estimate the input amplitude. (F) Correlation between input signal and decoded input signal. The amplitude of the input is 400 pA, its frequency goes from 0 to 100 Hz. (G) Correlation between input signal and decoded input signal. The amplitude of the input goes from 0 to 10000 pA, its frequency goes from 0 to 100 Hz. (H) Example of 1 s of colored noise input signal (A = 800 pA, mean = 400 pA, τ_filter = 100 ms) in red and decoded input in purple (correlation 0.8). (I) Pearson’s correlation coefficient between input and decoded input for static (plastic) network in black (gray) for different values of the mean initial GJ coupling strength, as function of the number of shared GJs. The simulation is repeated for 10 different inputs. (J) Pearson’s correlation coefficient between input and decoded input for static (resp. plastic) network in black (resp. gray) as function of the proportion of GJs removed. The simulation is repeated for 10 different inputs. (K) Mean gap junction change between the steady-state value obtained with all the gap junctions, and the steady-state value obtained after gap junction removal. The remaining gap junctions compensate for the missing ones as they become stronger in strength.

https://doi.org/10.1371/journal.pcbi.1006025.g005

As gap junction plasticity can regulate oscillations, we tested whether the plasticity can make this synchrony code robust to parameter variations or potential gap junction loss. First, as previously shown, gap junction plasticity enhances the ability of networks to synchronize. If initialized in the AI regime and with static gap junctions, there is no information transfer via frequency modulation (Fig 5J, left panel). However, with plasticity and fast gLTP, the oscillations are regulated and the network synchrony is recovered which results in successful information transfer (Fig 5J left panel). A critical amount of oscillation power and a critical number of shared gap junctions are required for information transfer, after which increasing each of them does not yield significant improvement (Fig 5J). Furthermore, we studied whether gap junction plasticity could restore information transfer if gap junctions were deleted. While there is loss in the quality of the transfer when static gap junctions are removed, plastic gap junctions maintain the quality of the transfer by increasing the strength of the remaining gap junctions. This mechanism compensates for the missing gap junctions (Fig 5J and 5K).

To summarize, gap junction plasticity expands the necessary conditions for information transfer. It regulates oscillations, and by promoting phase-locking of oscillations, it contributes to the propagation of information to downstream networks. Finally, if some gap junctions are failing, due to protein turnover perhaps, the remaining ones can increase their strength through plasticity. This helps to maintain accurate information transfer.

Gap junction coupling between interneurons affects network synchrony

Despite being less common than chemical synapses, gap junctions are ubiquitous in the central nervous system. Example includes the inferior olivary nucleus [71–73], the thalamic reticular nucleus [74, 75], the hippocampus [36, 76], the retina [52, 77], the olfactory bulb [78], the locus coeruleus [79], or also the neocortex [80, 81]. Moreover, they drastically alter the firing activity of their connecting neurons [82, 83], as well as the network dynamics [20–24]. Furthermore, gap junctions between inhibitory interneurons are reported in many cortical regions where global oscillations of neural activity are observed [21, 27, 84, 85]. These inhibitory neurons exhibit sub-threshold resonance that amplifies a specific frequency range [33]. Therefore, gap junction induced synchrony and inhibitory neurons frequency preference are a possible substrate for global oscillations in these cortical regions. Our work is consistent with recent results showing that together gap junction strength and sub-threshold resonance of inhibitory neuron promote oscillations of neuronal activity [24, 70].

Previous models of gap junction plasticity

There has been a recent interest in modelling gap junction plasticity. Snipas et al. [86] developed of model of gap junction coupling that would exhibit short-term plasticity. By combining a 36-state model of gap junction channel gating with Hodgkin-Huxley equations [87], they show that gap junction channel gating, induced by bursting activity, could lead to short term depression. In future work, it would be interesting to combine this model of gap junction short-term plasticity with our model. Chakravartula et al. [88] introduced a new type of adaptive diffusive coupling in a network of Hindmarsh-Rose neurons [89, 90]. They assumed that connections between pairs of neurons would follow a Hebb’s law [91], where neurons with simultaneous activity would strengthen their connection, while others with dissimilar activity would weaken their coupling. They observe the emergence of locally synchronized groups of neurons, whose synchronization could be transient or permanent. Their results are consistent with ours showing synchronization of subnetworks coupled with gap junctions.

Model of gap junction plasticity: Bursts induce gLTD, spikes induce gLTP

Recently, Haas et al. [55] reported the first experimental evidence of activity-dependent gLTD of gap junctions of interneurons in the thalamic reticular nucleus, even though the mechanism remains to be investigated [62]. Also Sevetson et al. [56] found that calcium-regulated mechanisms support gap junction gLTD in the thalamic reticular nucleus. The mechanisms are similar to those observed for the plasticity of chemical synapses. We designed a rule for activity-dependent gLTD consistent with those results. We assumed that a cortical fast-spiking interneuron would exhibit the same plasticity properties as a thalamic reticular neuron because gap junctions are mostly made from the connexin Cx36 throughout the central nervous system [74, 92]. To our knowledge, there is no study yet on activity-dependent gLTP of gap junctions. However recent studies suggest that gLTD and gLTP share a common pathway [48, 56]. Therefore, we propose a rule for activity dependent gLTP, assuming that low frequency spiking activity leads to gap junction potentiation. However, our results do not depend on the exact formulation of gLTP. As we have shown, an activity-independent rule yields similar behavior (supplementary material, S1 Fig). Moreover, we did not observe significant changes by modelling asymmetrical gap junctions (supplementary material, S2 and S3 Figs).

Gap junction plasticity regulates oscillations and propagates transient information

Our model demonstrates that the regulation of oscillations is mediated by gap junction plasticity. Fast potentiation leads to bursting activity while slow potentiation leads to asynchronous irregular activity. Our first hypothesis assumed that the potentiation is slow and the network is in the AI regime. Thus, at the steady-state, gamma power is weak or non-existent. Evidence from Tallon-Baudry et al. and Ray et al. [93, 94] is consistent with our results. When no stimulus is provided or task required, electroencephalogram recordings show that power in the gamma-band is weak. After the onset of a sensory stimulus, gamma oscillations can be detected in cortical areas. This has been reported for example with visual stimuli triggering gamma oscillations in the mouse visual cortex [95]. In our model, the neurons oscillate transiently when receiving a constant external stimulation. This mechanism operates by crossing the bifurcation boundary between the AI and SR regime. However, over time the mean gap junction strength decays due to the additional bursting activity. The gap junction depression leads to a loss of synchrony and the network returns to the AI regime. Therefore we predict a loss in gamma power for sustained stimulus. A similar mechanism may be involved in the reduction of gamma oscillation induced by slow smooth movements [96, 97].

We wondered what could be the functional role of this transient oscillatory regime. Projecting the excitatory activity of our network model to downstream neurons revealed that they fire sparsely, for a short duration after stimulus onset, and are quiescent otherwise. Thus, gap junction plasticity could efficiently encode the change in incoming stimuli. This could allow for energy conservation as oscillations are energetically expensive [65]. Moreover, [98] show that cortical circuits near the onset of oscillations could promote flexible information routing by transient synchrony.

Plastic gap junction coupling for robust information routing

The role of gamma oscillations is highly debated [94]. They could play no role and simply be a marker of the excitation-inhibition interaction. However others studies suggest they could be involved in information transfer. It is thought that retinal oscillations carry information to the visual cortex [99]. Moreover they could serve as inter-area communication by promoting coherence in neural assemblies which would align their windows of excitation. This would allow for effective spike transmission [68, 94, 100]. Furthermore, Roberts et al. [101] observed high gamma coherence between layers 1 and 2 of macaque’s visual cortex by dynamic frequency matching. Here, we demonstrate one potential mechanism for information transmission through gamma oscillations. Our networks make use of gamma frequency modulation to transmit information in a robust manner, similar to the principle used for FM radio broadcasting. The amplitude of the input signal modulates the oscillation frequency, which increases almost linearly with the amplitude. Our model demonstrates that gap junction plasticity robustly mediates network oscillations and cross-network synchronization. If some gap junctions are removed, the remaining gap junctions become stronger and compensate for the missing ones. Thus, gap junction plasticity insures the phase-locking of the coupled network and it allows for information routing. In particular, there is evidence suggesting that gap junctions could promote long-distance signaling by implementing frequency modulation of calcium waves in astrocytes [102]. Moreover, correlation was found during gamma activity between amplitude and frequency modulation of local field potential of CA3 pyramidal neurons of anesthetized rats [103]. In addition, our network models could also represent the subnetworks of the TRN, with each connected to a separate excitatory neuron of thalamus [104]. However, TRN inhibitory neurons exhibit longer bursts than those of cortical fast-spiking neurons, due to long lasting T-current (about 50ms) and further work is necessary to make predictions on this brain region behaviour [105].

Failure to regulate oscillations, could be the origin of several cognitive pathologies. Disruption of brain synchrony in the inferior olive is thought to contribute to autism due to the loss of coherence in brain rhythms [106]. Excess of high frequency network wide oscillations in the cortex have been observed to also correlate with autism in young boys [12]. The inferior olive differs for its density of gap junction being the highest in the adult brain [71, 72]. It may be involved in the generation of tremors in Parkinson’s disease, however the severity of induced tremors in Cx36 knockout mice remained the same as in wild-type mice [107, 108]. This could be due to gap junctions made from other connexins (such as Cx43) taking over for the knocked-out ones.

Recent studies highlight the critical role of gap junctions and their plasticity in efficient cognitive processing [109]. As experimental and computational techniques improve, new efforts can further unveil their properties and expand our understanding of cortical functions. Our computational model shows that gap junction activity-dependent plasticity may play an important role in network-wide synchrony regulation.

Neuron model

We model Fast Spiking (FS) interneurons with Izhikevich type neuron models [63]. This model offers the advantage to reproduce different firing patterns as well as a low computational cost [112]. The voltage v follows (1) (2) combined with the spiking conditions, (3) where τ_v is the membrane time constant, v_ra is the membrane resting potential, v_rb is the membrane threshold potential, k_u is the coupling parameter to the adaptation variable u, R is the resistance and I is the current. The adaptation variable u represents a membrane recovery variable, accounting for the activation of K⁺ ionic currents and inactivation of Na⁺ ionic currents. It increases by a discrete amount b every time the neuron is spiking and its membrane potential crosses the threshold v_threshFS. It provides a negative feedback to the voltage v. The recovery time constant is τ_u, a is a coupling parameter, v_resetFS, and v_rc are voltage constants and b is a current constant.

For the FS neurons, we chose the membrane potential reset v_resetFS and the spike-triggered adaptation variable b to account for the onset bursting activity observed in vivo. Modifying k_u, v_ra, v_rb and v_rc was sufficient to observe the emergence of a resonance frequency. We set the time constant τ_u to obtain a resonance frequency of 45 Hz, which is in the same range as observed in vivo by [33] (Fig 1B). To measure the sub-threshold resonant property (Figs 1B and 3B and 3D), we recorded the amplitude of the neuronal membrane potential V_E in response to different oscillation frequencies f of low level sinusoidal currents I(t) = I₀ cos(2πft) (with I₀ = 0.01 pA). We then normalized the amplitude response as follow (4) for frequencies between 0 and 1 kHz. The || || denotes the maximum absolute value observed over time.

To model regular spiking excitatory neurons, we chose a leaky integrate-and-fire model, (5) where τ_m is the membrane time constant, v the membrane potential, I the current and R_m the resistance. Spikes are characterized by a firing time t_f which corresponds to the time when v reaches the threshold v_threshRS. Immediately after a spike, the potential is reset to the reset potential v_resetRS.

Network

In the single network model (Figs 1 and 2), each neuron is connected to all others by chemical synapses, but in addition, inhibitory neurons are connected via electrical synapses to all other inhibitory neurons, as in [24]. Thus, the current each individual neuron i receives can be decomposed in four components (6) where is the current coming from the transmission of a spike via electrical (i.e. spikelet) and chemical synapses, is the sub-threshold current from electrical synapses (for inhibitory neurons only), is the noisy background current and characterizes the external current. The current due to spiking on excitatory neurons is given by (7)

The current into inhibitory neurons are (8) where W^αβ is the coupling strength from population α to population β with {α, β} = {E, I}. Finally, is the inhibitory to inhibitory coupling between neuron i and j, consisting of the chemical synaptic strength W^II,c and the electrical coupling for supra-threshold current, also called the spikelet. There is no experimental data yet on the change of the spikelet as function of the strength of the gap junctions. We hypothesize that the contribution of the spikelet is proportional to the gap junction coupling , where γ_ij is the gap junction coupling between neurons i and j. This spikelet term is necessary due to the fact that our neuron model does not explicitly have a spike kernel in the voltage dynamics [24]. Note that W^EE, W^EI, W^IE, W^II,c are identical among neurons, but varies as the spikelet contribution depends on the coupling strengths γ_ij, which can be plastic. We also modeled the network with chemical weights following a log-normal distribution, which yielded similar results (data not shown).

We represent the post-synaptic potential response to a chemical or electrical spike with an exponential of the form for t > t_jk. The excitatory and inhibitory synaptic time constants are τ_E and τ_I respectively and t_jk represents the k^th firing time of neuron j.

In between spikes, for every pair of inhibitory neurons i, j, the gap junction mediated sub-threshold current is characterized by (9) where γ_ij is the gap junction coupling between inhibitory neurons i and j of respective membrane potential V_i and V_j. In our model, we suppose that gap junctions are symmetric with γ_ij = γ_ji. Gap junctions are initialized following a log-normal distribution with the location parameter μ_gap = 1 + ln(γ/N_I) and the scale parameter σ_gap = 1.

Neurons also receive the current I_noise which is a colored Gaussian noise with mean ν, standard deviation σ and τ_noise the time constant of the low-pass filtering (10) and (11) with ξ is drawn from a Gaussian distribution with unit standard deviation and zero mean.

Plasticity model of gap junctions

Our plasticity model is decomposed into a depression γ⁻ and a potentiation term γ⁺.

gLTD: Depression of the electrical synapses for high frequency activity

Haas et al. [55] showed that bursting activity of both neurons or one of the two neurons leads to long-term depression (gLTD) of the electrical synapses. To capture this effect in our model, we first defined a variable b_i which is a low-pass filter of the spikes of neuron i (12) where δ is the Dirac function and τ_b = 8 ms is the time constant. When b_i reaches a value of θ_burst = 1.3, this indicates that two or more spikes happened within a short time interval. Therefore, the burstiness of neuron i is characterized by H(b_i − θ_burst) where H is the Heaviside function that returns 1 for positive arguments and 0 otherwise.

In our simplified model, we consider that the individual electrical coupling coefficient γ between neurons are non-directional. Every time the interneurons burst, the gap junctions undergo depression, (13) where α_gLTD is the depression learning rate.

We fit α_gLTD to the data by implementing the stimulation protocol used in [55]. We applied a constant current injection of 300 pA for 50 ms every 0.5 s (2 Hz) and of -80 pA the rest of the time, to maintain the membrane potential at -70 mV. This protocol lasts for 5 minutes. We estimate α_gLTD = 15.69 nS ⋅ ms⁻¹ by such that it leads to a depression of 13% of the gap junction strength at the end of the stimulation protocol, as reported by Haas et al.

gLTP: Potentiation of the electrical synapses for low frequency activity

If gap junctions were only depressed, they would decay to zero after some time. Therefore, there is a need for gap junction potentiation. However, no activity dependent mechanisms was reported yet in the experimental literature, but several studies suggest that the calcium-regulated mechanisms leading to long-term depression could be involved in potentiation as well [48, 53, 56, 113].

We consider two gLTP rules. The first has a soft bound, i.e. the magnitude of modification is proportional to the difference between the gap junction value and a baseline coupling strength γ_b (14) where α_gLTP is the potentiation learning rate and sp_i(t) = ∑_tik<t δ(t − t_ik).

The second gLTP rule we consider has no maximum bounds (15)

Moreover, to show that our results do not depend on the specific gLTP rule, we also consider a different gLTP rule where the update is passive and therefore does not depend on neural activity. This alternative rule yields similar results (supplementary material, S1 Fig).

Coupling coefficient

The coupling coefficient is the ratio of voltage deflections when a step current was injected to one neuron of a coupled pair, which were maintained at a baseline voltage of -69 mV. During current injection, the injected neuron is hyperpolarized at -75 mV (16) when 1 is the index of the injected neuron. The gap junction conductance used for measuring the coupling coefficient was obtained from the mean value of the gap junction coupling at steady-state. The coupling coefficient is about 5% for a network of 200 inhibitory neurons. Please note that the gap junction conductance and the coupling coefficient scale inversely to the network size in our model. We chose to use the mean value, as there is very little variance (4 orders of magnitude lower that the mean value) in the gap junction coupling strength at steady state. For reference, the coupling coefficient was measure around 12%±8% averaged for 313 pairs in the TRN [55]. Moreover, for adult rats, for 91 paired recordings of adjacent IO neurons, the coupling coefficient varies from 1% to 8% [114], and for 14 pairs of fast-spiking cortical neurons, the coupling coefficient was around 1.5% [115].

Quantification of network spiking activity

To estimate the plasticity direction for different value of external input ν and gap junction strength γ, we observe the activity of the network (without plasticity) in a steady state over a duration T = 6 s. For a chosen tuple (ν;γ), we average over time and over neurons the bursting and spiking activity (17) and (18) Then, we explore the values of the ratio of bursting over spiking activity (19) as function of the coupling coefficient γ and of the mean external input ν over the parameter space .

Quantification of oscillation power and frequency

To quantify the frequency and the power of the oscillations in the neuronal activity, we perform a Fourier analysis of the population activity r which we define as the sum of neuron spikes within a population, during the time step dt (20)

We compute a Discrete Time Fourier Transform (DFT) and extract the power and the frequency of the most represented frequency in the Fourier domain. The formula defining the DFT is (21) where the r_n sequence represents N uniformly spaced time-samples of the population activities. We measure the amplitude of the Fourier components for k = 1..N/2 (because the Fourier signal is symmetric from N/2 to N). We identify the maximal one, its associated frequency and its power .

Downstream read-out neurons

To simulate the projection of a cortical layer onto another layer, we model downstream read-out neurons with the same regular spiking neuron model as the first cortical layer. The input I_j received by each downstream neuron is the projected activity of all excitatory and inhibitory neurons of the first cortical layer, multiplied by the coefficients W^ERON and W^IRON respectively: (22)

When delivering the step current I_step to the network (Fig 2D), the time at which the neurons receive I_step follows a normal distribution centered on the transition time, with variance 10 ms. This variability avoids the confound of transient and unstable synchronization of the network due to a strong input delivered to all neurons simultaneously.

Cross-network synchronization

We investigate the role of gap junction coupling and its plasticity in synchronizing networks having different oscillation frequency preferences. We design a network consisting of two subnetworks having the same topology as described in Network: Each subnetworks has 800 excitatory neurons and 200 inhibitory neurons. There are all-to-all chemical synapses within each subnetworks (their strengths are reported in Table 1). There are no cross-network chemical synapses. The intra-network gap junctions are all-to-all. In addition, we vary the number of sparse cross-network gap junctions from 0 to 40. The gap junction strengths are initialized following a log-normal distribution as described in Network. We take γ = 3 which yields AI behavior in the network and we take γ = 5.5 which yields bursting behavior in the SR regime.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1.

https://doi.org/10.1371/journal.pcbi.1006025.t001

One of the subnetworks is called the Slow Network (SN) and we change the value of the membrane time constant of its inhibitory neurons τ_v from 17 ms to 55 ms. This decreases the neuron sub-threshold resonance frequency, which also lowers the frequency of the subnetwork oscillation when it is in the synchronous regime. The second subnetwork has its neuron membrane time constant fixed at 17 ms and is called the gamma-network because it oscillates at gamma frequency. The simulations last 10 seconds, which is long enough for the gap junction coupling to reach its steady state when the gap junctions are plastic.

To quantify the similarity between population activities from both subnetworks, we evaluate the Pearson’s correlation coefficient between their population activities r_GN and r_SN from the gamma- and slow-network respectively. The firing rates, r_GN and r_SN are defined as in Eq (20).

We measure the mutual information between the mean currents from SN and GN with (23) where p(x, y) is the joint probability function of X and Y, and p(x) and p(y) are the marginal probability distribution functions of X and Y respectively. Time bins of 10 ms are used to estimate the probability functions.

For each subnetwork, we evaluate the frequency and power of their oscillations as described in the section Quantification of oscillation power and frequency. When the difference of oscillation frequency between both subnetworks is less than 1 Hz, we measure the cross-correlation of their population activities r_GN and r_SN (24) The phase difference is measured as the time delay relative to the oscillation period (25) where ⋆ is the convolution operator and T_period is the oscillation period.

Information routing

We investigate whether gap junction coupling and its plasticity play a role in routing information between networks. We consider the same system as described in the previous section, with two subnetworks coupled with gap junctions, except here all the inhibitory neurons have the same membrane time constant τ_v = 17 ms (e.g. corresponding to resonance frequency at gamma). The first network, called the Input Network (IN) receives an input projected to its N_IN neurons (N_IN = 1000) by N_IN weights drawn from a uniform distribution between 0 and 1. The second network is called the Output Network (ON, N_ON = 1000).

To examine if there is successful transfer of information between both networks, we attempt to reconstruct the input signal from the ON’s population activity r_ON. First, we obtain the low-pass filtered population activity of ON, r_filt, with (26) with τ_r = 3 ms. Then we detect the rising and falling times of the filtered population activity by detecting when it crosses a threshold θ_r = 2. This gives us rising times , when it crosses the threshold from below, and falling times, when it crosses the threshold from above. We obtain the peak intervals T_k by measuring the time difference between consecutive rising times .

For Fig 5D, we plot , the mean values of the input signal x between the rising times and as function of their corresponding peak intervals T_k (27)

We reconstruct the network input (Fig 5E and 5H) by doing a linear interpolation of the inverse of those peak intervals T_k, so that the input signal and reconstructed input have the same length. (28) Finally to estimate the quality of the reconstruction, we measure the Pearson’s correlation coefficient (which is invariant by affine transformation) between the input and the reconstructed input.

In order to test the robustness of the system we measure the quality of the reconstruction for an oscillatory input signal of which we vary the frequency f (Fig 5F) and amplitude A (Fig 5G). (29)

Then we measure the routing of random signals x(t) = ν_IN + σ_IN ⋅ η_IN(t), where ν_IN is the signal mean, σ_IN is the signal standard deviation, η_IN is an Ornstein Uhlenbeck fluctuation with correlation time τ_x = 100 ms and unit variance. We build a dataset of 10 input signals and then we measure the Pearson’s correlation coefficients between the input x(t) and the reconstructed input for those 10 inputs respectively. For Fig 5I, we scale the log-normal distribution of the gap junction strength (see Network) with γ = 3 to set the network in the asynchronous, with γ = 5.5 to set the gap junction near their plasticity fix point, and γ = 8 for a regime with strong oscillations.

To study the robustness of the information routing to gap junction deletion, we randomly delete an increasing number of gap junctions and measure the evolution of the Pearson’s correlation between x and . We also measure the change in the mean gap junction coupling, if there is plasticity, between the initialization (with γ = 5.5) and the steady-state (after 6 s of simulation).

All parameters are listed in Table 1 unless otherwise specified in a figure.

Parameters

We list in Table 1 the parameters used for our simulations.