Dual Mechanism for the Emergence of Synchronization in Inhibitory Neural Networks

During cognitive tasks cortical microcircuits synchronize to bind stimuli into unified perception. The emergence of coherent rhythmic activity is thought to be inhibition-driven and stimulation-dependent. However, the exact mechanisms of synchronization remain unknown. Recent optogenetic experiments have identified two neuron sub-types as the likely inhibitory vectors of synchronization. Here, we show that local networks mimicking the soma-targeting properties observed in fast-spiking interneurons and the dendrite-projecting properties observed in somatostatin interneurons synchronize through different mechanisms which may provide adaptive advantages by combining flexibility and robustness. We probed the synchronization phase diagrams of small all-to-all inhibitory networks in-silico as a function of inhibition delay, neurotransmitter kinetics, timings and intensity of stimulation. Inhibition delay is found to induce coherent oscillations over a broader range of experimental conditions than high-frequency entrainment. Inhibition delay boosts network capacity (ln2)−N-fold by stabilizing locally coherent oscillations. This work may inform novel therapeutic strategies for moderating pathological cortical oscillations.

change the rise and decay time of the postsynaptic current, and vary the synaptic conductance (Supplementary Methods I,II). Accordingly, individual synapses have a tuneable inhibition delay d which we vary from 20 μs to model the latency time of neurotransmitter release 35,36 , to 800 μs to model the transmission line delay of inhibitory signals as they diffuse along the dendrites towards to the axon hillock of dendrite projecting interneurons 37 . These inhibition delays are chosen to match the transit time of action potentials across the 200 μm-700 μm long dendrites of somatostatin interneurons 37 at an average speed of 1-100 m/s ( Fig. 1(a)). The decay (resp. rise) time of the postsynaptic current was set by the undocking (resp. docking) time of neurotransmitters on neuroreceptors (GABA), τ u (resp. τ d ). τ u was tuned over 0-8 ms, a range comparable to the period of neuron oscillations: 5-20 ms 38 (Fig. 1(b)).
Synaptic kinetics of the half-center oscillator. We began to study the emergence of synchronization by probing the synchronization phase diagram of a pair of mutually inhibitory neurons as a function of synaptic kinetics in all connections (d, τ u ) and current stimulation applied to all neurons (I stim ). When inhibition delay is small ( μ < d 150 s), three modes of synchronized oscillations are observed as I stim increases ( Fig. 1(c)). Above the depolarization threshold (I th = 8 μA), neurons oscillate out-of-phase (antiphasic synchronization). They suddenly lock in phase (phasic synchronization) at I s = 14 μA. Higher current stimulation ( > I I stim s ) increases the frequency of neuron oscillations and makes inhibition increasingly tonic. As a result neurons decouple gradually. This loose coupling regime is characterized by higher order phase locking where a neuron entrain the other at a frequency which is a rational multiple of its own ( Fig. 1(c)).
Longer inhibition delays ( μ > d 150 s) broaden the synchronization current I s to a window of finite width [I L , I H ] ( Fig. 1(d)) which increases and eventually diverges at μ > d 300 s. The observation of phasic synchronization at longer inhibition delay concurs with similar results obtained by Van Vreeswijk et al. 11 when the synaptic response time becomes slower. Antiphasic, phasic, and loose coupling regimes form 3 domains in the d − I stim phase diagram of Fig. 1(e) showing that phasic synchronization may be induced either by delaying inhibition or by applying a stimulation current close to I s . Delayed inhibition gives each neuron in the pair the time to depolarize prior to receiving inhibition from its partner. This condition is necessary but not sufficient to explain phasic synchronization. Inhibition delay also decreases the slope of the phase response curve of the post-synaptic neuron near the origin (Supplementary Methods II). This reduces the phase correction that mutual inhibition applies to the early and late firing neurons which has the effect of stabilizing synchronous oscillations.
For shorter inhibition delays ( μ < d 150 s), the synchronization current (I s ) and frequency (f s ) decrease when τ u increases. This dependency is well explained by calculating the frequency of phase synchronized oscillations f p (Supplementary Discussion I) and its intercept with the excitatory response curve of a neuron ( Fig. 1(d)). We find τ ∼ − f s u 1/3 ( Fig. 1(f)). This result concurs with the onset of γ-oscillations shifting to lower frequency (current stimulation) following pharmacological manipulations that increase the recovery time of the postsynaptic current 7,8 .

