Minimal Size of Cell Assemblies Coordinated by Gamma Oscillations

In networks of excitatory and inhibitory neurons with mutual synaptic coupling, specific drive to sub-ensembles of cells often leads to gamma-frequency (25–100 Hz) oscillations. When the number of driven cells is too small, however, the synaptic interactions may not be strong or homogeneous enough to support the mechanism underlying the rhythm. Using a combination of computational simulation and mathematical analysis, we study the breakdown of gamma rhythms as the driven ensembles become too small, or the synaptic interactions become too weak and heterogeneous. Heterogeneities in drives or synaptic strengths play an important role in the breakdown of the rhythms; nonetheless, we find that the analysis of homogeneous networks yields insight into the breakdown of rhythms in heterogeneous networks. In particular, if parameter values are such that in a homogeneous network, it takes several gamma cycles to converge to synchrony, then in a similar, but realistically heterogeneous network, synchrony breaks down altogether. This leads to the surprising conclusion that in a network with realistic heterogeneity, gamma rhythms based on the interaction of excitatory and inhibitory cell populations must arise either rapidly, or not at all. For given synaptic strengths and heterogeneities, there is a (soft) lower bound on the possible number of cells in an ensemble oscillating at gamma frequency, based simply on the requirement that synaptic interactions between the two cell populations be strong enough. This observation suggests explanations for recent experimental results concerning the modulation of gamma oscillations in macaque primary visual cortex by varying spatial stimulus size or attention level, and for our own experimental results, reported here, concerning the optogenetic modulation of gamma oscillations in kainate-activated hippocampal slices. We make specific predictions about the behavior of pyramidal cells and fast-spiking interneurons in these experiments.

