Emergent bimodal firing patterns implement different encoding strategies during gamma-band oscillations

Upon sensory stimulation, primary cortical areas readily engage in narrow-band rhythmic activity between 30 and 90 Hz, the so-called gamma oscillations. Here we show that, when embedded in a balanced network, type-I excitable neurons entrained to the collective rhythm show a discontinuity in their firing-rates between a slow and a fast spiking mode. This jump in the spiking frequencies is characteristic to type II neurons, but is not present in the frequency-current curve (f-I curve) of isolated type I neurons. Therefore, this rate bimodality arises as an emerging network property in type I population models. We have studied the mechanisms underlying the generation of these two firing modes, in order to reproduce the spiking activity of in vivo cortical recordings, which is known to be highly irregular and sparse. We have also analyzed the relation between afferent inputs and the single unit activity, and between the latter and the local field potential (LFP) phase, in order to establish how the collective dynamics modulates the spiking activity of the individual neurons. Our results reveal that the inhibitory-excitatory balance allows two encoding mechanisms, for input rate variations and LFP phase, to coexist within the network.


INTRODUCTION
The type and distribution of ionic channels across the membrane of a neuron determine its excitability behavior. This behavior can be experimentally tested by injecting a pulsed, constant, or ramp current. In the latter case one observes how the membrane potential changes as the input amplitude increases, eventually switching from a resting to an oscillating (tonic firing) state. Mathematically, this transition is described as a bifurcation of the membrane potential from a stable fixed point to a periodic orbit.
Depending on how the firing frequency within the oscillatory regime behaves near the bifurcation, neural excitability can be classified into type I or type II (Rinzel and Ermentrout, 1989). In type I excitability, the frequency increases continuously from zero as the bifurcation is crossed. Type II neurons, on the other hand, exhibit a discontinuous jump in frequency as the tonic regime is entered, and their range of firing frequencies is quite narrow compared to type I neurons, which can achieve arbitrarily low frequencies (Izhikevich, 2007). These different behaviors are associated with distinct bifurcations. In type I excitability, the periodic orbit emerges due to the collision of a stable fixed point (the resting state) and an unstable equilibrium point, occurring on top of an invariant circle (SNIC bifurcation). Type II excitability, on the other hand, can arise in three different ways: via a subcritical or a supercritical Hopf bifurcation, or through a saddle-node bifurcation outside the invariant circle (Rué and Garcia-Ojalvo, 2011). The integrate-and-fire and conductance-based models used in the literature to describe cortical networks are usually of type I (Wang and Buzsáki, 1996;El Boustani et al., 2007).
The phase response can also be used as a criterion to distinguish between excitability classes: when operating in a tonic regime, type I neurons always advance their phase (defined with respect to their spiking period) when perturbed by a brief depolarizing pulse, while type II neurons can either advance or delay their phase depending on the instant of perturbation relative to the period of oscillation (Hansel et al., 1995). Within a network, the input pulses come from other neurons, and therefore synchronized periodic spiking might lead to lower or higher frequency rhythms depending on the type of neuronal excitability. However, high spike synchrony and periodic firing patterns are rarely seen in vivo; in fact the coefficient of variation of interspike intervals in the cortical neurons is typically larger than 1 (Shinomoto et al., 1999(Shinomoto et al., , 2005. Consequently, the rhythm of the neural population must emerge from the recurrent synaptic activity in the network, rather than from spike-to-spike synchrony. It has been experimentally observed (Softky and Koch, 1993) that the high synaptic bombardment acting upon cortical neurons is far from being constant, leading to stochastic fluctuations that affect these neurons in vivo. Thus, neither the single-neuron frequency-current curve (f -I curve) nor the phase response can perfectly describe the behavior of single neurons when they are embedded in a network.
Cortical oscillations thus arise as a collective phenomenon, which does not require individual-neuron firings to be oscillatory themselves. These global oscillations are discernible in averaged population activities such as the local field potential (LFP), whose troughs (or peaks, depending on the relative position of the recording electrode and the generating current dipoles) correspond to the minima of the synaptic inhibitory flow. In these temporal windows neurons are more likely to spike, producing an increase in the excitatory synaptic current, followed by a new burst of inhibition. These oscillations, at frequencies in the beta (12-30 Hz) and gamma (30-90 Hz) ranges, are experimentally seen in the cortex upon sensory stimulation (Buzsáki and Draguhn, 2004). Specifically, gamma-band synchronization has received special consideration, as it is hypothesized to be a mechanism for the dynamic generation of functional cell assemblies and for the flexible communication between brain areas [for a review, see Singer (1999) and Fries (2009)].
Here we study the generation of gamma-band (30-90 Hz) oscillations that occur without spike-to-spike synchrony [known as synchronous irregular dynamics (van Vreeswijk and Sompolinsky, 1996;Brunel, 2000;Hansel and Mato, 2003)], and how these rhythms affect single unit activity. Specifically, we consider networks of type I excitatory and inhibitory neurons described by conductance-based models and stimulated by slowly varying excitatory inputs, and study the emergent firing patterns of the neurons once embedded in a balanced network. Our results show that the synchronous irregular state giving rise to gamma-band oscillations is composed of two firing modes: a high-frequency (or bursty) regime and a low-frequency regime separated by a gap of quasi-forbidden firing frequencies. We also study the relationship between this firing behavior and the input: the fast spiking mode encodes for the stimulus rate, whereas the existence of a slow spiking mode allows for a phase code to appear.

DESCRIPTION OF THE NEURONAL MODELS
The dynamical equation for the neuronal membrane voltage is given by a conductance-based model: where g K , g Na , and g L are the maximum conductances for the potassium, sodium and the leak current, respectively, and I syn is the synaptic current coming from the neighboring neurons impinging on one neuronal cell. The dynamics of the potassium and sodium channels is represented by the time-varying probabilities that a channel is open: where x stands for n in the case of the potassium current, and for m and h in the case of the sodium current. α(V) and β(V) are voltage-dependent rate constants, and φ is the temperature factor, defined by φ = 3 (T−6.3)/10 , where T is measured in degrees Celsius.
The parameter values used throughout this study are those of Ref. Gutfreund et al. (1995): g K = 4.74 μS, g Na = 12.5 μS, and g L = 0.025 μS. The reversal potentials of the different channels are V K = −80 mV, V Na = 40 mV, and V L = −65 mV, and the membrane capacitance is C m = 0.25 nF (0.125 nF) for the excitatory (inhibitory)neurons. The leak conductance defines an effective membrane time constant for the isolated neuron according to the expression τ = C m /g L , which is taken to be 10 ms for the excitatory neurons and 5 ms for the inhibitory neurons. The temperature factor φ is set to 21, which corresponds to T = 34 o C. These parameter values lead to type I excitability.
The rate functions α and β for each gating variable are: for the gating variable n, for the gating variable m, and for the gating variable h. Due to the rapid activation of m we replace it by its steady-state value m ∞ = α m α m +β m .

DESCRIPTION OF THE NETWORK MODEL
We consider a network composed of 2000 neurons, 80% of which are excitatory while the remaining 20% are inhibitory (Soriano et al., 2008). All connections between cells are chemical synapses-no gap junctions are considered-and each neuron connects randomly with 200 other neurons, on average, belonging to both populations. Therefore, no architecture is imposed on the connectivity. We have also introduced a synaptic delay in the transmission of the action potential between neighboring neurons, taken from a gamma distribution of mean 2 ms and variance 4 ms 2 . The synaptic current is described using again a conductance-based formalism: where g syn is the synaptic conductance and E syn is the reversal potential of the synapse. For E syn greater than the resting potential V rest the synapse is depolarizing, i.e., excitatory, otherwise it is hyperpolarizing, i.e., inhibitory. We consider two temporal time constants, τ d and τ r (decay and rise synaptic time, respectively, The equilibrium potential at I syn = 0 for Equation (1) is V rest = −65 mV.
where g syn , shown in Table 2, is tuned in order to obtain a balance between excitation and inhibition, given the f -I relation.
We have chosen the maximal conductances, g syn , to maintain the postsynaptic potential (PSP) amplitudes within physiological ranges: the excitatory PSP in the range from 0.42 to 0.83 mV, and the inhibitory PSP from 1.54 to 1.88 mV.
Additionally, all neurons receive an heterogeneous Poisson train of excitatory presynaptic potentials with a mean event rate that varies following an Ornstein-Uhlenbeck process. This incoming external current mimics the direct input from any other area external to the network considered here. The instantaneous rate, λ(t), of the external excitatory train of spikes is generated according to an Ornstein-Uhlenbeck process as considered in Mazzoni et al. (2008): where σ(t) is the standard deviation of the noisy process and is set to 0.6 spikes/s. τ is set to 16 ms, leading to a 1/f power spectrum for the λ time series that is flat up to a cut-off frequency f = 1 2πτ = 9.9 Hz. η(t) is a Gaussian white noise. Excitatory synapses outnumber inhibitory ones, and yet the brain avoids epileptic states because inhibition is able to balance excitation, and thus neurons remain below threshold, firing only occasionally. Inhibitory neurons have higher firing-rates than excitatory neurons for any given input current, as can be seen in the f -I curve of Figure 1. Additionally, GABAergic (inhibitory) synapses are stronger than glutamatergic AMPA (excitatory) synapses to compensate their relative small number (Markram et al., 2004).

NUMERICAL SIMULATIONS
The model has been integrated using the Heun algorithm (Garcia-Ojalvo and Sancho, 1999), with a time step of 0.05 ms. All the simulations represent 3 s of activity and the connectivity, initial conditions, and noise realization were varied from trial to trial.

MODEL OF LFP
We quantify the activity of the network in different ways. At the single-neuron level we consider the instantaneous firing-rate as a measure of the individual spiking dynamics. At the population level we use two observables: the time-resolved average firing-rate of the whole neuronal population (defined as the total number of spikes per unit time in the population divided by the number of neurons), and the LFP, computed as the sum of the absolute values of the excitatory and inhibitory synaptic currents acting upon the excitatory neurons, averaged over this population (Mazzoni et al., 2008;Buzsáki et al., 2012): Here · · · denotes an average over all excitatory neurons (Berens et al., 2010). The term I AMPA accounts for both the external excitatory heterogeneous Poisson spike train and the recurrent excitatory synaptic current due to network connectivity, while I GABA corresponds to the recurrent inhibitory synaptic current. R e represents the resistance of a typical electrode used for extracellular measurements, here chosen to be 1 M . In Figure 9 the LFP was filtered with a 4th order Butterworth bandpass filter using MATLAB function filter.m.

SPIKE TRIGGERED AVERAGE
We have calculated the spike triggered average (STA) of the LFP and of the inhibitory synaptic current impinging on the neurons. For the considered spikes we have registered these signals during a time window starting 50 ms prior to the spike and ending 20 ms after it, and computed the mean across the total number of action potentials.

COMPUTATION OF POWER SPECTRA
The LFP power spectrum was estimated using the multitaper method (Thomson, 1982) commonly used to reduce the variance of the spectra of recorded signals, which are usually very noisy. This estimator was implemented in Chronux 2.10 (Bokil et al., 2010). The multitapered power spectrum, S(f ), is the average of the power spectrum of the LFP signal multiplied by K orthogonal Slepian functions (in our case K = 5), and further averaged over N = 20 trials: Here LFP n,k (f ) is the discrete Fourier transform of the LFP(t) signal of the n-th trial, multiplied by the k-th Slepian function (or taper). We have considered data segments within a 500-ms sliding time window with an overlap of 50 ms, padded with zeros up to a length of 512 in order to obtain an increased sampling rate in the frequency domain. The resolution bandwidth is thus ±6 Hz. The firing-rate power spectra are also obtained by the multitaper algorithm, with the same sliding time window, overlap, and padding. In this case the average instantaneous firing-rate is obtained from an histogram of the spiking times, with a 1-ms bin. All histograms and power spectra are averaged over 20 trials.

RESULTS
An adequate neural network model should reproduce several experimentally observed features of in vivo cortical oscillations: a prominent peak of the LFP power spectrum in the gamma-band resulting from external stimulation, irregular individual firing at frequencies lower than the gamma rhythm, and partial phase locking of individual spikes to the gamma cycle. We show in what follows that the neural network described above is able to reproduce these features.

FIRING-RATE OF NETWORK-EMBEDDED TYPE I NEURONS
In this section we study how the network shapes the firing patterns of the single neuron components. Before studying the firing modes of the cells when embedded in a network, we first characterize their behavior as single cells. When isolated, a neuron displays a characteristic transition to the tonic firing regime as the input current is gradually increased (solid lines in Figure 1). Type I integrator cells can achieve arbitrary low spiking frequencies.
This behavior is a consequence of the saddle-node bifurcation on an invariant circle that this type of neuron undergoes with increasing injected current. For the same current, inhibitory cells (with a smaller membrane time constant) fire with a higher rate. However, both excitatory and inhibitory cells have the same spiking threshold. The behavior of the firing-rate under constant injection current, shown in Figure 1A for adiabatically increasing values of the current, is not suited for characterizing the behavior of a neuron within a network. When embedded in a population, neurons are subject to synaptic input that fluctuates in time. Even when the firing of the presynaptic neurons projecting to a given postsynaptic cell is uncorrelated between them, the total input current is Gaussian distributed (Roxin et al., 2011). This input is the sum of the synaptic currents coming from all possible sources, which include recurrent excitation and inhibition (generated by the network itself), and afferent pathways terminating onto the population (represented in the model by a Poisson-distributed external train of spikes). The net synaptic current impinging on a neuron has a nearly stationary average value and fluctuates rhythmically around it (due to the network-induced gamma oscillations), crossing randomly the threshold and leading to an irregular individual firing-rate. Here, we capitalize on the natural variability of such recurrent and external synaptic drive to explore the effective f -I curve of neurons once embedded in a balanced recurrent network. Figure 1B (circles) represents the f -I response curve of type I neurons when embedded within a population. Given the above-mentioned fluctuating character of the synaptic current received by the neurons within the network, the f -I response is quantified here in terms of the interval between two consecutive spikes (interspike interval, or ISI) as a function of the mean synaptic current received by the neuron within that interval. This quantity is plotted in Figure 1B (circles), comparing its behavior with the corresponding f -I curves of the isolated neurons (solid lines). Note, however, the limitations of this comparison, as the mean synaptic current over an ISI differs from the actual fluctuating input received by a cell, specially for long periods. It is precisely in the range of low frequencies where both curves differ most. For instance, the firing-rate in the resting state region departs from zero close to the bifurcation, due to the appearance of noise-induced spikes throughout the network.
Differences between the firing-rate of a network-embedded neuron and its characteristic f -I curve in isolation also exist in the tonic regime, as can be seen in Figure 1B. We can understand this  (Ermentrout, 2002). (B) Equivalent f -I curve for neurons embedded in the network (circles) representing the instantaneous firing-rate versus the net synaptic current averaged over the corresponding inter-spike interval, for an excitatory (black) and inhibitory (red) neuron. In that panel the single-neuron f -I curves are also shown in solid lines for comparison. The network is excited by a train of spikes of rate 8500 spikes/s. from the fact that, due to the synaptic bombardment, neurons within a network are persistently perturbed, and consequently their firing-rates are altered. The characteristic phase response curve (PRC) of type I neurons (Figure 2) reveals a phase advance for any perturbation time within a period. Thus, when a type I neuron has just spiked and is excited again by a presynaptic potential its period is reduced (and thus its firing-rate is increased). This explains why the firing-rate of type I neurons embedded in a network, and thus continually perturbed by presynaptic activation, is above the f -I curve of an isolated neuron with tonic firing.

Frontiers in
It is worth noting that the firing-rate discontinuity shown in Figure 1B (circles) is an emergent network property. Within a network, type I neurons display two clearly defined groups of inter-spike intervals, in a sort of population-level bimodality that does not exist for individual neurons in isolation.

NETWORK OSCILLATIONS
When the induced spiking activity is large enough to percolate the network, the internal synaptic currents exhibit a fast cyclic behavior, alternating epochs of high excitation followed by high inhibition, due to the rhythmic synchronization of GABAergic cells via their recurrent connections. The net oscillatory synaptic current leads to a rhythmic behavior in averaged population measures such as the LFP and the time-resolved firing-rate of the population, as shown in Figures 3A,D, respectively. The LFP power spectrum is depicted by solid lines in Figure 3B, revealing a frequency peak whose precise position is determined by the GABAergic synaptic time constants and the synaptic strengths, as excitatory inhibitory FIGURE 2 | Phase response curve for a type I excitatory and inhibitory neuron representing the phase change in the spiking response of the neuron to a depolarizing pulse, as a function of the phase of the input (normalized to the spiking period T 0 of the neuron). The excitatory (inhibitory) neuron is spiking tonically at an unperturbed period of 8.09 ms (6.00 ms), and receives the injection of a pulse with amplitude 1.0 nA and duration 0.2 ms. The PRCs, obtained with the XPPAUT software (Ermentrout, 2002), are defined as 1 − T (τ)/T o , where T (τ) is the period after perturbing at time τ.
well as by the characteristics of the input (Henrie and Shapley, 2005). For a range of biologically plausible values, the frequency of such peak is in the gamma range (30-90 Hz). In particular, for the parameters given in Tables 1, 2 (see section 2), the frequency peak of the LFP spectrum is located around 45 Hz for an external stimulation of mean rate 8500 spikes/s. The low-frequency components are caused by the external incoming spike trains, whose mean rate is an Ornstein-Uhlenbeck process. When the effective coupling between neurons is low, because either the induced spiking activity or the synaptic conductances are weak, the collective rhythm disappears from the LFP power spectrum, as shown by the dashed line in Figure 3B.
Similarly to the LFP, the global firing-rate exhibits a marked peak at ∼45 Hz in its spectrum ( Figure 3E). This rhythmicity reflects epochs of synchronization between subsets of neurons. Despite the collective rhythmic dynamics (Figure 3D), single neurons display strongly irregular individual spiking (Figure 3C), firing mostly in a sparse and single-spike mode occasionally accompanied by some high-frequency tonic firing, which nevertheless can be compatible with population rhythmic activity at gamma frequencies.

FIRING-RATE DISTRIBUTIONS
As described in the last sections there is a striking difference in the spiking activity exhibited by the type I cells once embedded in a network compared to their intrinsic f -I curve. To better resolve the firing mode differences and their origin, we first characterize the distribution of instantaneous firing-rates, and then we describe the relation of firing behavior to the input and LFP dynamics.
To characterize the distribution of instantaneous firing-rates of individual neurons within the network described in section 2.2, we compute their histogram for the parameter values given in that section, which lead to synchronous irregular firing (Figures 3C,D) with global oscillations in the gamma range (Figures 3B,E). The distribution of instantaneous firing-rates is bimodal (Figure 4A), indicating the prevalence of two firing modes, with high-firing-rate events (short ISIs) reflecting bursts of spiking activity and short-firing-rate events (long ISIs). This effect arises from the network dynamics and cannot be predicted from the characteristic f -I curves of type I isolated neurons, which are continuous and hence do not forbid any particular range of firing-rates (Figure 1B; note the gap in the circles, which occurs in network-embedded neurons). This bimodal response was also observed in a network of LIF neurons (Roxin et al., 2011), when the inhibitory population firing-rate exceeded the excitatory population rate. We will address the origin of this bimodal behavior at the end of this section.
The high-firing-rate peak in Figure 4A arises from the spikes of the network-induced bursts (see upper cloud in Figure 1B). The low-firing-rate events are more frequent than the highfiring-rate ones (64% of the interspike intervals of the excitatory population across 20 trials lie below ∼58.31 spikes/s-see arrows in Figure 4A). In other words, more neurons are simultaneously found in a silent state or low-firing mode.
In order to understand the genesis of this bimodal firing behavior, we now examine the time evolution of the synaptic  currents acting upon the excitatory population (which by definition determines the dynamics of the LFP, see section 2.4). Figure 4B shows the time traces of the external excitatory input (brown lines) and the net synaptic current (green lines). This latter current accounts for both the excitatory component (arising from the external and the recurrent excitatory spikes) and the inhibitory component (arising from the recurrent inhibitory spikes). The thick lines are population averages over the excitatory neurons, whereas the thin lines represent the current values impinging locally on a given neuron. As shown in this figure, the average external synaptic input (thick brown line) lies well above the spiking threshold (horizontal dashed line), whereas the average net synaptic current (thick green line) is below the spiking threshold, due to the inhibitory flow that counteracts the strong external excitation. Therefore, the current impinging on a neuron [thin green line, corresponding to the I syn term of Equation (1)] is typically below threshold. Thus neurons spike rarely, only when the excitation-inhibition balance is lost during a certain time window, in which the external input brings the neuron above threshold (see asterisks in Figure 4B). From time to time these intervals are long enough for several spikes to occur in quick succession, giving rise to periods of high-firing-rate (for instance around ∼1325 ms). On the other hand, if the excitationinhibition balance is briefly lost, an isolated spike is elicited only if inhibition is low enough (as we will see in detail in section 3.5), giving rise to low-firing-rate events (for instance at ∼1230 ms). Therefore the combined dynamics of excitation and inhibition is the basis of the bimodal distribution of instantaneous firing-rates. According to the previous discussion, the ratio of fast to slow firing events is determined by the characteristic timescale of the periods in which the excitation-inhibition balance is lost. In order to verify this reasoning, we now analyze the behavior of the network for a longer value of the decay time constant τ d of the GABAergic synapses. Increasing τ d to 30 ms leads to a disappearance of the gamma rhythm (compare the thick green lines in Figures 4B, 5B) (Fisahn et al., 1998;Heistek et al., 2010).
Since the inhibitory currents are slower, the periods in which the excitation-inhibition balance is lost are longer, giving rise to an increase in the number of fast spiking events (see asterisks in Figure 5B). This is quantified in Figure 5A, which shows that high firing-rates are much more frequent than low firing ones. Thus the inhibition decay time, τ d , determines the principal individual firing-rate mode (mostly single-spike for fast inhibition and bursty for slow inhibition).

CODING IMPLICATIONS OF THE FIRING-RATE BIMODALITY
The type I network oscillating in the gamma range, Figure 4A, clearly has a sparse activity, in agreement with the sparse coding of sensory inputs (Olshausen and Field, 2004;Wolfe et al., 2010). We now ask how the bimodality between slow and fast spiking regimes, with a band of forbidden firing-rates, affects the coding capabilities of the neuronal network. To address this question, we establish how the individual firing-rate depends on the external input to the population, rather than on the internal synaptic current (as shown in Figure 1 above). In our case, the external input is modeled by a set of spike trains perturbing each neuron, mimicking either an external sensory input representing a stimulus or neuronal activity arriving from other areas. As mentioned above, the time course of the instantaneous rate is an Ornstein-Uhlenbeck process, equal for all the spike trains (but the specific realization of the Poisson process is different for each neuron). Its mean value is a measure of the input intensity. Figure 6 shows the instantaneous individual firing-rate of the excitatory neurons composing the network as a function of the mean external firing-rate, averaged over the corresponding ISI for which the instantaneous firing-rate is calculated, for three different time scales of the GABAergic synaptic dynamics, all of them generating LFP oscillations in the gamma range. All three plots clearly reveal the bimodal character of the firing-rate distribution described above. Two clouds of data points are clearly discernible, corresponding to distinct ISIs. The cloud at the top exhibits a clear correlation with the external rate, which provides for a standard mechanism of rate coding: the higher the external rate is, the faster the resulting firing-rate of the neurons in the network. The lower cloud, on the other hand, is associated with longer firing periods. Naturally, the instantaneous external rate of that low-firing-rate state converges to the mean value of the external rate (marked by a vertical dashed line in the three plots), due to the large periods over which that mean is calculated. Accordingly, for instantaneous external rates clearly above the mean external rate the network only responds with high-frequency firing, following the f -I curve of the individual neurons (i.e., operating in a rate coding mode). For moderate external rates, on the other hand, the distribution of firing-rates is bimodal, with low-firing and high-firing events coexisting for the same value of the external rate. This allows the network to encode for the phase of the global oscillations, since in each one of the two modes the locking of the individual spikes to the gamma rhythm is different, as described in section 3.5 below.
The scatter plots shown in Figure 6 also allow us to further investigate the origin of the firing-rate gap evident in those plots (and in Figure 1B above). We can anticipate that the basis of that gap is the fact that due to the recovery time of the network following an increase of inhibition, firing-rates that match the duration of network depression are forbidden. To verify that expectation, we compare the behavior of the network for different values of the decay time constant τ d of the GABAergic synapses, while keeping the IPSP amplitude and AMPA synapses constant. Our results show that longer inhibition [growing from panel (A) to panel (C) in Figure 6] increases the range of forbidden firing periods, because the probability of spiking after an inhibitory barrage is lower during a longer time interval. Accordingly, as τ d increases the low-firing mode reduces its area because the slower inhibitory synapses forbid the shortest ISIs within this mode. Therefore, the low-firing state is the result of a competition between the external slow fluctuations and the recurrent inhibition, whose oscillating frequency decreases with increasing GABA synaptic decay time. This leads to the counterintuitive discontinuity in the firing-rate of type I neurons described above. In summary, neurons behave like all-or-none detectors of rapid stimulus fluctuations, with instantaneous firing-rates faster than the LFP gamma peak encoding for the fast dynamics of the stimulus.

PHASE LOCKING TO THE GAMMA CYCLE
We have seen that, even though the firing-rate of the individual neurons is far from being tonic, the probability of firing across the network varies rhythmically in time with a frequency around 45 Hz. This rhythm is generated by recurrent excitatory and inhibitory connections, and is revealed in the LFP and firingrate dynamics. Both the synaptic flow (apparent in the LFP) and the spiking activity are mutually interacting, given that a decrease in synaptic inhibition triggers an increase in the mean firing activity of the network. This causal relationship implies that the peaks of the population firing activity and the troughs of the LFP are displaced in time within the millisecond range, controlled by the synaptic delay.
To further characterize the effect of the global rhythm on the firing activity, we have computed STA across trials of both the LFP and the inhibitory synaptic current, I GABA , acting upon each neuron. The LFP carries information about the global synaptic activity affecting the excitatory neurons of the population (see Equation 2), and is therefore a measure of the global activity of the network. In contrast, I GABA only accounts for the inhibitory synaptic current impinging on an neuron from the firing activity of its presynaptic inhibitory neurons, and is thus a local measure. We have considered the STA in a time window of 70 ms around a spike (50 ms pre-spike and 20 ms post-spike). Furthermore, to avoid having a previous spike from the same neuron fall within the time window being considered, we have only taken into account spikes at least 50 ms apart. When computing the STA we have classified the spikes, according to the bimodal distribution of instantaneous firing-rates shown in Figure 4A, as slow or fast firing events, the latter corresponding to network-induced bursts. Within the fast firing mode we selected the first spike of each burst, since we are interested in the events leading to burst initiation. We now discuss separately the mechanisms underlying the generation of the slow and fast firing modes and their relationship with the phase of the global oscillations.

Fast firing mode
As shown by the thin lines in Figure 7A, the high-frequency regime appears when the inhibitory synaptic current I GABA impinging on a neuron (i.e., its local inhibition) is close to zero for a relatively large time interval right after t = 0. In other words, the presynaptic inhibitory neurons must be silent for a long enough period of time after the neuron fires, in order for a burst to be initiated. In those conditions, it is more probable that the excitatory external input brings the neuron above the spiking threshold for a sustained amount of time, giving rise to a fast spiking period.
For the excitatory neurons this occurs preferentially at a high level of global inhibition, i.e., at the maxima of the LFP (thin line in the top panel of Figure 7B). These two events, a maximum of the global inhibition (i.e., of the LFP) and a minimum of the local inhibition (i.e., of |I GABA |), occur simultaneously, since it is when global inhibition is strong that some neurons can be locally surrounded by strongly inactivated inhibitory neurons. During this time interval, those neurons fire in bursts due to the external excitation, before all the inhibitory neurons become excited and prevent the firing of the rest of the population. A scheme of the interplay between local inhibition and global activation leading to the fast firing mode is shown in Figure 8B. In contrast with the excitatory neurons, in inhibitory neurons the initiation of the fast firing mode does not depend so strongly on the level of global inhibition (thin line in the bottom panel of Figure 7B). We hypothesize that these neurons, which have a smaller membrane time constant than the excitatory ones (see section 2), react rapidly to the external excitatory fluctuations and spike regardless of the level of global inhibition, provided their presynaptic inhibitory current is zero (thin line in the bottom panel of Figure 7A). Vertical thick ticks represent spikes. The black and red circles represent excitatory and inhibitory neurons, respectively. For simplicity, only contacts to the gray neuron are represented (black arrows). In the top panel of (A) (slow firing mode), the reference neuron (in gray) is surrounded by presynaptic inhibitory neurons that are starting to fire (i.e., the level of I GABA , represented by a horizontal red bar, is moderate). The global activity of the network (bottom panel, represented by a white horizontal bar) is low and the earliest inhibitory spikes prevent the cells from bursting. In the top panel of (B) (fast firing mode) the neuron is surrounded by inactive inhibitory neurons (i.e., the level of I GABA is low), whose silence is driven by the high activity of the rest of the inhibitory population (reflected in a high level of the LFP).

Slow firing mode
The low-frequency regime, on the other hand, takes place when the local inhibitory current, I GABA , impinging on the neuron is momentarily zero, as shown by the thick line in Figure 7A. The rapid increase of inhibition following a spike (at t = 0, dashed line) prevents the neuron from spiking again. This happens at low levels of global inhibition, i.e., close to the minima of the LFP (thick line in Figure 7B). Immediately after a spike the firing activity of the population increases, and thus the local inhibition grows as well, increasing the Frontiers in Computational Neuroscience www.frontiersin.org distance of individual neurons to threshold (see the scheme of Figure 8A).
The STA analysis discussed above shows that the simultaneous occurrence of slow instantaneous rates at both the excitatory and inhibitory populations (thick lines in Figure 7B) is only possible near the minima of the LFP, whereas fast instantaneous rates can only be simultaneously present in both populations at the maxima of the LFP (thin lines in Figure 7B). Given this interaction between the individual firing modes and the rhythmic network dynamics, one can ask whether a partial representation of the input is coded in the spiking timing of cells relative to the phase of the population oscillation. We now address this question by proceeding to filter the LFP signal around the gamma frequency peak (46.88 ± 5 Hz) and assign an instantaneous phase to the LFP time series via the analytical signal approach (Hilbert phase) (Le Van Quyen et al., 2001). A histogram of the LFP phase values at which the neurons spike is shown in Figure 9, with π corresponding to the LFP troughs (vertical dashed line in the figure).
Spike firings are widely distributed over the 2π cycle of the LFP phase because, as seen from the raster plot of Figure 3C, the firing activity is noisy and the correlation between neurons is weak, due to their low coupling and the absence of a common drive. We have performed the spike-LFP phase locking analysis for both high-firing (thin lines) and low-firing (thick lines) events, taking only into account the first spikes of short and long ISIs, respectively. Figure 9 shows that the slow firing regime of the inhibitory neurons (red thick line) is more sensitive to the LFP phase than the fast firing regime, resulting in a more pronounced locking to the troughs of the LFP, with a delay ≥2 ms due to the mean synaptic delay. The high-firing mode of the excitatory neurons, on the other hand, is phase-locked with no delay to the LFP peaks 0 π/2 π 3π/2 0 500 1000 phase (rad) Number of spikes excitatory inhibitory FIGURE 9 | Distribution of LFP phase values at which neurons spike. The plot distinguishes between spikes flanking ISIs with fast (thin lines) and slow (thick lines) firing modes, for both excitatory (black) and inhibitory (red) neurons. Measures were performed over the same data used in Figure 7. The dashed line marks the π phase that corresponds to the LFP minimum.
around 0 (note that no synaptic delay is considered in the external train of spikes, which as mentioned above controls the high-firing regime).
In section 3.4 we have shown that the fast fluctuations of the excitatory external rate induce high instantaneous firing-rates, and thus this mode encodes the fast dynamics of the input. In this section we have additionally shown that the occurrence of spikes relative to the LFP phase depends on the firing mode. In particular, the high-frequency mode appears with higher probability when inhibition is maximal (the LFP phase is 0) and the low-frequency mode appears when inhibition is minimal (at the rising LFP phases near π). This constitutes a mechanism of phase coding that complements the rate-coding mechanism depicted in Figure 6, according to which fast external fluctuations were unambiguously encoded by a high-frequency bursting activity of the neurons.

DISCUSSION
In this paper we have studied the firing characteristics of a neural network composed of type I neurons. From a single-neuron perspective, both the f -I and PRC have been widely studied and related to the kind of bifurcation that takes the cells from rest to tonic firing. On the other hand, it is known that cortical neurons, wired among them in a complex manner, do not spike periodically but rather sparsely and randomly. The f -I curve of the excitatory and inhibitory neural types helps to determine how the synaptic connection strength has to be tuned in order to balance excitation with inhibition, in agreement with the firing-rates of each population. Moreover, according to their phase response, type I neurons tend to fire faster within the network than in isolation.
The recurrent connections between the excitatory and inhibitory population lead to rhythmicity in the synaptic current impinging on a neuron, composed by a succession of a flood of excitation and a flood of inhibition. These oscillations are a collective phenomenon arising from the interactions among neurons, and their frequency can be slowed down with increasing GABA synaptic duration. In that way, gamma oscillations in the firingrate activity and LFP signal arise naturally from the recurrent connections between excitatory and inhibitory neurons, when the synaptic conductances are adjusted to balance excitation with inhibition.
As long as inhibitory neurons fire more intensely than the excitatory ones to compensate for their relative small number, the network settles in a non-periodic firing state with a significant gamma-band component. At frequencies closer or below the gamma LFP peak the external drive is balanced by inhibition, which leads to a low-firing mode that arises separated from the more "natural" fast-firing mode (given by the single-neuron f -I curve) by a gap of quasi-forbidden instantaneous firing-rates. This leads to a discontinuity in the firing-rates of the neurons, similar to what happens in single type II neurons but that does not occur in type I neurons in isolation. Thus, this firing-rate discontinuity constitutes an emerging dynamical property of the network.
Although the code strategy used by cortical neurons is still a matter of debate, bimodal distributions of instantaneous firing-rates are found experimentally in both the auditory cortex (Shih et al., 2011) and the visual cortex (DeBusk et al., 1997;Snider et al., 1998). Other modeling works have also shown a bimodal distribution of ISIs (Wang, 1998) using intrinsically bursting neurons, in contrast with our case in which the bimodality appears naturally from the balanced network as in (Roxin et al., 2011). Moreover, in Bereshpolova et al. (2011) the recorded single unit activity through different awake brain states, shows that during alert periods some particular neural types reduce their bursting with respect to the non-alert periods. In agreement with our results, the alert state corresponds to higher power in the gamma range.
In conclusion, the balance of an excitatory synaptic current by a strong inhibitory current yields to a discontinuity in the firing-rate of individual neurons forming a neural network. The highest instantaneous rates encode fast fluctuations of the external stimulus, while the spiking times of the network occurring at moderate fluctuations of the input with respect to the mean encode the phase of the LFP. Therefore neural networks can efficiently implement two coding strategies: (1) a rate code for the fast bursting mode, sensitive to rapid changes in the processing of stimuli and (2) a phase code for the slower input fluctuations, according to which isolated spikes occur at the troughs of the LFP, whereas bursts of frequency higher than the gamma peak of the global oscillations arise at the peaks of the LFP. This second feature might contribute to an internal cortical representation of the input. In summary, these results show that cortical population activity depends non-trivially on the dynamical properties of the underlying neurons, and that global population measures shape the firing dynamics of the constituent cells, allowing for multiple encoding mechanisms to be implemented in networks with balanced excitation and inhibition.