3-cell mutually inhibitory network.
Larger inhibitory networks ( ≥ N 3) generally have chaotic dynamics which makes network oscillations highly dependent on the timings of current stimuli. We defined the state of the system using the phase lags of individual neurons relative to a reference (neuron 1) and obtained the state trajectories by measuring the temporal evolution of these phase lags ΔΦ { } i p 1 ( ) , i = 2, 3 ...N over consecutive periods p = 1-50. The phase lag map of a 3-neuron network with 300 μs inhibition delay shows state trajectories converging towards 6 point attractors ( Fig. 2(a)). These attractors are sub-divided into 3 categories according to the duration of their interspike intervals (ISI): T/3, T/2 and T where T is the period of synchronized oscillations ( Fig. 2(b)). Two attractors (circle symbols) correspond to three neurons discharging in the clockwise and anticlockwise sequences, 1 → 2 → 3 and 1 → 3 → 2 (ISI = T/3). Three attractors (square symbols) correspond to 3 modes of partially synchronized oscillations including the sequence → 1 2 3 and its 2 permutations (ISI = T/2). The single coherent attractor (diamond symbol) corresponds to all 3 neurons discharging in phase (ISI = T). The 3-neuron map shows the basins of attraction becoming smaller as oscillations become more coherent. This demonstrates the greater fragility of coherent states relative to the oscillations of sequentially discharging neurons. We find that for μ > d 300 s, coherent and partially coherent oscillations become stable over the entire range of current stimulation. If μ < d 150 s however, the network only supports the oscillations of sequentially discharging neurons, as we shall see below. We find that substituting non-delayed inhibitory synapses (d = 0) with gap junctions 29 produces qualitatively similar phase portraits in that they only support sequentially discharging neurons ( Fig. 2(c)). For completeness, we also considered gap junctions between excitatory neurons. We find that the excitatory network hosts a single state of collective oscillations ( Fig. 2(d)). This expected result validates the correct operation of our analogue network. Returning to the 3-neuron network connected by non-delayed inhibitory synapses, and varying current stimulation applied to all neurons, we find that partially coherent oscillations vanish except in a very narrow range of current stimulation centered on I s -as in the neuron pair.
In the 3-neuron and 4-neuron networks, the synchronization current I s is the current that maximises the size of the coherent basin of attraction and stabilizes the coherent attractor with respect to noise (Fig. 3). For long inhibition delays (d = 350 μs), the network supports coherent oscillations over the entire range of current stimulation. When μ < d 150 s, coherent oscillations only form in a narrow range of current stimulation about I s . These observations generalize the d − I stim phase diagram of Fig. 1(e) to larger networks and demonstrate that synchronization may be achieved either through increases in inhibition delay or current stimulation.
Emergence of synchronization in all-to-all inhibitory networks. We next demonstrate the emergence of synchronization in larger networks (N = 3, 4, 5) and the critical importance of inhibition delay in stabilizing locally coherent oscillations. The maximum number of attractors in a N-neuron network was calculated by counting the number of cyclically invariant discharge patterns allowing partial synchronization (Supplementary Discussion II). We find that the maximum network capacity increases as 39 . The minimum capacity, allowing sequential discharges only, is L N = (N − 1)! Experimental results show that the capacity of an inhibitory network to encode information about its environment lies between L N and T N , depending on inhibition delay (Fig. 4). Longer inhibition delays (d = 400 μs) stabilize oscillations which range from purely phasic (Fig. 4: (a) diamond, (b) triangle, (c) hexagon) to purely sequential ( Fig. 4(a-c) circles). In between, all intermediate states of partial synchronization are observed ( Fig. 4(a-c)). For example, the 4-neuron map in Fig. 4(b) has 6 sequential attractors with 1 spike per ISI giving ISI occupancies (1, 1, 1, 1) (circle symbols), 12 partially synchronized attractors with ISI occupancies (2, 1, 1, 0) (square symbols), 4 + 3 partially synchronized attractors with (3, 1, 0, 0) and (2, 2, 0, 0) occupancies respectively (diamond symbols), and the coherent attractor (4, 0, 0, 0) (triangle symbol). Therefore the 4-neuron network hosts 26 attractors in total.
When inhibition delay is reduced further (d = 100 μs), the only attractors left are sequential oscillations ( Fig. 4(g-i)). The network capacity then scales as: 2 (N = 3), 6 (N = 4), 24 (N = 5) which matches the L N sequence above. These results demonstrate that, provided the inhibition delay is sufficiently large, the number of attractors increases according to sequence T N . For this, the inhibition delay needs to be at least 1/3 of the duration of the action potential ( > d W/3). The network capacity was found to be less sensitive to neurotransmitter kinetics. Increasing τ u from 1.5 ms to 3.5 ms marginally increased the number of attractors. No further change was observed beyond τ > .