: The I-cell is a Wang-Buzsáki model neuron [7]: with for x = m, h, or n, (S4) Here v denotes voltage in mV, t denotes time in ms, and I denotes current density in µA/cm 2 .
The parameters are C = 1 µF/cm 2 , g Na = 35 mS/cm 2 , g K = 9 mS/cm 2 , g L = 0.1 mS/cm 2 , v Na = 55 mV, v K = −90 mV, and v L = −65 mV. The excitatory pulse is assumed to arrive at time 0, modeled by adding the term −G ei e −t/τ e v to the right-hand side of Eq. S1, with τ e = 3 ms and G ei > 0. (The added term drives v towards 0 mV, hence it is excitatory.) We denote by T s the time at which the I-cell spikes; throughout the paper, we take the "spike time" to be the time when the membrane potential v rises above 0 mV. Since T s is a function of G ei and the external drive I i , we also write T s = T s (G ei , I i ). With this notation, panel A of Fig. 2 shows T s (G ei , −0.15) − T s (G ei , 0.15), as a function of G ei , and panel B shows T s (0.7G ei , 0) − T s (G ei , 0), as a function of G ei . Figure 3: The I-cells are as in Fig. 2. The E-cells are as in ref. [5]. Eqs. (S1)-(S4) are as in the I-cell model, with I i replaced by I e in Eq. (S1). Eq. (S5) is replaced by The rate functions α x and β x , x = m, h, and n, are The parameter values of the E-cell model are C = 1 µF/cm 2 , g Na = 100 mS/cm 2 , g K = 80 mS/cm 2 , g L = 0.1 mS/cm 2 , v Na = 50 mV, v K = −100 mV, and v L = −67 mV. We denote in general the number of E-cells in our network by N e , and the number of I-cells by N i . We adopt the synaptic model of ref. [2]. Each synapse is characterized by a synaptic gating variable s associated with the presynaptic neuron, with 0 ≤ s ≤ 1. This variable obeys where H denotes a smoothed Heaviside function: and τ R and τ D are the rise and decay time constants, respectively. To model the synaptic input from neuron j to neuron k, we add to the right-hand side of the equation governing the membrane potential v k of neuron k a term of the form where g( j, k) denotes the maximal conductance density associated with the synapse, s j denotes the gating variable associated with neuron j, and v rev denotes the synaptic reversal potential. For AMPA-receptor-mediated synapses, we use τ R = 0.1 ms, τ D = 3 ms, and v rev = 0 mV; for GABA A -receptor-mediated synapses, τ R = 0.3 ms, τ D = 9 ms, and v rev = −80 mV. We also refer to the maximal conductance g( j, k) as the "strength" of the synapse. The strength of the synaptic connection from the j-th E-cell to the k-th I-cell is random, chosen prior to the beginning of the simulation and fixed from then on: where γ ei and p ei are independent of j and k, γ ei > 0 and 0 < p ei ≤ 1. We define to be the sum of all excitatory synaptic conductance densities impinging upon the k-th I-cell, and Analogous formulas define g ie , G ie , G ie , g ii , G ii , G ii , g ee , G ee , and G ee . We do not include E-to-E-synapses throughout most of this paper (so usually g ee ( j, k) = 0 for all j and k), but see Supporting Information C. The strength of the external input to the j-th E-cell is also random, chosen prior to the beginning of the simulation and fixed throughout: where I e is independent of j, the Z j are independent standard Gaussians, and r e ≥ 0. We say then that the heterogeneity in the external drives to the E-cells is r e × 100%. For instance, when r e = 0.15, we say that there is 15% heterogeneity in the drives to the E-cells. The strength of the external input to the k-th I-cell is defined by where the U k are independent random variables uniformly distributed in [−1, 1], and r i ≥ 0. We do not use a formula exactly analogous to Eq. (S8) for I i,k because drives to the I-cells will often be taken close to zero, and it can therefore be misleading to specify relative fluctuations of those drives. In Fig. 3, N e = 80 and N i = 20. In panels A-I, G ie = 0.3 mS/cm 2 and G ii = 0.05 mS/cm 2 . In panels A, D, and G, G ei = 0.12 mS/cm 2 . In panels B, E, and H, G ei = 0.08 mS/cm 2 . In panels C, F, and I, G ei = 0.04 mS/cm 2 . In panels A-C, p ei = p ie = p ii = 1. In panels D-F, p ei = 0.5 and p ie = p ii = 1. In panels G-I, p ei = p ie = p ii = 0.5. In all panels, I e = 1.5 µA/cm 2 and I i = 0 µA/cm 2 . In panels A-F, external drives are homogeneous: r e = r i = 0. In panels G-I, they are heterogeneous: r e = 0.15 and r i = 0.2.
The parameters in panel J are the same as in panel I, except I i = 0.4 µA/cm 2 and G ii = 0 mS/cm 2 . Figure 4: The rhythmicity measure ρ associated with a network is defined as follows. We simulate the network starting at time t = T 0 = −100 ms and ending at time t = T 1 = 1, 000 ms. To avoid initialization effects, we disregard the initial 100 ms, analyzing the results for t between 0 and T 1 = 1, 000 ms only. We denote the average of the synaptic gating variables associated with the E-cells by s E . Thus our simulation generates s E (t), t = k∆t, with k = 0, 1, 2, ..., T 1 /∆t. We write M = T 1 /∆t. In the simulations of this paper, we always use ∆t = 0.02 ms, and therefore M = 50, 000. We expand s E into a discrete Fourier series of period T 1 : Note that the period of the Fourier mode e iν2πt/T 1 equals T 1 /ν = 1, 000/ν ms, and therefore its frequency is ν Hz. We define (S10) By Parseval's identity, the denominator of the right-hand side of Eq. (S10) equals, up to a factor of ∆t/T 1 , the energy (i.e., discrete L 2 -norm) of s E = s E (t), 0 ≤ t < T 1 . Similarly, the numerator, up to the same factor, equals the energy of the component of s E in the gamma range, defined for the purposes of this figure to be the 30-50 Hz range. (In the simulations of Fig. 4, there are no significant rhythms outside the 30-50 Hz range.) Thus ρ is the fraction of the energy of s E that lies in the gamma range. Figure 5: The E-cell is defined as described earlier in the supporting information about Fig.  3. The drive to the E-cell is I e = 1.6 µA/cm 2 , far above the spiking threshold. We assume that at time zero, the point (v, n, h) lies on the limit cycle, and v = 0 mV, dv/dt > 0. (These conditions determine (v, n, h) uniquely). We add to the right-hand side of the equation describing the evolution of the membrane potential v of the E-cell a term of the form where v inh = −80 mV denotes the reversal potential of inhibitory synapses, and τ i = 9 ms is the decay time constant of inhibition. We denote byT the time of the next spike, where "time of spike", as before, means the time at which v rises above 0 mV, i.e., v = 0 mV and dv/dt > 0. We define δ =T − t * . Figure Fig. 7 as well, except G ie is ten times greater: G ie = 2.0 mS/cm 2 . Figure 9: The quantity ρ plotted in Fig. 9 is essentially that plotted in Fig. 4 (see Eq. (S10)), with one modification: We define the gamma range to be the range from 30 to 60 Hz here. The reason is that the rhythms become fairly fast as G ie is weakened, making a wider definition of the gamma range more appropriate. Thus the definition of ρ used in Fig. 9 is (S11) Figure 10. This is an E/I-network like those in the previous figures, but with N e = 20 and N i = 1. In all three panels of the figure, G ei = 0.12 mS/cm 2 and G ii = 0.05 mS/cm 2 , p ei = p ie = p ii = 1, the drive to the j-th E-cell is (1.38 + 0.02 j) µA/cm 2 , and the drive to the I-cell is zero. In panel A, G ie = 0.24 mS/cm 2 ; in panel B, G ie = 0.12 mS/cm 2 ; in panel C, G ie = 0.06 mS/cm 2 .
The expectation of the total number of points that fall into the unit disk is therefore We would like to choose M and d such that π .
We therefore choose a positive number d, as small as possible but ≥ 1, such that 4d 2 N/π is the square of an integer, then define M = 4d 2 N/π. The number of points in the unit disk generated in this way is random, with expected value N. To obtain exactly N locations, we simply repeat the procedure, with different initial seeds for the random number generator, until, by chance, the exact number of selected locations becomes N. (This never takes very long, and is not an important part of the cost of our simulations.) Figures 15, 16, and 18 show results of simulations in which the neurons are placed in this way in the unit disk, and connection probabilities decay with distance. In these simulations, the E-cells in a smaller disk at the center of radius R are driven strongly; an example is indicated in yellow in Fig. 14. Figure 15. There are N e = 320 E-cells and N i = 80 I-cells, with random spatial locations in a unit disk as described in the supporting information on Fig. 14. We let the probability (not the strength) of E-to-I-connections decay as follows: where d jk is the (non-dimensionalized) spatial distance between the j-th E-cell and the k-th I-cell, δ syn > 0 is the (non-dimensionalized) length scale characterizing the decay of the probability of synaptic connections, and γ ei > 0 is the strength of an individual E-to-I-synapse. The I-to-E and I-to-I-synapses are described by analogous formulas, with the same decay length scale δ syn .
The strengths of individual synapses are γ ei = 0.006 mS/cm 2 and γ ie = γ ii = 0.024 mS/cm 2 . The length scale δ syn characterizing the decay of connection probabilities is taken to be 0.25. The E-cells in the disk centered around the origin with radius R receive constant heterogeneous drive with I e = 1.5 µA/cm 2 , r e = 0.   Fig. 18A are precisely as in Fig. 16A. In panels B and C of Fig.  18, the parameters are the same except for δ syn , which is reduced to 0.20 in panels B and C, and the synaptic strengths, which are tripled in panel C: The strengths of individual synapses are γ ei = 0.009 mS/cm 2 and γ ie = γ ii = 0.036 mS/cm 2 in panel C.
Numerics. All differential equations were solved using the explicit midpoint method with time step ∆t = 0.02 ms.
C. Effect of E-to-E-synapses on PING. We have left out E-to-E synapses in the model networks used in this paper. Here we present numerical experiments showing that to first approximation, E-to-E-synapses simply add excitation to the E-cells, raising the PING frequency but not affecting strong PING in a qualitative way. The additional excitation can of course be crucial, as illustrated for instance by Spencer [6, Fig. 1]. However, in this paper, we provide the needed excitation of the E-cells via external drive. Figure S1A is identical with Fig. 7A. In panels B and C of Fig. S1, we have added rapidly (panel B) or slowly (panel C) decaying E-to-E-synapses; for the details, see the figure caption. The E-to-E-synapses raise the frequency from 44 Hz (panel A) to 60 Hz (panel B) or 68 Hz (panel C). Figure S1 supports the assertion that the effect of E-to-E-synapses is somewhat akin to the effect of raising external drive to the E-cells. We think, however, that the role of E-to-E-synapses in PING rhythms is worth further study in the future. Both slow and fast E-to-E-synapses do affect weak PING: Not surprisingly, they can greatly raise the rate at which the E-cells participate in the rhythm.

