**FIGURE 6 | Detection can be improved by including past activity and weighting neurons depending on how many inputs they receive from directly stimulated neurons. (A,B)** Density representation of the firing rate return map, wherein the probability of obtaining consecutive rate values (r*t*, r*t*+1) is represented by a color scale, with red indicating the highest probability and blue indicating a near zero probability. The results are shown for **(A)** spontaneous activity and **(B)** stimulated activity. The plusses indicate the location of the peak in panel **(B)**. **(C)** Analysis of factors that contribute to a neuron's weight in detection decision that is outputted by the perceptron procedure. There was significant but small correlation between weight and (top left) in-degree or (top right) out-degree. There was a correlation between the weight and the number of direct inputs from stimulated neurons (bottom left), but only a weak correlation with the number of inputs from cells that received direct input from the stimulated cells (bottom right). The network was comprised of *N* = 2000 neurons, coupling constant *J* = 18, baseline firing rate *r*0 = 1 Hz. In the stimulated network, *np* = 8 neurons were stimulated for a duration of *T*stim = 6 time bins.

cells that are not directly connected to the stimulated neurons would display firing rate fluctuations that are unrelated to the stimulation, hence act as noise that reduces probability of detection. The signal to noise of the firing rate fluctuations could be improved by weighing those neurons less. To explore this hypothesis we applied a perceptron procedure (see Methods) to learn the optimal weights for classifying the network state vectors (Duda et al., 2001). An equal number of network states for spontaneous activity and for stimulated networks were supplied to the perceptron routine together with the corresponding class labels. The output was a weight for each neuron. As before the activity of the directly stimulated neurons was not included in this analysis. To determine what features contributed to the weight we plotted the weight vs. feature value in a scatter plot and calculated the corresponding Pearson correlation. There was a small, but significant correlation between the weight and the in-degree (**Figure 6C**, top left, correlation 0.073 ± 0.03, *p* = 0.0011) and with the out-degree (**Figure 6C**, top right, correlation−0.052 ± 0.027, *p* = 0.018). There was a strong correlation between the weight and the number of direct inputs the neuron received from stimulated neurons (**Figure 6C**, bottom left, correlation 0.410 ± 0.15, *p* = 0.0). The number of indirect inputs from stimulated neurons was less relevant (**Figure 6C**, bottom right, correlation 0.072 ± 0.032, *p* = 0.012). We calculated this by determining the number of inputs from cells that received direct inputs.

Taken together, these analyses show that our estimates for the detection of stimulation based on overall firing rate are underestimates and can be improved by taking into account network history and by selecting which neurons to listen to. The latter of which may be achieved through synaptic plasticity and the appropriate learning rules.

## **DETECTING ANTI-CORRELATION IN THE DEGREE DISTRIBUTION WITH LIMITED DATA**

The results here establish that anti-correlation between in- and out-degrees results in more stable, but equally sensitive networks, compared to networks without correlations between in- and outdegree, or positive correlations between them. Hence, learning to detect a stimulation could proceed by altering the correlation between in- and out-degree. To demonstrate such a learning effect, the in- and out-degree of a number of neurons needs to be sampled. Classical tracing techniques are not appropriate because they involve the connections to or from multiple nearby neurons (Lanciego and Wouterlood, 2011). For instance, when the retrograde tracer horseradish peroxidase is injected, it is absorbed by multiple axon terminals and transported to their respective cell bodies. These axon terminals do not necessarily synapse on one and the same neuron near the injection site. Hence, the data cannot be used to determine the in-degree of a neuron near the injection site.

New viral-based techniques could help, because they work by infecting a few cells in the neighborhood where the virus is injected (Wickersham et al., 2007; Osakada et al., 2011). The virus will then retrogradely label the cells presynaptic to these cells by crossing one synapse and one synapse only. In the presynaptic cells the infection stops because the virus misses the proteins necessary to cross another synapse. The challenge with this method is to infect only one cell, with both an anterogradely and retrogradely crossing virus.

Currently, the gold standard is to simultaneously record multiple cells *in vitro* and assess connections by inducing action potentials in one neuron at a time and recording the post-synaptic responses in the other cells. The current record is 12 cells recorded simultaneously (Song et al., 2005; Perin et al., 2011). This means that the anti-correlation in the degree distribution will have to be assessed indirectly, by sampling from sub-networks.

Motifs represent patterns in the connectivity that occur more often than expected if the connections were made random (Milo et al., 2002). For instance, consider a network for which the average probability of a connection is *p*. For two neurons, if these connections are made randomly, the probability of having no connection is (1 − *p*)2, for having one connection 2*p*(1 − *p*) and for having a bidirectional connection *p*2. When it is found that bidirectional connections occur significantly more than the expected *p*2 then there is additional, non-random, structure in the network (Song et al., 2005). Motifs most often refer to triplets of neurons and the patterns of connectivity between them that occur more often than expected in a random network (Milo et al., 2002). A motif distribution is the number of times each motif occurs in a network and a motif is considered present when it occurs more often than in a control network. Motif distributions are affected by many network properties such as, for instance, the degree distribution. The networks studied here, even when uncorrelated, have a different degree distribution than the ER network, which means that ER random networks are not a good control. Hence, we have to numerically generate the control distributions rather than having access to the analytical expression for the expected rate of each motif. In addition, in experimental settings we do not have access to the whole network from which to determine the motif distribution, we have to do with sub-networks. These subnetworks do not come from the same network, rather they come from networks sampled from an ensemble of networks with similar properties. To obtain estimates for how to observe evidence for anti-correlation in the degree distribution we need to deal with each of these issues.

The overall goal is to distinguish between pairs of networks with anti-correlated, uncorrelated and correlated degree distribution with the same marginal distribution for in- and out-degree.

We considered 13 different motifs that consisted of three connected neurons and gave each motif a numerical label as shown in **Figure 7A**. We determined the number of motifs in each realization of a network with correlated, anti-correlated or uncorrelated degree distribution and took the average. This was done for the full network (here reduced to *N* = 200) as well as for subnetworks (size *N*sub). The complexity of a motif corresponds to the number of edges in the pattern, ranging from 2 to 6, which determines how often it is counted in a network. We normalized the counts such that they took values on the order of unity in order to better compare them across motifs. The mean count as a function of *N*sub converged to a constant for network sizes between 50–100 neurons (**Figure 7B**), with more complex motifs requiring larger *N*sub. The width of the count distribution, quantified as the standard deviation, decreased with *N*sub as the-3/2 power (**Figure 7C**). Hence, for large enough networks the differences in mean counts across network type can be detected with certainty. This power law behavior is consistent with the results for a Binomial process with probability *p* and on the order of *n* ∼ *N*3 sub trials, for which the mean is *np* and the variance is *np*(1 − *p*). In that case the normalized mean is *p*, and its variance (*1-p*)*/n* (see also Equations 2, 3), leading to a standard deviation varying as *n*−1/2 = *N*−3/2 sub .

We are looking for motifs whose counts are different between the analyzed network types. In **Figure 7D** we show the count as a function of motif for the three network types, with the highest difference occurring for the complex motifs 110 and 238. However, the counts for these motifs, which have the largest number of edges, are characterized by a large standard deviation. When we plot the motifs in an order based on the ratio of count difference over standard deviation, motif 98 comes up as winner instead (**Figure 7E**). **Figure 7D** shows that there are fewer motif 98 in