Discussion
Our results suggest that inhibitory networks may synchronize via two mechanisms that exploit the distinct neurophysiological properties of fast-spiking interneurons 36 and the inhibition delay introduced by dendrite projecting synapses 37 . This study considers the primary effect of dendrite targeting synapses to be the introduction of a transmission line delay because the network frequency covers a very narrow range set by the constant step amplitude of current stimuli. The complex spectral response of dendrites is however known to be important and would need to be considered if the amplitude of current stimulation was varied. Dendrite projecting somatostatin interneurons introduce transmission line delays of the order of 0 -800 μs by projecting synapses on the 200-700 μm long dendrites of the mammalian visual cortex 37 . Transmission line delays of this magnitude postpone the onset of inhibition sufficiently to stabilize the coherent oscillations of inhibitory neurons ( Fig. 4(a-c)). The anatomical properties of somatostatin neurons would thus warrant robust phasic synchronization which is weakly dependent on current stimulation or postsynaptic kinetics but is strongly dependent on the timings of stimulation. This result is consistent with the rapid attenuation of visually induced γ-oscillations observed when visual stimuli become uncorrelated 6 . The coherent attractor is unique and its basin occupies a very small volume of phase space (triangle symbol, Fig. 4). As a result, the state of collective synchronization is the least robust of all states with respect to noise and structural inhomogeneity. In contrast, the bulk of the phase space is filled with partially coherent attractors whose proportion increases very rapidly according to 1 − (ln2) N as the network size increases. Using this expression, one calculates that partially coherent attractors form > .
98 7% of all attractors for the typical neuronal population, > N 12, excited during optogenetic experiments 6 . Besides being more numerous, partially coherent states also have wider basins which offer protection from decoherence by noise and structural heterogeneities (Fig. 4(a-c)). Accordingly, partially coherent states are the most thermodynamically stable with respect to coherent and sequential states and are the most likely to support synchronized electrical activity in the noisy environment of real cortical networks. Within partially coherent states, however, the neurons which oscillate in phase may distribute differently over the volume of the network. A subset of L neurons ( < L N) may oscillate in phase at different locations of the network, producing spatially homogeneous firing akin to the fully synchronized state. Two partially coherent states with identical L-number differ through the permutations of stimuli. The equivalence of these states is demonstrated by the six-fold symmetry of phase maps of the 4-neuron network ( Fig. 4(b)).
Our results suggest that spatially homogeneous firing within partially coherent states may be promoted by local repulsion through gap junctions 41 . These junctions are known to predominantly couple neighbouring inhibitory cells of the same population 42,43 . As we have seen in Figs 2(c) and 4(g), gap junctions and fast inhibitory synapses share the property of supporting sequential neuronal oscillations. Electrical synapses thus have a destabilizing effect on local neural synchronization as reported in earlier numerical simulations 21,23,44 . At the same time, Fig. 4 show that transmission line delays promotes synchrony. An inhibitory network can thus achieve a homogeneous distribution of phasic neurons 45 by breaking local coherence using gap junctions. Homogeneous firing is established from the long range attraction of delayed inhibition and the short range repulsion of electrical synapses. Note that many physical systems achieve long range order through short range repulsion. For example, the Wigner crystal arises from Coulomb repulsion between electrons 46 and vortex-to-vortex repulsion is responsible for the Abrikosov lattice in type II superconductors 47 . The effect of introducing heterogeneity in the network is seen in Fig. 4(a-c) where residual imbalance in network conductance breaks the symmetry of phase lag maps. Introducing a range of inhibition delays or mixing gap junctions with chemical synapses would similarly increase the volume of some basins -those associated with spatially homogeneous firing -to the detriment of others 26 . In contrast to somatostatin neurons, the wiring of parvalbumin neurons introduces delays which are too short to warrant automatic synchronization. Instead parvalbumin neurons may achieve synchronization through high frequency entrainment. This corresponds to the current induced synchronization which we observe at small d ( Fig. 1(c,e)). Because frequency f s is dependent on neurotransmitter kinetics ( Fig. 1(f)), this synchronization mechanism allows the onset of synchronized oscillations to be tuned using pharmacological manipulations targeting GABA receptors 4,5,7,9,48 .
Our study leads us to propose that local cortical circuits may have adapted to exploit the robustness of synchronization by delayed inhibition versus the tunability of synchronization by fast-spiking interneurons ( Fig. 1(e)). These synchronization mechanisms suggest strategies to reduce pathological cortical oscillations which include: inactivating dendrite targeting synapses, blocking GABA B receptors to accelerate the recovery of the postsynaptic potential, and applying visual stimuli lacking spatial coherence at frequencies in the γ band. This study has focussed on purely inhibitory networks (ING) which have intrinsically chaotic dynamics. The consideration of excitatory neurons and feed-forward processes within the pyramidal-interneuron-gamma (PING) mechanism invokes regular dynamics which has been treated elsewhere 8 .