D. Effect of heterogeneity on synchronization of I-cells by E-cells when excitation is
weak: Analysis of a model problem. For the purposes of the following discussion, we model an I-cell as the integrate-and-fire neuron defined by When using integrate-and-fire neurons, we non-dimensionalize most physical variables. In particular, v denotes non-dimensionalized membrane potential, shifted and scaled in such a way that the re-set value is 0 and the spiking threshold is 1. However, since firing frequency is crucial to our study, and since we are accustomed to measuring it in Hz, we leave the time t dimensional, measuring it, as everywhere else, in ms. As a result, τ is also a time (also measured in ms), and I is a reciprocal time (measured in ms −1 ). We assume that the neuron does not spike intrinsically: τI < 1. (Note that τI is non-dimensional.) Suppose that at time t = 0, v is at its equilibrium: v(0) = v eq = τI, and then an excitatory input pulse arrives, altering the sub-threshold equation as follows: for t ≥ 0, as long as v < 1, where w > 0 is, like I, a reciprocal time, ϕ ≥ 0, and ∞ 0 ϕ(t)dt = 1. (For simplicity, we model the excitatory pulse as current, not synaptic input. Note that due to the partial non-dimensionalization explained in the previous paragraph, I and w are not, dimensionally, currents here, but reciprocal times.) We are interested in how variations in I or w affect the time that it takes to reach the spiking threshold. We show that if the excitatory pulse can only just barely elicit a response, the dependence of the time to spike on the parameters I and w becomes great.
Note that v(t) = τIe −t/τ + t 0 (I + wϕ(u)) e (u−t)/τ du for t ≥ 0, as long as v < 1. Let T s denote the time at which v reaches the threshold 1, assuming that it ever does. Then We note for later reference that this equation can be re-written as (S17) Note that T s is a function of w and I. We are interested in understanding how sensitively T s depends on w and I, i.e., in the effects of heterogeneity in w and I. We therefore analyze the derivatives ∂T s /∂w and ∂T s /∂I. Differentiating Eq. (S16) with respect to w, we find Using Eq. (S17), Eq. (S18) simplifies as follows: We solve this equation for w∂T s /∂w: Eq. (S15), together with v(T s ) = 1, implies that the denominator of Eq. (S20) is v (T s ): We now assume that I is near 0, so that the model I-cell is driven significantly below the spiking threshold. For I = 0, Eq. (S21) becomes We interpret Eq. (S22) as follows. The left-hand side is the change in T s , namely ∂T s , divided by the relative change in w, namely ∂w/w. When the input pulse is strong, the time that v takes to rise from its resting value to the threshold value should be on the order of a few milliseconds at most, and v (T s ) should not be very much smaller than 1. On the other hand, when the input pulse is only just barely able to trigger a response in the I-cell, then v (T s ) is nearly zero; in this case, and only in this case, the right-hand side of Eq. (S21) is large.
Similarly, differentiating Eq. (S16) with respect to I, we find Using Eq. (S17), Eq. (S23) becomes Evaluating the integral in this equation yields Thus ∂T s /∂I can get large only if v (T s ) gets small, and as discussed earlier, this means that the input pulse is only just barely strong enough to elicit a response in the I-cell. Note that Eq. (S24) implies that T s depends sensitively on I if v (T s ) is small provided that τ is not small as well. We argue in Supporting Information F that τ should not be chosen much smaller than 10, since otherwise the intrinsic period of the integrate-and-fire neuron would become too sensitively dependent on I.
In Eq. (S24), the left-hand side is the change in T s divided by the absolute change in I, not the relative change in I. By contrast, in Eq. (S22), the left-hand side is the change in T s divided by the relative change in w. It is not natural to consider the relative change in I here, since we are primarily interested in values of I near zero; by contrast, w, the strength of the excitatory pulse, is not assumed to be near zero, and it is therefore natural to consider relative, not absolute changes in w in Eq. (S22).