**FIGURE 7 | Motif 98 is the most sensitive to degree correlations. (A)** There are 13 motifs that involve 3 connected nodes. Below the graphical representation we plot the numbering used here, which follows Itzkovitz et al. (2003). The expected number of motifs depends on network size, hence we normalize the count by *N*3(*k*/*N*)*e*, with *N* the number of nodes, *k* the expected number of edges per node and *e* the number of edges in the motif. In addition, we include a numerical factor representing the equivalent permutations [listed in Table 3 in Itzkovitz et al. (2003)]. **(B)** The normalized counts, averaged across a thousand realizations, converge to constant values for sub-networks larger than 50–100 nodes, with the precise value depending on the complexity of the motif involved. **(C)** The standard deviation of the normalized counts fall off as *N*−3/2. We illustrate the results for the anti-correlated network, which are typical for the correlated and uncorrelated network also. In addition, we omitted motif 238 because it occurs at such a low probability that it makes the statistics noisy. **(D)** The normalized counts for each motif for the (red) correlated, (blue) anti-correlated and (green) uncorrelated networks. We used the counts for the full network, rather than sub-networks. Network size in panel **(D)** and **(E)** was *N* = 200. We used a bivariate Gaussian degree distribution with a mean number of nodes equal to 10, a standard deviation along the long axis of σ*y* = 3.33 and along the short axis of σ*x* = 1.0. **(E)** The maximum difference in mean count between all three possible comparisons (black bars), relative to the mean standard deviation of these counts across the three network types. The motifs are ordered on the count over standard deviation ratio, starting with the largest. According to this analysis motif 98 should be used to best distinguish between different network correlation structures.

anti-correlated networks compared to correlated networks. This can be understood intuitively by noting that in "ring" motif 98 each neuron has the same number of inputs as outputs, namely 1, which is more representative for correlated networks (**Figure 1D**, inset) than for anti-correlated networks (**Figure 1C**, inset). Furthermore, this means there is a lower probability of closing the ring, because in an anti-correlated network a neuron with many inputs has fewer outputs to get to the next neuron in the ring.

The count distributions are not Gaussian for small subnetworks. **Figure 8A** shows the count distribution for motif 98 for networks with *N* = 200. Each network gives rise to a symmetric appearing distribution, with the peak at a different location depending on the network type. The distribution for the anticorrelated and correlated network were farthest apart, with that of the uncorrelated distribution situated in the middle. For *N*sub = 30 (**Figure 8B**), the corresponding distributions fell on top of each other and are asymmetric because the counts are always positive. To compare the distributions we therefore performed an ROC analysis. As expected based on the reduced overlap between distributions, the AUC increases with sub-network size, and motif 98 comes out on top with the highest AUC (**Figure 8C**). Furthermore, given the lower overlap between the anti-correlated and correlated distribution (**Figure 8A**), the AUC values for the comparison between anti-correlated and correlated network is higher (**Figure 8C**) than for either the comparison between anti-correlated and uncorrelated (**Figure 8D**) or correlated with uncorrelated (not shown).

In experiments only relatively small networks can be mapped, up to 12 cells using paired recordings and a few tens to hundreds using population calcium imaging. For these numbers the degree correlations cannot be reliably distinguished based on a single measurement. We therefore pooled measurements to see if this improved discriminability for more experimentally accessible smaller sub-networks. This procedure (pooling motif counts across *Nav* = 50 network realizations) indeed reduced overlap between distributions (**Figure 8E**, compare to **Figure 8B**). The more motif counts were pooled, the higher the AUC was (**Figure 8F**). Furthermore, the value of unity, corresponding to perfect discriminability is reached for smaller sub-network sizes. For *Nav* = 50, *N*sub = 30 networks are perfectly discriminable and the AUC transitions from values just above 0.7 to unity between *Nav* = 30 and 50 (**Figure 8F**). Taken together, subnetworks of a few tens of neurons could be used to test our hypothesis experimentally.

The question is whether this result can be improved by including counts for multiple different motifs (**Figure 9A**). Without pooling, motif 98 by itself outperforms any pair of motif counts, according to the AUC value (**Figure 9B**). To determine the AUC value for pairs of motif counts we used the FCM procedure as outlined in the methods section. When counts are pooled (**Figure 9C**), some motif pairs outperform motif 98 by a small margin. The pairs are highlighted in **Figure 9D**, and involve motif 98 itself. The more separated the cloud of points corresponding to different network types is, the better the FCM procedure classifies the networks, compare the plusses (correct discrimination) and dots (incorrect) in **Figure 9A**.

## **DISCUSSION**

The overall firing rates in barrel cortex (de Kock et al., 2007; Greenberg et al., 2008; Barth and Poulet, 2012) are much lower than might be expected based on the classic experiments in macaque visual cortex (Hubel and Wiesel, 1968). Neural activity is also variable, which can be characterized as across trial reliability, or in terms of the coefficient of variation and Fano factor of spontaneous activity (Shadlen and Newsome, 1998). These measures reveal that the activity is similar to that of a Poisson process, in which the occurrence of a spike in a time bin is uncorrelated with whether or not a spike occurred in previous bins. The mean

**FIGURE 8 | Degree correlations can be distinguished by pooling fifty measurements of networks at least thirty neurons in size. (A)** Motif count varies across network realizations, but degree correlations can be distinguished when the corresponding distribution show little overlap. We show the distribution of the normalized counts of motif 98 for (red) correlated, (blue) anti-correlated, and (green) uncorrelated networks with 200 neurons. **(B)** For smaller sub-networks (*N*sub = 30), the distributions overlap. Furthermore, these distributions are not Gaussian as they are skewed because counts are always positive. Hence, a more general procedure, such as the ROC analysis needs to be used instead of looking at the differences in mean count relative to the standard deviation. **(C,D)** The area under the ROC curve (AUC) as a function of sub-network size *N*sub for the comparison **(C)** between correlated and anti-correlated networks and **(D)** between anti-correlated and uncorrelated networks. Each motif is labeled with a line style and color as indicated in the legend. Motif 98 is most sensitive in both cases (as well as for the correlated vs. uncorrelated comparison that is not shown). It is more difficult to distinguish an anti-correlated network from an uncorrelated one than to distinguish it from a correlated network. The average AUC values were determined based on the AUC value for each of twenty different motif distributions of 500 network realizations, which were sampled randomly with replacement out of 1000 realizations. **(E)** The motif distribution for *N*sub = 30 can be pooled across *Nav* = 50 network realizations in order to shrink the width of the distribution, so that the differences in mean counts become clearer [(compare to panel **(B)]**. **(F)** The AUC for larger *Nav* values reaches unity (distributions are perfectly distinguishable) for smaller sub-network sizes. We show (green) no pooling, (blue) pooling across *Nav* = 5 realizations and (red) pooling across *Nav* = 50 realizations. The AUC goes from 0.7 to 1.0 between *N*sub = 30 and 50 when pooled across *Nav =* 50 realizations, indicating that networks of size 30 can be used to determine degree correlation structure.

firing rate is maintained by the intrinsic excitability of neurons and their synaptic inputs, including recurrent excitation. High variability together with a low firing rate implies that the network dynamics should be stable against fluctuations in the mean activity in the sense that these fluctuations do not generate states

**FIGURE 9 | The gain in discriminability by using the joint distribution of motif counts rather than the marginal distribution is limited. (A)** The outcome of the FCM procedure to find two clusters. Red points indicate counts obtained from a correlated network and blue points are those for the anti-correlated network. The plusses indicate points correctly classified by FCM and the dots represent incorrectly classified network realizations. Each motif-pair ROC curve is obtained by applying the ROC analysis to the FCM generated probability of each network to belong to cluster one (see methods). **(B)** The resulting AUC as a function of sub-network size for motif 98 (red) is higher than for all possible pairs of motifs (black curves). The AUC shown is for comparing anti-correlated and correlated networks. **(C)** When the motif count is pooled across *Nav* = 20 realizations some two-motif curves exceed the single motif curve, which are shown separately in panel **(D)**. This suggests that for a specific size of the sub-network distinguishability can be improved by considering pairs of motif counts.

with networks bursts in which all neurons in the network are active at the same time.

Experiments show that rodents can detect single-cell stimulation in barrel cortex, in which a single neuron is electrically stimulated to produce a high-frequency train of action potentials (Houweling and Brecht, 2008). This may mean that single-cell stimulation can cause an increase (or decrease) in the firing rate of the local network that is significantly different from that occurring during spontaneous activity. Taken together, this means that cortical networks with a low firing rate should at the same time be stable against fluctuations in firing rate and sensitive to weak stimulation. The overall goal of this paper was to a find a potential explanation for how the contrasting demands of sensitivity and stability can be realized. To achieve this we examined the dynamics of binary neural networks with correlation between the inand out-degrees of neurons. In the following we summarize the main results with the aim of linking the detection performance of the network to experimentally obtained behavioral results, the mechanism by which sensitivity and stability can be achieved, and predicting the anatomical signatures of the hypothesized network. We also discuss the role of other biophysical factors, such as inhibition, not taken into account in the present study.

## **GENERATING THE NETWORK CONNECTIVITY UNDERLYING ENHANCED STABILITY**

Our guiding hypothesis is that networks with an anti-correlation between the in- and out-degree of neurons are more stable and equally sensitive as other networks with comparable marginal degree distributions. Network sensitivity is generated by neurons with a high out-degree, because these would amplify the effect of nanostimulation the most. This amplification would also destabilize the network, so these cells should not be activated during spontaneous activity. As the input to neurons is proportional to the mean firing rate and their in-degree, this can be achieved by making sure that high out-degree neurons have low in-degrees. To maintain the average degree, both in and out, there then also need to be neurons with a low out-degree and a high in-degree. We implemented this hypothesis as an anti-correlation in the inand out-degree.

In the standard Erdos-Renyi networks, the relative variance in the degree distribution for large networks becomes too small to have out-degrees that are much larger than the mean degree, which is needed to reach the desired sensitivity. Hence, we needed to broaden the degree distribution artificially by using a truncated bivariate Gaussian distribution. Networks with this sampled degree distribution were generated via the configuration model (Newman, 2010). This configuration model generates networks with self-edges and multi-edges. Analytical calculations show that the probability for obtaining a network with one or more of these edges is close to one for the large mean degrees we consider (Blitzstein and Diaconis, 2006). Nevertheless, the number of these edges is low and their impact on the dynamics was limited.

There are a number of ways to address the multi and self-edge problem in a more principled approach that differ in computational efficiency and ease of implementation. First, one can use the configuration model procedure, but reject an invalid edge and find a valid replacement. This carries the risk that the algorithm stops when there are no valid edges available, which means that the whole procedure has to be restarted. Alternatively, as mentioned before, one can identify the invalid edges when the network construction has been completed and remove them or replace them by valid ones. See Blitzstein and Diaconis (2006) for a review. Second, one can find one graph that satisfied the degree distribution using the Havel-Hakim procedure (Viger and Latapy, 2005; Erdos et al., 2010; Chatterjee et al., 2011) and generate samples from the overall graph distribution by swapping links (Blitzstein and Diaconis, 2006). Swapping links refers to the procedure where randomly chosen existing links *i* → *j* and *k* → *l* are swapped into *i* → *k* and *j* → *l* when this yields a simple graph without self-edges and multi-edges. This requires careful calibration of the number of swaps and also introduces bias because these swaps do not change the number of triangles in the network (Roberts and Coolen, 2012). Third, a sequential method can be defined that produces all possible graphs, by randomly selecting amongst the allowed edges that keep the residual degree distribution graphical (Del Genio et al., 2010; Kim et al., 2012). A degree distribution is graphical when there exists a simple graph with that distribution, after each step the degree distribution is lowered to account for the connections realized, and this is referred to as the residual degree distribution. This method does not produce the graphs with the correct probability. Hence, averages based on these graphs have to be reweighted to take this into account. Furthermore, in our hands, an implementation of this method produces graphs with a correlation between the in- and/or out-degrees between different nodes, which is referred to as assortativity. This necessitates a number of link swaps to remove these correlations. Fourth, edges can be sampled according to a Boltzman function (Park and Newman, 2004), where the expectation value of the degree of a node is fixed through a Lagrange multiplier, for which the appropriate value has to be picked, which can be achieved, for instance, through a maximum likelihood approach or iterative rescaling (Chatterjee et al., 2011). Taken together, we opted to use the simplest method here, because these alternative methods for network generation were computationally more intensive and also suffered from aforementioned additional drawbacks, such as graphs that were not sampled according to a uniform probability (Del Genio et al., 2010) or other biases in the network statistics (Roberts and Coolen, 2012). Recently developed methods for generating networks with degree correlations, both in a single neuron as well as between pairs of neurons look very promising (Roberts and Coolen, 2012).

## **STABILITY IS ENHANCED WHEN THE IN- AND OUT-DEGREE ARE ANTI-CORRELATED**

Our aim was to find stable networks, by which we mean that fluctuations do not cause a cascade of recurrent excitation resulting in all cells being active at the same time. One solution would be to have inhibitory neurons, but this does not affect the stability of the LFS, it just changes the ultimate level of activity reached (Avermann et al., 2012). Stability can be assessed in a number of different ways. First, stability in the nonlinear system sense: is the LFS a fixed point of a noise-less, infinite size system? We determined that there was a range of coupling strengths *J*, below *Jc*, for which such a LFS exists. The higher the baseline firing rate, the smaller that range is. Finite-size systems have a smaller range of stable coupling strength, because there is heterogeneity, not every neuron has the same in-degree. For instance, the uncorrelated network had a higher variance in the degree distribution than the ER network, and also had a smaller *Jc*. Interestingly, networks with a positive correlation between in- and out-degrees reduced stability even more, leading to a lower *Jc*, whereas for networks with a degree anti-correlation, *Jc* was higher, even exceeding the value for the ER network of the same size.

These calculations ignore the effects of fluctuations, which we subsequently introduced by making the dynamics stochastic. This did not alter the stability as determined before in terms of the existence of the LFS, but introduced other features. The LFS has a BOA with a fuzzy boundary due to the stochastic dynamics. A network can then be unstable when the fluctuations are large enough to leave the BOA when you wait long enough. This is primarily a concern for *J* values close to (and below) *Jc*. We determined the fraction of trials during which the network left the LFS BOA during the simulated time interval. As expected the anti-correlated network is more stable, because *Jc* is larger. For coupling constants away from *Jc*, this way of characterizing the BOA does not work. Hence, we started the network in states with many more neurons active than would be expected as a result of any normal fluctuation, and determined whether it converged to the LFS or HFS. This revealed that the BOA was larger for the anti-correlated network even away from *Jc*.

Taken together, these results clearly show that anti-correlated networks are more stable than uncorrelated ones, which means they can operate stably at higher coupling strengths and baseline firing rates, which confers advantages when the sensitivity is higher for higher coupling strengths and baselines rates. Furthermore, their sensitivity is enhanced compared to ER networks with the same connection probability, because of a subset of neurons with a high out-degree.

Recent experiments summarized in Barth and Poulet (2012) show that the average firing rate in sensory cortex is low, especially in superficial layers. This holds for spontaneous as well as evoked activity, and for both anesthetized animals and awake animals and is the basis for the parameter settings in the model. Nevertheless, there is a small subset of cells that display high firing rates. Cells that have recently been active express the immediateearly gene *c-fos*. When the *c-fos* promoter is used to express the fluorescent marker GFP, the recently active cells can be targeted for recording *in vivo* and *in vitro*. The so called fosGFP+ cells had a higher firing rate both *in vivo* and *in vitro* and received more excitatory inputs and less inhibitory inputs (Yassin et al., 2010). Furthermore, these cells are more likely to be connected amongst themselves. In the anti-correlated networks, there are neurons with a high in-degree but a low out-degree which make the network more stable, and neurons with high out-degree but low in-degree that make the network more sensitive. The fosGFP+ neurons could correspond to the former group, which form the backbone for the spontaneous activity. We did not explicitly build in assortativity in the network to preferentially connect high indegree neurons to each other as suggested by Yassin et al. (2010). We take from this result that the prevailing homeostatic processes create networks with more strongly connected sub-networks and produce cell-to-cell heterogeneity in the balance between excitation and inhibition. Training to detect electrical stimulation should thus be able to induce similar changes in network structure.

## **THE SENSITIVITY ESTIMATED USING DIFFERENT MEASURES OF NETWORK ACTIVITY**

Rodents were able to distinguish between patterns of neural activity during spontaneous activity and those caused by single-cell nanostimulation. Nevertheless, this distinction was small, given the effect size measured experimentally (Houweling and Brecht, 2008). One hypothesis is that the total amount of activity (firing rate) due to nanostimulation significantly exceeds that expected of a typical fluctuation. For a stationary network dynamics, this implies a fixed threshold above which a fluctuation is more likely caused by nanostimulation, whereas fluctuations below the threshold are more likely due to spontaneous activity. This can be quantified using a ROC curve, and the area under it, the AUC. The ROC is the curve traced out by varying this threshold and plotting the true positive rate (nanostimulation above threshold) vs. false positive (spontaneous fluctuations above threshold). When both distributions for the fluctuations are Gaussians, the AUC corresponds to the difference in means divided by the (common) standard deviation (Kingdom and Prins, 2010). Hence it is a measure of the difference in response relative to the size fluctuations around it. We found that the main determinant of the AUC is the out-degree of the stimulated neurons, independent of the correlation between in- and out-degree in the network. The AUC increases with coupling strength and baseline firing rate. The anti-correlated network has an advantage because it allows for a broader range of *J* and *r*0 values. It thus has an increased stability at equal sensitivity.

The above represents an underestimate of the sensitivity, because it assumes that the activity of each neuron contributes equally to the detection (decision) and that the temporal signature of the firing rate fluctuation is not informative. Our further analysis shows that each of these factors would improve detection performance and makes it therefore likely that state-of-the-art classification approaches such as support vector machines would even further improve performance. Taken together this means that as a system the rodent brain could reach a much higher sensitivity than predicted here, when it could utilize all the information available in the network activity. Model simulations of spike pattern detection by cortical networks (Haeusler and Maass, 2007) suggests that laminar models with plastic synapses allow for more accurate estimates of the detection capability compared to neural networks that do not take into account the layered structure of cortex.

## **DETECTING SIGNATURES OF ANTI-CORRELATED DEGREE DISTRIBUTIONS**

The model makes the prediction that anti-correlated networks would be more appropriate for the detection of nanostimulation in stable networks. To test this prediction we need to be able to distinguish correlations in the degree structure of the network without having access to all the inputs and all the outputs of a subset of neurons. We find that anti-correlations change the frequency of specific network motifs in a way that is independent of the network size, which means that it can be determined by averaging across many smaller sub-networks. A "ring" motif, number 98, which was a projection from neuron 1 to 2, from 2 to 3 and from 3 to 1, discriminated best between correlated and anticorrelated networks (**Figure 7**). Pairs of motif counts increased discriminability to a small extent, and only when the counts were pooled. This shows that these networks can be detected experimentally based on sampling sub-networks comprised of 30 neurons, when enough samples are available.

## **FUTURE STUDIES SHOULD INCORPORATE MULTIPLE TYPES OF INTERNEURONS**

The model was highly simplified so that we could focus on the connectivity structure. Having established the advantages of anticorrelation, our goal is to study the effects in more realistic networks. There are many other biophysical features that could be included in the model that would change the results quantitatively or, in some cases, even qualitatively. Here we highlight a small selection of the most relevant ones.

The first issue is inhibition. Experimental evidence shows that two types of inhibitory neurons, those expressing parvalbumin (PV) and somatostatin (SOM), are relevant in determining the gain of the response of pyramidal cells to whisker stimulation, visual stimulation or current injection (Gentet et al., 2010; Kwan and Dan, 2012; Lee et al., 2012). Avermann and coworkers (Avermann et al., 2012) constructed a model of L2/3 in barrel cortex constrained by *in vitro* measurements and studied the effect of stimulating varying amounts of pyramidal cells expressing channelrhodopsin by light pulses. In this model the strongest projection, in terms of the connection probability and synaptic strength, was from pyramidal cells to fast spiking (FS) interneurons (corresponding to PV neurons). Even when a relatively small fraction of the pyramidal cells were stimulated, almost all FS cells were recruited. For higher fractions of stimulated pyramidal cells, the non-fast spiking (NFS) interneurons (such as SOM interneurons) would become gradually activated. As a result the pyramidal cell activity remained low despite strong stimulation. The authors hypothesize that the strong inhibition is a mechanism to maintain sparse spiking in the pyramidal cells, with the NFS cells providing a back-up inhibitory mechanism. It is not clear how this computational model would be applicable to *in vivo* dynamics where FS cells are already spontaneously active. Furthermore, the level of activity in the different interneurons depends on brain state (Gentet et al., 2010). We have simulated binary networks with inhibitory neurons and find that anti-correlated degree distributions in the E–E sub-network improve stability (and yield the same sensitivity).

Detection could also take place by a state change in the network. The network has a LFS and a not too biologically plausible HFS, which in the context of a network with inhibition would perhaps correspond to something like an upstate. The true positive rate would correspond to how often single-cell stimulation would drive the network out of the BOA for the LFS, whereas the false positive rate would correspond to how often this would happen in the spontaneous state. The latter is given by the fraction of trials the system goes to the HFS state (**Figure 3**). The former can be tuned by changing the number of neurons and the duration of stimulation. A proper examination of this issue would require a network with a population of inhibitory neurons (Avermann et al., 2012).

A second issue is the effect of including spike timing. Synapses are sensitive through short-term depression and facilitation to the temporal patterns of stimulation (Abbott and Regehr, 2004), which could thereby affect the postsynaptic response in a nonlinear fashion, thereby preferentially activating specific populations of neurons. Dendritic nonlinearities also affect the impact of synaptic inputs based on their temporal coincidence and whether they arrive on the same part of the dendrite (Gasparini et al., 2004; Major et al., 2008; Polsky et al., 2009; Lavzin et al., 2012). Either of these effects could increase the sensitivity to external stimulation, while not appreciably changing the stability, thereby strengthening the results reported in this paper. However, to fully quantify these effects would require new and more extensive simulations that fall outside the scope of this paper.

## **PERCEPTUAL RELEVANCE OF ELECTRICAL OR OPTICAL STIMULATION IN EXPERIMENT**

Our study explores a hypothesis for how to achieve detection of an electrical stimulation by quick recurrent excitation that escapes before being shut down by inhibition, without destabilizing the spontaneous state. We now review the relevant literature focusing on the difference between electrical and sensory stimulation and the role of inhibition.

The barrel cortex normally processes thalamic activity generated in response to whisker stimulation. According to the canonical cortical circuit (Douglas and Martin, 2004; Lefort et al., 2009; Petersen and Crochet, 2013) this activity arrives first in layer 4 (L4) of the barrel column representing the stimulated whisker and then goes to L2/3 and subsequently to L5. It stands to reason that when during a task an animal needs to make a decision based on whisker stimulation, this is based on activity in L2/3 or L5 that came there by way of L4. The path taken by activity induced by optical, micro- or nanostimulation does not necessarily directly involve L4 and improving detection could thus require altering the underlying cortical circuit.

When monkeys were trained to detect microstimulation at a location in the visual cortex corresponding to a specific retinotopic location, the stimulation threshold for detection was reduced from about 50μA to 5μA over a few thousands of trials (Ni and Maunsell, 2010). At the same time the contrast threshold needed to detect real visual stimuli at the same retinotopic location increased from 4–8% to 8–60%. When the monkeys were subsequently retrained on detecting visual stimuli, the sensitivity was recovered in another few thousand trials, but the sensitivity to electrical stimulation was reduced. One possible interpretation is that learning to detect electrical stimulation reorganizes the recurrent circuits in L2/3 to become more sensitive at the expense of the L4 to L2/3 feedforward connection.

The animal improves its performance when learning to detect microstimulation, which could also be the case for single-cell nanostimulation modeled here. This improvement could occur because of one or more of the following reasons. First, the network could become more anti-correlated by changing the in-degrees. This means that the stability of the network would improve over time and perhaps that the number of false positives would reduce. Second, the out-degree of the stimulated neurons could increase, so that the nanostimulation signal becomes louder, hence the true positives should increase. Third, the neurons involved in the detection process become more sensitive to neurons directly downstream of the stimulated cells.

The threshold for detecting microstimulation in monkey visual cortex matches the strength necessary to elicit action potentials in mouse and cat cortex in the neighborhood of the electrode, 5–10μA (Histed et al., 2009) and in rat barrel cortex 2–5μA (Houweling and Brecht, 2008). These numbers did not depend on whether metal or glass pipette electrodes were used. Stimulation close above this threshold activated a set of widely dispersed neurons within a few hundred microns from the electrode, through antidromic action potentials in axons that are close to the electrode. As a result, the spatial pattern of activation was very sensitive to small changes in the location of the electrode.

Similar stimulus strength, 10μA for 0.1 to 0.5 ms, yielding charge transfers on the order of 1 nC, applied in the infragranular layers could be detected in rats (Butovas and Schwarz, 2007). The authors (Butovas and Schwarz, 2003) estimate that this corresponds to activating 80% of the pyramidal cells within 450 micron of the electrode, yielding an increase in their firing rate of 25% corresponding to about 0.5 excess spike per neuron. Interestingly, trains of electrical stimulation were more effective, indicating that temporal correlation may be necessary to distinguish stimulation from spontaneous activity. Physiological measurements indicated that synapses of pyramidal cells on fast spiking interneurons depress more than the pyramidal to pyramidal synapses, which means that pulse trains could lead to more a prominent increase in activity than single stimuli (Holmgren et al., 2003).

Optogenetics was used to determine how many neurons in L2/3 would be required to generate a change in activity that would be detectable by a mouse (Huber et al., 2008). The authors' estimate of 300 neurons producing one action potential was based on a measured distribution of light intensity thresholds necessary to elicit an action potential, the number of neurons expressing the light-sensitive channelrhodopsin (ChR2) channels and the spatial fall off of the light intensity, and represents according to these authors an overestimate. The number of 300 neurons corresponds to about 5% of the approximately 6500 neurons present in a mouse barrel column (Lefort et al., 2009).

Nanostimulation refers to electrical activation of an individual neuron with a glass pipette in the juxtacellular configuration. Nanostimulation in rat barrel cortex must have led to behaviorally relevant changes in network activity, as the animal was able to detect nanostimulation, but the average effect size was rather small (Houweling and Brecht, 2008). The nature of this activity could not be assessed, but experiments in mouse visual cortex may shed some light on this. Single-cell stimulation led to spikes in the stimulated neuron and calcium transients in some of the surrounding neurons that could be detected using two-photon microscopy (Kwan and Dan, 2012). Such stimulation induced postsynaptic activity in very few other pyramidal cells, 20 out of 1152 measured. SOM interneurons [corresponding to the NFS of Avermann et al. (2012)] were most strongly activated, 5 out of 17 measured. PV expressing cells did not respond to this stimulation, but their calcium transients were most strongly correlated to the network activity produced by the rest of the measured cells. This indicates that in this state the SOM cells would be required to damp the increase in activity generated by the recurrently connected pyramidal cell network.

#### **SUMMARY**

Taken together, experimental results suggest that detection of single-cell stimulation requires a quick propagation of excitatory cell activity, before the various types of inhibition kick in. Our studies indicate that anti-correlated degree distributions could be an important strategy for increasing sensitivity while maintaining stability.

#### **AUTHOR CONTRIBUTIONS**

Paul Tiesinga and Arthur R. Houweling designed the research project, Paul Tiesinga and Juan C. Vasquez wrote the code, Juan C. Vasquez and Paul Tiesinga performed the simulations, Paul Tiesinga wrote the manuscript together with Arthur R. Houweling.

#### **ACKNOWLEDGMENTS**

This work was supported by the Netherlands Organization for Scientific Research (NWO), through a grant entitled "Reverse physiology of the cortical microcircuit," grant number 635.100.023.

#### **REFERENCES**

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 10 July 2013; accepted: 17 October 2013; published online: 07 November 2013.*

*Citation: Vasquez JC, Houweling AR and Tiesinga P (2013) Simultaneous stability and sensitivity in model cortical networks is achieved through anti-correlations between the in- and out-degree of connectivity. Front. Comput. Neurosci. 7:156. doi: 10.3389/fncom.2013.00156*

*This article was submitted to the journal Frontiers in Computational Neuroscience. Copyright © 2013 Vasquez, Houweling and Tiesinga. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Statistical evaluation of synchronous spike patterns extracted by frequent item set mining

#### *Emiliano Torre1 \*, David Picado-Muiño2, Michael Denker 1, Christian Borgelt 2 and Sonja Grün1,3*

*1 Institute of Neuroscience and Medicine (INM-6) and Institute for Advanced Simulation (IAS-6), Jülich Research Centre and JARA, Jülich, Germany*

*2 European Centre for Soft Computing, Mieres, Spain*

*3 Theoretical Systems Neurobiology, RWTH Aachen University, Aachen, Germany*

#### *Edited by:*

*Ruben Moreno-Bote, Foundation Sant Joan de Deu, Spain*

#### *Reviewed by:*

*Shigeru Shinomoto, Kyoto University, Japan Srdjan Ostojic, Ecole Normale Superieure, France*

#### *\*Correspondence:*

*Emiliano Torre, Institute of Neuroscience and Medicine (INM-6) and Institute for Advanced Simulation (IAS-6), Jülich Research Centre, Wilhelm-Johnen-Strasse, 52425 Jülich, Germany e-mail<EMAIL_ADDRESS>
We recently proposed frequent itemset mining (FIM) as a method to perform an optimized search for patterns of synchronous spikes (*item sets)* in massively parallel spike trains. This search outputs the occurrence count (*support*) of individual patterns that are not trivially explained by the counts of any superset (*closed frequent item sets*). The number of patterns found by FIM makes direct statistical tests infeasible due to severe multiple testing. To overcome this issue, we proposed to test the significance not of individual patterns, but instead of their signatures, defined as the pairs of pattern size *z* and support *c*. Here, we derive in detail a statistical test for the significance of the signatures under the null hypothesis of full independence (*pattern spectrum filtering*, PSF) by means of surrogate data. As a result, injected spike patterns that mimic assembly activity are well detected, yielding a low false negative rate. However, this approach is prone to additionally classify patterns resulting from chance overlap of real assembly activity and background spiking as significant. These patterns represent false positives with respect to the null hypothesis of having one assembly of given signature embedded in otherwise independent spiking activity. We propose the additional method of *pattern set reduction* (PSR) to remove these false positives by conditional filtering. By employing stochastic simulations of parallel spike trains with correlated activity in form of injected spike synchrony in subsets of the neurons, we demonstrate for a range of parameter settings that the analysis scheme composed of FIM, PSF and PSR allows to reliably detect active assemblies in massively parallel spike trains.

**Keywords: higher-order correlations, neuronal cell assemblies, spike patterns, spike synchrony, multiple testing, data mining**

## **1. INTRODUCTION**

The cortex is comprised of a highly interconnected network of neurons and thus one may speculate that information processing in the brain may only be understood on the basis of the concerted activity of the neuronal population. Hebb (1949) suggested that neurons coordinate their activities by organizing in functional groups, termed cell assemblies. Synchronous spike input to receiving neurons is known to be more effective in generating output spikes (Abeles, 1982; König et al., 1996), which leads to the hypothesis that temporal coordination of spiking activity or correlational processing is the defining expression of an active cell assembly (Singer et al., 1997; Harris, 2005). As excitatory postsynaptic potentials are small in amplitude compared to the gap between the resting potential and the neuronal firing threshold, it is expected that a cell assembly is composed of many neurons firing in a correlated fashion. This observation is the basis for the assumption that higher-order synchronous spiking activity serves as a signature expression of an active assembly (Riehle et al., 1997; Berger et al., 2010; Staude et al., 2010b; Shimazaki et al., 2012).

In order to observe and detect such signatures in the brain, the spiking activities of many neurons must be recorded simultaneously. Fortunately, in recent years considerable progress has been made in the development of multi-electrode recording techniques [e.g., Nicolelis, 1998; Buzsaki, 2004; Hatsopoulos et al., 2007; Riehle et al., 2013], which enable to record the activity of hundred(s) of neurons. Such massively parallel spike train data pose statistical challenges due to the inherent complexity of the required multivariate approaches. Most notably, increasing the number of observed neurons leads to a combinatorial explosion of the number of potential spike patterns that need to be detected and tested. Based on pairwise correlation analyses only, the existence and functional relevance of neuronal correlations could be demonstrated in various cortical systems and behavioral paradigms [e.g., Gerstein and Aertsen, 1985; Riehle et al., 1997; Kohn and Smith, 2005; Berger et al., 2007; Fujisawa et al., 2008; Feldt et al., 2009; Humphries, 2011; Masud and Borisyuk, 2011]. Nevertheless, a correlation analysis considering the complete set of simultaneously recorded spike trains is required to uncover also higher-order correlations among neurons. In recent years several such approaches were developed, each of which focuses on different aspects: (i) methods to determine the presence of higher-order spike correlations with a minimum order without explicitly identifying the participating neurons [e.g., Louis et al., 2010a; Staude et al., 2010a,b]; (ii) methods that test whether individual neurons participate in synchronous spiking activity without identifying the groups of correlated neurons [e.g., Berger et al., 2010]; (iii) methods that test for the presence of correlation as predicted by a specific correlation model such as a synfire chain (Abeles, 1991), that is, spatio-temporal spike patterns or propagation of synchronous spiking activity [e.g., Abeles and Gerstein, 1988; Schrader et al., 2008; Gerstein et al., 2012; Gansel and Singer, 2012]; (iv) methods that directly identify the members of cell assemblies on the basis of the patterns of synchronous spiking activity [e.g., Gerstein et al., 1978; Pipa et al., 2008; Feldt et al., 2009; Gansel and Singer, 2012; Shimazaki et al., 2012; Picado-Muiño et al., 2013].

In Picado-Muiño et al. (2013) we presented the basic approach and relevant statistics to employ frequent item set mining (FIM) to identify significant patterns of spike synchrony in massively parallel spike trains. FIM enables fast and efficient counting of synchronous spike patterns by pruning the tree of all possible patterns. To address the problem of multiple testing, statistics are not computed for individual patterns, but on the pattern spectrum that collects the number of observed patterns based on their signature. A signature is defined as the pair (*z*,*c*) of pattern size *z* (i.e., number of participating neurons) and *support c* (i.e., number of pattern occurrences). In *pattern spectrum filtering* (PSF) those identified sets of neurons for which patterns with the same signature (*z*,*c*) occur also in appropriate surrogate data are then marked as chance patterns and discarded.

Here, we extend the approach of Picado-Muiño et al. (2013) in three ways that will enable the application of the method to biological data. First, we refine the statistical test employed in pattern spectrum filtering for reporting significant patterns of a given signature (Section 2). Then, we introduce a subsequent analysis step, termed *pattern set reduction* (PSR), to additionally filter out those patterns that are detected as significant, but are compositions of chance spikes or patterns and the actual cell assembly pattern (Section 3). Finally, we report on the performance of our method related to features describing the data (e.g., coincidence rate, assembly pattern size, firing rate heterogeneity or non-stationarity) and analysis parameters (Section 4). The discussion (Section 5) includes a step-by-step instruction on how to utilize the proposed method in the context of massively parallel spike trains obtained from electrophysiological recordings.

## **2. SPIKE PATTERN DETECTION AND STATISTICAL TESTING**

In this section we introduce our approach to detect frequent synchronous spike patterns in massively parallel spike trains (MPST). We first briefly review frequent item set mining (FIM) and related terminology and definitions as proposed in Picado-Muiño et al. (2013) as a tool to efficiently detect and count synchronous spike patterns in MPST. Then we derive a modified version of the FIM-based statistics proposed in Picado-Muiño et al. (2013) for assessing pattern significance.

### **2.1. FREQUENT ITEMSET MINING**

Given *N* parallel spike trains with neuron ids 1, 2,..., *N*, observed in the time window [0, *T*), we partition [0, *T*) into *b* exclusive bins {*bi*} *b i* = 1 of identical width *w* = *T*/*b* (typically chosen as a few ms): *bi* = [(*i* − 1) · *w*, *i* · *w*). If one or more spikes of one neuron fall into a bin, we consider the bin occupied and reduce the entry to 1 (*clipping*), so that each time bin contains at most one spike per neuron. Spikes from different neurons falling into the same time bin are defined as *synchronous* (see **Figure 1A**). Borrowing terminology from FIM, we define each neuron id as an *item*, the set *Ti* of all items spiking in *bi* as the *i*-th *transaction* in the binned data, and {*Ti*} *b i* = 1 as the *transaction list*. Given a *minimum pattern size z*0, each set of *z* ≥ *z*0 items in *Ti* constitutes a *pattern of synchronous spikes*, or *item set* (see **Figure 1B**). Here we

**FIGURE 1 | From spike data to closed frequent itemsets. (A)** Sketch of a raster plot of 4 neurons firing in parallel. Shaded colors separate adjacent bins. Red spikes mark the occurrences of the synchronous pattern composed of neurons 1, 3, 4. **(B)** Transaction list derived from the spike data in **(A)** after binning. **(C)** List of item sets obtained from **(B)**, together with their occurrence counts. Black boxes mark non-frequent item sets (support set to 2), blue boxes mark non-closed frequent item sets, red boxes mark CFISs.

**(D)** Average number of item sets (dashed black line), frequent item sets (dashed blue line) and CFISs (dashed red line) obtained from 100 simulations of 100 parallel independent spike trains with a firing rate of 20 Hz, as a function of the simulation time. Other parameters are bin width *w* = 3 ms and minimum pattern size *z*0 = 2. Bars mark ±1 std. dev. The solid line indicates the number of time bins (and thus transactions) as a function of the simulation time.

set *z*0 to 2. Due to clipping, each item set occurs at most once per transaction. The number of occurrences of an item set in the transaction list is the *support* of that item set.

A transaction that contains *K* items yields 2*K* − *K* − 1 different (but possibly overlapping) item sets of size *z* ≥ 2, that is, all 2*K* possible subsets without the empty set and the *K* singletons. The total number of different item sets in a transaction list can thus largely exceed the number of transactions (i.e., time bins). This number grows with the duration of the data set (see **Figure 1D**) and with the number of parallel spike trains (not shown).

In order to limit the data to potentially interesting and nontrivial item sets, we select only item sets whose support *c* is larger than or equal to a *minimum support c*0 (*c*0 ≥ 1) as introduced by Picado-Muiño et al. (2013). Here we set *c*0 to 2. An item set whose support equals or exceeds the minimum support is called *frequent item set*. For *c*0 > 1, frequent item sets are usually a small fraction of all item sets (**Figure 1D**, compare black dashed line to blue dashed line). Furthermore, we discard any frequent item set occurring as many times as any of its supersets. These patterns are trivially explained by the occurrences of their supersets, which are more significant due to the larger number of neurons involved. Non-trivial frequent item sets are called *closed frequent item sets* (CFISs; see **Figure 1C**). Discarding non-closed frequent item sets does not yield any loss of information. Indeed, the set *F* of all frequent item sets can be reconstructed from the set *C* of CFISs by

$$\mathcal{F} = \bigcup\_{I \in \mathcal{C}} \bigcup\_{J \subset I, \ |J| \ge z\_0} J.$$

The support *s*(*I*) of a non-closed frequent item set *I* ∈ *F* can be computed as *s*(*I*) = max*J*∈*C*, *J*⊃*Is*(*J*).

If *A* and *B* are two CFISs such that *B* - *A*, and *cA*, *cB* their respective supports, it follows from the definition of CFISs that *cB* > *cA* (*a priori* property). We refer to the (non-empty) set *A*\*B* as the *excess items* of *A* with respect to *B*, and to the difference *cB* − *cA* as the *excess occurrences* of *B* with respect to *A*.

Following Picado-Muiño et al. (2013), we make use of frequent itemset mining [FIM; for a review, see Goethals (2010), Borgelt (2012)] to extract CFISs and their support from an MPST transaction list. FIM performs a non-redundant search for spike patterns, starting from those of size *z*0 and then moving on to supersets of increasing size. Starting at lowest-size patterns, the search is organized in a search tree in layers of increasing pattern size. A branch connects two patterns if one is a subset of the other. Each pattern is visited at most once. FIM exploits the *apriori* property to stop the search at infrequent patterns, as no supersets of an infrequent item set can be frequent. The output of FIM is a list of all CFISs with their support (**Figure 1C**).

#### **2.2. PATTERN SPECTRUM FILTERING**

Direct statistical tests of all individual patterns occurring in MPST are not suitable, as they cause a severe multiple testing problem yielding large occurrences of false positives (FPs), or enhanced levels of false negatives (FNs) after statistical corrections. Therefore Picado-Muiño et al. (2013) proposed to pool CFISs according to their size *z* (number of neurons involved) and their support *c* (number of occurrences) in a two-dimensional histogram (*pattern spectrum*) and to evaluate patterns of the same signature (*z*,*c*) for significance by a Monte-Carlo approach using surrogate data. Here we present a refinement of this original approach, named *pattern spectrum filtering* (PSF), that bases the test for a specific signature (*z*,*c*) also on patterns of higher size and support than specified by the signature.

In order to implement the null hypothesis *H*0 of independent spiking, and to approximate the *p*-values of the signatures (*z*,*c*), from the original data (**Figure 2A**) we repeatedly generate surrogate data (**Figure 2B**), collect from each one its CFISs through FIM as done for the original data, and compute the corresponding surrogate pattern spectrum (**Figure 2C**). The surrogates are generated from the original data by intentionally destroying correlations while keeping other features, such as firing rates, intact [e.g., by spike randomization or spike dithering, Louis et al. (2010b)].

Let be the partial ordering on the real plane, that is, (*x*∗, *y*∗) (*x*, *y*) if *x*∗ ≥ *x* and *y*∗ ≥ *y*, where holds if at least one inequality is strict. From each surrogate pattern spectrum we compute a binary spectrum which takes value 1 at each signature (*z*,*c*) such that at least one signature (*z*∗,*c*∗) (*z*,*c*) is occupied, and value 0 otherwise [in contrast to Picado-Muiño et al. (2013) where only the occupation of signature (*z*,*c*) is checked]. Formally, we define the *signature operator* sgt(·) such that, given a CFIS *A* with size *zA* = |*A*| and occurrence count *cA*, sgt(*A*) := (*zA*,*cA*). For each list *Si* of CFISs from one surrogate data set, let *P*ˆ*i* be the *binary pattern spectrum*, defined for each *z*,*c* ≥ 2 by:

$$\hat{P}\_i(z,c) := \begin{cases} 1 \text{ if } \exists A \in \mathcal{S}\_i \text{ : } \text{sgt}(A) \succeq (z,c) \\ 0 \text{ otherwise} \end{cases}.$$

Averaging the binary spectra at each signature, we get the *p-value spectrum P*ˆ:

$$\hat{P}(z,\mathcal{c}) := \frac{1}{K} \# \left( \mathcal{S}\_i : \exists A \in \mathcal{S}\_i : \text{sgt}(A) \succeq (z,\mathcal{c}) \right).$$

*P*ˆ(*z*,*c*) yields an estimate of the probability to observe (one or more) patterns with signature (*z*∗,*c*∗) (*z*,*c*) under *H*0 (see **Figure 2D**).

We then classify any signature (*z*,*c*) whose *p*-value is lower than the significance level α∗ as significant. Given the desired overall significance level α for PSF, we derive α∗ from α by Bonferroni correction for the number *m* of tests, i.e., the number of signatures in the data to test for: α∗ = α/m. Any signature (*z*,*c*) for which *P*ˆ(*z*,*c*) < α∗ is classified as significant. Formally, we introduce the *significance spectrum S*ˆ defined at each (*z*,*c*) by

$$
\hat{S}(z,c) := \begin{cases}
1 & \text{if } (z,c) \text{ is significant} \\
0 & \text{otherwise}
\end{cases}
$$

In **Figure 2E** *S*ˆ(*z*,*c*) = 1 is marked in white, *S*ˆ(*z*,*c*) = 0 in gray. The border between the two is the *detection border*, on the left of which signatures in the original data are classified as not significant and rejected. Signatures to its right (*S*ˆ(*z*,*c*) = 1) are considered as significant (marked in red in **Figure 2E**). The corresponding patterns and their supports are listed in **Figure 2F**.

**FIGURE 2 | PSF on artificial data. (A)** Raster plot of 100 parallel simulated spike trains consisting of independent Poisson activity plus 6 injections of one pattern of synchronous spikes (highlighted in red) from neurons 1 to 10, occurring at random times (see Section 4 for details). The total firing rate of each neuron is 20 Hz, the simulation time is 3 s. **(B)** Same as in **(A)**, but without injection of synchronous patterns. The spike trains are therefore completely independent. **(C)** Pattern spectrum of CFISs extracted from the data in **(A)** by FIM (*z*0 = 2, *c*0 = 2, *w* = 5 ms). Counts are color-coded (logarithmic scale). **(D)** *P*-value spectrum drawn from 5000 surrogate, independent data sets of the type shown in **(B)**. *P*-values are color-coded (logarithmic scale). **(E)** Significance spectrum (overall significance α = 0.01, Bonferroni-corrected for *m* = 50 tests yielding α∗ = 2 · 10−4). Gray squares indicate signatures that are not significant, white squares mark potentially significant signatures. Red squares mark significant signatures of the pattern spectrum shown in **(C)**, i.e., which fall into white squares of the significance spectrum. **(F)** List of patterns detected by PSF. Besides the injected pattern *A* = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, PSF also classifies additional patterns as significant, all being subsets or supersets of *A*.

## **3. PATTERN SET REDUCTION**

PSF tests the significance of patterns under the null hypothesis *H*0 of fully uncorrelated spike trains. However, PSF might fail in rejecting patterns that result from combinations of chance spikes or chance patterns with the assembly pattern (see list of detected patterns in **Figure 2F** besides the injected one). These patterns are a specific kind of false positive, not resulting from merely independent data. They may be subsets or supersets of the assembly pattern, or patterns that partially overlap with it (**Figures 3A–C**). In this section we define the type of FPs that may occur, investigate why PSF is prone to return such FPs, and propose an additional statistical analysis, termed *pattern set reduction* (PSR), to remove them.

## **3.1. TYPES OF FPs**

#### *3.1.1. Chance subsets*

If a CFIS *A* repeats *cA* times and a subset *B* of *A* (with |*B*| ≥ *z*0) has *c* additional chance occurrences, *B* represents a CFIS repeating *cB* = *cA* + *c* total times. We call *B* a *chance subset* of *A*, having *c* excess occurrences (**Figure 3A**). PSF is designed to test the significance of signature (|*B*|, *cB*) under *H*0 (complete independence), thus disregarding the fact that *cA* occurrences are due to pattern *A*. As a result it classifies *B* as a significant pattern, thus yielding an FP outcome. This is illustrated in **Figure 2F**, where e.g., pattern {4, 6, 10} occurs twice by chance plus 6 times as a subset of pattern {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}. The corresponding signature (3, 8) is significant compared to the surrogates (**Figure 2E**), so that PSF does not reject it.

## *3.1.2. Chance supersets*

If a CFIS *B* occurs *cB* times and another set *C* of neurons fires by chance synchronously with *B* in *c* of those *cB* transactions (with *c* ≥ *c*0), then the pattern *A* = *B* ∪ *C* represents a CFIS repeating *cA* = *c* times. We call *A* a *chance superset* of *B*, with |*C*| excess neurons (**Figure 3B**). PSF tests the significance of signature (|*A*|, *cA*) under *H*0, disregarding the fact that |*B*| of the |*A*| neurons of *A* are due to the presence of pattern *B*. The test is therefore prone to classify *A* as significant. This is the case for patterns {1, 2,..., 10, 80}, {1, 2,..., 10, 28} and {1, 2,..., 10, 24} in **Figure 2F**, each of which occurs twice as a superset of {1, 2,..., 10}. The corresponding signature (11, 2) is significant compared to the surrogates (**Figure 2E**), so that PSF classifies these patterns as significant.

**FIGURE 3 | Excess occurrences and excess items.** Sketch of the possible relationship between a reference pattern and patterns sharing neuron identities and/or time occurrences with it. In each panel, ticks represent individual spikes. Rows correspond to neurons and columns to transactions, i.e., time bins. Spikes forming a pattern are grouped by an ellipse. The reference pattern of each panel is shown by black ticks and is indicated by a

solid ellipse. **(A)** *B* is a subset of *A* with excess occurrences (red). **(B)** *A* is a superset of *B* with excess items (blue). **(C)** *B* is a subset of *A* with excess occurrences (red). Neurons in *C* (blue) additionally fire synchronously to *A* and to excess occurrences of *B*. Thus pattern *D* = *B* ∪ *C* forms a CFIS, which partially overlaps with *A*. **(D)** Patterns *A* and *B* are disjoint: they are composed of different neuron identities and occur at different time bins.

## *3.1.3. Chance overlapping sets*

The simultaneous presence of excess items and excess occurrences can yield yet another type of FP outcome, namely patterns that overlap with the actual assembly. Given an assembly *A*, assume that a subset *B* of *A* has additional chance occurrences. If an additional set *C* of neurons disjoint from *A* fires synchronously to *A and* to an excess occurrence of *B* for a total of *c* ≥ *c*0 chance times, then the set *D* = *B* ∪ *C* represents a CFIS which partially overlaps with *A* (**Figure 3C**). PSF is prone to classify *D* as significant.

## *3.1.4. Disjoint patterns*

Two patterns which do not have items in common are *disjoint* (**Figure 3D**). In contrast to the previous classes of chance patterns, the presence of an active assembly does not enhance chance patterns disjoint from it. PSF therefore correctly estimates their significance and manages to filter out almost all of them, as shown in 4.

## **3.2. PSR STATISTICS**

Let *P* be the class of CFISs reported as significant by PSF. Given a pair (*A*, *B*) ∈ *P* × *P* such that *B* ⊂ *A* (therefore *cB* > *cA* by definition of CFIS, and |*B*| < |*A*|), we propose statistical tests to assess the conditional significance of either *A given B* (*A*|*B*) or *B given A* (*B*|*A*), i.e., of one pattern given that the other represents an assembly pattern. These tests can be applied, using different strategies, to the class of all such (*A*, *B*) pairs, reducing *P* to a subclass *Q* of patterns which are mutually significant given each other.

## *3.2.1. Subset filtering*

This procedure aims at rejecting FPs that are chance subsets of other CFISs. For each pair (*A*, *B*) ∈ *P* × *P* such that *B* ⊂ *A* (so that *cB* > *cA*), *B* has *cB* − *cA* excess occurrences with respect to *A*. Subset filtering tests *B*|*A*, i.e., the null hypothesis *HB*|*A* 0 that *B* is a chance subset of the actual assembly *A*, by assessing the significance of the excess occurrences of *B*. Equivalently, *HB*|*A* 0 states that the pattern *B* defined by the same items as *B* but its excess occurrences only (red spikes in **Figure 3A**) is a chance pattern. If *HB*|*A* 0 is rejected, *B* is kept and *A* discarded, otherwise *A* is kept and *B* discarded. Thus, the procedure keeps either *A* or *B* and discards the other (*exclusive*). We present two alternatives to test *HB*|*A* 0 .

*<IP_ADDRESS>. Exact test.* This test computes the *p*-value of the signature (|*B*|,*cB* − *cA*) of *B* . If *cB* − *cA* < *c*0, *B* is classified as a chance subset of *A*. Otherwise, let *T A* be the transaction list obtained from *T* by discarding the transactions where *A* occurred, and keeping in the remaining transactions only the items composing *A*. All the excess occurrences of subsets of *A* must be contained in *T A*. *B* itself is a CFIS in this transaction list: it is an item set because |*B* |=|*B*| ≥ *z*0, it is frequent because *cB* − *cA* ≥ *c*0, it is closed because otherwise *B* itself would be non-closed. To test the significance of *B* , one can therefore run FIM and PSF on surrogates of *T A* to estimate the significance of its signature (|*B*|,*cB* − *cA*)**.** If (|*B*|,*cB* − *cA*) is significant, *B* is significant in *T A* and *B* is classified as significant in *T* (given *A*). Otherwise, *B* is classified as non-significant.

*<IP_ADDRESS>. Approximate test.* This test approximates the *p*-value of the signature (|*B*|,*cB* − *cA*) in *T A* by the *p*-value of the signature (|*B*|,*cB* − *cA* + *h*), *h* ≥ 1, in *T*, already obtained when performing PSF. In contrast to *T A*, *T* is composed of more neurons than those which can actually form chance subsets of *A* (because it does not contain the items of *A* only), and more transactions than those where such subsets could actually display excess occurrences (because it also contains the transaction where *A* is already present). Therefore, the *p*-value of (|*B*|,*cB* − *cA*) would be underestimated if computed over *T* instead of *T A*. Parameter *h* heuristically corrects for this by substituting it with the *p*value of a signature with the same size but higher support. The lower *h*, the higher the probability to reject *B*. If *h* ≥ *cA*, then (|*B*|,*cB* − *cA* + *h*) (|*B*|, *cB*) and *B* is necessarily reported as significant. This test avoids to run FIM and PSF on *T A* and is therefore computationally more efficient.

## *3.2.2. Superset filtering*

This procedure aims at rejecting FPs that are chance supersets of other CFISs. For each pair (*A*, *B*) ∈ *P* × *P* such that *B* ⊂ *A* (so that |*B*| < |*A*|), *A* has |*A*|−|*B*| excess items with respect to *B*. Subset filtering tests *A*|*B*, i.e., the null hypothesis *HA*|*B* 0 that *A* is a chance superset of the actual assembly *B*, by assessing the significance of the excess items of *A*. Equivalently, *HA*|*B* 0 states that the pattern *A* defined by the same transactions as *A* but containing its excess items only (blue spikes in **Figure 3B**), is a chance pattern. If *HA*|*B* 0 is rejected, *A* is kept and *B* discarded from *P*, otherwise *B* is kept and *A* discarded from *P*. Thus, the procedure keeps either *A* or *B* and discards the other (*exclusive*). We present two alternatives to test *HA*|*B* 0 .

*<IP_ADDRESS>. Exact test.* This test computes the significance of the signature (|*A*|−|*B*|, *cA*) of *A* . If |*A*|−|*B*| < *z*0, *A* is classified as a chance superset of *B*. Otherwise, let *T B*¯ be the transaction list obtained from *T* by keeping only the transaction where *B* occurred, and discarding from them the items constituting *B*. All groups of excess items of *B* (i.e., neurons that fire synchronously to *B*) must be contained in *T B*¯ . *A* itself is a CFIS of this transaction list: it is an item set because |*A* |=|*A*|−|*B*| ≥ *z*0, it is frequent because *cA* ≥ *c*0, it is closed because otherwise *A* itself would be non-closed. To test the significance of *A* , one can therefore run FIM and PSF on surrogates of *T B*¯ to estimate the *p*-value of its signature (|*A*|−|*B*|, *cA*)**.** If (|*A*|−|*B*|, *cA*) is significant, *A* is significant in *T B*¯ and *A* is classified as significant in *T* (given *B*). Otherwise, *A* is classified as non-significant.

*<IP_ADDRESS>. Approximate test.* This test approximates the *p*-value of the signature of *A* in *T B*¯ by the *p*-value of signature (|*A*|−|*B*| + *k*, *cA*), *k* ≥ 1, in *T*, already obtained when performing PSF. In contrast to *T B*¯ , *T* is composed of more neurons than those that can actually form excess items of *B* (because it contains the items of *B*, too), and more transactions than those where supersets of *B* could actually occur (because it contains also transactions where *B* does not occur). Therefore, the *p*-value of (|*A*|−|*B*|,*cA*) would be underestimated if computed over *T* instead of *T B*¯ . Parameter *k* heuristically corrects for this by substituting it with the *p*-value of a signature with the same support but higher size. The lower *k*, the higher the probability to reject *A*. Note that if *k* ≥ |*B*| then (|*A*|−|*B*| + *k*, *cA*) (|*A*|, *cA*) and *A* is necessarily reported as significant. This test allows to avoid running FIM and PSF on *T B*¯ for each *B*.

## *3.2.3. Covered-spikes criterion*

This simple selection strategy consists of taking all pairs (*A*, *B*) ∈ *P* × *P* for which *B* ⊂ *A*, and keeping for each pair the pattern covering the largest number of spikes, while rejecting the other. Specifically, the criterion prefers *A* to *B* if *zA* · *cA* ≥ *zB* · *cB*, *B* to *A* otherwise. It does not involve significance tests, but is based on the observation that, given the probability *p* for a neuron to spike in a time bin, the probability for *z* neurons to fire synchronously in a bin is approximately *pz*, so that the probability that this pattern occurrs *c* times is binomially distributed and approximately proportional to *pz*·*c* . The larger the *z* · *c* score, the less likely a pattern of such size and support. This matches the finding that the detection border separating non-significant signatures (marked gray in **Figure 2E**) from significant ones (marked white in **Figure 2E**) in the significance spectrum exhibits a hyperbolic shape. The criterion thus keeps the less likely of the two patterns.

A variant consists in keeping the pattern with the largest (*z* − 1) · *c* score. This choice is motivated by the observation that a pattern of size *z* and support *c* can be seen as *z* − 1 spike trains which synchronize their spikes to another train *c* times. Thus, (*z* − 1) · *c* spikes are coincident to spikes in another spike train. Keeping the pattern with the largest (*z* − 1) · *c* score amounts to keeping the pattern which covers more coincident spikes. Geometrically, penalizing the pattern size corrects for the fact that the hyperbolic shape of the detection border in **Figure 2E** is elongated toward the pattern support (*y*-axis) rather than being equilateral.

## *3.2.4. Combined filtering*

Subset filtering, superset filtering and covered-spikes criterion can be combined into a filtering procedure which tests for both excess coincidences and excess items. Combined filtering tests for each pair (*A*, *B*) ∈ *P* × *P* both the null hypotheses *HB*|*A* 0 (i.e., that *B* is a chance subset of *A*) *and HA*|*B* 0 (i.e., that *A* is a chance superset of *B*). If one of the null hypotheses is rejected, the corresponding pattern is retained as significant. Thus, if both hypotheses are rejected, both patterns are retained (*inclusive*). Accepting one null hypothesis does not necessarily lead to the rejection of the corresponding pattern (in contrast to subset or superset filtering): the pattern is rejected only if the other pattern is accepted, i.e., if the other null hypothesis is rejected. If both *HB*|*A* 0 and *HA*|*B* 0 are accepted, one of the two patterns is kept based on the covered-spikes criterion.

## **4. CALIBRATION ON ARTIFICIAL DATA**

In this section we compare the performance (in terms of FPs and FNs) of PSF to PSF followed by PSR to illustrate the advantages yielded by the latter. For the sake of computational efficiency we employ the approximate versions of the tests for the subset and superset filtering with parameters *h* = 1 and *k* = 2, respectively. We test different types of artificial data that involve typical features of experimental data. After studying the general behavior of the analysis method for stationary, homogeneous data, we study data sets with heterogeneous firing rates across neurons, and with non-stationary firing rates in time.

## **4.1. CORRELATED DATA**

As a model for data containing assembly activity, we generate correlated spike trains by a modified version of the singleinteraction-process [SIP; Kuhn et al. (2003); Berger et al. (2010)], which we keep calling SIP for convenience. First, we simulate *N* = 100 parallel independent Poisson spike trains as background activity. Then we model assembly activity by inserting synchronous spike events in a subset of *z* of the *N* neurons (the *SIP neurons*, with ids 1 to *z*). This is done by generating a hidden Poisson process with the desired number *c* of pattern occurrences, from which spikes are copied into each of the *z* spike trains of the SIP neurons. Thus, as compared to the model proposed by Kuhn et al. (2003) we insert correlated firing only in a specific subset of the parallel processes. Before insertion of the synchronous patterns, the background firing rate of the SIP neurons is reduced by the rate of the hidden process to ensure the same firing rate for all neurons. In the simplest scenario, the firing rates and the pattern occurrence rate are stationary over time and homogeneous across neurons. More complicated cases will include either nonstationarity or heterogeneity of rates. The purpose of the analysis of such data is to test under controlled conditions if the simulated assembly is indeed detected and can be distinguished from background activity.

### **4.2. INDEPENDENT DATA**

To implement the null-hypothesis *H*0 of complete independence needed to derive the significance of signatures of the correlated data, we generate independent Poisson processes of the same rates as the data to be tested, thus keeping the same marginal statistics. This is one way of implementing the null-hypothesis. However, in the context of analyzing real experimental data, one may want to keep more statistical features of the experimental data (e.g., nonstationary and heterogenous firing rates, deviation from Poisson, and so on). This can be realized by the use of more complex surrogates derived by manipulation of the original data, e.g., spike dithering (Grün, 2009; Louis et al., 2010b).

#### **4.3. ASSESSING SIGNIFICANCE**

We evaluate the performance of our analysis in terms of the average number of FPs and FNs obtained with PSF and PSR in *R* = 1000 iterations on the same model of correlated data (SIP of size *z* in *N* = 100 parallel spike trains). To study the performance of our analysis, we investigate 243 models differing in the size of the injected assembly *z* = 2,..., 10, its injection count *c* = 2,..., 10, and the firing rates *r* = 5, 10 or 20 Hz (here: homogeneous for all neurons). We analyse each model with a bin width *w* = 3 ms and *w* = 5 ms for the detection of synchronous spike patterns. See **Table 1** for an overview of the parameter combinations. For the significance estimation we generate surrogate data, i.e., independent Poisson processes with the same firing rates as the correlated data, and analyse them with FIM as done for the correlated data. This procedure is repeated for *K* = 5000 times to derive the *p*-value spectrum and then the significance spectrum by employing an overall significance level of α = 0.01, Bonferroni-corrected for the number of signatures tested. The latter is given by the number of signatures existent in the correlated data, which never exceeded *m* = 50. In order to have the same corrected significance level for each of the 1000 iterations of each SIP model, we always correct for *m* = 50 tests, instead of correcting for the individual number *m* < *m* of signatures found in each data set. This yields the corrected significance level α∗ = 2 · 10−4, which is typically more conservative than correcting individually for *m* tests. This procedure allows us to use a single significance spectrum for all 81 SIP models with the same firing rates, differing by parameters *z* and *c* only, and for all 1000 realizations of each model. To obtain the *p*-values with precision α∗ we generate *K* = 1/α∗ = 5.000 surrogates, compute their binary spectra and average them to draw the *p*-value spectrum (see Section 2.2).

**Figure 4** shows significance spectra obtained from surrogate data for models differing by the firing rate *r* (5, 10 or 20 Hz) analysed with different bin widths *w* (dark gray for *w* = 3 ms, light gray for *w* = 5 ms; α∗ = 2 · 10−4). The set of non-significant signatures shows a hyperbolic shape, which grows with both *r* and *w* to higher *z* and higher *c*. Both factors, higher firing rates and

*Parameters for the background activity: N: number of neurons, r: firing rate, T : simulation time. Parameters for correlated data: z: number of neurons in correlated activity (size of SIP), c: SIP occurrences. Analysis parameters: w: bin width, c0: minimum item set support, z0: minimum item set size. Statistical parameters:* α*\*: Bonferroni-corrected significance level (for m* = *50 tests), K: number of surrogates, R: number of simulation runs per SIP model.*

**FIGURE 4 | Significance spectra for different parameter sets.** Independent Poisson spike trains (*N* = 100; *T* = 3 s) of different firing rates (*r* = 5, 10 or 20 Hz) serve as surrogates for the computation of three significance spectra (from left to right). Each square represents a (*z*, *c*) signature. Dark-shaded gray squares mark non-significant signatures obtained with *w* = 3 ms. Light-shaded squares represent further non-significant signatures for *w* = 5 ms. White squares indicate significant signatures for both choices of the bin width. Other parameters: *z*0 = 2, *c*0 = 2, α∗ = 2 · 10−4, *K* = 5000.

larger bin width, cause more spikes per bin, and therefore larger and more frequent chance patterns.

## **4.4. PERFORMANCE, HOMOGENEOUS FIRING RATES**

For each SIP parameter set we simulate the corresponding model *R* = 1000 times, and evaluate FPs and FNs of each realization. Their averages measure the performance of the analysis for each parameter constellation.

As previously discussed (Section 3), in the presence of correlations PSF tends to classify chance subsets, supersets or overlapping sets as significant, thus yielding FPs. **Figure 5**, top row, shows this effect on simulations of SIP models differing by SIP size (*x*-axis of each panel) and injection count (*y*-axis). For each model, the FP level is computed as an average over 1000 stochastic simulations. The total amount of FPs increases as the SIP size and/or the number of injections get larger. The contribution of FP supersets (green) and FP subsets (blue) is about the same, while in comparison FP overlapping sets (yellow) occur only at higher values for *z* and *c*, and FP disjoint patterns (purple) are almost never observed. As shown in **Figure 5**, bottom row, PSR (here, combined filtering) largely reduces the amount of FPs. Although the PSR statistical tests apply to chance subsets (blue) and supersets (green) only (Section 3.2), they successfully remove most of the overlapping patterns (yellow) as well. The reason is that, if there is a CFIS *D* overlapping with the actual assembly *A* by *z*0 or more items, their intersection *B* is a CFIS as well (**Figure 3C**). In most cases PSF classifies *B* as significant together with *A* and *D*. If so, PSR likely rejects *D* when testing *HD*|*B* 0 , and rejects *B* when testing *HB*|*A* 0 .

A reduction of the amount of FPs typically comes at the expense of enhanced FNs. In particular, FNs may occur if the real pattern is rejected in favor of one of its subsets or supersets. **Figure 6** shows, for a range of combinations of SIP size and injection count, the resulting level of FPs, FNs, and the maximum of the two (as a measure of overall performance) after performing each of the proposed PSR strategies. The significance spectrum used to determine significance for all realizations of the SIP models is the one for *w* = 3 ms shown in **Figure 4** (top right, dark-shaded entries). For the FPs shown in **Figure 6**, top row, the color-coded level refers to the fraction of simulations (out of 1000) containing one or more FPs. This measure takes values between 0 and 1, unlike the average FP counts shown in **Figure 5**. This representation simplifies the comparison with the average FN level, which ranges between 0 to 1 since here only a single spike pattern is injected in every simulation. To aid the comparison between the performances of PSF and PSR, gray dots mark those squares that correspond to models where the error rates exceed 5%. PSF on its own never performs well in terms of FNs and FPs simultaneously, while all PSR strategies yield a range of models for which both quantities are low. In summary, the relative improvement of PSR versus PSF shows that any PSR strategy reduces the FP rate considerably, while causing only a minor increase in the FN rate.

## **4.5. PERFORMANCE, HETEROGENEOUS FIRING RATES**

If neurons have the same spiking statistics, the spike pattern statistics depends on the pattern size only. Thus, the *p*-value of

**FIGURE 5 | Average number of FPs, distinguished by type, after PSF and PSR.** Average number of FPs obtained for different SIP models on *R* = 1000 model simulations. FPs are shown after performing PSF (top) and then PSR with combined filtering (bottom), and are distinguished by type (columns from left to right: FP supersets, FP subsets, FP overlapping, FP disjoint

patterns). Each panel shows the average number of FPs obtained for different SIP models, each corresponding to a square in the grid: the models differ by the SIP size (from 2 to 10; *x*-axis) and its injection count (from 2 to 10; *y*-axis). Other parameters (same for all simulations): *N* = 100, *T* = 3 s, *r* = 20 Hz, *w* = 3 ms, *K* = 5000, α∗ = 2 · 10−4.

(third row) indicates the combined error rate. Each matrix shows the performance for 81 different SIP models varying by SIP size (from 2 to 10, *x*-axis) and number of SIP injections (from 2 to 10, *y*-axis), of stationary and applied after PSF, from left to right: no filtering, subset filtering, superset filtering, covered-spikes criterion, combined filtering. Other parameters (same for all panels): *N* = 100, *T* = 3 s, *w* = 3 ms, *K* = 5.000, α∗ = 2 · 10−4.

each pattern is fully determined by the pattern signature. This does not hold when neurons have different spiking statistics, and in particular different firing rates. Here we discuss the case of heterogeneous firing rates across neurons, which are often present in electrophysiological data. Higher firing rates lead to a higher spiking probability per time bin. Patterns composed of neurons with higher firing rate are more likely to occur by chance, and are therefore less significant than patterns composed of neurons with lower rates. Thus, the *p*-values of patterns with the same signature (*z*,*c*) differ for different compositions of the firing rates. Pooling patterns by size and support in the pattern spectrum does not take into account the heterogeneity of firing rates across neurons and thus may lead to a biased statistics.

To investigate the robustness of our method against firing rate heterogeneity, we first simulate independent data consisting of 100 neurons, with a small population of neurons (2 to 10) firing at a higher rate (20 Hz) than the rest of the neurons (5 Hz). We simulate 1000 data sets of this type, and evaluate FPs in each of them by means of FIM and PSF (*K* = 5000 surrogates). In none of the simulations we detect significant signatures, i.e., FPs. The opposite scenario, where 2 to 10 neurons fire at 5 Hz and the others at 20 Hz, does not yield FPs as well. Thus, employing rate-preserving surrogates allows PSF to correctly estimate the significance of signatures under *H*0, also when rates are heterogeneous across neurons.

Next we study correlated data characterized by heterogeneous background firing rates. We investigate two cases based on a SIP model. In scenario S1, a pattern is injected in a set of neurons firing with lower firing rate (*rS* = 5 Hz) than the independent neurons firing at rate *rI* = 20 Hz (**Figure 7**, left column). In contrast, in scenario S2 the pattern is injected in neurons with higher firing rates (*rS* = 20 Hz, *rI* = 5 Hz; **Figure 7**, right column). In comparison to the homogeneous case where all neurons fire at 5 Hz (data not shown), the overall performance drops significantly, but does not so compared to the 20 Hz homogeneous case

**FIGURE 7 | Performance of PSR with heterogeneous firing rates.** Performance of PSR (combined filtering with parameters *h* = 1, *k* = 2) in terms of FP rates (*top* row), FN rates (middle row) and the combined error rates (maximum of FP and FN rates) (*bottom* row) of data with heterogeneous rates. Left column: SIP neurons fire at *rS* = 5 Hz, independent neurons fire at *rI* = 20 Hz. Gray dots mark entries where the error rate is above 5%. Right column: SIP neurons fire at *rS* = 20 Hz and independent neurons at *rI* = 5 Hz. Other parameters (same for all panels): *N* = 100, *T* = 3 s, *w* = 3 ms, *K* = 5.000, α∗ = 2 · 10−4.

(see **Figure 6**, right column). This is consistent with the previous finding that higher rates worsen the performance by shifting the detection border in the significance spectrum to the right (**Figure 4**, left vs. right). This also explains why FP and FN rates in scenario S1 are higher than in scenario S2: the average firing rate in the former ranges (depending on the SIP model) from 18.5 to 19.7 Hz, in the latter from 5.3 to 7 Hz. Our choice of using PSR with combined filtering leads to a better performance in this scenario than the covered spikes criterion (not shown). Taken together, these results indicate that the method can deal well with heterogeneity of firing rates without severe performance loss.

#### **4.6. PERFORMANCE, NON-STATIONARY FIRING RATES**

Now we want to consider the case when the firing rates of the neurons are not stationary in time. To explore the sensitivity of our method to non-stationarities we employ simulated data, again consisting of 100 parallel spike trains, which fire in two consecutive epochs of length *T*1 and *T*2 (the total simulation time *T* = *T*1 + *T*2 is 3 s, as in the data previously analysed) at different rates (*r*1 = 5 Hz and *r*2 = 20 Hz; or vice versa), homogeneously across the neurons in both epochs. In the first epoch, correlated activity is inserted by the SIP model. SIP of size 2 to 10, injected 2 to 7 times, amount to a coincidence rate of 1.33 to 4.66 Hz in the first epoch. The background rate is reduced correspondingly. For comparison, we also study the stationary case, where all neurons fire at *r* = 10 Hz. The performance for the three scenarios is shown in **Figure 8** (first column: *r*1 = 5 Hz, *r*2 = 20 Hz; second column: *r*1 = 20 Hz, *r*2 = 5 Hz; third column: *r*1,2 = 10 Hz). Although our analysis performs better (detection border more to the left) in the stationary case (*r* = 10 Hz; third column), it can still recover SIP activity with no FPs in a large portion of the parameter space, provided that rate-preserving surrogates are employed. As in the heterogeneous case, FPs increase when the SIP neurons have higher firing rates and thus more FP subsets occur. As apparent from **Figure 8**, bottom row, the method can correctly detect significant patterns in a wide range of models also in the presence of non-stationary rates. To study whether short transients in the firing rates tend to generate FPs, we repeated the analysis for *T*1 = 0.5 s, *T*2 = 2.5 s, setting first *r*1 = 5 Hz, *r*2 = 20 Hz and then *r*1 = 20 Hz, *r*2 = 5 Hz. In all cases we do not find enhanced FPs (data not shown), indicating that employing ratepreserving surrogates suffices to correct for rate non-stationarity in independent data.

## **5. DISCUSSION**

In this study we have presented a method to detect significant patterns of synchronous spiking in a subset of massively parallel spike trains in the presence of background activity. Our work is rooted in Picado-Muiño et al. (2013), where we demonstrated how to efficiently detect spike patterns in such data, and assess their significance under the null hypothesis of independent firing. Here we refined this significance test, which evaluates the significance of patterns using PSF on the basis of the pattern signature (size and support). PSF is prone to report FP patterns that arise due to the activation of an actual assembly mixed with chance synchrony because of background activity. To identify and remove these FP

comparison, the third column shows the performance for the stationary case with all neurons firing at rate *r* = 10 Hz, and a duration of *T* = 3 s. Other parameters (same for all panels): *N* = 100, *w* = 3 ms, *K* = 5.000, α∗ = 2 · 10−4. Gray dots mark entries where the error rate is above 5%.

detections, we introduced here PSR as an additional statistical testing step. As shown in **Figure 6** (second to last columns), PSR succeeds in eliminating FPs for a wide range of parameters, at the expense of a minor increase in FNs. A series of calibrations demonstrates the effectiveness of our approach under conditions of heterogeneous and non-stationary firing rates.

The relevance of higher-order correlations for information processing in the nervous system is hotly debated. Approaches based on maximum entropy models, such as Schneidman et al. (2006), suggest that higher-oder correlations contribute by a negligible fraction to the total network correlation, which appears to be dominated by pairwise correlations. However, it is important to stress that for correlations of a specific order, maximum entropy models estimate the overall magnitude of that correlation order, and are not sensitive to individual correlation structures of that order. Thus, the presence of a single group of correlated neurons with a certain size in the data is not enough for maximum entropy models to report significant correlation of the corresponding order. The study by Shlens et al. (2006) addresses this point, discussing that maximum entropy models may miss higher-order correlations because they overall contribute only by a negligible fraction to the total correlation. Besides, Roudi et al. (2009) showed that the statistical power of maximum entropy models describing spike correlations in heavily undersampled biological systems (such as parallel recordings with electrode arrays) is low. Despite these challenges, Ohiorhenuan et al. (2010) have shown using a maximum entropy model approach that in visual cortex local microcircuits exhibit evidence of higher-order interactions, whereas correlation statistics across long-range connections are explained on the basis of pair-wise interactions. However, methods designed to investigate individual spike patterns are needed to investigate the detailed structure of correlation in groups of spiking neurons.

A majority of current methods for spike correlation analysis limit themselves to fully synchronous patterns or to patterns with a specific size of typically low order [e.g., Grün et al., 2002a,b; Berger et al., 2007, 2010; Shimazaki et al., 2012]. Other approaches, such as CuBIC (Staude et al., 2010b), conclude on the presence of higher order correlations based on the statistics of the population activity without identifying the specific units engaged in such correlations. While Gansel and Singer (2012) presented a method for the detection of higher-order patterns, they identify pattern subsets by a purely heuristic procedure that is not accessible by analytic treatment, and that tests patterns directly, which requires a number of statistical corrections to avoid FPs (at the expense of FNs). Our proposed method instead first tests the significance of pattern signatures. PSF eliminates non-significant signatures based on surrogate data through the significance spectrum (see **Figure 4**), and determines the class *P* of associated significant patterns. Testing patterns on the basis of their signature rather than testing individual patterns reduces the number of required statistical tests to the number of signatures found in the data. We have shown that the composition of assembly and background spikes typically leads to the identification of additional significant patterns (i.e., FPs). In order to remove this type of FPs, we introduced here the PSR procedure that is based on conditional pairwise tests.

We have tested the performance of our analysis on artificial data where we embedded groups of synchronously spiking neurons in background activity of independent Poisson spike trains [SIP, cf. Kuhn et al. (2003)]. We studied the rate of FP and FN detections for occurrence rates of the synchronous pattern varying from 0.66 to 3.33 Hz, which reflect plausible values for the activation frequency of the assumed assemblies (Grün et al., 1999; Denker et al., 2010). The analysis shows in particular that by introducing PSR, assembly detection becomes possible with near perfect reliability and precision for a large range of SIP parameters. The transition shifts toward higher support and assembly size as the bin width or the firing rates increase (cf., **Figure 4**). Nevertheless, for physiologically realistic parameters only for very small or very infrequent SIP injections these patterns cannot be distinguished from chance synchrony. Moreover, evaluating patterns obtained from a larger set of simultaneously recorded neurons will have only minor impact on our findings due to a slight increase in the average size of observed patterns.

Non-stationarities of the firing rate in time or across neurons are a common concern faced by correlation analysis methods. The effect of non-stationary firing rates on PSF is two-fold. First, the surrogates used to calculate the significance estimates on pattern signatures should adequately reproduce the experimental rate profiles. Even if the underlying rate profile is not known, a variety of suitable approaches for surrogate generation is available for this task (Grün, 2009; Louis et al., 2010b). However, the sensitivity of detecting assembly activations is further affected by where these occur with respect to the rate non-stationarity. In this respect we tested the performance of PSF and PSR in a scenario of step-wise non-stationary firing rates where spike patterns were injected at selected rate levels only. Compared to the stationary case, the method retains a high performance for large parameter regimes (**Figure 8**), and shows only a slight increase in the number of FNs. For very large rate non-stationarities, a time-resolved analysis may be used to additionally aid the detection, as done, e.g., in the Unitary Events analysis (Grün et al., 2002b). In a similar framework, we found that also heterogeneous firing rates across neurons (**Figure 7**) exhibit a performance similar to the stationary case. While we see minor increases in the number of FPs, we remark that to a large extent these are indeed supersets of the injected pattern due to the high probability of gaining an additional coincident spike by chance from the set of neurons spiking at high rates.

In this study we assumed that assemblies occur at the time resolution of the data, i.e., that spike times of the assemblies are not jittered in time. In electrophysiological data this is a rare scenario, and instead spike synchrony typically occurs with a temporal jitter of up to several milliseconds [Grün et al., 1999; Pazienti et al., 2008]. In order to capture such slightly imprecise synchrony, exclusive binning is typically applied (Grün et al., 1999), where the bin width is chosen large enough to capture the jittered spike pattern. However, the spikes of the pattern may be split into adjacent bins with a probability that depends on the jitter, bin size, and pattern size. Therefore, the original synchronous events are destroyed, leading to increased FN rates (Grün et al., 1999). In **Figure 9** we show how this effect can have a substantial impact on the performance of the method. We applied PSF followed by PSR (combined filtering) on data where synchronous patterns are injected with a jitter of ±1 ms, and analysed with a bin width of *w* = 3 ms (left column) and *w* = 5 ms (right column). The performance drops considerably due to an increase of the FP rate for higher *z* and *c*, and an overall increase of the FN rate. The performance is slightly better for a bin width of 5 ms. Consistent with these findings, Grün et al. (1999) showed that for two parallel spike trains about 60% of the synchronous events are lost if the bin width corresponds to the jitter width. An earlier modification of exclusive time binning [multiple shift method, Grün et al., 1999] that avoids the splitting of jittered synchrony was not trivially applicable to large numbers of parallel spike trains. In Picado-Muiño et al. (submitted) we demonstrate how to implement a method for pattern detection based on the inter-spike distances rather than discrete time binning. This approach successfully detects jittered spike patterns and therefore trivially exhibits a performance in the context of PSF that is similar to that achieved in the absence of jitter (see Picado-Muiño et al., submitted, for details). Thus, it also complements the PSR framework presented in this study. Therefore, we suggest to detect jittered synchrony by the continuous detection method and perform the analysis by the proposed sequence of FIM, PSF, and PSR.

A further scenario that remains to be addressed in the future is unreliability in spiking activity that causes neurons to selectively skip participation in assembly activations. This scenario

was discussed in the context of the synfire chain model, where it was shown that stable propagation of synchronous spike packages through the network happens reliably although the probability that individual neurons participate in each activation of the synfire chain is lower than 1 (Diesmann et al., 1999). Selective participation may arise as a consequence of synaptic failure. The multiple interaction process [MIP; Kuhn et al., 2003] was proposed as a stochastic model implementing such a behavior. Our method would interpret the variable composition of spikes in a single MIP event as occurrences of multiple SIP events of lower support.

We conclude with a discussion of the practical implementation of the proposed analysis on data from electrophysiological recordings. Given a set of parallel spike recordings obtained at a resolution (i.e., binning) *w*, we choose the minimum pattern size *z*0 and the minimum pattern support *c*0 of the analysis. First, the spike data is binned and, using FIM, the CFISs and the corresponding pattern signatures are obtained from the transaction list. While this approach is feasible for the experimental data available today, with several hundreds of parallel recordings the computational effort may become too large. In this scenario, we suggest to pre-filter the data entering the analysis as suggested by Berger et al. (2010) before applying FIM on the reduced set of neurons. To monitor dynamic changes in the correlation structure of the activity, e.g., if assemblies are time locked to a particular behavioral event, one may choose to additionally perform the analysis in sliding windows.

Next, the significance of the observed patterns is evaluated by PSF under the null-hypothesis of full independence implemented by uncorrelated surrogate data. For experimental data, several techniques for surrogate generation based on stochastic sampling have been proposed in the past [for a review, see Grün, 2009]. Surrogates that preserve the firing rate profiles, such as spike dithering, seem most appropriate since PSF determines pattern significance based on the firing rates. Given the significance level α and *m* detected pattern signatures, a minimum of *K* = *m*/α surrogates are required to achieve the Bonferronicorrected significance level α∗ = α/*m* . Once the surrogates have been generated, we follow the procedure described for the simulated data. CFISs, pattern signatures and the resulting binary pattern spectrum are obtained for each surrogate run. Next, the *p*-value spectrum is obtained as an average of the binary spectra (see Section 2.2). The signatures whose *p*-values do not exceed the Bonferroni-corrected significance level α∗ are marked as significant, and the CFISs of significant signatures are collected into the class *P* of potential assemblies. Finally, PSR with combined filtering is performed to reduce *P* to a subclass *Q* of patterns which are mutually significant with respect to each other.

In summary, the use of FIM combined with the statistical tests described in this study and in Picado-Muiño et al. (submitted)

#### **REFERENCES**

Buzsaki, G. (2004). Large-scale recording of neuronal ensembles. *Nat. Neurosci.* 7, 446–451. doi: 10.1038/ nn1233

represents a powerful tool to extract candidate assemblies from experimental data. The method is statistically rigid, computationally feasible, robust against heterogeneity in the data, and powerful enough to deal with the limited amount of data typically available from electrophysiological experiments. We expect that our approach will help to reveal how precise spike synchronization observed by pairwise analysis in relation to behavior (Riehle et al., 1997) is manifested at the level of neuronal populations.

## **ACKNOWLEDGMENTS**

We thank Günter Palm for stimulating discussions. This work was partially supported by the Helmholtz Alliance on Systems Biology, BrainScales (EU Grant 269912), the Helmholtz Portfolio Supercomputing and Modelling for the Human Brain (SMHB), and the Spanish Ministry for Economy and Competitiveness (MINECO Grant TIN2012-31372).

## **SOFTWARE AND SUPPLEMENTAL MATERIAL**

The FIM library underlying the Python scripts with which we carried out our experiments is available at http://www.borgelt. net/pyfim.html. Python and shell scripts for related experiments as well as more extensive result diagrams are available at http://www.borgelt.net/accfim.html and http://www.borgelt.net/ cocofim.html. Please also consult http://www.spiketrain-analysis. org for these codes and further information on the analysis of parallel spike trains.

Staude, B., Rotter, S., and Grün, S. (2010b). Cubic: cumulant based inference of higher-order correlations in massively parallel spike trains. *J. Comput. Neurosci.* 29, 327–350. doi: 10.1007/s10827-009- 0195-x

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 31 May 2013; accepted: 11 September 2013; published online: 23 October 2013.*

*Citation: Torre E, Picado-Muiño D, Denker M, Borgelt C and Grün S (2013) Statistical evaluation of synchronous spike patterns extracted by frequent item set mining. Front. Comput. Neurosci. 7:132. doi: 10.3389/fncom.2013.00132*

*This article was submitted to the journal Frontiers in Computational Neuroscience.*

*Copyright © 2013 Torre, Picado-Muiño, Denker, Borgelt and Grün. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Correlations in background activity control persistent state stability and allow execution of working memory tasks

## *Mario Dipoppa1,2\*† and Boris S. Gutkin1,3\**

*1 Departement d'Etudes Cognitives, Ecole Normale Superieure, Group for Neural Theory, Laboratoire des Neurosciences Cognitives INSERM U960, Paris, France*

*2 Ecole Doctorale Cerveau Cognition Comportement, Université Pierre et Marie Curie, Paris, France*

*3 Centre national de la recherche scientifique, Paris, France*

#### *Edited by:*

*Robert Rosenbaum, University of Pittsburgh, USA*

#### *Reviewed by:*

*Carson C. Chow, National Institutes of Health, USA Zachary P. Kilpatrick, University of Houston, USA*

#### *\*Correspondence:*

*Mario Dipoppa and Boris S. Gutkin, Departement d'Etudes Cognitives, Ecole Normale Superieure, Group for Neural Theory, Laboratoire des Neurosciences Cognitives INSERM U960, 29 rue d'Ulm, 75005 Paris, France e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
#### *†Present address:*

*Mario Dipoppa, University College London, 21 University Street, London WC1E 6DE, UK*

Working memory (WM) requires selective information gating, active information maintenance, and rapid active updating. Hence performing a WM task needs rapid and controlled transitions between neural persistent activity and the resting state. We propose that changes in correlations in neural activity provides a mechanism for the required WM operations. As a proof of principle, we implement sustained activity and WM in recurrently coupled spiking networks with neurons receiving excitatory random background activity where background correlations are induced by a common noise source. We first characterize how the level of background correlations controls the stability of the persistent state. With sufficiently high correlations, the sustained state becomes practically unstable, so it cannot be initiated by a transient stimulus. We exploit this in WM models implementing the delay match to sample task by modulating flexibly in time the correlation level at different phases of the task. The modulation sets the network in different working regimes: more prompt to gate in a signal or clear the memory. We examine how the correlations affect the ability of the network to perform the task when distractors are present. We show that in a winner-take-all version of the model, where two populations cross-inhibit, correlations make the distractor blocking robust. In a version of the mode where no cross inhibition is present, we show that appropriate modulation of correlation levels is sufficient to also block the distractor access while leaving the relevant memory trace in tact. The findings presented in this manuscript can form the basis for a new paradigm about how correlations are flexibly controlled by the cortical circuits to execute WM operations.

#### **Keywords: correlations, background activity, working memory, spiking neural network, persistent activity**

## **INTRODUCTION**

Working memory (WM), defined as short term storage of information that is actively used on-line to carry out actions and decisions and drive learning, is one of the key processes that underpins our cognitive abilities. WM is characterized by an information bottleneck with resources restricting its "on-line" capacity to a relatively limited number of items at high levels of performance (Miller, 1956; Luck and Vogel, 1997; Cowan, 2001; Vogel et al., 2001) and a rapid decrease in performance with item number due to limited resource allocation (Wilken and Ma, 2004; Bays and Husain, 2008; van den Berg et al., 2012) as suggested by the recent experiments. Furthermore by its very nature, WM is characterized by the need to operate on the stored information rapidly. Such limitations and rapid operations of WM create the need for selective gating and rapid updating as well as active information maintenance to enable its immediate use (Frank et al., 2001). One of the central unresolved issues is how the multiple requirements for WM are carried out by the brain circuits: whether the maintenance, read-in, gating, and read-out are implemented by separated systems (e.g., as suggested by Baddeley, 2003) or by operations within the same neural circuit (e.g., as recently put forward by Machens et al., 2005).

Electrophysiological data from primate performing delayedresponse tasks show that persistent neuronal activity in prefrontal cortex (PFC) underlies the maintenance of WM: during the delay period between the stimulus presentation and the read-out, neurons selective to the memorized stimulus fire spikes at an elevated rate with respect to the resting state (Fuster and Alexander, 1971; Fuster and Jervey, 1981; Funahashi et al., 1989; Miller et al., 1996; Romo et al., 1999).

In order to highlight the unique requirements of the WM as a neural process let us focus on the DMS task with distractors as a prototypical example (Miller et al., 1996). In this task the subject must remember the identity of an item briefly shown (the sample) and respond correctly only when the item is shown again (match) all the while ignoring other items flashed (distractors). To execute correctly this task, the neural circuitry needs to perform three operations (**Figure 1A**): first, encode and maintain in memory the sensory stimulus during the delay period; second, robustly maintain the memory face to distractors presentation; third, erase the memory trace at task completion to make the store available again, given the limited WM capacity. These operations are translated in terms of neural activity as follows: item-related activity is turned on rapidly and selectively by the sample-stimulus, is

**FIGURE 1 | Outline of the models. (A)** Time sequence of the delay match-to-sample task for the working memory network. Active neurons are represented in full colors. Successively: (I) both populations are in a quiescent state, (II) sample stimulus (blue arrow) activates blue population, (III) the network prevent a distracting stimulus (red arrow) to activate the red population, (IV) match stimulus allows the read-out of the encoded memory in the blue population, and (V) persistent activity is erased in the blue population. **(B)** Correlations in external

background activity generated by a common source of noise, in addition to independent sources of noise. **(C)** Single unit network receiving shared and independent sources of noise. **(D)** Winner-take-all network with two competing excitatory populations coupled through one inhibitory population. In addition to independent sources the excitatory populations receive background activity by two different common noise sources. **(E)** Two-unit network with two excitatory populations receiving shared noise.

protected from distractors during the delay period, and is rapidly turned off on response by the match.

A number of spiking network models have been conceived to describe the neural substrate for WM where persistent activity is maintained by recurrent connections that allow for co-existing attractor memory states and a ground non-memory state (Amit and Brunel, 1997; Compte et al., 2000; Brunel and Wang, 2001; Gutkin et al., 2001; Laing and Chow, 2001; Machens et al., 2005; Miller and Wang, 2006; Ardid et al., 2010). In some of these models, protection from distractors and memory clearance are performed through the recruitment of inhibition (Compte et al., 2000; Brunel and Wang, 2001; Machens et al., 2005). As an alternative to the erasing-by-inhibition paradigm, it has been shown, in a spatial WM model, that a transient excitatory stimulus matching the memory trace "location" on the network extinguishes the persistent state by transiently synchronizing the spike-times of the neurons (Gutkin et al., 2001; Laing and Chow, 2001). This work, along with Machens et al. (2005) showed how the read-out and clear-out can be merged into a single operation. However, in these alternative frameworks, protection from distractors, or selective gating, was not addressed. Here we propose that the gating is obtained by flexibly controlling the spike-time structure of the WM network activity. In support of this idea, it has been shown that spike-time synchronization is modulated in association with cognitive processing (Abeles et al., 1993; Riehle et al., 1997; Funahashi and Inoue, 2000) and in particular in WM (Sakurai and Takahashi, 2006; Pipa and Munk, 2011).

Critically, WM trace appears in the context of on-going background activity. While background activity is not related to task parameters, this is not without structure. Correlations have been found broadly in spontaneous neural activity in the cortex (Tsodyks et al., 1999). In particular, it has been shown that nearby neurons receive common inputs from afferent neurons making their voltages correlated (Lampl et al., 1999). Effects of correlations have been widely studied for their effect on population code (Salinas and Sejnowski, 2001), to measure network connectivity (Aertsen et al., 1989; Cocco et al., 2009), on neural dynamics for coupled neurons (Ly and Ermentrout, 2009), and for multiple independent neurons (Galán et al., 2006; Moreno-Bote et al., 2008).

In computational models of WM the background activity has been largely seen as problematic for memory maintenance. For example one of the more sensitive technical issues addressed by several computational proposals is how to stabilize the WM trace in face of random background activity (Compte et al., 2000). The benefits of external input correlations on persistent activity in recurrent networks have only recently started to be addressed theoretically (Buice et al., 2010; Polk et al., 2012). For the specific case of line-attractor networks (modeling parametric WM) Polk et al. (2012) showed in a detailed analysis how properly tuned input noise correlations can promote stability of the persistent firing rate. This was further noted in Lim and Goldman (2012) who also showed that the correlation structure of background noise can suggest the optimal architecture of neural networks for short term memory performance.

Finally, in this article we examined the influence of input correlations on recurrent spiking networks, finding that the correlation level in fact may destabilize the persistent activity state, rendering it a slow transient state. Buice et al. (2010) used a path integral approach to integrate the effects of correlations and synchronization into a rate model of recurrent networks and examined the stability of the persistent state. For a bistable firing-rate network they noted that transient increases in input correlations (synchronizing noise input) can lead to a turnoff of the persistent activity. This approach may in fact provide an analytical framework of the observations we make in the present manuscript for recurrent spiking networks and the correlation-based control of the persistent state lifetime. In this manuscript we also go beyond noting that input correlations defined the lifetime of persistent activity; we show that input correlations can effectively control the access to WM by disallowing transient stimuli to initiate persistent activity. The functional consequences of these two effects are the central topic of this work.

To demonstrate that, by controlling the correlation-driven synchronization of the background activity it is possible to control the lifetime of the persistent state and to manipulate selectively the transitions in sustained activity and consequently to perform the required operations of the WM task, we first consider a minimal recurrent network. In this recurrent network the neurons receive an excitatory random background noise, and background correlations are induced by a common noise source. Then we implement a discrete item WM model where the modulation of the background correlation level sets the network into different regimes allowing for loading of memory, protection from distractors and memory persistence. In addition we show the possibility to merge the read-out and the clearance in a single operation since the presentation of the match stimulus can directly quench the persistent activity.

## **MATERIALS AND METHODS**

#### **NEURAL MODELS**

In this work we study recurrent spiking networks that show bistability between a ground state and an active persistent spiking state. Our goal is to construct and analyze a minimal network capable of showing the required bistability. Hence we consider networks of recurrently connected excitatory pyramidal neurons. The elements of the network are represented by non-linear "point" neurons that are sparsely connected by instantaneous excitatory recurrent synapses. The dynamics of a neuron's membrane potential *v* is described by the Quadratic Integrate and Fire (QIF) equation, which represents the normal form of type 1 spike generating dynamics (Ermentrout, 1996):

$$
\pi \frac{d\nu}{dt} = \nu^2 - b^2 + I\_{\text{syn}}(t) \tag{1}
$$

$$\nu(t) = V\_t \to V\_r \tag{2}$$

where τ represents the membrane time constant, −*b* is the resting potential, *I*(*t*) the input current, *Vt* a spike threshold, and *Vr* the reset membrane potential. The voltage of the neuron is scaled such that *v* is a non-dimensional variable. When the membrane potential neuron attains the threshold value *v* = *Vt*, a spike is emitted and a post-synaptic current (PSC) is transmitted to an output neuron. We set the parameters as follows: *Vr* = −20, *Vt* = 20, *b* = 1 and τ = 20 ms.

The input current to a given cell in the network is decomposed into three different components:

$$I\_{\rm syn}(t) = I\_r(t) + I\_s(t) + I\_{\rm ba}(t) \tag{3}$$

where *Ir*(*t*) represents the recurrent input due to other neurons in the network, *Is*(*t*) represents the input from external stimuli directed to the network, and *I*ba(*t*) represents a non-specific background activity. Each of the three currents corresponds to a sum of PSCs originating from synaptic inputs generated by the presynaptic neurons at times *tn*. The PSCs are modeled with delta pulses:

$$I(t) = \sum\_{a} \sum\_{\{t\_n\}} J\_a \pi \delta(t - t\_n) \tag{4}$$

where *Ja* represents the synaptic strength for a given connection and could be positive (corresponding to an AMPA synapse) or negative (corresponding to a GABA synapse).

#### **BACKGROUND ACTIVITY AND CORRELATIONS MEASURES**

Ample data shows that cortical neurons receive a large amount of non-specific cortical and subcortical inputs whose structure is not directly related to the specific task and stimulus [e.g., see Shadlen and Newsome (1994) and summary of data in Amit and Brunel (1997)]. We refer to this type of input as an external background activity. It is taken to be composed of sequences of excitatory PSCs of synaptic strength *J*0 and with the synaptic times generated by a Poisson process. The synaptic currents are depolarizing in accordance with the notion that cortical neurons receive inputs from long-range excitatory glutamatergic projections.

In our model, this background activity can be either unstructured (uncorrelated) or structured (correlated). The correlation level, between two spike trains *Si*(*t*) and *Sj*(*t*) is given by:

$$
\lambda\_{\vec{\imath}\vec{\jmath}} = \frac{1}{\langle S\_{\vec{\imath}}(t) \rangle} \int \text{CCV} \mathbf{F}\_{\vec{\imath}\vec{\jmath}}(s) ds \tag{5}
$$

where CCVF corresponds to the cross-covariance function (Brette, 2009). This function is normalized to zero if *Si*(*t*) and *Sj*(*t*) are generated by independent Poisson processes.

We consider two ways for constructing the background activity: