Oscillations emerging from noise-driven steady state in networks with electrical synapses and subthreshold resonance

Oscillations play a critical role in cognitive phenomena and have been observed in many brain regions. Experimental evidence indicates that classes of neurons exhibit properties that could promote oscillations, such as subthreshold resonance and electrical gap junctions. Typically, these two properties are studied separately but it is not clear which is the dominant determinant of global network rhythms. Our aim is to provide an analytical understanding of how these two effects destabilize the fluctuation-driven state, in which neurons fire irregularly, and lead to an emergence of global synchronous oscillations. Here we show how the oscillation frequency is shaped by single neuron resonance, electrical and chemical synapses.The presence of both gap junctions and subthreshold resonance are necessary for the emergence of oscillations. Our results are in agreement with several experimental observations such as network responses to oscillatory inputs and offer a much-needed conceptual link connecting a collection of disparate effects observed in networks.

C ognitive phenomena such as conscious perception and attention are associated with cortical oscillatory activity in the frequency range from lower than one to a few hundred Hertz [1][2][3][4] . A number of theories have been proposed for the generation of oscillatory rhythms that emphasize the importance of a specific subtype of neurons and their synaptic connections. Inhibitory interneurons are implicated in oscillatory cortical activity in many brain regions. An important feature of inhibitory interneurons is the presence of a subthreshold resonance that selectively amplifies firing rate response to select frequencies already present in isolated neurons [5][6][7][8][9] . This alone could suffice for the inhibitory interneurons to be the driving force behind global network oscillations. In addition, those same inhibitory interneurons are coupled by gap junctions in many cortical regions where global network oscillations are reported [10][11][12][13] . It is therefore conceivable that gap-junction-induced synchrony and subthreshold oscillations provide a substrate for global oscillations in any regions where they co-occur. Although a large number of theoretical studies have tackled the effects of chemical synapses on synchrony and global oscillations, the role of electrical synapses has received less attention and has been mostly studied in networks with identical neurons in the low noise limit and largely without considering the frequency preferences of single neurons [13][14][15][16][17] . In physiological conditions in vivo, neurons receive noisy input and fire irregularly a situation where global oscillations could be difficult to achieve and the synchronizing effect of gap junctions and single neurons could be cancelled out. It is therefore important to understand how the single neural properties interact with chemical and electrical connectivity in the noisy fluctuation-driven regime, in which neurons fire irregularly because they are driven by noisy inputs, and how these networks can transition to global oscillations and synchrony. The current lack of theoretical approaches so far is largely due to the technical challenges of solving multiple coupled Fokker Planck equations that are necessary for the inclusion of gap junctions and single-neuron frequency preference. Therefore, most of the computational studies consider either gap junctions or resonant single neurons, never both simultaneously and largely focus on numerical simulations [18][19][20][21][22] .
To overcome these difficulties and address analytically the interplay between the single neuron biophysics and the chemical and electrical connectivity, we develop an adaptive threshold model framework, where we consider neurons in the fluctuationdriven regime, that is driven by largely fluctuating inputs, that results in irregular firing, while neglecting their reset. Inhibitory neurons have a subthreshold resonance, they are connected electrically among themselves and chemically to excitatory neurons. Excitatory neurons do not exhibit subthreshold resonance and are chemically coupled. Within this framework we compute the sub-and superthreshold resonant properties of single neurons and use this framework to set up a self-consistent mean-field description. Using the mean-field approach, we compute the network dynamics in response to an incoming oscillatory current. Notably, we show that the resonance frequency of the network reflects directly the subthreshold preference in individual neurons of a subset of neurons, while chemical synapses act as a multiplicative gain that regulates the amplitude of this response but not its shape.
We then derive the oscillatory instability condition for the network to transition from irregular spiking to global selfsustained oscillations. In the oscillatory state, the neurons fire synchronously and therefore are no longer in the fluctuationdriven regime. We demonstrated analytically two requirements for the emergence of self-sustained global oscillations from the fluctuation-driven state, (i) frequency preference in the inhibitory neurons and (ii) inhibitory neurons connected with gap junctions, where the gap-junction-induced spikelet is larger than the g-aminobutyric acid (GABAergic) inhibitory postsynaptic current (IPSC), and leads to global excitatory coupling. It is important to emphasize that in our model synchronous global oscillations are characterized by one spike per neuron and cycle, and emerge from the irregular fluctuation-driven network state. The frequency of global oscillations is largely defined by single neuron subthreshold oscillations. Notably, this is a fundamentally different mechanism compared with previous studies that describe the emergence of sparsely synchronized fast oscillations 23 , or gap-junction-mediated synchronous oscillations in networks without subthreshold frequency preference 24 or global oscillations in network-coupled oscillators 15 .
Here we demonstrate that the theoretical predictions we derived are consistent with the recent physiological evidence from cell-specific light stimulation in the barrel cortex 22 . Adapting our framework to the cell properties of the barrel cortex, we show that the excitatory neurons amplify low frequencies at the network level, whereas the inhibitory neurons alone can support a resonance in the g-range, which is consistent with ref. 22. Our analytical adaptive framework is not limited to the case of the barrel cortex, but can be applied to any neuronal network in the fluctuation-driven state with subthreshold resonance and gap junctions.

