Gap junction plasticity as a mechanism to regulate network-wide oscillations

Cortical oscillations are thought to be involved in many cognitive functions and processes. Several mechanisms have been proposed to regulate oscillations. One prominent but understudied mechanism is gap junction coupling. Gap junctions are ubiquitous in cortex between GABAergic interneurons. Moreover, recent experiments indicate their strength can be modified in an activity-dependent manner, similar to chemical synapses. We hypothesized that activity-dependent gap junction plasticity acts as a mechanism to regulate oscillations in the cortex. We developed a computational model of gap junction plasticity in a recurrent cortical network based on recent experimental findings. We showed that gap junction plasticity can serve as a homeostatic mechanism for oscillations by maintaining a tight balance between two network states: asynchronous irregular activity and synchronized oscillations. This homeostatic mechanism allows for robust communication between neuronal assemblies through two different mechanisms: transient oscillations and frequency modulation. This implies a direct functional role for gap junction plasticity in information transmission in cortex.


Introduction
Oscillatory patterns of neuronal activity are reported in many brains regions with frequencies ranging from less than one Hertz to hundreds of Hertz.These oscillations are often associated with cognitive phenomena such as sleep or attention.Local field potential measurements in the neocortex and thalamus show the prevalence of delta oscillations (0.5-4Hz) and spindle oscillations  during sleep [1].Theta oscillations (4-10Hz) are also reported in hippocampus and other brain regions [2].Gamma oscillations  observed in the cortex are thought to be involved in attention [3][4][5][6], perception [7,8] and coordinated motor output [9,10].Thus, at the minimum, oscillations are present during the normal functioning of neural circuits.
However, oscillations are also associated with pathological circuit dynamics, such as hypersynchronous activity during epileptic seizures [11].Altered gamma-frequency synchronizations may also be involved in cognitive abnormalities such as autism [12] or schizophrenia [13].Thus, given both the functional and pathological effects of oscillations, a homeostatic mechanism is necessary to regulate oscillatory behavior.
Several mechanisms can lead to the emergence of oscillations.They can arise in homogeneous population of excitatory neurons, where the positive feedback loop of excitation is only limited by the refractoriness of the neurons [14].Alternatively, oscillations can also arise in a coupled network of excitatory and inhibitory neurons, where the excitatory and inhibitory neurons burst in opposing phase.[15][16][17][18][19]. Finally, gap junctions between inhibitory neurons promote synchronous oscillatory patterns [20][21][22][23][24].
The inhibitory network oscillations primarily involve fast-spiking interneurons.These neurons represent a large proportion of GABAergic interneurons [25].They are the main cells targeted by thalamocortical synapses transmitting sensory information to the cortex [26].They are coupled via chemical synapses and gap junctions.Gap junctions are mostly found between neurons of the same class [26][27][28] but they can also connect different subtypes, such as fastspiking and regular spiking cells [26,29,30].Moreover, there is evidence of the critical role of fast-spiking parvalbulmin (FS) interneurons in the emergence of cortical gamma activity in the cortex of rodents in response to sensory stimuli [31][32][33][34].
Recently, gap junction plasticity has been experimentally demonstrated [47][48][49][50][51].For example, the gap junctions between rod cells in the retina can vary their conductance during day and night cycles [52].Moreover, they can experience bidirectional long-term plasticity in an activity-dependent manner [49,53,54].High frequency stimulation of a coupled pair of thalamic reticular nucleus (TRN) neurons induces gap junction long-term depression (gLTD) [55].This occurs only when the TRN neurons burst.There is no data yet on the long-term potentiation of cortical gap junctions.However, [56] show that the pathways leading to gLTD are calcium-dependent which suggest that gap junction long-term potentiation (gLTP) could also be the result of an activity-dependent mechanism.Other passive mechanisms, such as gap junction connexin turnover could compensate for long-term depression as well [57][58][59][60][61][62].
Given the existence of gap junction plasticity and the omnipresence of oscillations in cortex, we wondered whether gap junction plasticity can regulate network-wide gamma oscillations in cortex.To that end, we developed a computational model of a network of excitatory and FS inhibitory neurons.As demonstrated analytically by [24], we observed two different network behaviors depending on the gap junction strength.For weak gap junction strength, the network exhibits an asynchronous regime, whereas for strong gap junctions, the network synchronizes into coherent gamma oscillations with bursting activity.We then modelled the gap junction plasticity observed by [55] showing that bursting activity leads to gLTD.The plastic network sets itself at the transition between the asynchronous regime, where sparse spiking dominates, and the synchronous regime, where network oscillations dominate and burst firing prevails.Thus, our model shows that gap junction plasticity maintains the balance between the asynchronous and synchronous network states.This is robust to different possible gLTP rules.We then show that the network allows for transient oscillations driven by external drive.This demonstrates that transient, plasticity regulated oscillations can efficiently transfer information to downstream networks.Finally we show that gap junction plasticity mediates crossnetwork synchronization and allows for robust information transfer trough frequency modulation.Critically, gap junction plasticity allows for the recovery of oscillation mediated information transfer in the event of partial gap junction loss.

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 subthreshold resonance in the gamma range [42][43][44][45][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 gammafrequencies [33].
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).
Activity-dependent gap junction long-term potentiation (gLTP) has not been reported experimentally yet in the mammalian brain.There is evidence for activity dependent shortterm 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 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][67][68][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 gammanetwork.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 (  For panels E, F, H, I, the x-axis represents the number of crossnetwork 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.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.
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 ).
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.

Discussion
Our modelling study tested whether gap junction plasticity can regulate gamma oscillations in cortical network models.Our findings suggest that gap junction plasticity can maintain a balance between synchronous regular and asynchronous irregular regimes.For strong electrical coupling, the network is in the oscillatory regime.The oscillations consist of synchronized bursting mediated by the inhibitory neuron network.These bursts trigger depression of the gap junctions [55] allowing the network to leave the oscillatory regime and spike asynchronously.However, the irregular asynchronous regime is dominated by sparse firing.Either this sparse firing, or constant protein connexin turnover may be a source of gap junction potentiation [48,[56][57][58][59][60][61].Thus, the asynchronous irregular regime tends to potentiate gap junctions.Therefore, the network behavior critically depends on the plasticity learning rate.Fast gLTP leads to synchronous activity while slow gLTP leads to asynchronous states.We demonstrate the functional role of plasticity in both cases.In the AI regime, the network can respond to changes in input drives through transient oscillations.Those transient oscillations could serve as an energetically efficient way to transfer information to a downstream neuron.In the SR regime, the network oscillations can serve as the substrate for information routing between networks.These results demonstrate how gap junction plasticity can regulate oscillations to mediate information transfer between cortical populations of neurons.

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][72][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][21][22][23][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.

Methods
We consider a network with N I inhibitory neurons (20%) and N E excitatory neurons (80%) with all-to-all connectivity (Fig 1A).Inhibitory neurons are modelled by an Izhikevich model and excitatory neurons by a leaky integrated-and-fire model (LIF) [63,110].The simulation time-step is dt = 0.1 ms.Inhibitory neurons are connected by both electrical and chemical synapses, whereas excitatory neurons have only chemical synapses.We designed a novel plasticity model for activity dependent plasticity of gap junctions and we investigated its impact on network dynamics and function.We then investigated the dynamics of two networks coupled by chemical and electrical synapses.We use a decoder to quantify the effects of gap junction plasticity on information transfer.The model is written in Python and takes advantage of the tensorflow library that leverages GPU parallel processing capabilities [111].It is available on ModelDB (http://modeldb.yale.edu/230324).

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 combined with the spiking conditions, 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 0 cos(2πft) (with I 0 = 0.01 pA).We then normalized the amplitude response as follow 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, 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 where is the current coming from the transmission of a spike via electrical (i.e.spikelet) and chemical synapses, I gap i is the sub-threshold current from electrical synapses (for inhibitory neurons only), I noise i is the noisy background current and I ext i characterizes the external current.The current due to spiking I spike i on excitatory neurons is given by The current I spike i into inhibitory neurons are where W αβ is the coupling strength from population α to population β with {α, β} = {E, I}.
e ij is the inhibitory to inhibitory coupling between neuron i and j, consisting of the chemical synaptic strength W II,c and W II;e ij the electrical coupling for suprathreshold 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 W II;e ij ¼ k spikelet Ã g ij , 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 W II ij 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 exp À tÀ t jk t a 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 subthreshold current I gap i is characterized by 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 and 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 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, 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 −1 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 where α gLTP is the potentiation learning rate and sp i (t) = ∑ t ik <t δ(t − t ik ).
The second gLTP rule we consider has no maximum bounds 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 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 and Then, we explore the values of the ratio of bursting over spiking activity ratio ¼ A bursting A spiking ð19Þ as function of the coupling coefficient γ and of the mean external input ν over the parameter space P 1 ¼ ½0; g max Â ½0; n max .

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 rðtÞ ¼ 1 dt 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 where the r n sequence represents N uniformly spaced time-samples of the population activities.
We measure the amplitude of the Fourier components rk for k = 1..N/2 (because the Fourier signal is symmetric from N/2 to N).We identify the maximal one, its associated frequency f max ¼ k N and its power P ¼ ðjr k j=NÞ 2 .

Downstream read-out neurons
To simulate the projection of a cortical layer onto another layer, we model downstream readout 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: 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.
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 IðX; YÞ ¼   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 ðr GN ?r SN ÞðtÞ ¼ 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 t r _ r filt ðtÞ ¼ À r filt ðtÞ þ r ON ðtÞ; ð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 t Ã k , 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 T k ¼ t Ã kþ1 À t Ã k .For Fig 5D , we plot " x k , the mean values of the input signal x between the rising times t Ã k and t Ã kþ1 as function of their corresponding peak intervals T k 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.
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 ( 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 x.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.

Fig 1 .
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

Fig 2 .
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 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.

Fig 3 .
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 crossnetwork 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.
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 crossnetwork 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

Fig 4 .
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

Fig 5 .
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

X
x;y pðx; yÞ log pðx; yÞ pðxÞpðyÞ ; ð23Þwhere p(x, y) is the joint probability function of X and Y, and p(x) and p(y) are the marginal

def Z 1 À 1 r
GN ðtÞ r SN ðt þ tÞ dt: ð24ÞThe phase difference is measured as the time delay relative to the oscillation periodD0 ¼ arg max t ððr GN ?r SN ÞðtÞÞ T period ;ð25Þ Fig 5F) and amplitude A (Fig 5G).xðtÞ¼ A½cosð2pftÞ þ 1ð29ÞThen we measure the routing of random signalsx(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 xðtÞ 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.