## Abstract

The striatum is composed of GABAergic medium spiny neurons with inhibitory collaterals forming a sparse random asymmetric network and receiving an excitatory glutamatergic cortical projection. Because the inhibitory collaterals are sparse and weak, their role in striatal network dynamics is puzzling. However, here we show by simulation of a striatal inhibitory network model composed of spiking neurons that cells form assemblies that fire in sequential coherent episodes and display complex identity–temporal spiking patterns even when cortical excitation is simply constant or fluctuating noisily. Strongly correlated large-scale firing rate fluctuations on slow behaviorally relevant timescales of hundreds of milliseconds are shown by members of the same assembly whereas members of different assemblies show strong negative correlation, and we show how randomly connected spiking networks can generate this activity. Cells display highly irregular spiking with high coefficients of variation, broadly distributed low firing rates, and interspike interval distributions that are consistent with exponentially tailed power laws. Although firing rates vary coherently on slow timescales, precise spiking synchronization is absent in general. Our model only requires the minimal but striatally realistic assumptions of sparse to intermediate random connectivity, weak inhibitory synapses, and sufficient cortical excitation so that some cells are depolarized above the firing threshold during up states. Our results are in good qualitative agreement with experimental studies, consistent with recently determined striatal anatomy and physiology, and support a new view of endogenously generated metastable state switching dynamics of the striatal network underlying its information processing operations.

## Introduction

The striatum is the main input structure of the basal ganglia (BG), receiving excitatory glutamatergic inputs from the entire cerebral cortex (McGeorge and Faull, 1989). Medium spiny neurons (MSNs), which account for >90% of striatal neurons (Oorschot, 1996), form inhibitory synapses with each other. Many authors have interpreted this anatomy as a winner-take-all (WTA) network, the most strongly activated neuron-suppressing activity in its neighbors (Groves, 1983; Wickens et al., 1991; Beiser and Houk, 1998; Fukai, 1999; Suri and Schultz, 1999; Bar-Gad and Bergman, 2001). However, several experimental findings argue against this interpretation (Czubayko and Plenz, 2002; Tunstall et al., 2002; Koos et al., 2004; Taverna et al., 2004). First, direct testing for interactions among nearby MSNs showed sparse connectivity, the probability of a connection being ∼10%. Second, individual connections typically involve one or only a few synapses, suggesting that weak interactions predominate. Third, reciprocal interactions are rare, with the majority of MSN pairs involved in only one-way connections. Finally, experimental observations show that highly irregular firing predominates (Wilson, 1993) and that cells form assemblies that fire in coherent episodes *in vivo* (Miller et al., 2008) and *in vitro* (Carrillo-Reid et al., 2008).

These findings are inconsistent with the proposal that WTA operations are a fundamental aspect of striatal information processing, but the absence of a good alternative model raises the question of what contribution lateral inhibition between MSNs makes to striatal network dynamics. To address this question and as a first step to understand the real function of lateral inhibition between MSNs in the striatum, it is important to characterize the dynamical behavior generated by the MSN network under the realistic assumptions of sparse, weak and random connectivity using a biophysically minimal spiking cell model representing only relevant definitive *in vivo* up-state MSN characteristics.

We ask what quantitative predictions of up-state MSN firing activity can be made based on this minimal network model and how these predictions depend on connectivity. In particular, we investigate whether the activity generated by the network at striatal connectivities can account for the puzzling conundrum of observations of broad but low MSN firing rate distributions, highly irregular episodic firing but with only single peaked broad interspike interval (ISI) distributions, and slowly varying cross-correlation and autocorrelation but without precisely synchronized spiking.

We first describe the network model of deterministically interacting up-state MSNs. We describe spiking activity in large 500-cell networks demonstrating highly irregular episodic firing and complex identity–temporal patterns. We find that the striatal network itself endogenously generates episodically firing sequentially switching coherent cell assemblies without the need for patterned cortical excitation. We show how this activity arises from dynamical metastable state switching and depends on network structure details. Finally, we investigate quantitatively how assembly formation and MSN firing activity varies with network connectivity. We suggest that the striatal network may be adapted to generate episodically firing cell assemblies producing large slow sequential fluctuations on behaviorally relevant timescales for exploratory sequence generation.

## Materials and Methods