Results
Subthreshold resonance leads to firing rate resonance. We model inhibitory neurons via a leaky voltage dynamics coupled with an adaptive variable. This model exhibits subthreshold oscillations that depend on the relative time constants of the voltage (t V ) and of the adaptive variable (t w ), as shown in Fig. 1a and mathematically derived in the Methods section. Next, we show the effect of this subthreshold resonance on the firing rate susceptibility to periodic inputs. Figure 1b demonstrates that the Q-value (peak/width) of the resonance peak is highest in the parameter regime where subthreshold resonance is most salient. The preferred frequency for the rate response is very closely related to the preferred subthreshold resonance frequency (Fig. 1b,c). Albeit sizable differences between rate and subthreshold resonance can be observed in Fig. 1c close to the edge of resonant parameter space. Importantly, this tight correspondence between subthreshold and rate resonance we observe for the threshold model framework is consistent with the previous observation by Richardson et al. 21 They found that fluctuationdriven but not mean activity can unmask the subthreshold oscillations in the spike firing patterns of leaky integrator neurons. Finally, the predicted steady-state firing rate we calculated analytically is consistent with numerical simulations, Fig. 1d. The firing rate increases with the current drive and decreases with the adaptive coupling a, consistent with previous experimental and theoretical reports 21 .
Response of a single neuron to an oscillatory current. To address how networks respond to dynamical stimuli, the first step is to understand how the firing rate of single neurons responds to weak sinusoidal current drive of frequency f (refs 19,21). We therefore derive the linear response functions R j (o) of a neuron receiving an oscillatory current (see Methods). Figure 2 demonstrates the linear response functions for inhibitory and excitatory neurons.
The salient feature of the inhibitory linear response is the resonance emerging from the adaptive voltage coupling a. Figure 2a shows that increases in a result in a higher resonance peak and a modest shift to higher frequencies. Figure 1b,c also demonstrate the firing rate frequency response and its sharpness in the frequency domain. Q-value that we use to quantify sharpness is defined as Df/f R where Df is the full width at half maximum and f R is the peak frequency. The excitatory neurons, however, show amplification of lower frequencies only Fig. 2b because, by assumption, they lack subthreshold adaption (a ¼ 0). Their firing rate response drops at higher frequencies due to the membrane time constant of the leaky integrator.
To show the robustness of our findings across different spiking models, we compare the threshold model with the adaptive exponential integrate-and-fire (aEIF) model, which is shown to reproduce different firing patterns 25 . In Fig. 2c, we show that the inhibitory response function in both models are very similar even at high rates, firing rate around 35 Hz in both models. A more detailed response functions comparison across models can be found in Figs 1 and 2 in the Supplementary Methods.
Network response -Link to experiments. Studying the network response to external periodic drive in Fig. 3, we discovered that stimulation of excitatory and inhibitory neurons yield qualitatively different results. If the oscillatory current is delivered to excitatory neurons, the network only amplifies lower frequencies (Fig. 3b). If the oscillatory current is delivered to inhibitory neurons, the network shows a resonance peak that resembles the rate resonance of single inhibitory neurons (Fig. 3c). This indicates that the location and width of the inhibitory subthreshold resonance rather than recurrent chemical synaptic connectivity is a major determinant of network dynamics in response to external stimuli. We note that the form of rate response functions does not depend on the details of the spike generation mechanism such as voltage reset and is preserved for the adaptive exponential and leaky integrate-and-fire models (see Fig. 3c the black and grey line, as well as Supplementary Fig. 1b). We therefore expect these results to generalize to other classes of networks where subpopulations of neurons exhibit differences in subthreshold frequency preferences.
Next, we explore whether parameter regimes exist where stimulation of the excitatory neurons could show signatures of the inhibitory subthreshold resonance. To investigate this hypothesis, we increased threefold the E to I connectivity strength G IE and shortened the voltage time constant of the excitatory population t V,E to enable excitatory response function to cover the g-frequency range. Supplementary Fig. 3 shows that these two interventions are sufficient to induce g-band resonance on excitatory stimulation. However, it is hard to find experimental evidence supporting G IE 4 4G EI , G EE such that we consider the parameter set of similar recurrent connectivities in Fig. 3 to be a more plausible case in cortical networks.
We find that results in Fig. 3 are consistent with the experimental finding of Cardin et al. 22 in the barrel cortex in vivo. These experiments show that inhibitory neurons rather  Table 1.   Table 1.
NATURE COMMUNICATIONS | DOI: 10.1038/ncomms6512 ARTICLE than excitatory neurons amplify g-range frequencies. The external stimulation at varying stimulation frequencies in Cardin et al. 22 corresponds to the network linear response we have studied in the previous two sections. In Fig. 3b-d, we show that our model can reproduce the high spike probability per cycle (around 0.9 for 40 Hz) and recover the frequency response dynamics found in neural networks under excitatory and inhibitory stimulation in Cardin et al. 22 This framework offers us the opportunity to investigate the contribution of the single neuron properties and connectivity to the network response dynamics. Cardin et al. 22 addressed this question by blocking successively excitatory and inhibitory connections. They showed that stimulating the inhibitory neurons optogenetically at 40 Hz lead to a large local field potential (LFP) response, which decreased by B70% when the excitatory connections are blocked, and decreased to almost baseline when all the chemical connections were blocked (see  Fig. 3f, we mimic this experiment and reproduce quantitatively the experimental results. However, in our model, blockade of connectivity is a combination of removing recurrent excitatory synapses and reducing external excitatory drive. Similarly, blockage of GABAergic connections correspond to a removing of inhibitory synapses in our model, increasing the network drive and reducing the variance of the noise. This is due to the fact that in the irregular activity of cortical neurons, the mean and the variance of the firing rate co-vary 26 .
Notably, other models such as PING/ING models have been tried to explain the findings of Cardin et al. 22  Predictions for network's response to an oscillatory drive. The response of a network subject to external stimuli is influenced by a number of effects such as external drive and recurrent chemical and electrical coupling. To consider each effect separately we turn to Fig. 3g. In Fig. 3g, we show the separate contribution of the mean drive, variance drive, recurrent connectivity and gap junctions on the network response to an oscillatory current at 40 Hz. We find that the inhibitory subthreshold resonance rather than recurrent chemical synaptic connectivity is a major determinant of network dynamics in response to external stimuli. We find the largest drop in response amplitude when excitatory external drive's mean or variance is reduced. On the other hand, reductions in recurrent conductivities by and large introduce only modest effects. This could be tested experimentally by opto-    Table 1 genetically silencing the thalamic neurons providing inputs to the cortical neurons. This manipulation should show a drastic lowering of LFP g-power. However, silencing the local cortical excitatory neurons should only show a subtle drop in g-power.
In summary, our findings support the hypothesis that intrinsic biophysics of excitatory cells is the primary cause for the reduction of excitatory neurons ability to respond primarily to lowfrequency stimulation rather than synaptic recurrent connectivity, or even the recruitment of a distinct subgroup of low-threshold spiking neurons that selectively enhance lower-frequency oscillations 20 . Therefore, our study offers a complementary explanation for the observed experimental findings and supplements previous multi-compartment large-scale simulations by offering a compact analytical treatment that combines both electrical, chemical and single-neuron subthreshold contributions.
Predictions for self-sustained network oscillations. In this section, we characterize the emergence of self-sustained global oscillations occurring when the asynchronous irregular network state becomes unstable with respect to spontaneous oscillatory perturbations. We consider the stability conditions outlined in refs 23,29 and apply them to our network. Starting with an infinitesimal time-dependent perturbation in input current, we obtain changes firing rate via R j (o) (equation (3)) and propagate these changes sequentially back to the current level via S j (o) (equation (4)). The stability of this transformation is given in equation (6) and its two population analogon can be found in equation (47) in the Supplementary Methods. It is important to emphasize that these stability conditions can predict where this oscillatory transition boundary should occur in parameter space and what oscillatory frequency emerges at this transition.
We also address the stability of the irregular steady state in numerical simulations and contrast these predictions with the above theory. Figure 4a,b illustrates the global network coherence (Cor, as defined in section 'Quantification of global network coherence') for a network as a function of gap-junction strength g c and normalized current drive m ext,I /V thI . The dotted line shows the theoretical boundary, matching the very onset of global coherence computed numerically. We find good correspondence between the theoretically predicted stability boundary and numerical simulations, both for the threshold based model in Fig. 4a and the corresponding leaky and integrate-and-fire (EIF) models with reset ( Supplementary Figs 1 and 2). As individual single-neuron response functions and network firing rates are not affected by the presence of voltage reset (see Supplementary  Fig. 1b-d), we expect the oscillatory boundary derived here for the reset-free model to correspond to that of leaky integrate-andfire and models with reset. Indeed, we find that the transition boundary to the global oscillatory state remains unchanged for the adaptive integrate-and-fire model as well as for the aEIF model 30 , which is known to reproduce different neuronal firing patterns 31,32 . Figure 4b shows the corresponding spike rasters in the asynchronous irregular and in the oscillatory state. Note that in the oscillatory state, the neurons fire synchronously and therefore are not in the fluctuation-driven regime anymore. Increasing the excitatory drive m ext,E /V thE however does not lead to selfsustained oscillations (Fig. 4a, bottom), as shown numerically and theoretically. This is due to the assumption that excitatory neurons do not have a subthreshold resonance (a ¼ 0). In summary, we observe that inhibitory current is a potent drive of global self-sustained oscillations, while excitatory drive is less effective. In Fig. 4c, we show that the frequency of these global self-sustained oscillations corresponds closely to the subthreshold resonance as well as the rate preference ( Fig. 1) of inhibitory neurons. The grey region in Fig. 4c delimits the frequency range of emerging self-sustained oscillations. Importantly, we show 1,000 Self-sustained global Firing rate Subthreshold   Table 1). Current input to excitatory neurons is less effective in eliciting global oscillations. (b) Spike rasters of the inhibitory neurons in the oscillatory (top, parameter choice indicated by # in a) and in the asynchronous irregular regime (bottom, parameter choice indicated by * in a). (c) Membrane constant as a function of oscillation frequency: for self-sustained global oscillation (grey area), firing rate resonance (dark grey) and subthreshold resonance (black). In all three conditions, oscillatory frequencies are closely related. NATURE COMMUNICATIONS | DOI: 10.1038/ncomms6512 ARTICLE mathematically that the self-sustained oscillations can occur only if (1) the net effect of the GABAergic synaptic current and gapjunction-mediated spikelet is positive (G II 40) and (2) at least one subpopulation of the neurons expresses subthreshold oscillations (see Methods). We would like to emphasize that these requirements hold for a network transitioning from the fluctuationdriven regime. Of course, oscillations can occur through different mechanisms, such as through coupling of mean-driven neurons that effectively act as oscillators 1,13-15,33 , or through synaptic delays 23,29 .

Discussion
We presented a comprehensive network study of the effects arising from a combination of subthreshold resonance, electrical and chemical synaptic coupling. In our analysis we focused on the fluctuation-driven network state expressed in many cortical regions 26,28,34 . Our results demonstrate that in a situation where neuronal networks are subject to external oscillatory stimuli, the subthreshold frequency preference of a subpopulation defines the global network frequency preference, while electrical and chemical connectivity provide multiplicative gain factors. When an external stimulus of varying frequency impinges on a subpopulation, not only the response amplitude of the stimulated population but also the complete network is closely related to the subthreshold frequency preference of the stimulated population. We found that the electrical synapses are effective at synchronizing the neuronal network if the spike triggered coupling, composed of the GABAergic current and the gapjunction-mediated current, is net positive and a subthreshold resonance is present in at least one subpopulation. The absence of a subthreshold resonance (a ¼ 0) stabilizes the irregular steady state and precludes the generation of a sustained global oscillation, regardless of electrical or chemical coupling. Notably, we find that subthreshold resonant frequency closely corresponds not only to the firing rate preference in independent neurons but also determines the frequencies of the global network oscillations emerging at the instability boundary of the irregular steady state.
Electrical synapses as well as subthreshold oscillations have been previously recognized as key ingredients in the generation of global rhythms 1,10 . Yet, typically, these two properties are studied separately in model settings. Here we briefly recapitulate previous results on the role of electrical synapses as well as subthreshold frequency preference and highlight their relation to our results and model. We start by considering the studies with a focus on neurons with subthreshold frequency preference and those describing the irregular steady state of neural networks. The irregular asynchronous steady state of cortical networks and its oscillatory instability has first been quantified by Brunel and Hakim 23,29 , and Brunel and Wang 35 in a series of landmark papers. Using the Fokker-Planck formalism, they have shown that sparse synchronized fast oscillations can emerge from instabilities of an irregular network dynamics through a Hopf-bifurcation. The emergent oscillation frequency is then inversely proportional to the synaptic delays, which results in a fast rhythm Z200 Hz. This could be lowered by considering synaptic dynamics. Here we considered the same oscillatory instability condition in equation (6) to derive the stability boundary of the irregular steady state, but developed an alternative framework, the adaptive threshold framework, rather than the Fokker-Planck treatment, to explicitly calculate the necessary transfer functions. One of the transfer functions, rate response function, has been calculated in the presence of subthreshold oscillations for the integrate-and-fire model driven by white noise by Richardson et al. 21 A key result of this study is that fluctuation-driven but not mean activity can unmask the subthreshold oscillations in the spike firing patterns of leaky integrator neurons. However, this study considered only isolated neurons rather than an interconnected network and limited its noise statistics to the white noise case. Both of these assumptions could significantly affect the response characteristics as it has been shown for leaky integrate-and-fire neurons 36 . However, extending the Fokker-Planck framework in the presence of coloured noise and subthreshold resonance is analytically hard and has so far not been attempted to the best of our knowledge. Recent studies have started addressing networks of neurons with a subthreshold frequency preference and rate adaptation; however, the analytical treatment has so far relied on time scale separation between t w and t V , low noise activity and numerical simulations 18 . Yet, our study shows that the regime where t w and t V interact is particularly interesting, because it is here that rate and subthreshold resonances are the most pronounced. Despite these differences, we can confirm the observation by Augustin et al. 18 that in the absence of adaptive variable w even excitation-dominated networks lack resonances at any frequency.
Next, we discuss how our results complement previous studies addressing the effect of gap junctions on network dynamics. Ever since gap junctions were first observed in higher brain regions, theorists began to address the similarities and differences between chemical and electrical synaptic communication. Mostly, these studies considered a low-noise regime where all neurons fire regularly and can be described by phase-coupled oscillators [13][14][15] . Counterintuitively, these studies have shown that direct voltage coupling through gap junctions does not necessarily lead to network-wide synchronization but, depending on parameters, could also be a desynchronizing force. However, although the regular firing regime assumed in these studies can occur in the periphery, this is not a typical situation for the cortex where coefficient of variation (CV) and Fano values are close to 1 and fluctuations occur rather than the mean drive spiking 26,37 . This is a situation where global network oscillations could be difficult to achieve and the synchronizing effect of gap junctions could be cancelled. One of the few studies addressing the stability of the irregular spiking regime in the presence of gap junctions is by Ostojic et al. 24 Studying an inhibitory leaky integrate-and-fire network coupled by gap junctions, Ostojic et al. found that network with gap junctions can indeed synchronize in the presence of noise and heterogeneities, and the period of oscillation is determined by the intrinsic firing rate of the neurons. This is in contrast to inhibitory chemically coupled networks where the oscillatory period is determined by the time course of the inhibitory currents following a spike. Importantly, the authors did not consider the subthreshold frequency preference of inhibitory neurons and neglected the synaptic time constants by assuming a white noise drive. Both of these phenomena can significantly alter the transfer rate function R i (o) and thereby affect the stability and oscillatory properties of networks 38 .
Here we proposed an analytical framework to incorporate these phenomena in a network of neurons connected by electrical and chemical synapses, where inhibitory neurons express subthreshold resonance. We show that in the fluctuation-driven regime, inhibitory subthreshold resonance is the substrate for firing resonance of single neurons as well as global oscillations. From a fluctuation-driven regime, self-sustained oscillations only emerge if the inhibitory neurons have subthreshold resonance and are connected by gap junctions, leading to a net spike-induced depolarizating current. Our model therefore indicates that excitatory neurons that typically lack subthreshold resonance are not essential for global rhythm generation but rather amplify self-sustained oscillations. Our theory predicts that this resonance is amplified by chemical connectivity in the network and gap-junction-mediated subthreshold current, but the resonance peak location is not affected by the connectivity. Our framework nicely accounts for experimental evidences 22 showing amplification of g-frequencies by periodic stimulation of the inhibitory neurons and of low frequencies by stimulation of excitatory neurons. Thus, our study provides a very-muchneeded mechanistic understanding of experimental evidence showing the special role of subthreshold resonance for defining the preferred global oscillation frequency of noise-driven networks as well as in defining the response profile of networks to external stimulation.

Methods
We consider a network of N ¼ N E þ N I neurons, where N I is the number of inhibitory and N E is the number of excitatory neurons. Each neuron is modelled by an integrate-and-fire model without reset and an additional adaptation variable, and simulated for a time T sim with a resolution of dt ¼ 0.01 ms. The code used to simulate the model will be posted on ModelDB (https://senselab.med.yale.edu/ modeldb).
Dynamics of a single neuron. We consider a single neuron that exhibits a subthreshold frequency preference and describe the voltage at cell j by the following two dimensional differential equations, t V;j _ V j ðtÞ ¼ À V j ðtÞ þ a j w j ðtÞ þ X j ðtÞ ð 1Þ where t V,j and t w,j are the time constant of the voltage and the adaptation variable respectively, X j (t) is the current, a j and b j are coupling constants. jA{E, I} where E and I represents the excitatory and inhibitory neurons, respectively. The parameters for both populations are summarized in Table 1. We note that this simple leaky integrator dynamics with one adaptation variable w is a linearized representation of a physiologically more accurate conductance-based Hodgkin-Huxley 21 . A spike is emitted by the neuron j, whenever V j (t) crosses the voltage threshold V th,j from below [39][40][41][42] . This reset-free spike implementation is compatible to the classical integrate-and-fire model for low rates and finite time constants 41 , and has been shown to capture essential features of spike correlations and response dynamics of cortical neurons in the noise-driven regime 39,43 . The subthreshold voltage dynamics in equation (1) is inspired by four recent publications [39][40][41]44 . It can be viewed as a minimal model that due to the adaptation variable w i can exhibit subthreshold resonance. Although there are a number of cellular mechanisms other than adaptation that can contribute to a subthreshold frequency preference 45,46 , for reasons of mathematical tractability we opted for the above implementation. The subthreshold resonance emerges in equation (1) if the eigenvalues of the twodimensional transformation are complex, see Fig. 1a. The subthreshold frequency is determined by the imaginary part of the eigenvalues and grows when the membrane time constant t V decreases, see Fig. 1c.
To show the robustness of our results across spiking models, we compare the threshold model with the aEIF model 25 , where the voltage V j at cell j and the adaptation variable w j are described by the following two-dimensional differential equations t V;j _ V j ðtÞ ¼ À V j ðtÞ þ D T expððV j À V th Þ=D T Þ þ a j w j ðtÞ þ X j ðtÞ and We find that all our model predictions can also be observed in the aEIF model and features such as firing rate response, response dynamics and network stability regimes are preserved across models, see Single-neuron response to dynamical stimuli. To address how interconnected networks respond to dynamical stimuli, the first step is to understand how the firing rate of a single neuron responds to weak sinusoidal current drive of frequency f (refs 19,21). We therefore derive the linear response function R j (o) for a neuron jA[E, I] where o ¼ 2pf, V 0,j is the mean voltage of the neuron j, L ¼ exp À s V,j is the variance of V j (derived in equations (15)(16)(17)(18)(19)(20)(21) in the Supplementary Methods). All other parameters as in equation (1) and their values for excitatory and inhibitory neurons are summarized in Table 1. We choose the parameter set for both neurons to match the adaptation and membrane time constants, rate resonance peak and the reliably encoded frequency range that are observed experimentally in cortical neurons 6,19,43 . The reverse characteristic transfer function of a single neuron j is the current response function S j (o) that governs how the average input current of single neuron responds to weak sinusoidal rate drive of frequency f.
Detailed derivation of R j (o) and S j (o), as well as their generalization for electrically connected neurons, can be found in the Supplementary Methods. It it is important to emphasize that the linear response function describes the rate response at the applied stimulus frequency o. Even though this linear response is derived mathematically for small inputs, it is also valid for larger input stimuli. Intuitively, we know that strong stimuli at frequency o excite responses at o and higher-order harmonics at no. As the amplitude grows, the amplitude of no components grows too. However, as we are interested in the response at o, the higher-order harmonics, no, are negligible relative to the linear response component at o.
Gap junctions. Inhibitory interneurons are often connected by both gap junctions and chemical GABAergic synapses 11 . Between the spikes, the gap-junctionmediated coupling, g c , is proportional to the difference of voltages between the neurons, see equation (5). Consider two inhibitory neurons chemically and electrically coupled. When a spike arrives at a GABAergic synapse, a brief inhibitory current pulse (IPSC) is triggered in the postsynaptic neuron (see Fig. 3a, bottom) for a schematic depiction. This is the well-known spike transmission at chemical synapses, which we model here via an exponential IPSC G I y(t À t j ) exp( À t/t I ) at the time of spike t j . On the other hand, the same spike is also transmitted through the gap junction resulting in a so-called 'spikelet' or electrical postsynaptic potential. The spikelet depends on the shape of the presynaptic spike: spikes characterized by a brief depolarization and large hyperpolarization tend to evoke predominantly inhibitory spikelets, while spikes with a prolonged depolarization lead to a primarily excitatory spikelets. A number of studies have investigated the spikelet shape and reported a great diversity 47,48 . For example, parvalbumin-positive interneurons in the somatosensory cortex in L2/3, the gap junction evoked spikelets are purely excitatory 47 , while in the amygdala the spikelet of parvalbumin-positive interneurons can have a late hyperpolarizing transient 49 .
In this paper, we assume that the spikelet contribution is excitatory, although our framework in not limited to that. Mathematically, we model the spikelet as an exponential at the time of spike t j with the same time constant as the IPSC. We can therefore sum the contribution of the GABAergic IPSC and the spikelet, leading to a net exponential postsynaptic potential G II y(t À t j ) exp( À t/t I ). This single exponential approximation serves to keep the analytical complexity at bay while considering the dominant finite time scale of synaptic interactions.
Dynamics of the network. Interested in the temporal spike dynamics of a network with electrically and chemically coupled inhibitory neurons and only chemically coupled excitatory neurons, we model the input current X j (t) to a neuron in the irregular steady state as a sum of chemical and electrical coupling. The topology of the network is all-to-all with the connectivity matrix G ij . The all-to-all connectivity is not critical and we observed equivalent results with sparse connectivity, as long as the product between the connectivity probability and the connectivity strength ARTICLE remains the same. Note that the inhibitory spike-mediated coupling G II ¼ G E,gap þ G I consists of a spikelet-mediated excitatory contribution (G E,gap ) and the contribution of GABAergic chemical synapses (G I ). We represent each spike triggered electrical and chemical synaptic potential by an exponential of the form G ij y(t) exp( À t/t I ) such that the external fluctuations Z ext (t) are of the Ornstein Uhlenbeck type with time constant t I (ref. 50). We neglect synaptic delays and assume that both chemical and electrical spike interactions are instantaneous. We operate in the irregular steady state where the mean recurrent current drive is proportional to the firing rate 34 . For the sake of tractability, we assume that the fluctuations in the external input s I dominate over network-generated noise and are sufficient to describe the firing rate of the network (see Fig. 1d for a comparison between theory and numerical simulations). Explicitly, we model the input current X j (t) to a neuron in the irregular steady state as X j ðtÞ ¼ g c hV I;k ðtÞi k;2inh þ G jI tn I ðtÞ þ G jE tn E ðtÞ þ m ext þ s I Z ext ðtÞ: Here, /V I,k (t)S k,Ainh is the average voltage of the inhibitory population, v E (t) is the firing rate of the excitatory population, v I (t) is the firing rate of the inhibitory population, m ext is the constant external drive and sm ext (t) the fluctuating external input. g c is the subthreshold coupling contributed by gap junctions. The strength of excitatory and inhibitory spike triggered interactions for the population j are G jE and G jI , respectively. The parameters are given in Table 1.
Steady-state firing rate. First, we study the firing rate as a function of current drive in a single, isolated neuron, driven by a mean current X 0 ¼ m ext þ s I Z ext where Z ext is an Ornstein Uhlenbeck fluctuation with correlation time t I and unit variance. The firing rate v i is given in the Supplementary Methods and is shown on Fig. 1d. The theoretically derived firing rate v i as a function of m ext /V th is in good agreement with numerical simulations. In a mixed irregular, asynchronous network with both excitation and inhibition, the steady-state firing rate of inhibitory neurons v I and excitatory neurons v E are derived in the Supplementary Methods and shown in the Supplementary Fig. 1c,d, again matching numerical simulations.
Condition for self-sustained oscillations. To assess the stability of the asynchronous irregular network state with respect to oscillatory perturbations, we consider the stability conditions outlined in refs 23,29. Starting with an infinitesimal time-dependent perturbation d I (o) exp(iot N ) in current applied at time t N , we obtain changes in firing rate and resulting changes in current induced by as change in firing rate. For a network consisting of only inhibitory neurons, the stability condition reads G II R I ðoÞS I ðoÞ ¼ 1; where G II is the recurrent connectivity strength, and R I and S I are response functions given in equation (3) and equation (4), respectively. For the general twodimensional network with an excitatory and inhibitory subpopulations, the stability condition takes a matrix form and is provided in the Supplementary Methods. Studying equation (6) and its real part, we find that for all G II o0 no real o solution can be obtained implying that the irregular regime is stable for inhibitory interactions. Another condition for the stability of the irregular steady state is obtained if subthreshold resonance is lacking (a ¼ 0). In this case we find Solving for o we obtain only imaginary solutions implying that the oscillatory instabilities can not emerge and the irregular steady state is stable for any G II under these conditions. Although this condition is derived for the threshold neuron model, we show in the Supplementary Fig. 2 that it also holds for the adaptive leaky and aEIF neurons with reset.
Quantification of global network coherence. In simulations, we assess the global network coherence using the rate-normalized correlation r. For each excitatory (E) and inhibitory (I) population, we compute the correlation according to r I;E ¼ ðhs s I;E ðtÞ Á s s I;E ðtÞi=n I;E À n I;E Þ=n I;E . Here, / Á S denotes the ensemble average, v I,E are the rates of excitatory and inhibitory neurons, respectively. s I,E are the compound spike trains of the respective populations. To compute the spike correlation hs s I;E ðtÞs s I;E ðtÞi, we convolve each spike with a Gaussian kernel with s ¼ 3 ms, leading to s s I;E ðtÞ. If the network is in an asynchronous, irregular state, then we expect r ¼ 0. In contrast, if the network is in an oscillatory state with each spike aligned across the population then lim s-0 r ¼ N. We therefore use r I,E as an index of population synchrony. In a mixed excitation-inhibition network we use the average Cor ¼ (r E þ r I )/2 as a measure of network coherence 51 , for example in Fig. 4a. Note that in practice in the asynchronous irregular regime r will fluctuate around zero, with the amplitude of fluctuations decreasing for growing network size N.