Methods
Electronic models. We synthesized two VLSI networks interconnecting 6 Mahowald-Douglas neurons 30 with either inhibitory synapses or gap junctions (Supplementary Methods I). VLSI neurons modelled the dependence of the membrane voltage V on current stimulus I stim using the analogue electrical equivalent circuit of the neuron membrane. Its equation was stim where E Na and E K are the sodium and potassium reversal potentials and C is the membrane capacitance. The sodium and potassium conductances, g Na and g K , are modelled by the transconductances of p− and n− type field effect transistors respectively 30 . The gate variables m, h and n of the Hodgkin-Huxley model are represented in the analogue circuit by currents ι which are either activated or inactivated according to: where x ≡ {m, h, n}, V x is the threshold voltage of each ion gate, and dV x is the width of the transition from the closed to the open state of that gate. The V τ,x variables follow a first order dynamics which describes the recovery of each gate variable and is characterized by recovery time τ x 29 . Chemical synapses were implemented using a differential pair integrator 31 (Supplementary Methods II). As our transistors functioned with above threshold currents as opposed to below threshold 31 , the postsynaptic current was approximately given by I post (t) = gS(t)(V post (t) − V rev ) where V rev = 7 V was the reversal potential, V post (t) the membrane voltage of the postsynaptic neuron, g the maximum conductance and S(t) was the fraction of docked neurotransmitters at time t. The neurotransmitter docking rate was given by: The empirical inhibition delay d, decay time τ u and synaptic conductance g were controlled by 3 gate voltage parameters: V th , V W and V τ in the circuit (Supplementary Methods II). The synaptic conductance varied in the range g = 1-3 μs.
We implemented gap junctions electronically using a differential transconductance amplifier to model electrical coupling between GABAergic-like interneurons 41 . Their current-voltage transfer characteristics has been measured by Zhao and Nogaret 29 . The gap junction current varies linearly as post post pre near the balance point of the pre-synaptic and post synaptic membrane potentials 41 . The transconductance ′ g is tuneable in the range 24 μs μ < ′ < g 45 S using the gate bias V M of the current source transistor (Fig. S7). Away from the balance point, saturation effects reduce the rate of current injection 29 . We were able to change the sign of the injected current by swapping the voltage inputs and in this way obtain either an inhibitory or an excitatory link ( Fig. 2(d)).
Circuits were built from VLSI current mirrors (ALD1116, ALD117). The depolarization threshold of neurons was adjusted to match the range of synaptic currents. This was done by adjusting the leakage conductance of the neuron membrane. The current thresholds were I th = 8 μA (synaptic coupling) and 86 μA (gap junction coupling). The duration of an action potential was W = 1 ms.
Data acquisition and analysis. Individual neurons were stimulated by timed current steps of constant amplitude I stim . These stimuli were generated by the analogue outputs of two DAQ cards (NI PCI6259) and a bank of 6 voltage-to-current converters. Labview code was written to vary the timings of current stimuli in a systematic manner so that initial conditions meshed the (N − 1)-dimensional phase space with a grid size of T/20. The Labview/DAQ card recorded the membrane voltage time series of individual neurons during each current protocol. The sampling frequency was 20 kHz. Between the end of one protocol and the beginning of the next, a 200 ms long time window was inserted during which no stimulation was applied to let the system return to its steady state.
The where θ 4 = π/12, and: where θ 5 = π/4. The state trajectories pertaining to the same basin were regrouped using Matlab code which calculated the coordinates of experimental attractors and their total number.