The network (Fig. 1) is composed of model MSNs randomly connected. In the simulations reported here, we use random networks in which cells *i* and *j* are connected with fixed probability *p*, and there are no self-connections. This produces random networks with binomial degree distributions for both incoming connections and outgoing connections, which is the simplest type of network topology. Furthermore, if a pair of cells *i* and *j* are connected, the synaptic conductance *k*_{ij}^{syn} has a random strength given by *k*_{ij}^{syn} = (*k*^{syn}/*p*)ε* _{ij}*, where ε

*is a uniform quenched random variable on [0.8, 1.2] independent in*

_{ij}*i*and

*j*, and

*k*

^{syn}is a fixed parameter that is rescaled by the connection probability

*p*so that, when the network connectivity is varied, the average total inhibition on each cell is constant independent of

*p*.

Striatal MSNs are likely to be contacted by of the order of ∼500 other cells and similarly to contact ∼500 (Oorschot et al., 2002; Plenz, 2003; Koos et al., 2004; Tepper et al., 2004; Wickens et al., 2007). Furthermore, the probability of a connection between nearby MSNs is estimated to be fairly sparse, between *p* = 0.1 and *p* = 0.2 (Oorschot et al., 2002; Plenz, 2003; Koos et al., 2004; Tepper et al., 2004; Wickens et al., 2007). To simulate a striatal network that respects these two figures would require a network of ∼500/0.15 = 3300 cells. However, this argument neglects the fact that not all cells are cortically excited into the up state. Such cells have no effect on striatal dynamics and can therefore be left out of network simulations. We assume that ∼10–30% of MSNs are cortically excited at any time and accordingly simulate networks of 500 up-state MSNs with connectivities of approximately *p* = 0.15. This translates into each up-state MSN being inhibited by, and inhibiting around, 50–150 cortically excited potentially spiking MSNs.

To describe the MSN cells, we use the *I*_{Na,p} + *I _{k}* model described by Izhikevich (2005). The

*I*

_{Na,p}+

*I*cell model is two-dimensional and described by the following: having leak current

_{k}*I*, persistent Na

_{L}^{+}current

*I*

_{Na,p}with instantaneous activation kinetic, and a relatively slower persistent K

^{+}current

*I*(

_{K}. V_{i}*t*) is the membrane potential of the

*i*th cell,

*C*the membrane capacitance,

*E*

_{L,}_{Na,k}are the channel reversal potentials, and

*g*

_{L,}_{Na,k}are the maximal conductances.

*n*(

_{i}*t*) is K

^{+}channel activation variable of the

*i*th cell. The steady-state activation curves

*m*

_{∞}and

*n*

_{∞}are both described by the following: where

*x*denotes

*m*or

*n*, and

*V*

_{∞}

^{x}and

*k*

_{∞}

^{x}are fixed parameters. τ

*is the fixed timescale of the K*

_{n}^{+}activation variable. The term

*I*(

_{i}*t*) is the input current to the

*i*th cell.

In this study, we ask how far network dynamics can account for observations of up-state *in vivo* MSN firing behavior. We have therefore chosen a simple cell model and made no attempt to make a biophysically detailed model including all the MSN ion channels or to reproduce resting state membrane potential values and other MSN characteristics such as up–down-state transitions or long latencies to first spike. The important point of the cell model is that the cell parameters are chosen so that the cell is the vicinity of a saddle-node on invariant circle (SNIC) bifurcation. As the current *I _{i}*(

*t*) in Equation 1 increases through the bifurcation point, the stable node fixed point and the unstable saddle fixed point annihilate each other, and a limit cycle having zero frequency is formed (Izhikevich, 2005). Increasing current further increases the frequency of the limit cycle. This models the dynamics when the cells are in the up state receiving excitatory drive to firing threshold levels of depolarization.

The response of the model MSN cell to excitatory current is illustrated in Figure 2*a* for four different values of injected current. For injected currents below the SNIC bifurcation current of 4.51 μA/cm^{2}, the cell depolarizes to a fixed membrane potential level that depends on the level of the injected current. At injected current levels above the bifurcation, the cell fires regularly at arbitrarily low frequency without spike frequency adaptation.

The SNIC bifurcation is an appropriate bifurcation to use for a model of an MSN in the up state, because its dynamics are in good qualitative agreement with studies of MSN firing (Hikosaka et al., 1989; Wilson, 1993; Wilson and Kawaguchi, 1996; Wickens and Wilson, 1998; Kitano et al., 2002). First, the SNIC bifurcation allows firing at arbitrarily low frequencies (Izhikevich, 2005), which is important because MSNs are known to fire with very low frequencies (Hikosaka et al., 1989; Wilson, 1993). Second, MSNs do not show subthreshold oscillations (Nisenbaum and Wilson, 1995; Wilson and Kawaguchi, 1996) (as models based on Hopf bifurcations do). Third, during the spike repolarization, the membrane potential for the SNIC bifurcation can decrease below the spike threshold potential, appropriately for MSNs (Plenz and Aertsen, 1996; Wickens and Wilson, 1998). Fourth, the SNIC bifurcation does not allow bistability between a spiking state and a quiescent state (Izhikevich, 2005), in agreement with studies of up-state MSNs (Nisenbaum and Wilson, 1995; Wilson and Kawaguchi, 1996; Wickens and Wilson, 1998). Finally, the membrane potential in the vicinity of an SNIC bifurcation shows a point of inflection between successive spikes as observed in MSNs (Wickens and Wilson, 1998). In contrast to integrate-and-fire models, the cell can easily switch between limit cycle spiking and resting behavior without the need for a dynamically discontinuous reset after firing yet has only two variables making large-scale simulations tractable. Neurocomputationally, the SNIC bifurcation is appropriate for cells that may act as integrators and coincidence detectors of synaptic inputs, which do not display resonance to inputs of a certain preferred frequency, which fire well defined all-or-nothing spikes, and which can fire arbitrarily slowly so that firing frequency varies gradually with input current without a sudden transition from a quiescent state to a high-frequency firing state (Izhikevich, 2005).

It is important to point out, however, that the particular form of the bifurcation is not crucial for the formation of episodically firing cell assemblies. As described below, the important factor the cell model must have is simply proximity to an appropriate bifurcation between firing and quiescent states. The SNIC bifurcation has been chosen for the above described similarity to *in vivo* up-state MSN activity, in particular the low firing rates and absence of subthreshold oscillations. However, other bifurcations, such as a saddle-homoclinic bifurcation (Izhikevich, 2005), may be used if one desires other cell properties, such as bistability between a strongly hyperpolarized resting state and a depolarized firing state, under application of NMDA for example. As will be described, it is the appropriate network structure and attributes that allow the formation of switching episodically firing assemblies, not the particular cell model. However, the appropriate choice of bifurcation does allow comparison with experimental studies of up-state MSN spiking activity, such as ISI distributions.

The input current *I _{i}*(

*t*) in Equation 1 is composed of two parts and is given by the following: The second term represents the inhibitory input from the MSN network (see below). The first part, described by the term

*I*

_{i}

^{c}represents all other input from sources other than the MSN network. This predominantly includes excitatory input from the cortex and thalamus but also inhibitory input from the striatal interneurons. In general, this term would include complex time-varying fluctuations patterned across cells, generated internally in the nervous system or by changes in the external environment that would modulate MSN network dynamics. Here, however, we want to understand the dynamics endogenously generated by the MSN network itself before considering how the dynamical variation of inputs further affect the network dynamics. In this paper, we therefore restrict ourselves to the two simplest forms this input can take: (1) constant and fixed or (2) fluctuating randomly. In the fixed input condition,

*I*

_{i}

^{c}has a fixed magnitude for the duration of a simulation but varying across cells. This can also be considered an input varying very slowly on the timescale of MSN internally generated network dynamics. In the simulations reported here, the

*I*

_{i}

^{c}are quenched random variables drawn uniformly randomly once at the start of the simulation. They are drawn from the interval [

*I*

_{bif},

*I*

_{bif}+ 1], where

*I*

_{bif}= 4.51μA/cm

^{2}is the current at the firing threshold so that all cells receive net excitatory input. Indeed, suppose that

*I*

_{i}

^{c}for cell

*i*were below the firing threshold

*I*

_{bif}(attributable, for example, to strong input from inhibitory interneurons), then, because the MSN cell would never fire, it can be left out of the simulation altogether. In the fluctuating input condition, the

*I*

_{i}

^{c}are drawn randomly in exactly the same way, but this random drawing is performed anew every 10 ms throughout the duration of the simulation. Notice that this second condition is a very stringent test of the tendency of the network to form assemblies because, in this case, all cells receive the same average levels of excitation. In both cases of fixed and fluctuating excitation, these values of input current mean that all cells would be on limit cycles just slightly above the firing threshold and firing with low rates if the MSN network inhibition were not present. In the case of fixed input, all cells would fire perfectly regularly with their own specific quenched random rates, whereas, in the case of fluctuating input, all cells would fire with the same average rate but not perfectly regularly. In fact, the MSN inhibitory network causes some cells to become quiescent by reducing the total input current to below the bifurcation point.

The assumption that the neurons would be depolarized above threshold were it not for the synaptic inhibition is compatible with the behavior of real striatal neurons *in vitro* and *in vivo*. Bernardi et al. (1976), using intracellular recording in anesthetized whole animals, showed that intravenous injections of the GABA antagonist bicuculline led to repetitive action potential firing in response to stimulation of the cerebral cortex. Also in anesthetized animals, West et al. (2002) showed that local perfusion of bicuculline rapidly depolarized neurons and increased spike activity. In cortex–striatum cocultures, Plenz and Aertsen (1996) showed that local injection of bicuculline caused individual MSNs to enter sustained periods of strong firing. In the more artificial conditions of the NMDA-treated slice, GABA_{A} antagonists increased the duration of the up states (Vergara et al., 2003). These findings support the assumption that striatal cells in the up state fire if GABA input is removed.

Because the inhibitory current part is provided by the GABAergic collaterals of the MSN striatal network, it is dynamically variable. These synapses are described by Rall-type synapses (Rall, 1967; Destexhe et al., 1968) in Equation 3, in which the current into postsynaptic neuron *i* is summed over all inhibitory presynaptic neurons *j*, and *V*_{syn} and *k*_{ij}^{syn} are channel parameters. *g _{j}*(

*t*) is the quantity of postsynaptically bound neurotransmitter given by the following: for the

*j*th presynaptic cell. Here,

*V*

_{th}is a threshold, and θ(

*x*) is the Heaviside function.

*g*is essentially a low-pass filter of presynaptic firing. The timescale τ

_{j}*should be set relatively large so that the postsynaptic conductance follows the exponentially decaying time average of many preceding presynaptic high-frequency spikes.*

_{g}Figure 2*b* illustrates how IPSPs elicited by presynaptic spikes depend on the membrane potential of the postsynaptic cell. The synaptic strength parameter *k*^{syn}/*p* here has been set to the same low value as used in the 500-cell network simulations at striatally relevant connectivities of *p* = 0.2 connections per cell, described below. The presynaptic cell fires as a result of current injected at 4.52 μA/cm^{2} (as in one of the time series in Fig. 2*a*) from *t* = 400. Its spike times are shown as dashed vertical lines. Postsynaptic cells are depolarized to different membrane potentials by excitatory current, and IPSPs are generated by neurotransmitter release when the presynaptic membrane potential exceeds *V*_{th} = −40 mV (Eq. 4). The bottom-most panel shows the PSPs for current injection of 0. IPSPs are depolarizing in this case because the synaptic reversal potential of *V*_{syn} = −65 mV (Eq. 3) exceeds the holding potential. At the larger injected currents (but still below the firing threshold) of 4.1 and 4.5, hyperpolarizing IPSPs are generated. Near the firing threshold at ∼4.51, the effect of an inhibitory presynaptic spike is to very slightly depress the membrane potential of the postsynaptic cell, in good qualitative agreement with previous studies (Czubayko and Plenz, 2002; Tunstall et al., 2002; Plenz, 2003; Koos et al., 2004; Taverna et al., 2004), and IPSPs have amplitudes of ∼200 μV, very similar to real striatal IPSPs.

The top two panels of Figure 2*b* show examples in which the excitatory current was high enough for the postsynaptic cell to be in the firing state before activation of the presynaptic cell. The spiking of the presynaptic cell after *t* = 400 ms is able to completely quench the firing of the postsynaptic cell with current injection of 4.53 μA/cm^{2}, although the current injected to the postsynaptic cell is greater than that to the presynaptic one. When postsynaptic cell current injection is increased to 4.55, however, the presynaptic spikes cannot quench the postsynaptic activity because of the weak strength of the inhibitory connection. However, the presynaptic spikes do lower the average postsynaptic firing rate and also cause the postsynaptic spiking to become irregular. This is because the presynaptic spike can sometimes alter the timing after the postsynaptic spike (for example, the second, fourth, and fifth spikes in Fig. 2*b*, top panel), whereas at other times, a presynaptic spike has little effect. Indeed, depending on the phase relationship of the membrane potentials of presynaptic and postsynaptic cells, presynaptic spikes delay or even hasten the next postsynaptic spike (Rinzel and Ermentrout, 1989; Van Vreeswijk et al., 1994; Ermentrout, 1998; Izhikevich, 2005), which is consistent with observations of MSNs (Plenz, 2003, his Fig. 5*a*).

It should be noted that the effect an inhibitory connection has on the postsynaptic cell depends on both the connection strength and the proximity of the postsynaptic cell to the bifurcation. The results do not require an unrealistically strong synaptic current.

The network model we present here is minimal and contains only the minimal requirements needed to produce the desired dynamical activity of coherent episodically firing cell assemblies. This choice of minimal requirements allows us to maximize model generality and robustness while elucidating pertinent and poorly understood basic facts of striatal anatomy. Providing the basic assumptions of the model, i.e., proximity to a bifurcation (of appropriate type) produced by weak excitation and weak lateral inhibition and sparse to intermediate random connectivity, are not violated, the model may be applicable to understand the formation of coherent episodically firing assemblies in several different types of *in vivo* and *in vitro* striatal preparations, although specific details of the spiking may differ.

All simulations were performed with fourth-order Runge–Kutta.

## Results

Below in *Episodic firing in striatal network*, we describe the typical episodic firing of MSN cells in the network. In *Sequentially firing cell assemblies and metastable state switching*, we describe sequentially switching cell assembly formation and show how it arises from metastable state switching depending on network structure details. In *MSN firing statistics predicted by the network model with striatal connectivity*, we use model simulations at striatal connectivities to derive quantitative predictions of MSN firing and also to address how closely our model accounts for known MSN firing patterns. Finally, in *Connectivity variation*, we study how predictions of MSN firing derived from network simulations vary with network connectivity.

### Episodic firing in striatal network

Sparsely connected networks produced irregular switching between firing and quiescent states. This occurs even when excitatory input is fixed and network dynamics is therefore completely deterministic. Figure 3*a* shows a time series segment of membrane potentials *V _{i}*(

*t*) for some randomly selected cells from such an

*N*= 500-cell network with connectivity 0.2 under the fixed excitatory input condition. Irregular switching between firing and quiescent states can clearly be seen. Even when cells are in the firing state, they do not fire spikes regularly, and average frequencies vary between cells and between firing episodes in the same cell. The periods of time cells spend in firing episodes or quiescent states also appears to vary randomly. However, the network simulation has no stochastic variables, and therefore this irregular episodic firing is caused by a type of deterministic chaos. The small chaotic fluctuations in membrane potential can be seen in the top two panels of Figure 3

*a*, which show cells that do not fire spikes at all during the period shown. Notice that the membrane potential of these quiescent cells fluctuates a few millivolts below spiking threshold in agreement with studies of up-state MSNs (Wilson and Groves, 1981; Wilson, 1993; Wilson and Kawaguchi, 1996; Stern et al., 1998; Wickens and Wilson, 1998).

The corresponding time series segment of *g _{j}*(

*t*), the quantity of postsynaptically bound inhibitory transmitter contributed to postsynaptic cell

*i*by presynaptic cell

*j*, is shown in Figure 3

*b*. During firing episodes,

*g*(

_{j}*t*) slowly rises to a stable value reflecting the presynaptic firing rate (always below 0.1 for the parameter settings here), decorated by oscillations that reflect the individual spikes of the presynaptic cell. Between firing episodes,

*g*(

_{j}*t*) decays exponentially to zero on a slow timescale of several hundred milliseconds. The inhibitory current (Eq. 3) to postsynaptic cell

*i*is made up of a sum over a subset of

*g*(

_{j}*t*), such as the ones shown in Figure 3

*b*.

MSNs show relatively low firing rates with long irregular periods of quiescence between episodes of irregular spike firing, despite the fact that individual MSN cells do not show burst firing intrinsically (DeLong, 1973; Hikosaka et al., 1989; Wilson, 1993; Jaeger et al., 1995; Jog et al., 1999; Barnes et al., 2005). Firing episodes last on the order of many hundreds of milliseconds (Plenz and Aertsen, 1996; Miller et al., 2008). This is usually attributed to episodic excitation from cerebral cortex, which gives rise to large-amplitude subthreshold membrane potential fluctuations. However, our results show that this activity may also arise as a network property in the absence of cortical fluctuations. This occurs for several reasons.

First, the proximity of the parameters to the bifurcation, as explained in Materials and Methods, in the presence of excitatory driving maintains the cells in the up state just above the firing threshold permitting low firing rates. The inhibitory MSN network connections are weak, and individual cells do not quench postsynaptic cells but produce irregular spiking. However, fluctuations in network inhibition as a result of simultaneous activation of many presynaptic MSN cells can produce prolonged periods of quiescence in the postsynaptic cell. Network inhibition is never strong enough to depress the membrane potential more than a few millivolts below the up-state firing threshold, however, and therefore episodes of spiking appear repeatedly from time to time.

Second, the slow synaptic timescale τ* _{g}* (Eq. 4) facilitates episodic firing. In the simulations reported here, we have set the synaptic timescale τ

*≈ 50, which means that postsynaptically bound transmitter exponentially decays to half its value in time τ*

_{g}*(2) ≈ 34 ms. While sufficient presynaptic cells remain in firing states, the postsynaptic cell will remain quiescent, and only when sufficient presynaptic cells cease to fire will the postsynaptic cell slowly revert to its firing state as a result of weak excitation and slow decay of postsynaptic neurotransmitter.*

_{g}lnThird, the other important ingredient to produce network burst firing is to have the appropriate sparse network connectivity, as will be described in detail below.

### Sequentially firing cell assemblies and metastable state switching

In *Episodic firing in striatal network*, we showed that cells fire in irregular episodes in striatal network simulations, but do cells tend to fire cooperatively in assemblies or more individually? Here we study raster plots of cell spiking that illustrate the formation of sequentially switching assemblies in which cells fire in coherent episodes with complex identity–temporal structure through metastable state switching. Subsequently, we describe how cell assemblies depend on the details of the network structure.

#### Complex identity–temporal patterns

Figure 4*a* shows a spike raster plot from a 500-cell network simulation with striatally realistic connectivity of 0.1 (so that each cell is inhibited by ∼50 others) under fluctuating excitatory input. Because our network model has no intrinsic cell ordering, the cells in the raster plot have been ordered in a special way to reveal the presence of cell assemblies varying on slow timescales. To produce this cell ordering, we first calculate firing rate time series *R _{i}*(

*t*) for each cell

*i.*The calculation of a firing rate depends on the bin size used to count spikes. Here we use a sliding bin of 2000 ms to extract assemblies on this long timescale. Next the zero time lag cross-correlation matrix

*C*is calculated from the rate time series for all pairs of cells

_{ij}*i*and

*j*(see Appendix). Finally, a clustering algorithm,

*k*-means, (see Appendix), is applied to the cross-correlation matrix to group the cells into clusters. The cells in the raster plot in Figure 4

*a*and the corresponding cross-correlation matrix in Figure 4

*b*are then ordered according to these clusters.

Although the excitatory input to the network is noisy and fluctuates rapidly, the spike raster plot (Fig. 4*a*) demonstrates slowly evolving dynamical structure with complex “identity–temporal” patterns. Identity–temporal patterns are the analog of spatial–temporal patterns in systems without a real spatial dimension. The patterns occur not in real space but in the space defined by the labeling of the cells ordered by the clustering algorithm. The complex patterns are produced by switching cell assemblies in which many cells switch between firing and quiescent states together.

To demonstrate the emergence of switching cell assemblies in more detail, it is instructive to use smaller networks. Figures 5 and 6 show examples of “metastable state switching” in randomly connected 100 cell networks with connectivity of 0.2 so that each cell is inhibited by, and inhibits, 20 cells on average. Figure 5 shows an example in which excitatory input is constant, and Figure 6 shows an example in which excitatory input is noisy and fluctuating rapidly. Again the cells in the spike raster plots (Figs. 5*a*, 6*a*) and cross-correlation matrices (Figs. 5*b*, 6*b*) have been ordered by the clustering algorithm. To visualize the assemblies further, the cells in the spike raster plots have also been colored according to their assigned clusters. We have also calculated “cluster firing rate time series,” shown in Figures 5*c* and 6*c*. These are produced by merging the spikes from all cells in a given cluster together, preserving the timing of each spike, and then calculating rate time series from these cluster spike time series.

The spike raster plots and the cluster firing rate time series demonstrate that sequentially switching episodically firing cell assemblies are the normal dynamical mode in these networks for both constant excitatory input (Fig. 5) and noisy fluctuating excitatory input (Fig. 6). This is particularly evident from the large oscillatory type switches between dominant assemblies in the cluster firing rate time series (Figs. 5*c*, 6*c*). The example with fixed excitatory input shows more precisely patterned repetitive activity than the example with fluctuating input, but noisy fluctuation in input does not destroy the tendency of the network to form assemblies. This is because the patterned activity is structured endogenously in the striatal MSN network by metastable state switching of the underlying deterministic dynamical system. In both examples, assemblies are made up of many cells firing episodes of multiple irregular spikes, not single spikes, so that the cells are not, in general, synchronized at the spiking level even when excitatory input is fixed and the system is therefore completely deterministic. Variability in the timing of individual spikes does not interfere with the tendency of the cells to form assemblies. Different cell assemblies sometimes seem to enter firing episodes simultaneously in a positively correlated way, but, at other times, this correlation seems to be broken (Fig. 6*a*,*c*, orange and green assemblies). Any particular assembly can go through spells of periodically switching between firing episodes and quiescence before becoming quiescent for longer periods, (Fig. 5*a*,*c*, orange assembly). Other cell assemblies fire very occasionally in isolated episodes (Fig. 5*a*,*c*, red assembly). The division of cells into assemblies is not perfect, however, and some cells can seem to fire with one assembly at some times and a different assembly at other times.

Sequential state switching between different active assemblies can be behaviorally relevant because it occurs repetitively on the typical behavioral timescales of seconds. Furthermore, multiple different timescales can be represented together. For example, in Figure 5, *a* and *c*, during the period from approximately *t* = 13,000 to *t* = 17,500, the pink and black assemblies are dominant and display approximately periodic firing episodes with a period of ∼650 ms, whereas the other cell assemblies are primarily quiescent. However, this approximate periodic state is interrupted from time to time (e.g., from *t* = 17,500 to *t* = 22,000) on a longer timescale by the orange and blue assemblies, which also fire episodes rhythmically, whereas the pink assembly is suppressed.

This assembly structure can be confirmed from the corresponding cross-correlation matrices shown in Figures 5*b* and 6*b*. Cells in the same cluster are positively correlated (see scale bars) with each other, which produces the violet and pink colored squares on the main diagonal. Cells in different clusters are negatively correlated, as can be seen by the cyan and yellow regions off the main diagonal. The scale bars show that both positive correlations within assemblies and negative correlations between assemblies can reach large highly significant values. The cross-correlation matrix for the case of fluctuating input (Fig. 6) exhibits similar structure to the case of constant excitatory input (Fig. 5), further demonstrating the resilience of the network generated assemblies to noisy input. When networks are larger, the number of cell assemblies and the complexity of the relationships between them increases, as is shown by the checkerboard appearance of the cross-correlation matrix for the 500-cell simulation at striatal connectivity shown in Figure 4*b*. Again, cells within an assembly on the main diagonal show strong positive correlation, but cells in different clusters may show strong positive or negative correlation. Although some cell assemblies show positive correlation between them, they differ in their correlations to other assemblies, so that no two assemblies have an identical pattern of positive and negative relationships.

To demonstrate metastable state switching in the network in a different and striking way, we also calculated network state transition matrices, *D*(*t*_{1},*t*_{2}) (Schreiber et al., 2003; Sasaki et al., 2006b, 2007; Carrillo-Reid et al., 2008). These are shown in Figures 5*d* and 6*d* for the same two network simulations. To construct these plots, the scalar products of vectors *R*(*t→*_{1}) and *R*(*t→*_{2}) made from the firing rates of all cells, *R _{i}*(

*t*) (where now a small 40 ms time bin size is used), are taken at all different times

*t*

_{1}and

*t*

_{2}(see Appendix). The square blocks of violet and pink colors persisting for several seconds indicate that the network is in a certain persistent firing state,

*R*(

*t→*) ≈

*R*(

*t→*+ 1) ≈

*R*(

*t→*+ 2) …. This persistent state repeats periodically separated by green and blue epochs in which the network visits a different dissimilar state. For example, the large solid block of violet from

*t*= 13 s to

*t*= 17 s in Figure 5

*d*indicates a persistent approximately fixed network state that is then revisited at later times, at approximately

*t*= 22 s and

*t*= 25 s. This provides clear demonstration of abrupt transitions between metastable states in the network that are then revisited later, and it occurs despite rapid noisy fluctuation in the excitatory input (Fig. 6

*d*).

The off-diagonal streaks of violet *R*(*t→* + *n*) ≈ *R*(*t→*), (*n* ≠ 1) (Fig. 5*d*, near *t*_{1} ∼ 20 s, *t*_{2} ∼ 12.5 s, for example) demonstrate another interesting dynamical phenomenon. They indicate that the network travels through a sequence of different metastable assembly states, in which the same sequence recurs later, demonstrating the sequential repetitive nature of assembly switching dynamics. Furthermore, these switching state sequences last for periods of seconds and may therefore be behaviorally relevant. Sequential assembly dynamics in this network may be consistent with the prototype of deterministic metastable state switching in inhibitory networks (Rabinovich et al., 2000; Nowotny and Rabinovich, 2007) known as winner-less competition (WLC) (see Discussion) and more generally is type of transient (not fixed point) dynamical activity suggested by theoretical and experimental studies to be broadly relevant in neural information processing (Buonomano and Merzenich, 1995; Maass et al., 2002; Mazor and Laurent, 2005; Rabinovich et al., 2008a).

Sequentially switching episodically firing assemblies are not a consequence simply of weak lateral inhibition, however; the network also requires the sparse to intermediate connectivity observed in the striatum. The episodically firing assembly behavior at striatally realistic connectivities of 0.1 in 500-cell networks under fluctuating excitatory input shown in Figure 4, *a* and *b*, is contrasted with the dynamical behavior of the 500-cell network simulation shown in Figure 4, *c* and *d*, at the much higher connectivity of 0.82, which is much higher than the striatally realistic connectivity regime. In this simulation, all cells can be seen to be firing very independently and randomly, and few assemblies or interesting identity–temporal patterns can be discerned in either the raster plot or the cross-correlation matrix. The absolute values of the cross-correlations at high connectivity in Figure 4*d* can also be seen to be generally smaller than those at lower connectivity in Figure 4*b*. Simulations at high connectivities under fixed excitatory input also show independent randomly firing cells without assembly formation.

There is evidence that cell assemblies fire in coherent episodes in the *in vivo* striatum during periods of behavior and rest. Miller et al. (2008) found episodic firing activity in MSNs corresponding to periods of high-frequency firing was common in striatum in mice, and firing was correlated between pairs of MSNs. Coincident firing episodes, which the authors defined as bursts between pairs of neurons that overlap in time, occurred more often in pairs of MSNs that exhibited correlated firing. The appearance of episodically firing cell assemblies in our model is in good agreement with this study.

#### Effect of detailed network structure

We have demonstrated that episodically firing cell assemblies exist in random networks of inhibitory neurons at appropriate connectivities, but how do the cell assemblies depend on the details of the network structure? The network structure is shown by the black circles superimposed on the two 100 cell network cross-correlation matrices in Figures 5*b* and 6*b*. A black circle is shown if two cells are connected regardless of the variable connection strength. The violet–pink squares of positive correlation have few connections between their members. For example, the cluster of cells with index numbers from 0 to 10 in Figure 5*b* corresponding to the black assembly in Figure 5*a* have only a few connections between themselves. It is this that allows this set of cells to fire together continuously as a result of the external excitation if the rest of the MSN network that inhibits them is quiescent for a period. Other areas in which there are strong positive correlations also display few network connections, whereas areas of weak and negative correlation show many network connections.

How switching cell assemblies depend on the details of the network structure is further illustrated using a simple small network example. Figure 7*a* shows a spike raster plot from a very small nine-cell network under fixed excitatory input, which has been ordered, as above, by the clustering algorithm. This network is similar to the example described by Rabinovich et al. (2000) in their study of WLC. The individual spikes can be seen in this raster plot, demonstrating that all cells fire many spikes in an firing episode and that the cells have different firing rates during the firing episodes. Figure 7*b* shows the correspondingly ordered cross-correlation matrix with network connections superimposed. Three clusters are shown, two of four cells and one of one cell. The clusters switch perfectly periodically in this small network. Between the member cells of each cluster, strong positive cross-correlation exists, whereas between cells of different clusters, strong negative correlation is seen. In fact, every member of each cluster is inhibited by at least one member of the cluster that inhibits it. The feedforward inhibitory loop this produces is depicted more simply in Figure 7*c*.

At higher striatally realistic connectivities, such as the example whose time series is shown Figure 4*a*, connections are much more numerous and much weaker than the nine-cell example shown in Figure 7, *a* and *b*, and “net excess” connection effects, rather than individual connections, are the relevant quantities. This is depicted schematically in Figure 7*d*. Cells can fire together as an assembly even if there are connections between them, as can be seen in the network connections superimposed in Figures 5*b* and 6*b*. This is because, as described above, connections have weak strength, and the presynaptic firing rate is low because of the low level of excitatory drive to the cell. Such a weak presynaptic cell only lowers the firing rate of the postsynaptic cell and causes it to fire irregularly without quenching it completely. Furthermore, switching cell assemblies are actually a multiple timescale phenomenon. The cross-correlation matrix calculation based on a 2000 ms time bin has extracted assemblies on this 2000 ms timescale. However, within an assembly at this timescale, there can be smaller subassemblies that switch on faster timescales, and inhibitory connections can exist between these smaller subassemblies. In general, assemblies are arranged in the hierarchical structure illustrated in Figure 7*d*, and the assembly scale we focus on depends on the timescale. Another factor is that, during metastable states, cells can also occasionally transiently entrain at the spiking level and become synchronized with fixed phase relationship for a short spell, thereby continuing to fire despite active inhibitory connections (Ponzi and Wickens, 2009). This is possible in cells with close firing rates because the effect an inhibitory spike has on a postsynaptic cell depends on the postsynaptic membrane potential (Van Vreeswijk et al., 1994; Ermentrout, 1998; Izhikevich, 2005). Therefore, in general sparse random networks sets of cells with sufficiently few or sufficiently weak connections between themselves exist, and these cells can fire together as an assembly if the inhibition exerted by the rest of the network is sufficiently weak for a period. Roughly speaking, feedforward chains of assemblies, such as those described in Figure 7*d*, occur when the total inhibition exerted by all members of an upstream assembly on individual members of a downstream assembly is greater than the total inhibition exerted by all members of the downstream assembly on individual members of the upstream assembly. When the upstream assembly ceases firing, the downstream assembly will start firing after the postsynaptically bound neurotransmitter has decayed. The member assemblies of such circuits slowly switch between quiescence and episodic firing in sequence. In contrast to the simple feedforward loop described in Figure 7*c* for the small nine-cell network, when networks are larger, assembly chains become interlocked so that the slow switching of one circuit interferes with the dynamical switching of another. Furthermore, as depicted in Figure 7*d*, any given cell can be a member of several such sets of weakly connected cells. This role for sets of unconnected cells in inhibitory network clustering has also been noticed by Assisi and Bazhenov (2007). Such interlocked assembly chains with partially shared members produce sequentially switching episodically firing assembly dynamics at striatally relevant connectivities, as exemplified by the complex identity–temporal patterns in the spike raster plot in Figure 4*a* and the checkerboard relationship of assemblies in the cross-correlation matrix in Figure 4*b*.

Statistically, networks with sparser connectivities are most favorable for the formation of cell assemblies. This is because, at high connectivities, the total inhibition on all cells takes a uniform value across cells as a result of the law of large numbers, and fluctuations in cell connectivities are weak. In effect, all cells are in the same environment because they are all connected to the same set of cells that becomes an inhibitory “pool” of constant inhibition. Active cells therefore have similar dynamics without large fluctuations, as shown in Figure 4*c* at high connectivity.

### MSN firing statistics predicted by the network model with striatal connectivity

Studies of *in vivo* MSN firing have displayed a puzzling conundrum of observations. In general, cells fire with very low rates, but firing rate distributions are very broad. Cells exhibit episodic firing, but ISI distributions are only single peaked. Firing is highly irregular, showing very high coefficients of variation (CVs) of the ISI distributions and also very broad CV distributions. Furthermore, cells show cross-correlations and autocorrelations slowly varying on the timescales of hundreds of milliseconds, but precisely synchronized spiking is absent. We asked what precise quantitative predictions our network model at striatal connectivity makes for these MSN statistical characteristics and how well our model can account for the existing observations. Could the striatal MSN network itself, based on only the simplest assumptions of weak sparse lateral inhibition driven by only weak incoherent excitation and without a biophysically detailed MSN cell model, produce these nontrivial statistical observables? We wondered whether this collection of MSN firing observations could in fact derive from the tendency of the network to endogenously generate episodically firing assemblies.

First, we study firing rate and ISI distributions. Next, we investigate irregularity. Subsequently, we address coherent firing between cells, and, finally, we directly measure the quantity of episodically firing assemblies. We concentrate on the case of randomly fluctuating excitatory input because it is the most biologically relevant but also describe results for fixed input when it is instructive.

#### Highly variable low firing rates

Figure 8*a* shows the distribution of mean firing rates for each active cell for three different 500-cell network simulations with striatally relevant connectivities of 0.06, 0.1, and 0.2 in the fixed excitatory input condition and in the inset for connectivities 0.1 and 0.15 in the fluctuating input condition. These distributions are calculated from time series with a long time period of 10 min. The distributions are shown in log-log scale because they are very broad and, especially in the fixed input condition, are consistent with a power law. This broad distribution implies that cells have very variable firing rates, and there are many more cells firing very slowly than would be expected from a normal (Gaussian) distribution. In a single network simulation, most of the cells fire very slowly at <0.1 Hz, but there are also significant amounts of cells firing at 1 Hz and a few firing as fast as 10 Hz.

The mean firing rates of ∼3.5 Hz here are quite realistic (Miller et al., 2008). Very broad distributions of MSN firing rates have been observed by Wilson (1993) and in WT mice by Miller et al. (2008). Broad firing rate distributions are also observed in many other brain regions (Lee et al., 1998; Nevet et al., 2007) and predicted by modeling studies of balanced random networks (Van Vreeswijk and Sompolinsky, 1996; Amit and Brunel, 1997).

ISI distributions for some individual cells from the 10 min network simulation with striatally relevant connectivity 0.15 under fluctuating excitatory input whose firing rate distribution is shown in Figure 8*a* are shown in Figure 8*b*. The two insets show different magnifications of the same results. Although individual cells have very broad power-law-distributed firing rates the ISI distributions at large ISIs are consistent with exponential distributions, as can be seen from the linear relationships when plotted in log-linear axes (Fig. 8*b*, main, left inset). This indicates that the cells in this simulation obey Poisson processes for large ISIs, albeit with very different rates. However, Figure 8*b* shows that the ISI distributions at small ISIs deviate strongly from exponential. First, ISIs cannot be shorter than a fixed value determined by the excitation level. Network inhibition will only decrease this maximal firing rate, and therefore we expect a minimum ISI. Second, above this minimum ISI, there are many more short ISIs than would be expected for a Poisson process. This reflects the fact that cells fire in short episodes of relatively high-frequency spikes, between longer and very variable periods of quiescence. The distribution of these short ISIs is consistent with a power law across two orders of magnitude as can be seen from their approximately linear relationships when plotted in log-log axes (Fig. 8*b*, right inset). Power-law ISI distributions (Usher et al., 1994; Baddeley et al., 1997; Tsubo et al., 2009) indicate self-similarity of spike trains and can be formed from summation over many exponential distributions with different rate parameters. These ISI distributions are obtained from network cells under fluctuating excitatory input, but ISI distributions obtained from network simulations under fixed excitatory input display very similar structure.

Although burst firing is typical of striatal MSNs (Wilson, 1993; Miller et al., 2008), ISI histograms of MSNs do not generally show the bimodal distribution one might expect for burst firing, i.e., one mode corresponding to intraburst intervals and one to intervals between bursts (Wilson, 1993). It has been suggested that this is because interburst intervals are much longer and very variable compared with ISIs within bursts and so make a much smaller contribution to the ISI histogram (Wilson, 1993). The ISI distributions shown here are consistent with these observations. Poisson-like exponential tailed ISI distributions are also consistent with many experimental and modeling studies of neuronal firing patterns (Tomko and Crapper, 1974; Softky and Koch, 1993; Holt et al., 1996; Van Vreeswijk and Sompolinsky, 1996; Amit and Brunel, 1997; Shadlen and Newsome, 1998; Compte et al., 2003; Miura et al., 2007; Renart et al., 2007; Barbieri and Brunel, 2008).

#### Highly irregular firing

The firing rate distributions for individual cells shown in Figure 8*a* are the mean rates calculated over a long time period. This masks the fact that firing rates of individual cells can vary greatly within the time period.

The CV of the ISI distribution of a cell is a useful quantity to characterize firing irregularity. This quantity is defined to be the ISI SD normalized by the mean ISI (see Appendix). It is unity for Poisson processes and greater than unity for processes more variable than Poisson. The distribution of this quantity across the cells for the connectivity 0.15 network simulation under fluctuating excitatory input studied above is plotted in Figure 9*a* (black circles), denoted CV_{cell}. The cells have CVs that are broadly distributed around 1.7. This is a high value, indicating a high degree of firing variability. The distribution is also very broad, with some cells acquiring CV values as high as 2.1. Miller et al. (2008) show that MSNs in WT mice display high ISI CVs, with broad CV distributions centered around 2 or 3, very similar to our CV_{cell} distribution at striatally relevant connectivity. Our model is therefore in good agreement with this study. Such CV values, somewhat larger than unity, are also often observed in experimental (Tomko and Crapper, 1974; Softky and Koch, 1993; Holt et al., 1996; Lee et al., 1998; Compte et al., 2003) and modeling (Van Vreeswijk and Sompolinsky, 1996; Amit and Brunel, 1997; Shadlen and Newsome, 1998; Renart et al., 2007; Barbieri and Brunel, 2008) studies in various brain regions.

The high cell coefficients of variation tell us that cells have highly irregular firing. To discover whether cells are actually firing in distinct episodes or not, it is useful to measure a related quantity, CV_{2}. This local coefficient of variation (see Appendix) is defined for a cell through CV_{2}^{n} = |ISI_{n+1} − ISI* _{n}*|/(ISI

_{n+1}+ ISI

*), where ISI*

_{n}*is the*

_{n}*n*th ISI in a spike train (Holt et al., 1996; Compte et al., 2003). For Poisson processes, the distribution of this quantity is uniform and flat across the interval from zero to one. Bursty processes, conversely, show CV

_{2}distributions with peaks near zero, produced by episodes of similar successive ISIs, and peaks near one, produced by the very different successive ISIs (occurring at the starts and ends of periods of quiescence), with depletion at intermediate values. The inset of Figure 9

*a*shows the distribution of CV

_{2}for all cells combined for the same connectivity 0.15, 10 min network simulation under fluctuating excitatory input. The distribution has a peak near zero, depletion at intermediate values and a peak near one. The distribution is characteristic of episodically firing cells, and the simulation therefore contains many cells firing episodically. Such CV

_{2}distributions are often observed in neuronal studies (Compte et al., 2003), but unfortunately there is no relevant striatal study and this result is a prediction of our model.

#### Slowly varying cell correlation

We might expect that the simultaneous cross-correlation between cells in a network would provide an indication of the amount of coherent assembly formation in the network. Wang and Buzsáki (1996) in their study of interneuron networks measure a similar quantity (Wang and Buzsáki, 1996; Buzsáki et al., 2004). In that study, the authors were mainly interested in spike synchronization and therefore used a small timescale for the calculation of cross-correlation. In this study, we focus on the formation of episodically firing assemblies, in which episodes can last several hundreds of milliseconds and use a longer bin size of 2000 ms for the extraction of assemblies. If cells form coherent assemblies, there should be some large positive cross-correlations between members of an assembly and possibly large negative correlations between members of different assemblies. That this is indeed the case can be seen for the cross-correlation matrix shown in Figure 4*b*.

Figure 8*c* shows time-lagged cross-correlation coefficients (see Appendix) between some pairs of randomly selected cells calculated from rate time series with a 1000 ms bin size for the same connectivity 0.15, 10 min simulation under fluctuating excitatory input studied in the preceding section.

The large positive or negative zero time lag cross-correlation coefficients confirm the results already shown in the matrix in Figure 4*b*, but, as can be seen, large highly significant slowly varying “oscillatory” fluctuations with magnitudes of approximately ±0.1 also extend for long time lags of many tens of seconds. These oscillatory fluctuations typically have a timescale of many hundreds of milliseconds, but cross-correlation can also be significantly persistently above or below baseline on longer timescales of many seconds. This demonstrates that the network displays coherent oscillations across cells on multiple long behaviorally relevant timescales even when excitatory input rapidly fluctuates.

The inset in Figure 8*c* shows the same results calculated with a 20 ms bin size. Cross-correlation is now much weaker with magnitudes of only ±0.01 and lacking large slow fluctuations. Of course, Figure 8*c* only shows a few cell pairs of the 499^{2} in the simulation, but the absence of narrow strong peaks or troughs at distinct time lags in the 20 ms bin results indicates that cells do not show strong spiking synchronization on average.

Because the cross-correlation results shown in Figure 8*c* are calculated from a network under fluctuating excitatory input, in which the input is randomly changed every 10 ms, it is perhaps unsurprising that the 20 ms bin results do not show strong cross-correlation. However, this result is not a consequence of fluctuating excitatory input, as can be seen in Figure 8*d*, which shows cross-correlation coefficients for some randomly selected cells from a long 10 min connectivity 0.2 simulation under constant excitatory input. The 1000 ms bin size results again show persistent large slow oscillatory fluctuations with magnitude of approximately ±0.1, whereas the 20 ms bin size results again show much lower magnitude, ±0.02, fluctuations. In fact, in completely deterministic simulations under fixed excitatory input, certain cell pairs do show transient phase-locked synchronization (as mentioned above), but this occurs in short transient bursts and the effect is averaged out over long time periods. The slowly varying cross-correlation justifies our use of the long 2000 ms timescale for the extraction of episodically firing assemblies.

Our results in good qualitative agreement with studies of cross-correlation in MSNs (Stern et al., 1998), which show that correlation is apparent between up-state MSNs on broad timescales of tens to hundreds of milliseconds but not between spikes on a millisecond timescale.

Autocorrelograms of MSNs have been observed to show rhythmic low-frequency (0.5–3 Hz) patterns on long timescales but exponentially decaying autocorrelograms on short timescales of hundreds of milliseconds (Wilson, 1993; Stern et al., 1998). Plenz and Aertsen (1996) show autocorrelograms decaying on a timescale of seconds. The black dashed lines in Figure 8, *c* and *d*, show the time-lagged autocorrelation coefficient of a single cell. Consistent with observations, the autocorrelation decays from unity to near zero on the order of 1 s and subsequently shows strong fluctuations of magnitude at approximately ±0.1 on slow timescales. The smaller bin size shows the same initial decay over the first several hundred milliseconds but much weaker correlation (±0.01) at longer time lags. Again, these effects are not dependent on the fluctuating excitatory input, as can be seen in Figure 8*d* (inset) for fixed excitatory input.

#### Episodically firing assemblies

Finally, we asked whether we could directly quantify the amount of episodically firing assemblies. Two quantities are important for the formation of slowly switching episodically firing assemblies. First, individual cells must fire episodically, and second, there must be sets of cells that show excess cross-correlation on long timescales. The above results show that, at striatally relevant connectivities, both these conditions occur. To demonstrate this more explicitly, we combine the coefficient of variation analysis with the clustering analysis based on cross-correlation matrices. We apply the clustering algorithm to the cross-correlation matrices to find cell assemblies and then merge the spike time series of the individual cells to form a cell assembly spike time series, preserving the timing of each spike (see Appendix). We then measure the coefficient of variation of the cell assembly ISI distribution.

Figure 9*a* shows the distribution of cell assembly coefficients of variation for the connectivity 0.15, long 10 min network simulation under fluctuating excitatory input studied above, denoted CV_{assem} (red squares). CV_{assem} shows a broad distribution similar to CV_{cell}, with a peak at relatively high CV near 1.3, indicating the presence of episodically firing cell assemblies. Some assemblies have CV values as high as 2.0 like the individual cells do.

Figure 9*a* also shows the distribution of two control measures, denoted CV_{rand} (green diamonds) and CV_{scram} (blue triangles) (see Appendix). Both measures are calculated exactly the same way as CV_{assem} except that, in CV_{rand}, the cells are randomly collected in clusters, whereas in CV_{scram}, the ISIs of the individual cell are all initially scrambled before the clustering algorithm is applied (this procedure does not change the coefficients of variation or ISI distributions of the individual cells but does destroy the cell–cell cross-correlation). Both control distributions are quite similar, peak near one, and are much narrower than the *k*-means algorithm-generated cluster distribution CV_{assem}. This indicates than randomly assigning cells to clusters or scrambling ISIs reduces the formation of episodically firing assemblies and produces assemblies that fire in a Poisson-like way. This shows that this long 10 min network simulation at striatally relevant connectivity produces significant episodically firing assemblies.

Coefficients of variation averaged over their distributions, 〈CV〉, provide a useful single valued measure of the quantity of episodically firing assemblies in a given network simulation. The time series from a simulation at realistic connectivity of 0.1 plotted in Figure 4*a* has an average cell coefficient of variation 〈CV_{cell}〉 of 1.7, confirming the observed highly variable firing. Its average cluster coefficient of variation 〈CV_{assem}〉 of 1.47 is significantly higher than unity and the two controls 〈CV_{scram}〉 of 1.06 and 〈V_{rand}〉 of 1.14, confirming the presence of significant episodically firing assembly formation. This is in contrast to the network simulation shown in Figure 4*c* at much higher connectivity of 0.82. Its 〈CV_{cell}〉 of 0.97 is close to unity, confirming the random Poisson-like firing patterns. It has a much lower 〈CV_{assem}〉 of 0.99, close to unity, not significantly higher than the two controls 〈CV_{rand}〉 of 0.97 and 〈CV_{scram}〉 of 0.96, confirming the absence of episodically firing assemblies at high connectivity.

### Connectivity variation

We now address the question asked above: why does the MSN network have the particular sparse connectivity it does? Can the striatal MSN network subserve a WTA function, or is this particular connectivity more appropriate for the formation of sequentially switching cell assemblies and promotion of episodic firing and dynamical network coherence? Does variation in connectivity strongly affect the predictions of MSN firing patterns we observed in the previous section?

To investigate this, we perform many 500-cell network simulations while varying the connectivity. For each simulated network, we observe the spiking for a period of 50,000 ms after discarding a transient of 10,000 ms and calculate relevant quantities.

Increasing connectivity with fixed synaptic strength *k _{i,j}*

^{syn}(Eq. 3) would also increase the total inhibition on a cell. As explained in the Materials and Methods section, to isolate effects strictly relating to network structure, i.e., connectivity, we keep the total inhibition on each cell fixed on average and scale the synaptic strength parameter. This synaptic rescaling is also performed in the network simulations of Wang and Buzsáki (1996) in their study of random interneuron networks.

First, we investigate how the quantity of cells firing varies with the network connectivity. We next study how firing rates depend on connectivity. Subsequently, we study how firing irregularity varies with connectivity. Finally, we investigate directly how the formation of episodically firing assemblies depends on the network connectivity. Again, we concentrate on the case of fluctuating excitatory input because it is the most biologically relevant, but we also describe results for fixed input when it is instructive.

#### Connectivity regimes and WTA behavior

The striatum network has often been considered to subserve a WTA function, with the most strongly activated neuron suppressing activity in its neighbors (Groves, 1983; Wickens et al., 1991; Beiser and Houk, 1998; Fukai, 1999; Suri and Schultz, 1999; Bar-Gad and Bergman, 2001). However, the sparse structure of the MSN network and the fact that the inhibition exerted by a single MSN is too weak to suppress other MSNs has led several workers to conclude that this notion cannot be feasible (Tunstall et al., 2002; Plenz, 2003; Koos et al., 2004; Tepper et al., 2004; Wickens et al., 2007). Accordingly, we first investigate this question by studying how the number of cells firing varies with the network connectivity.

In Figure 9*b* (main and left inset), we show that, under fluctuating excitatory input as the connectivity increases from zero, the number of cells that fire at least once during the period of observation decreases to a minimum of ∼200 at connectivity of ∼0.01 (five connections per cell) before increasing to a plateau level, where all cells in the network fire at least once, at higher connectivities. In fact, below the transition at connectivity 0.01, connectivity is so low that the cells are divided into two types: those that have no inhibitory input from active cells at all, which therefore fire continuously without quiescent periods, and those that have at least one strong inhibitory input from an active cell that causes them to be completely quiescent for the duration of the simulation. Above the transition, inhibitory connections increase in number but have weaker strength (because of the connection strength rescaling), and it is possible for cells to fire in irregular episodes from time to time even with inhibitions from other active cells.

In the following, we concentrate on the behavior of networks under fluctuating excitatory input, but in Figure 9*b* (right inset), we digress to show that, when excitatory input is fixed, the number of active cells varies in a more complex way with connectivity. The low connectivity minimum is preserved but, rather than plateau out at a high level as connectivity increases, as in the case of fluctuating excitatory input, the number of active cells decreases from its low connectivity peak, when almost all the network cells are firing, to reach a plateau of ∼75 active cells above connectivity of ∼0.4. The reason for the difference in behavior is that, when excitatory input fluctuates rapidly, all cells experience the same average excitation when averaged on the slower timescale of changes in network inhibition. In contrast, when excitatory input is fixed, some cells receive permanently higher levels of excitation and some permanently lower. At high connectivities, the network behavior is consistent with *k*-winners-take-all (Fukai and Tanaka, 1997) model in which the active cells will generally be the ones most strongly excited, suppressing the more weakly activated ones. This WTA-type behavior at high connectivities cannot occur when all cells receive the same average levels of excitation and input fluctuates rapidly.

This also demonstrates that, at striatally relevant sparse connectivities, WTA behavior does not occur under weak fluctuating or weak constant excitatory input because, in both cases, most of the network is active. Indeed, we have chosen the excitation to be random across cells within a fixed range [*I*_{bif}, *I*_{bif} + *I*_{top}], where the input gain *I*_{top} = 1. At striatal connectivities, these weak excitation levels do not allow the network to find a permanent state of regularly firing active cells and permanently quiescent cells, but almost all cells fire episodically from time to time regardless of whether the excitatory input is fixed or fluctuating randomly. However, if excitatory input is fixed and *I*_{top} is strongly increased (*I*_{top} ∼ 100) it is possible, even in sparse asymmetric networks of striatally realistic connectivity with weak lateral inhibition, to find permanent states with some permanently active cells and some permanently quiescent cells (data not shown). However, because of the sparse asymmetric structure, it is usually not true that the *k* cells with the strongest excitations will be the *k* winners [as in a *k*-winner-take-all model (Fukai and Tanaka, 1997)], and therefore the striatal network cannot perform an unbiased WTA function even when excitation is fixed but varying across cells and strong.

#### Mean firing rate weakly depends on connectivity

Figure 9*b* (main and left inset) also shows how the mean firing rate varies with connectivity for the active cells that fire at least one spike during the course of the simulation under fluctuating excitatory input. Below the transition at connectivity ∼0.01, the cells that fire all do so with high frequency, ∼65 Hz, because of the lack of inhibitory inputs. Above the transition, the mean firing rate drops off rapidly with connectivity increase, reaching a plateau level of ∼3 Hz. The mean firing rate is fairly constant throughout the striatally relevant connectivity regime and up to high connectivities.

#### Firing becomes Poisson-like at higher connectivity

We now investigate how this observed firing irregularity varies with network connectivity. The black line in Figure 9*c* shows the mean single-cell coefficient of variation 〈CV_{cell}〉 obtained by calculating the coefficient of variation for each active cell separately and then averaging, versus network connectivity, in 500-cell network simulations under fluctuating excitatory input. As can be seen in the inset, this quantity also shows the transition at connectivity ∼0.01 from a regime of low CV to a high CV regime. Below this connectivity, almost all of the active cells are firing fairly regularly because they do not have inhibitory input from other active network cells, and they therefore have CVs close to zero. Notice that, although excitatory input fluctuates noisily, this does not produce high mean CVs when network inhibition is absent. Just above the transition, however, mean CVs get very large, indicating that most cells are firing highly irregularly. As connectivities increase further, mean CVs decrease until they approach and drop slightly below unity when the network is almost fully connected. This implies that, at high connectivity, active cells are firing in a Poisson-like way without episodes of firing and quiescence. This variation of CV with connectivity is confirmed in Figure 9*d*, which plots the CV distributions for cells for various 500-cell network simulations with varying connectivity under fluctuating excitatory input. At low connectivity, distributions are broad and have high mean CV, whereas as connectivity increases, distributions become increasingly concentrated around unity. At connectivities relevant for the striatum, average CVs range from ∼2 to 1.

The distribution of the local coefficient of variation CV_{2} for all successive ISIs of all cells combined in several 500-cell network simulations of various connectivity under fluctuating excitatory input is also plotted in Figure 9*d* (inset). At lower striatally relevant connectivities, the CV_{2} distribution is bimodal, with a peak near zero and a peak near unity characteristic of episodically firing cells. For higher connectivities, however, the distribution changes over to one resembling that of Poisson processes with an absolute refractory period (Holt et al., 1996; Compte et al., 2003). As explained above, each cell is driven by an excitation level with fixed upper limit that limits the rate at which a cell can fire, and the network inhibition will decrease the rate in general. It is this limitation that enforces the absolute refractory period in these simulations.

#### Episodically firing assemblies occur at sparse to intermediate connectivities

Finally, we investigate quantitatively how episodically firing assemblies depend on connectivity. The mean values across clusters for the measure of episodically firing assemblies described above and the two controls are plotted versus network connectivity in Figure 9*c* for 500-cell network simulations under fluctuating excitatory input. All three measures show a similar profile that is also similar to the profile for 〈CV_{cell}〉. There is a peak around connectivities of ∼0.01 and then a decrease as connectivity increases. However, 〈CV_{assem}〉 is significantly higher than either 〈CV_{rand}〉 or 〈CV_{scram}〉, from low connectivity up until connectivities of ∼0.3. Above this connectivity, all measures are quite close. Furthermore, in the low connectivity region in which the measures are significantly different, all measures are below 〈CV_{cell}〉, although 〈CV_{assem}〉 approaches close to it. At high connectivities exceeding ∼0.7, however, all the measures similarly exceed 〈CV_{cell}〉, which drops below unity. The fact that 〈CV_{assem}〉 is significantly higher than unity and the two controls up to ∼150 connections per cell indicates that cells tend to fire in episodic assemblies throughout this neurophysiologically plausible region for the striatum. In this region, individual cells also show high CV values (〈CV_{cell}〉) far above unity, suggesting that the irregular firing often observed in neural studies may facilitate the formation of episodically firing assemblies.

The time series from a simulation at striatally realistic connectivity of 0.1 plotted in Figure 4*a* is indicated by the point *A* in Figure 9*c*, and the network simulation shown in Figure 4*c* is indicated by point *B* in Figure 9*c* at much higher connectivity.

## Discussion

We propose a new view of the dynamical behavior of the striatum based on simulation of networks with realistic connectivity. Recently, the GABAergic synaptic connectivity among MSNs has been shown to mediate functional inhibition. However, the inhibitory interactions are relatively sparse, weak, and predominantly asymmetrical (Stern et al., 1998; Czubayko and Plenz, 2002; Tunstall et al., 2002; Koos et al., 2004; Taverna et al., 2004). This has led to rejection of a simple WTA model of network dynamics, leaving the role of lateral inhibition between MSNs in striatal network dynamics unexplained. We here show that random inhibitory networks of deterministic spiking neurons with sparse weak connections appropriate for the striatum exhibit a rich dynamical behavior characterized by firing in irregular coherent episodes, the formation of slowly sequentially switching cell assemblies, and complex identity–temporal patterns.

The dynamically interesting behavior is caused by lateral inhibition, which generates chaotic switching between metastable states and occurs because the cells are just above the firing threshold and because network connections are weak. This allows even sets of cells that inhibit each other to fire simultaneously but irregularly, whereas stronger fluctuations in network inhibition can cause cells to become quiescent for extended periods. This occurs despite the fact that excitatory drive is simply constant or fluctuating randomly, providing only that it is weak.

This network model is not intended to simulate detailed cell properties. The SNIC bifurcation used here allows the spiking frequency to be arbitrarily low without subthreshold oscillations, appropriately for up-state MSNs, permitting comparison with experimental studies of up-state MSN activity. However, similar coherent assembly switching dynamics could be observed with other bifurcations if connections are appropriate to allow the generation of large slow network fluctuations.

We studied how network activity depends on connectivity and found interesting transitions around realistic connectivities. Below the low connectivity transition, the small proportion of active cells all fire fairly regularly with high rates, whereas above the transition, almost all the cells fire highly irregularly. Firing becomes more Poisson-like and incoherent across cells as connectivity increases. The quantity of coherent episodically firing assemblies in the network also shows a sudden increase at low connectivity and then gradually decreases as connectivity increases. The quantity of assemblies is large at striatally relevant connectivities, suggesting that the striatum may have adapted to be in this regime.

We found that, at striatally realistic connectivities, cell firing rates are very broadly distributed, consistent with a power law, and the distribution of individual cell ISIs demonstrates irregular episodic firing with exponential tailed distributions. Striatal MSNs show strong firing irregularity (Wilson, 1993) with long periods of quiescence and high coefficients of variation (Miller et al., 2008), and such ISI distributions are often observed throughout the brain. We found that cells do not show spiking synchronization, except for short transient spells, but firing rates fluctuate coherently on long timescales of tens to hundreds of milliseconds in agreement with studies of up-state MSNs (Stern et al., 1998; Miller et al., 2008).

There have been many studies of cell assemblies throughout the brain (Cossart et al., 2003; Ikegaya et al., 2004; Brown et al., 2005; Sasaki et al., 2006b, 2007; Jones et al., 2007; Plenz and Thiagarajan, 2007; Carrillo-Reid et al., 2008; Miller et al., 2008; Pastalkova et al., 2008). Jones et al. (2007) recently observed *in vivo* sequentially switching cell assemblies involved in the representation of taste. In some studies, cells were observed to show precise spiking relationships, but in others, including *in vivo* (Miller et al., 2008) and *in vitro* (Carrillo-Reid et al., 2008) striatal studies, cells did not synchronize on precise timescales but fired episodically in assemblies on slower timescales, as in the present model.

Many theoretical studies have addressed population synchronization of coincident spiking in random inhibitory networks throughout the brain under variation of connectivity (Golomb and Rinzel, 1993, 1994; Wang and Buzsáki, 1996; Tiesinga and Jose, 1999; Golomb and Hansel, 2000). Our study complements these works, but, rather than find spiking synchronization, we have found that, in the striatum, cell assemblies in which individual cells spike irregularly but in coherent episodes on slower timescales are common. In fact, cell assemblies in which many cells fire in episodes coherently rather than synchronizing precisely may be just as (or more) effective at driving downstream targets. Carrillo-Reid et al. (2008) recently investigated the formation of striatal cell assemblies in NMDA-treated slices. The neurons fired in recurrent sequential coherent episodes, and generated spatiotemporal patterns and network transition matrices displayed abrupt transitions between different active cell assemblies. Caution should be exercised when drawing comparisons with *in vitro* NMDA-treated preparations as a result of the effects on the dynamics of the individual cells. However, by virtue of the fact that our model is minimal, in which switching cell assembly formation depends only on the assumptions of sparse weak random connectivity, proximity of cells to an appropriate bifurcation between firing and quiescent states, and fairly slow synaptic timescale, useful comparisons with preparations that should preserve these characteristics can be made. Our network model generates dynamically switching assemblies, in good agreement with the observations of Carrillo-Reid et al. (2008).

The power-law distributions we observed may be a ubiquitous feature of balanced spiking networks in which the individual cells are near bifurcation points. Studies show that balanced networks with strong excitatory and inhibitory connections display chaotic irregular firing with broadly distributed rate (Van Vreeswijk and Sompolinsky, 1996; Amit and Brunel, 1997). Critical amplification of fluctuations and power laws in neural networks has been addressed by Usher et al. (1994, 1995). In the low connectivity regime, our network can be considered a critical system in the vicinity of a phase transition in which the fluctuations arise from the underlying deterministic chaos.

The origin of the assembly switching dynamics in our model may be consistent with WLC (Rabinovich et al., 2000, 2008a). In this paradigm, neural activity does not reach a fixed steady-state attractor as in WTA but is transient, continuously switching between the vicinities of metastable states. Rabinovich et al. (2000) demonstrate that small strongly asymmetric inhibitory networks generate stimulus-dependent switching between neural ensembles. In large random networks as studied here, strong asymmetry is likely at sparse to intermediate connectivities. WLC has several computational advantages over WTA. It represents the input sensory information using both space (neural identity) and time, depends on the stimulus, and, although based on transient dynamics (Rabinovich et al., 2008a), is reproducible and robust against noise. Because of these factors, WLC has been applied to modeling of sequential firing patterns in olfactory systems (Rabinovich et al., 2000, 2001; Laurent et al., 2001), sequential decision making (Huerta and Rabinovich, 2004; Ashwin and Borresen, 2005; Rabinovich et al., 2006a,b, 2008b), and central patterns generators (Selverston et al., 2000; Varona et al., 2002; Levi et al., 2004, 2005).

The striatum is the main input structure to the BG. Coherent episodic firing activity in cortico-BG microcircuits is important in the encoding of movement sequences (DeLong, 1973; Hikosaka et al., 1989, 2000; Jaeger et al., 1995; Kasanetz et al., 2006), and the execution of learned motor programs and sequence learning (Kimura, 1990; West et al., 1990; Brotchie et al., 1991; Gardiner and Kitai, 1992; Kimura et al., 1992; Kermadi and Joseph, 1995; Mushiake and Strick, 1995; Aldridge and Berridge, 1998; Jog et al., 1999; Barnes et al., 2005). Such activity could result from coherent patterned cortical activity. However, here we have shown that excitatory input does not have to be patterned and coherent sequential episodic firing is generic even under fixed or randomly fluctuating (but weak) excitation. Why could this be useful?

The BG are also closely involved in reinforcement learning (Sutton and Barto, 1998; Doya, 2000), which requires random exploratory switching between motor sequences (Doya and Sejnowski, 1996; Barnes et al., 2005; Kao et al., 2005; Olveczky et al., 2005). Animals need to generate different movement sequences precisely in response to the same sensory cortical signals. Striatal assemblies firing in coherent episodes may generate random sequences of “macroscopic scale” fluctuations on long behaviorally relevant timescales facilitating such exploration. This is not possible in incoherent systems far from criticality in which independently fluctuating cells average themselves out. Such large slow macroscopic fluctuations are visible in the time series shown by Miller et al. (2008) (Fig. 9) for WT mice during exploration, generated by episodically firing cell assemblies. Furthermore, the endogenously generated assembly sequences we observed require weak input. If the input is strong, network activity becomes locked into input-dependent states (data not shown). It is possible that excitatory input strength is controlled, for example by dopamine, to produce stereotyped sequences when reward-motivated behavior is fully learned but maintain the cells in the vicinity of the firing threshold to promote exploratory sequence variability in early learning stages. In fact, dopaminergic modulation of sequence variability has been observed in birdsong (Sasaki et al., 2006a).

Our study shows that switching assemblies form spontaneously as part of the intrinsic dynamics of the striatal network. In the future, it is necessary to study how varying cortical input affects assembly switching and how this can be used in behavioral tasks.

## Appendix

From a network simulation, we generate spike time series for a certain period of observation for each cell *i* = 1, …, *N*, where here the number of cells *N* = 500. The period of observation here is usually *t*_{start} = 10,000 ms to *t*_{end} = 60,000, although some longer time series with *t*_{end} = 600,000 ms are also used when necessary for more accurate distributions (see below). A spike time is defined as the time the cell *i* membrane potential crosses *V _{i}*(

*t*) = −40.0 from below. Each cell

*i*therefore has an associated series of spike times, {

*t*

_{n}

^{i}}, where

*n*= 1, …,

*S*, and

_{i}*S*is the total amount of spikes fired by cell

_{i}*i*in the period of observation,

*t*

_{start}<

*t*

_{n}

^{i}<

*t*

_{end}, ∀

*n*,

*i*.

The number of cells that fire at least one spike is given by *N*_{fire} = *N* − *N*_{quies}, where the number of completely quiescent cells is given by *N*_{quies} = Σ_{i=1}^{N} δ_{0,Si} where δ is Kronecker delta.

Rate time series *R _{i}*(

*t*) for each cell

*i*are constructed so that

*R*(

_{i}*t*) =

*S*(

_{i}*t*)/

*T*

_{win}, where

*T*

_{win}is the time window length and the number of spikes fired by cell

*i*,

*S*(

_{i}*t*) in the time window is given by

*S*(

_{i}*t*) =

*n*(

_{e}*t*) −

*n*(

_{s}*t*), where

*n*(

_{s}*t*) is defined such that

*t*<

*t*

_{ne}

^{i}

*(t)*<

*t*+

*T*

_{win}but

*t*

_{ns}

^{i}

*(t)*−1 <

*t*, and

*n*(

_{e}*t*) is defined such that

*t*<

*t*

_{ne}

^{i}

*(t)*<

*t*+

*T*

_{win}but

*t*

_{ne}

^{i}

*(t)*+1 >

*t*+

*T*

_{win}. Here,

*T*

_{win}is usually set to

*T*

_{win}= 2000 ms and to construct

*R*(

_{i}*t*),

*t*runs from

*t*=

*t*

_{start}to

*t*=

*t*

_{end}−

*T*

_{win}in increments of

*T*

_{inc}= 20 ms. Sometimes for calculation of time-lagged cross-correlations, a window of

*T*

_{win}= 1000 and, for comparison, the shorter window

*T*

_{win}= 20 ms is also used. Cells that are quiescent in the observation period have

*R*(

_{i}*t*) = 0, for

*t*

_{start}<

*t*<

*t*

_{end}.

Network state transition matrices *D*(*t*_{1}, *t*_{2}) at two different times *t*_{1} and *t*_{2} are constructed from dot products of normalized rate vectors, *D*(*t*_{1}, *t*_{2}) = *R*(*t→*_{1}) × *R*(*t→*_{2})/(|*R*(*t→*_{1})| |*R*(*t→*_{2})|), where *R*(*t→*_{1}) is the vector of cell rate time series *R _{i}*(

*t*), and |

*R*(

*t→*)| is the length of the vector. For these calculations, a window of

*T*

_{win}= 40 ms is used, incremented in steps of

*T*

_{inc}= 40.

*D*(

*t*

_{1},

*t*

_{2}) varies from 0, meaning the network state at times

*t*

_{1}and

*t*

_{2}is very dissimilar, to 1, meaning the network state at times

*t*

_{1}and

*t*

_{2}is identical. The color scale for plots in Figures 5 and 6 has been rescaled to show 2D-1 for presentation clarity and goes from −1 to 1 in 11 equal-sized steps.

Cell–cell cross-correlation matrices *C _{ij}* are constructed from the rate time series

*R*(

_{i}*t*) as follows: where 〈

*R*(

_{i}*t*)〉 denotes expected value, 〈

*R*(

_{i}*t*)〉 = (1/

*J*)Σ

_{j=0}

^{J}

*R*(

*t*

_{start}+

*jT*

_{inc}), where

*J*is defined such that

*t*

_{start}+

*JT*

_{inc}<

*t*

_{end}and

*t*

_{start}+ (

*J*+ 1)

*T*

_{inc}>

*t*

_{end}. Only cells that fire at least one spike are included in the cross-correlation matrix calculation. The color scale for plots in Figures 4⇑–6 goes from −1 to 1 in 11 equal-sized steps.

Time-lagged cross-correlations and autocorrelations are defined as in the correlation matrices:
where τ is the time lag. [The quantities 〈*R _{i}*(

*t*+ τ)〉 and 〈

*R*(

_{i}*t*)〉 may be slightly different at long time lags τ because of the finite size observation period.]

ISI time series are defined by *I*_{n}^{i} = *t*_{n+1}^{i} − *t*_{n}^{i}. Cell *i* mean ISI is defined by *m*_{I}^{i} = 〈*I*_{n}^{i}〉 = (1/(*S _{i}* − 1)) Σ

_{n=1}

^{Si−1}

*I*

_{n}

^{i}; it is only defined for cells that fire at least two spikes,

*S*> 1, during the observation period. Cell

_{i}*i*ISI SD is defined by σ

_{I}

^{i}=

*S*> 2, during the observation period. The coefficient of variation CV

_{i}*of a cell ISI distribution is defined by CV*

^{i}*= σ*

^{i}_{I}

^{i}/

*m*

_{I}

^{i}when cells fire at least three spikes. The average CV for a network is defined by CV = (1/

*N*

_{fire′}) Σ

_{i}

^{N}

_{fire′}CV

*, where the sum that runs over the*

^{i}*N*

_{fire′}is the cells that fire at least three spikes.

Local coefficients of variation CV_{2} are defined through the quantity
which varies between zero and one. The cell CV_{2}^{i} = (1/*S _{i}* − 1)) Σ

_{n=1}

^{Si−1}CV

_{n}

^{i}is defined by averaging over

*n.*The network CV

_{2}is defined by averaging over all active cells CV

_{2}= (1/

*N*

_{fire′}) Σ

_{ifire′}

^{N}CV

_{2}

^{i}as above. Distributions of CV

_{2}are formed by combining all cells

*i.*

To study the clustering, we apply the *k*-means algorithm to the cross-correlation matrix *C _{ij}* to order the cell indices. To each of the

*N*cells

_{fire}*i*, we associate the vector

*C→*=

_{i}*C*

_{i}_{1},

*C*

_{i}_{2}, …,

*C*

_{iN}_{fire}of cross-correlation coefficients. We chose

*N*centroids

_{K}*K→*, where

_{k}*k*= 1, …,

*N*, which are also

_{K}*N*

_{fire}dimensional vectors. These are initially chosen randomly from the

*N*

_{fire}vectors

*C→*. We associate each cell

_{i}*i*to a single centroid

*k*, which is the centroid nearest to the cell, i.e., min

_{k ϵ NK}||

*C→*−

_{i}*K→*||

_{k}^{2}. Denoting the centroid to which cell

*i*is associated by

*k*= 1, …,

_{i}*N*, we then form new centroids as the centers of their associated cluster of cells,

_{K}*K→*= Σ

_{k}*δ*

_{i}C→_{i}_{k1,k}/Σ

_{i}δ

_{k1,k}. We then repeat the operation by associating each cell

*i*to the nearest of the new centroids and calculate the new centroids again. The algorithm is repeated until it halts when no cells change their associated centroids. This minimizes the total distance Σ

*||(*

_{i}*C→*−

_{i}*K→*

_{k1})||

^{2}. The cells are then reordered according to their clusters, so that the cluster of cells

*i*with

*k*= 1 come first, followed by the cluster of cells with

_{i}*k*= 2, and so on up to the cluster of cells with

_{i}*k*=

_{i}*N*. Within each cluster, cells are further reorganized according to their within-cluster cross-correlation coefficients. That is, we first select the cell with the maximum cross-correlation coefficient in a given cluster, then find the cell with which it has maximum cross-correlation within that cluster, and then find the cell in the cluster with the maximum cross-correlation to them and so on. We start with cell

_{K}*i*, given by max

*, where*

_{i,j}C_{ij}*i*and

*j*index the cells in a given cluster. Cell

*i*receives index 1 and cell

*j*index 2, then we look for cell

*k*defined by max

*(*

_{k}*C*

_{1k}+

*C*

_{2k}), and this cell receives index 3 and so on. The cells are thereby reorganized.

Cluster spike time series {*t*_{n}^{k}}, where *k* labels the cluster, are defined by combining all the spike times of the member cells, {t_{n}^{k}}=U_{ijk1=k}{t_{n}^{i}}. Cluster rate time series are calculated from these cluster spike time series in the same way as for individual cells described above. CV_{assem} are calculated from cluster spike time series in the same way as the single-cell CVs, and their distributions can be studied in the same way as single-cell CVs. Network cluster coefficients of variation can also be calculated by averaging the individual cluster CV_{assem}.

To test for significance of various quantities, we compare their values with those obtained after randomly scrambling the ISIs. That is, we create the ISI time series defined by *I*_{n}^{i} = *t* _{n+1}^{i} − *t*_{n}^{i} and then scramble the indices *n* = 1, …, *S _{i}* − 1 to make a new set

*n** (

*n*) and then make the new spike time series

*t*

_{n}

^{i},* =

*t*

_{1}

^{i}+ Σ

_{j=1}

^{Si−1}

*I**

_{n}_{(j}

_{)}. We then calculate the rates

*R*

_{i}*(

*t*) and the cross-correlation matrix

*C*

_{ij}* exactly as described above. The scrambling operation produces randomly ordered spike time series that nevertheless preserves quantities such as the mean rates, 〈

*R*(

_{i}*t*)〉 = 〈

*R*

_{i}*(

*t*)〉 and does not alter the distribution of ISIs but destroys cross-correlation between cells. The

*k*-means algorithm can then be applied to the new cross-correlation matrix to produce clusters from scrambled ISI time series. Cluster spike time series are created as above, and their coefficients of variation are calculated. The resultant cluster coefficients of variation CV

_{scram}can be studied, and they can be averaged to produce the network scrambled ISI coefficient of variation. As another method of significance testing, we find clusters in the normal way by applying

*k*-means algorithm to the unscrambled spike time series, but, before creating cluster time series, we scramble the cell indices, randomly associating cells to clusters but preserving cluster sizes. That is, we use the

*k*-means algorithm to obtain appropriate cluster sizes but then associate the cells to the clusters randomly. It is important that we preserve the correct cluster sizes because larger clusters obviously produce denser cluster spike time series with shorter ISIs. From the resultant cluster spike time series, we obtain their coefficients of variation CV

_{rand}and then the average of those to obtain the network random coefficient of variation.

To refine these methodologies further, we note that the centroids found by the *k*-means algorithm are not unique and depend on the initial choice of centroids. As mentioned above, we chose these at random from the *N*_{fire} vectors *C→ _{i}*, and, similarly for the randomized calculation, we chose the initial centroids from

*C→*

_{i}*. The

*k*-means algorithm relaxes until it has found a minimum of the total distance as described above. However, this need not be a global minimum and is likely to be local minimum. Therefore, different choices of initial centroids lead to different final centroids. To control for this, we perform the calculation many times on the correlation coefficients calculated from both the nonscrambled and the scrambled spike time series. Each time, we chose a different set of initial centroids. Each time, we calculate the cluster spike time series and the cluster coefficients of variation. From the nonscrambled ISI time series, we obtain many samples of CV

_{assem}and CV

_{rand}, and, from the scrambled ISI time series, we obtain many samples of CV

_{scram}and find their distributions. Then we average the results to obtain 〈CV

_{assem}〉, 〈CV

_{rand}〉, and 〈CV

_{scram}〉, where 〈.〉 denotes the average over the many different initial centroids.

## Footnotes

This work was supported by the Okinawa Institute of Science and Technology Promotion Corporation.

- Correspondence should be addressed to Jeff Wickens, Neurobiology Research Unit, Okinawa Institute of Science and Technology, 1919-1 Tancha, Onna-Son, Okinawa, 904-0412, Japan. wickens{at}oist.jp