E. Effect of heterogeneity on synchronization of E-cells by I-cells when inhibition is
strong: Analysis of a model problem. We give here an analysis for integrate-and-fire neurons showing that the effects of heterogeneity in the strength of inhibition and in the external drive to the E-cells is largely independent of the strength of inhibition, as long as inhibition is strong enough to bring a homogeneous population to approximate synchrony with a single inhibitory pulse, and the external drive to the E-cells is strong enough to drive a gamma-frequency rhythm. For theta neurons and quadratic integrate-and-fire neurons, similar results can be found in [1] and [3].
We consider a single integrate-and-fire neuron defined by Eqs. (S13), (S14), but now think of it as an E-cell, and assume supra-threshold drive: τI > 1; we examine how a strong inhibitory pulse acts on such a model neuron, and study in particular the dependence of the response on external drive and strength of inhibition. If this dependence is weak, then heterogeneity effects will be weak when an inhibitory pulse synchronizes a heterogeneous population of E-cells.
Assume that at time t = 0, a synaptic pulse of inhibition sets in, decaying exponentially with time constant τ i > 0. For simplicity, we assume that the synaptic reversal potential is the same as the reset potential, namely 0. (Recall that when using the integrate-and-fire model, we always assume that the membrane potential is shifted and scaled so that the reset potential is 0 and the firing threshold is 1.) Thus the equations governing with 0 ≤ v 0 ≤ 1. (Here g, like I, is a reciprocal time, measured in ms −1 .) Let T P > 0 be the time at which v reaches the spiking threshold 1. (Since τI > 1, there is such a time.) If v 0 = 0, one can think of T P as the period of the strong PING rhythm in a highly idealized setting. Even in this very simple model, T P is a function of five variables: I, g, τ, τ i , and v 0 . To make a complete parameter study feasible, we fix τ i = 10 ms and τ = 10 ms. The choice τ i = 10 ms is motivated by the fact that the decay time constant of GABA A -receptormediated inhibition is about 10 ms. We argue in Supporting Information F that τ should not be chosen much smaller than 10 ms, since otherwise the intrinsic period of the integrateand-fire neuron would become too sensitively dependent on I. We have, however, tried other values of τ, and found that our conclusions are not substantially different for larger or (somewhat) smaller values of τ. Now T P depends on only three variables: T P = T P (I, g, v 0 ), I > 0.1, g > 0, v 0 ≤ 1. We will focus on values of g large enough to cause approximate synchronization of a population of model neurons governed by (S25) with different initial values of v 0 . We will assume that g > I − 1/τ. Then the right-hand side of Eq. (S25) is negative for t = 0 and v = 1, so the pulse of inhibition prevents spiking for some positive amount of time, regardless of v 0 . Consequently "T P (I, g, 1)", which we define to be the limit of T P (I, g, v 0 ) as v 0 tends to 1 from below, is positive.
Unfortunately, not even for this very simple model is there a simple analytic answer to the question "For given I, how large does g have to be for synchronization to be tight?" However, the question can easily be answered numerically. For a given I, we define g I to be the value of g for which T P (I, g, 1) differs from T P (I, g, 0) by exactly 10%: T P (I, g I , 0) − T P (I, g I , 1) T P (I, g I , 0) = 0.1.
(The choice of 10% is arbitrary here. Our results would be qualitatively very similar if we had chosen 5% or 50% here.) For g ≥ g I , reasonably tight synchronization is achieved by a single inhibitory pulse: All T P (I, g, v 0 ), 0 ≤ v 0 ≤ 1, are within 10% of T P (I, g, 0). We now write I = 1/τ + J, J > 0, and g = g I + h, h > 0. Thus J is the amount of drive above the spiking threshold 1/τ, and h is the inhibitory conductance above the minimal conductance g I needed for synchronization within 10% as a result of a single pulse. We plot T P (I, g, 0) = T P (1/τ + J, g I + h, 0) as a function of J and h. This is shown in Fig. S2, panel A. In panel B of the same figure, we show the quantity κ = I T P (I, g, 0) plotted as a function of J and h. Note that the definition of κ can be re-written (omitting, for brevity, the arguments (I, g, 0)) as κ = ∂T P /T P ∂I/I .
A relative change in I causes a κ times greater relative change in the "PING period" T P (I, g, 0). If κ is large, then even modest heterogeneity in I can have a large effect. Thus κ is a measure of sensitivity to perturbations in I; it is what is called a "condition number" in Numerical Analysis. Note that κ and the expressions on the left-hand sides of Eqs. (S22) and (S24) are of a similar nature. However, in Eqs. (S22) and (S24), we considered absolute, not relative changes in T s . In Eqs. (S22) and (S24), T s is expected to be small, since it is the delay between the excitatory spike volley and the response of the inhibitory cells; therefore relative changes in T s are less informative than absolute ones. Similarly, in Eq. (S24), we considered absolute, not relative changes in I because we thought of I as small. Here we think of neither T P nor I as small, and therefore consider relative changes in both T P and I.
Panel A of Fig. S2 indicates that the PING period never varies very rapidly with h (that is, with g). Thus as long as inhibition is strong enough to cause approximate synchronization in a single pulse (h ≥ 0), heterogeneity in g should be expected to have moderate effects only. Panel B shows that the sensitivity measure κ is nearly independent of h, and becomes somewhat larger only when J gets small, i.e., when drive to the E-cell is near the firing threshold. (Panel A shows that in that case, the PING period is around 40 ms, corresponding to a 25 Hz oscillation, at the lower edge of the gamma range.) F. The rhythm in the I-cells in Fig. 7C is an ING-rhythm. In Fig. 7C, there is a clear rhythm in the I-cells. We noted earlier that this is an ING-rhythm. Figure S3 confirms this point. Panel A of the figure shows a closeup of Fig. 7C (only spike times of 10 E-cells and 10 I-cells are shown), demonstrating that there is a rhythm in the I-cells. Panel B is the same as A, but with the I-to-I-synapses removed. Rhythmicity is clearly much reduced, although a remnant of rhythmic behavior still appears to be visible in the I-cells. We believe that this is an accidental finite-network-size effect. In fact, if we perform the same simulation in a twice larger network (Panel C), there appears to be yet less rhythmicity than in Panel B. fire model (S13), (S14) not be very small. One might think of τ as the "membrane time constant". However, τ also governs the nature of the dependence of the intrinsic firing frequency of the integrate-and-fire neuron on I. We will argue here that for τ 10 ms, that dependence becomes unreasonably sensitive. The intrinsic period of the integrate-and-fire neuron defined by Eqs. (S13) and (S14) is The absolute value of the right-hand side of Eq. (S31) grows rapidly as T /τ increases, as shown in Fig. S4. For instance, if T /τ = 5, a 1% increase in I causes a decrease in T by approximately 30%, according to this analysis. (The analysis is approximate because it is based on infinitesimal perturbations dI of I.) In neuroscience, one is typically interested in neurons that spike with intrinsic periods T of several tens of milliseconds. If τ 10 ms, then the dependence of T on I becomes extremely sensitive for values of I that yield intrinsic periods of this order of magnitude.
Our analysis reflects a flaw of the integrate-and-fire model. In a more realistic neuronal model, the membrane time "constant" (which typically is not a constant, but varies as conductances vary) can be much smaller than 10 ms without there being any hyper-sensitivity issue.