#### *Uncorrelated background activity*

All *N* neurons receive spike trains generated by *N* independent channels with rate ν0. This leads to CCVF(*s*) = 0 and thus the correlation level is λ*ij* = 0.

#### *Correlations induced by a common source of noise (Figure 1B)*

All the *N* neurons receive inputs both from independent channels, with frequency (1 − λ)ν0, and from a common channel, with frequency λν0 and 0 ≤ λ ≤ 1. Each channel generates a spike train with Poisson statistics. The average background input rate is ν0 for each neuron. The cross-covariance function is then CCVF(*s*) = λ*Si*(*t*)δ(*s*) and the correlation level is λ*ij* = λ. This gives purely spatial correlations.

We measure the correlation level of the synaptic input among cells in the network with the mean Pearson correlation coefficient. We first compute a running mean (averaged over a time window of 5 ms) of the synaptic input *Ii a*(*t*) for each cell during a certain interval of time. Then we compute the Pearson correlation between the synaptic input of two cells:

$$\rho\_{ij} = \frac{\text{cov}(I\_a^i, I\_a^j)}{\sigma(I\_a^i)\sigma(I\_a^j)} \tag{6}$$

Finally we compute the average over all the cell pairs of the network ρ = [2/*N*(*N* − 1)] *Ni* = 1 *N j* = *i* + 1 ρ*ij*. In particular, in **Figure 4**, we performed this measure for the recurrent input *a* = *r* and background input *a* = *ba* .

#### **FUNCTIONAL NETWORK STRUCTURES IMPLEMENTING WM TASKS**

In this work we study three different networks. We start out by studying a homogeneous network of recurrently coupled excitatory neurons. This network can be also thought of as a encoding a single item of WM: a "single-unit network". The second model consists of two homogeneous excitatory networks coupled together through a population of inhibitory neurons: a "winner-take-all network" of two discrete competing short-term memory items. The third model is made up of two recurrent excitatory populations without mutual connections: a "two-unit network".

#### *Single-unit network*

A homogeneous network with *N* = 100 identical sparsely coupled neurons is represented in **Figure 1C**. Each neuron in the network receives synaptic inputs from *cN* other excitatory neurons, where *c* = 0.2 is the probability of connection, and *J* = 0.26 is the recurrent synaptic strength [described in Equation (4)]. Neurons receive excitatory inputs also from external background activity, with synaptic strength *J*0 = 0.151 and firing rate ν0 = 106 Hz. Neurons also receive an excitatory input from external sensory stimuli with synaptic strength *J*1 = 1.5 and firing rate ν1 = 56 Hz for a duration of 50 ms, as will be described hereafter. Parameters of the network are chosen such that the network sustains a quiescent state, with low firing rate (*f* < 5 Hz), and a persistent state, with high firing rate (≈20 Hz).

#### *Winner-take-all network*

The second model is a reduced version of the network proposed by Amit and Brunel (1997) (**Figure 1D**). The network is composed of two excitatory populations and one inhibitory neural population. Each of the two excitatory populations has *NE* = 40 neurons, and the third population is made up of *NI* = 20 inhibitory neurons. An excitatory neuron receives synaptic inputs from *c*EE*NE* (*c*EE = 0.45) neurons of the same population, with synaptic strength *J*EE = 0.3, and from *c*EI*NI* (*c*EI = 0.35) inhibitory neurons with synaptic strength *J*EI = −0.25. An inhibitory neuron receives synaptic inputs from *c*IE*NE* (*c*IE = 0.34) excitatory neurons with synaptic strength *J*IE = 0.05 from each excitatory population. Other parameters of the network are: *J*0 = 0.4, *J*1 = 1.5, ν0 = 60 Hz, and ν1 = 17 Hz. In addition we augment the mutual inhibition network with the added feature to control the amount of correlated noise in each excitatory population. More precisely each excitatory population receives background activity by common noise sources in addition to independent sources. In such a way the correlation level λ is regulated independently in each excitatory population.

#### *Two-unit network*

We devised a third version of our network models that is made of two excitatory independent populations, each one making recurrent connections with itself. Both networks share a common excitatory noise source projecting simultaneously to all excitatory neurons in addition to the independent uncorrelated background noise (**Figure 1E**). Since the common noise source is shared between the two populations the correlation level λ varies equally in the two excitatory populations. The parameters of each excitatory population are those given for the single-unit network, with the only difference that we used here a larger network (*N* = 1000) and we scaled accordingly the recurrent synaptic strength (*J* = 0.026).

#### **DELAYED MATCH-TO-SAMPLE TASK**

We study the spike-timing based mechanisms able to implement the DMS task (**Figure 1A**). The sequence of operations and the neural dynamics aim to reproduce the experimental results of Miller et al. (1996). For illustrative purposes, the discrete items can be viewed as corresponding to colors. The activation of one excitatory population encodes color blue, that we define population *B*, while the other encodes color red, that we define population *R*. If the populations are both in a quiescent state, the state of the network represents the absence of color information. For simplicity we represent the spontaneous state as the quiescent state (average firing rate ≈0 Hz).

During the task, the animal has to maintain a memory of an item (a color) during the delay period. In terms of neural activity, the corresponding excitatory population should be activated and maintained in a persistent sate. Additionally the model should protect the memory from the presentation of a distractor stimulus. At task completion after the decision, the system should erase rapidly the memory, i.e., the persistent activity should be deactivated to its quiescent state.

To establish that a network performs a WM task correctly we require it to perform all the operations of the task. The first operation, *load*, corresponds to loading the memory by the sample signal, and corresponds to *B* that is activated in a persistent state while *R* is in a quiescent state. The second operation, *protect*, corresponds to the maintenance of the blue item memory in the face of the distractor presentation. In terms of activity it corresponds to *B* maintained in the persistent state and *R* that is not activated to the persistent state even when the red stimulus is presented during the delay period. In networks (**Figure 1E**) where population *B* and population *R* are not connected, the operation *protect* can be separated in two independent sub-operations: *maintain* (maintain item memory in population *B*) and *block* (prevent activation of population *R*). The third operation, *clear* corresponds to the clearance of the memory encoded in the network. This is equivalent to the erasing of the persistent activity in the network. Note that in this work we do not focus explicitly on the read-out mechanism following the presentation of the match stimulus.

In particular in the winner-take-all network (resp. the twounit network), operation *load* is executed with success if the sample stimulus activates population *B*. This is measured before distractor presentation during 350–450 ms (resp. 350–450 ms): ν*B* > 5 Hz and ν*R* < 5 Hz, where ν*B* < 5 and ν*R* denote the average population firing rates of populations blue and red, respectively. Operation *protect* is executed with success if population *B* maintains the persistent state and population *R* is not activated. This is measured before match presentation during 750–850 ms (resp. 700–800): ν*B* > 5 Hz and ν*R* < 5 Hz. Operation *clear* is executed with success if population *B* is deactivated at task completion. This is measured during an interval after match presentation, during 1150–1250 ms (resp. 1050–1150 ms): ν*B* < 5 Hz and ν*R* < 5 Hz.

#### **NUMERICAL ANALYSIS**

All the numerical results are obtained by algorithms run in Python. The differential equations are integrated with Euler steps of *dt* = 0.1 ms. The mean population firing rate *f* is computed over population average in 10 ms.

Data points for networks and associated error bars are computed by averaging over simulated individual network realizations. We generated random connectivity matrices such that every neuron receives the same number of input connections. For a fixed network connectivity matrix, we computed the average over 100 realizations of background activity and stimuli for each of 30 random realizations of the network connectivity matrix when not otherwise stated.

## **RESULTS**

## **EFFECTS OF CORRELATIONS ON PERSISTENT ACTIVITY STATE IN THE SINGLE-UNIT NETWORK: ERASING AND BLOCKING THE MEMORY TRACE**

We examine how correlations in the background activity control selective persistent activity in WM networks. Hence we start out by analyzing how background correlations affect the transitions between the quiescent and self-sustained states in our network model.

Correlations in background activity are generated by the addition of a common noise source to independent stochastic channels (see **Figure 1B**). By changing the relative firing rate of the common source with respect to the independent channels we control the correlation level λ. We set two different protocols represented in **Figure 2A**. In the first protocol, the correlation level is increased instantaneously from λ = 0 to some value λ > 0 at 500 ms. Therefore given that the stimulus activates the persistent state, this protocol allows us to test the effects of the correlations on the probability that the active state is erased and we refer to it as the erasing protocol. In the second protocol the correlation level is set λ > 0 for all the time, before the transient stimulus appears. In this way it is possible to see the effect of correlations on blocking the ability of the stimulus, presented during 50–100 ms, to initiate the persistent state and we refer to it as the blocking protocol.

We first demonstrate the prevalent effect of correlations: control of active memory state and control of access to the memory. In an example of the erasing protocol the excitatory stimulus activates the network into a persistent state; at 500 ms correlations are increased and the persistent state is disrupted (**Figure 2A**). In an example of the blocking protocol the excitatory stimulus is not able to activate the persistent state (see **Figure 2B**). In order to understand how these effects depend on the activity parameters we ran a large number of simulations where we injected background activity with different correlation levels for 0 ≤ λ ≤ 1 to networks with different connection probability *c* to measure how this effect spread thanks to the network architecture (**Figure 2C**). We compared networks with the same scaled synaptic strength *J* such that *cJN* = const. = 5.2. In the erasing protocol we estimated the erasing probability *Pe*(*c*, λ) defined as the probability for the network to have the firing rate ν < 5 Hz in the interval 800–900 ms. We discarded trials where the network is not in a persistent state (ν > 5 Hz during 400–500 ms). In the blocking protocol we estimated the probability that the correlations block the stimulus; the blocking probability *Pb*(*c*, λ) defined as the probability for the network to have the firing rate ν < 5 Hz in the interval 400–500 ms. This could also be seen as a gating of the persistent activity. We observe that for both protocols, the increase of both *c* and λ disrupts the persistent state: in the first case by erasing it and in the second case by blocking its activation.

We hence wanted to assess how the network size influences the stability of the persistent state under the various background activity regimes (**Figure 2D**). We compared networks with different size *N* with an equal average synaptic input *cJN* = const. = 5.2. We measured both erasing probability *Pe*(*N*, λ) and blocking probability *Pb*(*N*, λ) as a function of λ. We observe that both the erasing probability and the blocking probability (*Pe* and *Pb*) increase with the network size *N*. We observe that both for fixed *c* and for fixed *N Pb* > *Pe*. Finally we studied the probabilities *Pe* and *Pb* as function of *N*, fixing both *J* = 0.26 and the number of inputs that each neuron receives, i.e., *cN* = cost = 20. We computed these probabilities averaging over 500 trials. We found that with such a scaling both *Pe* and *Pb* are approximately constant (**Figure 2E**).

In order to determine whether the optimal stimulus parameters to load of a memory (or activation of a persistent state) depend on the correlation strength we measured the loading probability (1 − *Pb*) as function of the stimulus strength (ν1) and for different values of λ (**Figure 3A**), we computed the probabilities of **Figure 3** averaging over 300 trials. Different values of λ change the amplitude of (1 − *Pb*) but do not shift the tuning with respect to ν1. We also found that there are two peaks of (1 − *Pb*): one at about ν1 ≈ 20 Hz and another at about ν1 ≈ 50 Hz. To test whether the positions of the two peaks depend on the recurrent network properties we measured the loading probability as function of ν1 and for different values of the recurrent synaptic strength *J* (**Figure 3B**). Similarly to the previous results, different values of *J* change the amplitude of (1 − *Pb*) but do not shift the peaks of the curves with respect to ν1. In summary this indicates that indeed the strength of the stimulus required to active the persistent state with a set probability is dependent on the background correlations, and yet the tuning is rather broad.

To further investigate the effect of correlation on the stability of the persistent state we determined the lifetime of the sustained activity and the level of correlations prior to the erasing time. We defined the end of the persistent state *t*stop (magenta vertical line,

**Figure 4A**) as the first period of 10 ms (after the correlation onset) during which the firing rate of the network falls below 5 Hz. Noticing that in most of the trails a peak of activity was preceding the persistent state erasing, we defined the time of such a peak *t*peak (black vertical line, **Figure 4A**) as the last period of 10 ms before *t*stop that the firing rate attains a local maximum (in time) and that is beyond 20 Hz. We determined for each trial where the persistent state was not erased before the onset of correlations *t*corr = 800 ms (red vertical line, **Figure 4A**), the interval *t*c.p. = *t*peak − *t*corr. We determined the interval *t*p.s. = *t*stop − *t*peak. We performed this protocol for three different values of the correlation level: λ = {0.3, 0.6, 0.9} (**Figure 4B**). We found that the distribution of *t*c.p. decreases with time for all the values of

λ = 0.8 all the time (yellow shaded areas). The network receives an excitatory stimulus during 50–100 ms. The stimulus fails to activate the

> correlations. When the level of correlations is larger (**Figure 4B**, top) the probability of reaching the peak earlier in time slightly increases with λ. Furthermore we found that the interval between the peak of activity and the erasing of the activity in the network is narrowly distributed in time. Finally this interval is independent of the correlation level, meaning that the correlations do not have a strong effect on this timing (**Figure 4B**, bottom). We computed these distributions averaging over 500 trials.

fixed λ = 0.6 and *J* = 0.26 and with *c* scaled such that *cN* = *const* = 20. *Pe* and *Pb* remain approximately constant.

The mean Pearson correlation coefficient ρ (see Materials and Methods) of the synaptic input in the network during the interval with uncorrelated background activity was compared with that during the interval with correlated background activity just preceding the peak. Only trials where *t*peak − *t*corr > 100 ms were considered. The interval with uncorrelated background activity is defined as the 100 ms preceding *t*corr (gray shaded area, **Figure 4A**). The interval with correlated background activity is defined as the 100 ms preceding *t*peak (red shaded area, **Figure 4A**). We computed ρ for the background input (red lines, **Figure 4C**) and for the recurrent input (black lines, **Figure 4C**) both for the uncorrelated interval (dashed lines) and for the

strength (ν1) for different values of the correlation level (λ). **(B)** Probability of activation as function of the stimulus strength (ν1) for different values of the recurrent synaptic strength (*J*).

correlated input (continuous lines) when λ = 0.6. We found that ρ during the correlated input is smaller in the recurrent input with respect than in the background input. However, ρ of the recurrent input is larger during the correlated interval than during the uncorrelated interval. Interestingly we found that during the correlated interval, while the correlation coefficient of the background input increases with λ (**Figure 4D**, bottom), the correlation coefficient of the recurrent input instead remains approximately equally distributed when λ is changed (**Figure 4D**, top). This suggests that the network has reached the maximal amount of sustainable correlations before turning off.

To understand whether the persistent activity deactivation is caused by an increase of spike synchrony we tracked the synchrony of the spike times using the multivariate SPIKE-distance measure *S* (Kreuz et al., 2013) (**Figure 5A**). The spike synchrony is given by 1 − *S* spanning the values between 0 (no synchrony) and 1 (perfect synchrony). We compared the average spike synchrony during two intervals, similarly to **Figure 4**: the first interval corresponds to the 100 ms preceding the start of correlated background activity and the second interval corresponds to the 100 ms preceding the last peak of activity before the deactivation of the persistent activity (provided that the onset of this last interval does not precede the start of the correlated background activity). The distribution of the average value of 1 − *S* (computed

**FIGURE 4 | Persistent state suppression is preceded by an increase of recurrent correlation. (A)** Timing outline of erasing persistent activity. Mean recurrent input (blue trace), mean background input (red trace), and external input (black trace) are represented together with the time at which activity is erased (*t*stop, magenta vertical line), time of the last peak of activity before erasing (*t*peak, black vertical line), and time of the onset of the correlations in background activity (*t*corr, red vertical line). The interval of uncorrelated activity and the interval of correlated background activity during which the correlation coefficient is measured

in panels **(C)** and **(D)** are represented with gray and red shaded areas, respectively. **(B)** Distribution of the interval *t*c.p. = *t*peak − *t*corr (top) and distribution of the interval *t*p.s. = *t*stop − *t*peak (bottom) for different values of λ. **(C)** Mean Pearson correlation coefficient (ρ) of the recurrent input (black traces) and of the background input (red traces), computed during the uncorrelated interval (dashed traces) and the correlated interval (continuous traces). **(D)** Coefficient ρ of the recurrent input computed during the correlated interval for different values of λ (top). Same analysis for the background input (bottom).

over 2000 trials) during these two intervals shows that there is a weak increase of spike synchrony preceding the persistent activity deactivation with respect to the case of uncorrelated background activity (**Figure 5B**). This weak increase could indicate that few spike coincidences might be the cause of the persistent activity turning off.

## **EFFECTS OF BACKGROUND ACTIVITY CORRELATIONS IN A WINNER-TAKE-ALL NETWORK**

We show here that modulating appropriately in space and time the correlation level of the background activity in a network performing a WM task significantly improves correct execution of all the required operations: *load*, *protect*, and *clear*.

We compared two different versions of the winner-take-all network; each made of two excitatory populations *B* and *R*, representing respectively colors blue and red. The two populations interact via a third population of inhibitory neurons that creates a winner-take-all mechanism. The two versions differ in that the first one receives only uncorrelated background activity while in the second each excitatory population receives also background activity from a different common noise source (**Figure 1D**).

We fixed the stimuli sequence as follows: during 50–150 ms a sample blue stimulus excites population *B*; during 450–550 ms a distractor red stimulus excites population *R*; during 850–950 a match blue stimulus excites again population *B*. The operations that the network has to do are to load the blue item in memory, to protect the memory at red item presentation, and to clear the memory after the match presentation.

Brunel and Wang (2001) pointed out that in order to perform the DMS task correctly, the distractor stimulus strength needs to be controlled with care: above a certain strength persistent memory-trace is perturbed by the distractor. For our case we suppose that it is reasonable to assume that all sensory stimuli in the task are of the same strength. As a preliminary test we want to confirm that in absence of background correlations the network without common noise source does not perform efficiently when the stimuli are too strong, as was already stated in the reference network described by Brunel and Wang (2001). In the example shown in **Figure 6A** the distractor activates *R* and via the inhibitory population the persistent state in *B* is deactivated leading to a failure of the operation *protect*.

We then consider the network represented in **Figure 1D** that allows the modulation of the correlation level λ in each excitatory population independently. The network initially receives uncorrelated background activity to λ = 0.9. After the first item has been loaded, the system increases the correlation level in the background activity of the other non-activated population *R*. After the match stimulus has been presented, the correlation level is increased also in population *B* to λ = 0.9.

We show an example of the network executing the WM task where the correlation level is modulated independently in the excitatory populations (**Figure 6B**). In this particular example we illustrate a trial where the network performs the required operations of the WM task (compare with **Figure 1D** and see below for statistics across trials). The distractor excites *R* only transiently such that excitation does not last enough to disrupt the activity in *B*, in addition as shown below this happens also for strong distractor stimuli. Therefore the operation *protect* has been executed with success and the memory is maintained. At the end of the match stimulus the persistent activity is disrupted also in population *B* caused by the increase of λ in that population too. Therefore the operation *clear* is executed with success and the memory is erased in the network. This example illustrates that the success of the operations *protect* and *clear* in the network with correlations are not due to the presence of inhibition as was set in the model of Brunel and Wang (2001).

To get quantitative measures of performance for these two networks (with and without correlations in background activity), we analyzed the statistics of *load* and *protect* performance, as a function of the stimulus intensity ν1 (and thus its strength) (**Figure 7A**). We consistently find higher *protect* performance for correlated background activity than for uncorrelated background activity throughout the whole range for ν1. In fact the success of the *protect* operation depends only gradually on

**FIGURE 6 | Selective correlations implemented in a working memory task.** Two competing populations network, with two item-selective excitatory populations (blue and red) and one inhibitory non-selective population (black). **(A)** Without background correlations, the distracting stimulus activates population *R* and population *B* is deactivated. **(B)** After the activation of *B* at 150 ms, a common source of noise increases the correlations in background

activity (λ = 0.9) in *R*. The correlations block the activation of *R* and maintain the persistent state *B*. After the completion of the task at 950 ms the correlations erase persistent activity in *B*. (Top) Raster plot of the neural activity in the task. (Bottom) Successively: sample stimulus to *B* (50–150 ms), distracting stimulus to *R* (450–550 ms), and match stimulus to *B* (850–950 ms).

optimal value for the network with correlations. **(B)** The performance of the network is measured on four different operations for the optimal value of

probability of preventing memory disruption by a distracting stimulus in the protocol with correlations (protect), probability of erasing of the memory at the end of the task (clear). All probabilities have a high value showing that the network has good task performance, above chance.

the distractor strength. On the other hand in order to perform operation *protect* above chance level in the network with uncorrelated background activity distractors should be carefully adjusted to have intensity ν1 < 5 Hz. However, in this range the operation *load* is suboptimal. Hence the uncorrelated model fails in the task. This fact illustrates a recurrent problem in the protect-by-inhibition paradigm: it needs finetuning and achieves only low performance if the stimuli are too strong. Instead, using correlations as a mechanism to protect the activity does not need precise fine-tuning as can been seen in the large range in which both *load* and *protect* are well above chance level. We found the value ν1 = 4.8 Hz maximizes the joint probability of executing with success *load* and *protect* (**Figure 7A**, vertical dashed line). We show in **Figure 7B** probability of success of the three operations *load*, *protect*, and *clear* finding that all of them score a value higher than chance level.

## **IMPLEMENTING WORKING MEMORY TASK BY FLEXIBLE CORRELATIONS MODULATION**

We now go on to show that mutual inhibition is not a required mechanism for implementing the WM task. We show here that modulating appropriately the background activity correlations in time in a network without inhibitory population allows correct execution of all the required WM operations: *load*, *maintain*, *block*, and *clear* ( Please note that since the network studied here is made of two separated excitatory populations the component *maintain* and *block* of the operation *protect* can be treated separately).

### *Network operating regimes*

In order to characterize the network performance statistics during the task we need to track three probabilities. The first probability *P*g.o. = *PePb* corresponds to the joint probability of deactivating by correlations the network that is in the persistent state (erase) and to block activation of the network that is in the quiescent state and is excited by a stimulus. When *P*g.o. dominates over the other probabilities the system is in a gate-out regime, i.e., memory cannot neither be loaded nor maintained in the network. The second probability *P*g.i. = (1 − *Pe*)(1 − *Pb*) corresponds to the joint probability that, despite the correlations, the network maintains the persistent state, if previously activated, and that the stimulus activates the persistent state when the network is in the quiescent state. When *P*g.i. dominates, the system is in a gate-in regime, i.e., the memory can be loaded and maintained in the network. Finally the third probability *P*s.g. = (1 − *Pe*)*Pb* corresponds to the probability of maintaining the persistent activity in the presence of correlations while blocking the activation of a persistent state with correlated background activity when the system is in a quiescent state and is excited by a stimulus. When *P*s.g. dominates the system is in a selective-gate regime, i.e., the memory is maintained but cannot be loaded. In a sense we want to show that correlations in the background activity can selectively switch the network from the gate-in regime at the outset of the task, to the selective-gate regime during the memory period. We do not consider the probability *Pe*(1 − *Pb*).

To obtain the network performance on the DMS task we considered the statistical results presented in **Figure 2** for the single excitatory population (*N* = 1000 and *c* = 0.2). We note that there is a difference between the erasing probability *Pe*(λ) and the

Joint probability of (not-)erase and (not-)block: *P*g.o.(λ) (black curve), *P*g.i.(λ) (blue curve), and *P*s.g.(λ) (red curve). Gate-out regime corresponds to domination of *P*g.o.(λ) and falls in the range λ > 0.11. Gate-in regime corresponds to dominance of *P*g.i.(λ) and falls in the range λ < 0.04. Selective-gate corresponds to dominance of *P*s.g.(λ) and falls in the range 0.04 <λ< 0.11, with maximal value at λ = 0.07 (red star).

blocking probability *Pb*(λ) in function of the correlation level. In **Figure 8** we present the results for network we consider in this manuscript. We see that when *P*g.i.(λ) dominates (λ < 0.04), the system is in a gate-in regime, i.e., the memory can be loaded and maintained in the network (**Figure 8**). When *P*s.g.(λ) dominates (0.04 <λ< 0.11) the system is in a selective-gate regime. We do not consider here the gate-out regime that corresponds to *P*g.o.(λ) dominating over the other probabilities (λ > 0.11 ). We set the gate-in regime at λ = 0 and the selective-gate regime at λ = 0.07.

### *Modulation of correlation level in time*

We now show that correlations induced by a global common noise source to the whole network executes the DMS task efficiently by modulating the correlation level λ during the different phases of the task. We note that in the mutual inhibition model, at task completion, increase in the correlation level induces the gateout regime and erases the memory. We show here, in a two-unit model (**Figure 1E**), how the presentation of the match stimulus directly can erase the memory thereby implementing a direct match-based suppression without requiring inhibition. In this model each of the two excitatory populations receives background activity from sources independent to each neuron and from a noise source common to all neurons.

An example of the network performing the DMS task is represented in **Figure 9A**. The stimuli are presented in the following sequence: sample stimulus to population *B* at time during 100– 150 ms, distractor stimulus to population *R* during 450–500 ms, and match stimulus to population *B* at time 800–850. In the beginning the system is in the gate-in regime (λ = 0): the sample stimulus activates *B*. From 300 ms the network is set in a selective-gate regime (λ = 0.07): the distractor stimulus activates transiently population *R* while persistent activity is maintained in population *B*. At the end of the task the match stimulus, first, increases the activity in *B* and, then, destroys it.

We can then compute the task performance of the network, corresponding to the success rate that the operations *load*, *maintain*, *block*, and *clear* are executed successfully (**Figure 9B**). These measure are all above chance level. Notice the high level of performance of the *clear* operation.

### **DISCUSSION**

#### **RESULTS AND DATA DISCUSSION**

In this work we present a novel paradigm explaining how the persistent activity can be modulated on-line by the mean of both information-related signal and background activity. This paradigm is based on our result showing that background correlations influence the transition between the persistent state and the quiescent state in a bistable recurrent neural network. We call this phenomenon correlation-induced gating.

In order to implement a multi-unit network performing a WM task, we began by establishing the basis of the correlation-induced gating on a single-unit network. We show that background correlations block and erase a persistent state in a homogeneous recurrent neural network representing a single unit. We found that the transition rate from the persistent state to the quiescent state increases, with the network size and with the connection probability. For all situations the probabilities increase

with the correlation level. Increasing the network size, while fixing the connections probability and renormalizing the synaptic inputs to keep the average input strength constant scales up the probabilities. In other words, in larger size network with weaker but more numerous synapses, correlations appear to have a stronger effect. On the other hand, when we fix the total number of synapses each neurons receives, growing network size does not appear to have much effect on the correlation driven probabilities. These effects could be related to the fact that the amount of correlation between neurons sharing common input is mainly determined by pooling (Rosenbaum et al., 2010, 2011).

We implemented a winner-take-all network composed of two excitatory populations and one inhibitory population. Each of the excitatory populations receives background input from independent noise sources and a noise source common to the neurons of such population. The amount of correlation could be changed independently in the two excitatory populations. By increasing the level of correlations in the populations encoding an irrelevant information we prevented a distractor from loading a memory item in such population. In particular we showed that this model allows to prevent stronger distractors with respect to a model inspired by Brunel and Wang (2001) where the distractor is blocked only by the mutual inhibition. Our model could therefore explain how the response to the distractor stimulus in a WM task could be as strong as for the sample stimulus (Miller et al., 1996). This effect would be in fact not compatible with a model where a distractor is prevented by mutual inhibition.

We implemented a WM network differing by the same previous one in the construction of background correlations that are induced by a shared source. We show that modulating the correlation level in background activity we can set the system in the different regimes. This time instead of modulating the correlations level "in space" we modulate it in time. Depending on the strength of the correlations the system is set in different operating points, namely the gate-in, selective-gate, and gate-out regimes. The gate-in regime allows to load a memory in the WM store and to maintain it subsequently. The selective-gate regime maintains a previously loaded memory but blocks the load of any new memory. The gate-out regime blocks the network both to load and to maintain a memory. We can switch instantaneously from one dynamic regime to the other by tuning the strength of the background activity correlations. We further show that the projection of a strong match stimulus can be sufficient to clear the memory at task completion, thereby suggesting that correlations also play a role in match-suppression.

We must also note that in this work we considered spatial correlations and their effect on the persistent activity and WM task executions. In a companion paper we have shown that temporal structure also has an important effect: the gating modes are modulated by the oscillatory frequency content of the background activity (Dipoppa and Gutkin, 2013). While in the companion paper the block and erase probabilities, and thus the gating modes of the network, are modulated by the oscillation frequency, in this work they are modulated by the correlation level. As opposed to the non-monotone relationship between the oscillation frequency and the block and erase probabilities [Figure 3A of Dipoppa and Gutkin (2013)], there is a monotone relationship between the correlation level and the same measures (**Figures 2C,D**). Hence control of the WM through spatial correlations could be implemented by a simple increase or decrease of activity within a neural population furnishing connections common to the WM store, while the oscillatory control would require more complex task-dependent shifting between the frequency bands. We would like to speculate that the two mechanisms could represent two independent modes of control over the WM networks. Furthermore, in this work we uniquely examine the role of mutual inhibition and show that the spatial correlation structure alleviates the network sensitivity to stimulus strength.

Although we do not propose a mechanism for the readout of the memory information we note that the mechanism proposed by Brunel and Wang (2001) for their network would be compatible with our model. This mechanism corresponds to the fact that a match stimulus will elicit a stronger response with respect to a distractor stimulus in the first few tenths of milliseconds since the first will excite a network that is already in the persistent state. Hence we might speculate that a complementary population of neurons sensitive to rapid transients in the activity might be a way to signal read-out differentially.

#### **MODEL PREDICTIONS AND OPEN QUESTIONS**

The novel paradigm that we present here allows to manipulate persistent activity through background correlations. An advantage of the correlation-induced gating with respect to the inhibition-induced gating is that the gate can be rapidly and flexibly opened or closed depending on the correlation level, instead of being fixed by the network connectivity structure.

The effects that we find for the flexible changes in the correlation levels is the major prediction of the model. We predict that an examination of multi-unit electrophysiological recordings of animals performing a WM task will show the following modulation of correlation level: low level during loading and intermediate level during maintenance (as in the two-unit model of **Figure 1E**) or alternatively high level of correlations in the population of neuron selective for a non-memorized item during the delay period (as in the winner-take-all model of **Figure 1D**). To our knowledge experiments specifically analyzing how correlations change in the PFC as the delay-response task unfolds are still lacking.

At the same time, there are several lines of indirect evidence that lead us to believe that task dependent correlation modulation is indeed possible. First, it has been found that there is a modulation of spike coincidences during different phases of a motor task (Riehle et al., 1997). Riehle et al. (1997) found that at times during the delay when the animal was expecting to generate a response there were transients of synchronized spikes. Furthermore for successful trials there were more synchronized spikes during the delay period than for failure trails. This indeed suggests that spike coincidence is modulated in a functional way. The increase of excess synchrony at response (or expected response) times is compatible with the correlation based memory clearance discussed in this manuscript. Furthermore, it has been found that a change in representation during the delayresponse task leads to an increase of synchronization (Sakamoto et al., 2008). Pipa and Munk (2011) analyzed multi-unit activity during the delay period of a match-to-sample task and found that on correct vs. incorrect trials there is a modulation of spike synchronization and further, synchronous spike events are more prevalent at match presentation. This last point again suggests that increased correlations may be involved in erasing the memory trace.

In fact there is ample literature relating changes in oscillatory synchrony, coherence and frequency during WM tasks (Tallon-Baudry et al., 1998; Pesaran et al., 2002; Lee et al., 2005; Pipa et al., 2009). For example Pesaran et al. (2002) found that gamma-band spiking coherence is increased during the delay period in the lateral intraparietal cortex (LIP) in primates performing a delayed response task. Given that LIP is coupled to the PFC and is also involved in WM trace (Chafee and Goldman-Rakic, 1998), this is suggesting of increased input correlations to the PFC during the WM task. In the context of irregular poisson firing, oscillatory coherence is nothing other than correlations organized both in time (the frequency) and space. Oscillatory effects are beyond the scope of this paper and are a subject of the companion manuscript (Dipoppa and Gutkin, 2013).

The data reviewed above does show that there is a modulation of activity correlations during the WM task, yet it does not provide the mechanism. Here we propose that the mechanism is in the background input correlations generated by a common source. One hence might ask where such inputs may be coming from. As hinted above, one source could be coherent firing activity in the cortical regions coupled to the PFC and involved in WM processing (e.g., LIP). In addition, we propose that the source of shared background input generating spatial correlations can reside in the striatum, a subcortical area thought to be involved in WM. In fact the structure of the cortico-striatal loops as been longly seen as a disadvantage for the WM capacity if the memory store is located also in the striatum. Since the number of striatal neurons is much lower than the number of pyramidal neurons (Lange et al., 1976) and the loop is based on divergence (resp. convergence) in the striato-cortical (resp. cortico-striatal) direction, then the striatum could not have the same memory capacity of the cortex. It has been suggested that instead that divergent/convergent structure could be useful since the basal ganglia do not encode the individual information of WM but they control the gate of other region and decide when they can be updated (Frank et al., 2001). We also suggest that striatum plays a gating role since it could be the source of the common noise that creates the different regimes.

The correlation-induced gating is a robust effect to parameters variation. We propose the following explanation for this phenomenon: background correlations induce spike-times synchronization in the recurrent network, as was found similarly for independent neurons by Galán et al. (2006), and this leads to persistent activity erasing and block because of the refractory period of the neurons, as was found by Laing and Chow (2001) and Gutkin et al. (2001). Providing a proof of this assumption and a mathematical explanation of the correlation-induced gating will be the subject of future research.

## **ACKNOWLEDGMENTS**

The authors want to thank Ole Jensen, Romain Brette, Christian Machens, and Thomas Kreuz for constructive discussions. Mario Dipoppa was partially supported by MESR (France) (for Mario Dipoppa). Boris S. Gutkin was partially supported by CNRS, ANR-Blanc Grant Dopanic, CNRS Neuro IC grant, Neuropole Ile de France, Ecole de Neuroscience de Paris collaborative grant, and LABEX Institut des Etudes Cognitives, INSERM, and ENS.

## **REFERENCES**

*Comput. Neurosci.* 11, 121–134. doi: 10.1023/A:1012837415096

*Acad. Sci. U.S.A.* 103, 201–206. doi: 10.1073/pnas.0508072103

working-memory tasks. *J. Neurosci.* 26, 10141–10153. doi: 10.1523/ JNEUROSCI.2423-06.2006

*Hum. Percept. Perfor.* 27, 92–114. doi: 10.1037/0096-15<IP_ADDRESS>

Wilken, P., and Ma, W. J. (2004). A detection theory account of change detection. *J. Vis.* 4, 1120–1135. doi: 10.1167/4.12.11

**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: 25 September 2013; published online: October 2013. 21*

*Citation: Dipoppa M and Gutkin BS (2013) Correlations in background* *activity control persistent state stability and allow execution of working memory tasks. Front. Comput. Neurosci. 7:139. doi: 10.3389/fncom.2013.00139*

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

*Copyright © 2013 Dipoppa and Gutkin. 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.*

## Single-unit activities during epileptic discharges in the human hippocampal formation

## *Catalina Alvarado-Rojas 1,2, Katia Lehongre1,2, Juliana Bagdasaryan1,2, Anatol Bragin3, Richard Staba3, Jerome Engel Jr. 3, Vincent Navarro1,2,4 and Michel Le Van Quyen1,2\**

*1 Centre de Recherche de l'Institut du Cerveau et de la Moelle Epinière, INSERM UMRS 975 - CNRS UMR 7225, Hôpital de la Pitié-Salpêtrière, Paris, France*

*2 Université Pierre et Marie Curie, Paris, France*

*3 Department of Neurology, David Geffen School of Medicine at UCLA, Los Angeles, CA, USA*

*4 Epilepsy Unit, Groupe Hospitalier Pitié-Salpêtrière, Paris, France*

#### *Edited by:*

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

#### *Reviewed by:*

*Emili Balaguer-Ballester, Bournemouth University, UK Abdelmalik Moujahid, University of the Basque Country UPV/EHU, Spain*

#### *\*Correspondence:*

*Michel Le Van Quyen, Centre de Recherche de l'Institut du Cerveau et de la Moelle épinière, INSERM UMRS 975 - CNRS UMR 7225, Hôpital de la Pitié-Salpêtrière, 47 Bd de l'Hôpital, 75651 Paris, Cedex 13, France e-mail: quyen@ t-online.de*

Between seizures the brain of patients with epilepsy generates pathological patterns of synchronous activity, designated as interictal epileptiform discharges (ID). Using microelectrodes in the hippocampal formations of 8 patients with drug-resistant temporal lobe epilepsy, we studied ID by simultaneously analyzing action potentials from individual neurons and the local field potentials (LFPs) generated by the surrounding neuronal network. We found that ∼30% of the units increased their firing rate during ID and 40% showed a decrease during the post-ID period. Surprisingly, 30% of units showed either an increase or decrease in firing rates several hundred of milliseconds before the ID. In 4 patients, this pre-ID neuronal firing was correlated with field high-frequency oscillations at 40–120 Hz. Finally, we observed that only a very small subset of cells showed significant coincident firing before or during ID. Taken together, we suggested that, in contrast to traditional views, ID are generated by a sparse neuronal network and followed a heterogeneous synchronization process initiated over several hundreds of milliseconds before the paroxysmal discharges.

**Keywords: interictal epileptiform discharges, microelectrode recordings, multiunit activity, temporal lobe epilepsy, spike synchronization**

## **INTRODUCTION**

Synchronization of local and distributed neuronal assemblies is thought to underlie fundamental brain processes such as perception, learning, and cognition (Varela et al., 2001). In neurological diseases, neuronal synchrony can be altered and in epilepsy may play an important role in enhanced cellular excitability (Jasper and Penfield, 1954). Besides ictal events or seizures, interictal discharges (ID) are a typical signature of abnormal neuronal synchronization, seen spontaneously between seizures in scalp and intracranial EEG. They are used as a clinical indicator for the location of the epileptogenic zone, the region that generates seizures. Furthermore, it is believed that this region contains both, the seizure onset zone and the surrounding "irritative zone," which generates ID and limits with normal tissue (Talairach and Bancaud, 1966). These transient epileptic synchronization events are characterized by a large-amplitude, rapid component lasting 50–100 ms that is usually followed by a slow wave of 200–500 ms duration (de Curtis and Avanzini, 2001). In some cases, they are associated with an oscillation in the high frequency range greater than 40 Hz (Bragin et al., 1999; Jacobs et al., 2011; Le Van Quyen, 2012). Despite their fundamental importance in diagnosing and treating epilepsy, little is known about the neurophysiological mechanisms generating these events in the human brain. Experimental work on animals and human tissue propose the paroxysmal depolarization shift (PDS) as the cellular correlate of ID (Prince and Wong, 1981; Avoli and Williamson, 1996). This event is defined as a burst of action potentials on a large depolarization, followed by a longer hyperpolarization. However, *in vivo* human evidence is scarce, because of the limited opportunities to study the behavior of single neurons in human subjects. To overcome this difficulty, epilepsy patients suitable for surgical treatment are sometimes studied with intracranial depth electrodes in order to record EEG activity from deep cortical structures and accurately identify the regions originating seizures. Using depth electrodes specially adapted with microelectrodes (Fried et al., 1997; **Figure 1A**), ID can be studied by simultaneously recording action potentials from individual neurons and the local field potentials (LFP). Studies using microelectrode technology, have reported a variable and complex relation between ID and the activity of individual neurons, more heterogeneous than simple PDS (Babb et al., 1973; Wyler et al., 1982; Ulbert et al., 2004; Keller et al., 2010; Alarcon et al., 2012). In particular, a large diversity of neuronal response were found including an increase or decrease in their firing rates or even changes in firing that precede the defining interictal discharge itself. Most of these studies were performed on patients with neocortical epilepsy that exhibit a wide range of heterogeneity. In the present work, we recorded ID in the hippocampal formation of 8 patients with drug-resistant mesial temporal lobe epilepsy. Our objective is to describe firing patterns and neuronal synchronization of single-unit activities during spontaneous IDs.

**FIGURE 1 | (A)** Macro- and micro-electrodes superimposed on a magnetic resonance imaging scan. Nine microwires (40μm diameter) extend beyond the tip of each macro-electrode and record the hippocampal formation. **(B)** Interictal discharges (ID) recorded with microelectrode local field potentials from adjacent electrodes in the

hippocampus of a patient. **(C)** Example of wide-band recording of an ID event with the corresponding extracted single unit activities. **(D)** Raster plot and peri-event histogram (bin size, 10 ms) of the single unit activity shown above. Note the strong changes in the firing rate and instantaneous frequency (red) during the ID.

## **MATERIALS AND METHODS**

#### **DATABASE**

Subjects were 8 patients [two female, mean age ± standard deviation (SD) 36.3 ± 10.5 years] with pharmacologically intractable temporal lobe epilepsy who were implanted with 8–14 intracranial depth electrodes in order to localize epileptogenic regions for possible resection. The placement of the electrodes was determined exclusively by clinical criteria (Fried et al., 1999). Extending beyond the tip of each electrode were nine Pt-Ir microwires (40μm diameter) with inter-tip spacing of 500μm, eight active recording channels and one reference. Each microwire was sampled at 28 kHz (Cheetah recording system, Neuralynx Inc., Tucson, AZ). Spatial localizations were determined on the basis of postimplant computed tomography scans coregistered with preimplant 1.5T MRI scans. Our results are based on microelectrode recordings located in the anterior hippocampus (*n* = 40 channels in 5 patients) and entorhinal cortex (*n* = 24 channels in 3 patients). The recording states were quiet wakefulness and slow waves sleep (stages 1–4). All studies conformed to the guidelines of the Medical Institutional Review Board at University of California, Los Angeles.

#### **SPIKE SORTING**

In order to detect single-units, all channels were high-pass filtered at 300 Hz and were visually examined for the presence of unit activities. In those microwires with clear unit activities, we performed spike detection (>4:1 signal to noise ratio) to obtain multi-unit activities (MUA). Single-unit activities were extracted with spike sorting using KlustaKwik 1.7 program (Software: http://klustakwik.sourceforge.net/; Harris et al., 2000) which employs the 10 principal components of the spike shape and an unsupervised Conditional Expectation Maximization (CEM) clustering algorithm (Hazan et al., 2006). After automatic clustering, the clusters containing non-spike waveforms were visually deleted and then the units were further isolated using a manual cluster cutting method. Only units with clear boundaries and less than 0.5% of spike intervals within a 1 ms refractory period are included in the present analysis. Typically we isolated 1 or 2 distinct neurons from each microwire, but in several cases we observed up to 4 distinct neurons from a single microwire. The instantaneous spike frequency was measured by convolving the timing of each unit with a Gaussian function of standard deviation of 20 ms (*Ts* = 1 ms), set close to the modal interspike interval (Le Van Quyen et al., 2008, 2010). This operation leads to an analog trace of the instantaneous firing rate (Paulin, 1996).

#### **OSCILLATION ANALYSIS**

LFP are complementary to action potential information and have shown prominent oscillatory activity within the high-frequency frequency range from 40 to 300 Hz (Worrell et al., 2012). A wavelet time-frequency analysis was used to determine precisely the mean frequency, maximum amplitude and onset and offset of these LFP oscillations. The advantage of the wavelet analysis lies in the fact that the time resolution is variable with frequency, so that high frequencies have a sharper time resolution (Le Van Quyen and Bragin, 2007). The Complex Morlet wavelet was here applied that uses a wave-like scalable function that is well-localized in both time and frequency:

$$
\Psi\_{\mathfrak{r},f}(\mathfrak{u}) = \sqrt{f} \exp(j2\pi f(\mathfrak{u} - \mathfrak{r})) \exp\left(-\frac{(\mathfrak{u} - \mathfrak{r})^2}{2\sigma^2}\right).
$$

This wavelet represents the product of a sinusoidal wave at frequency *f*, with a Gaussian function centered at time τ, with a standard deviation σ proportional to the inverse of *f*. The wavelet coefficients of a signal *x(t)* as a function of time (τ) and frequency (*f*) are defined as: *W*(τ, *f*) = +∞ −∞ *x*(*u*)τ,*f*(*u*)*du*. It depends solely on σ, which sets the number of cycles of the wavelet: *nco* = 6*f* σ. The value *nco* determines the frequency resolution of the analysis by setting the width of the frequency interval for which phase are measured. Here, we chose *nco* = 5. For baseline correction, the average and SD of power were first computed at each frequency of the baseline period. Then, the average baseline power was subtracted from all time windows at each frequency, and the result scaled by 1/SD, yielding baseline-adjusted *Z* scores. Significant increases with respect to baseline activity showed up as positive *Z*-values and tabulated probability values indicate that, for absolute values of *Z* > 3.09, we have *P* < 0.001. The Kolmogorov– Smirnov test was applied to assess the distribution normality of the wavelet coefficients, using a 0.05 probability level.

#### **SPIKE SYNCHRONIZATION**

Different measures exist to detect and quantify synchronization between spike trains (Brown et al., 2004; Kreuz et al., 2007). In this study, we used two complementary techniques: (1) Crosscorrelation analysis was performed for cell pairs (Perkel et al., 1967; Amarasingham et al., 2012). To evaluate the significance of the correlation, we used a boot-strap method that accounts for the firing rate changes of the neurons (Hatsopoulos et al., 2003; Grün, 2009). Since the widths of the peaks in the original crosscorrelograms were typically in the range of 5–30 ms (Krüger and Mayer, 1990), the spikes were jittered by adding a random value from a normal distribution with a 50-ms SD and 0 mean to the spike times. For each cell pair, 1000 jittered spike trains were created, and the expected cross-correlogram (and 99% confidence interval) was estimated on 1 ms time bins. For any given cell pair where at least one bin in the [1.5 ms, 30 ms] interval exceeded the 99% confidence interval, the interaction was considered significant. (2) A method for identifying statistically conspicuous spike coincidences was implemented to detect the number of quasi-simultaneous appearances of spikes over small coincidence windows, here of 5 ms (Gütig et al., 2002; Quian Quiroga et al., 2002). Their occurrence was then studied in relation to surrogate data generated by dithering the individual, original spike times within a given time interval. Here, each spike in the original data set was randomly and independently jittered on a uniform interval of [−5, +5] ms to form a surrogate dataset. By repeating the procedure 1000 times, the 99.9% confidence interval for each bin (*p* = 0.001) was calculated.

## **RESULTS**

Microelectrode recordings were selected by an expert electroencephalographer to have very abundant and persistent ID in the hippocampus (5 patients) or entorhinal cortex (3 patients) during quiet wakefulness or slow-wave sleep (recording durations from 10 to 118 min; total recording time: 6 h). All ID were recorded in the epileptic zone and appeared as spatially synchronous patterns emerging at about the same time on the same bundle of microelectrodes (**Figure 1B**). A standard, threshold-based ID detector was performed to automatically detect, from the LFP, events showing a pointed peak with a large amplitude, large slope and duration of 20–100 ms, appearing at a frequency of 0.07 ± 0.30 Hz (range: 0.01–0.21 Hz). After expert visual confirmation, 862 ID were identified showing a large pattern of morphological characteristics typical for sharp waves, spikes and spike-wave discharges (Niedermeyer, 2005). Events were aligned by the sharpest peak of the discharge (**Figure 1C**). In order to analyze the patterns of neuronal activity around the discharge, we defined a baseline period (–600 to –300 ms), pre-ID period (–300 to –50 ms), the interictal discharge (–50 to 50 ms), the post-ID period (50– 400 ms). The activities of different neurons per microelectrode were identified with a spike sorting algorithm and a total of 75 single units were selected for analysis. To visualize the dischargerelated activity of single neurons, peri-stimulus raster plots and timing histograms were constructed for the period of 1 s before and after each event (**Figure 1D**).

During the ID period (–50 to 50 ms), we found that around 40% of the recorded units showed some change in firing, whereas 60% remain unchanged. About 32% increased their firing rate more than 2 times during ID relative to baseline epochs [**Figure 2B**; right-tail *t*-test: *T*(23) = 1.78; *p* = 0.04; an example of a cell can be seen in **Figure 1D**]. The firing rate of these cells showed a considerable degree of variability (range from 1.4 to 99 Hz) with a mean of 9.4 ± 19.7 Hz during ID (baseline: 2.7 ± 3.1 Hz). During the post-ictal period, 40% of units decreased firing by half [50–400 ms, mean firing rate: 1.8 ± 2.7 Hz and baseline: 7.0 ± 2.7 Hz, left-tail *t*-test: *T*(29) = −3.73; *p* = 4.1 · 10−4, **Figure 2C**]. In addition to this modulated single unit activity during ID, many units showed a significant change in firing preceding the interictal discharge. From 30% of single units that significantly changed during the pre-ID period, 12% increased [mean firing rate: 10.0 ± 13.5 Hz and baseline: 4.2 ± 5.8 Hz; *T*(8) = −3.45; *p* = 0.004] and 18% decreased [mean firing rate: 2.3 ± 7.0 Hz and baseline: 5.2 ± 11.6 Hz; *T*(13) = −1.64; *p* = 0.06] their firing rate (–300 to –50 ms, **Figure 2A**; examples are given in **Figure 2D**). On the corresponding channels, we were interested in the relationship between these pre-ID firing changes and

LFP (<300 Hz). Spectral power was performed by using Morlet wavelet analysis (20–300 Hz) and pre-ID changes in LFP were tested for significant increases/decreases from baseline of specific frequency bands (*p* < 0.001). In 4 subjects we observed that pre-ID neuronal firing pattern was correlated with an increase in high-frequency oscillations between 40 and 120 Hz (mean peak from baseline SD: *Z* = 6.1, range from 4.3 to 9.1). **Figure 3** shows average time-frequency representations around the ID for the two patients of **Figure 2D**. Main changes in spectral power can be seen in the LFP preceding the interictal discharge and correlate closely with the increase or decrease in neuronal firing.

Finally, we analyzed unit synchronization during ID between pairs of units simultaneously recorded in two different microelectrodes. Because of the inter-tip spacing of 500 μm, the units are assumed to reflect adjacent but different neuronal populations. Two complementary methods have been used to address the synchrony between spike trains. First, analysis of cross-correlograms between pairs of units was performed for each cell pairs that showed a sufficient number of spikes (>100) during ID. The significance of the correlation was obtained by jittering each pair of spike trains and by computing the 99% confidence interval. Of the 120 cross-correlograms constructed, only 5 cross-correlograms (about 4%) had a significant peak that occurred within ±25 ms around the origin, indicating that these neuronal pairs were discharging in a correlated way. **Figure 4** (top) illustrates examples of significant peaks in crosscorrelograms of two units. In addition to cross-correlation analysis, we also analyzed the overall level of synchronicity from the number of quasi simultaneous appearances of spikes. In order to not overestimate the number of random synchronous spikes due to the elevated firing rate, we used jitter techniques to infer millisecond-precise temporal synchrony (Hatsopoulos et al., 2003). Here, spikes of one of the pairs of neurons were time jittered by ±5 ms to generate jittered peri-stimulus raster plots of unit coincidence that could be used to assess the statistical significance of bin fluctuations in the non-jittered spike series. Because the jittered data sets preserve firing rates on timescales much broader than that of the jitter interval (in this case, 5 ms), the overall effect of the analysis is to identify those pairs that showed excessive co-firing at short latencies that cannot be accounted for by firing rates varying at timescales of tens of milliseconds. Despite the strong increase in about 30% of the recorded units during ID, only a very small subset of cells (18 of 120 analyzed pairs, about 15%) showed significant coincident firing before or during ID. For two patients, **Figure 4** (bottom) illustrates pairs of units that showed significant coincident firing (*p* = 0.001) during ID (A) and the pre-ID period (B).

## **DISCUSSION**

We found that a large subset of the recorded units showed significant changes in firing in or around ID in the hippocampal formation of patients with mesial temporal epilepsy. Around 30% of the unit increased their firing rate during ID while 40% showed a decrease during the post-ID period. This percentage of modulated neurons agrees with that described by Wyler et al. (1982), who found that 44% of recorded neurons showed primarily an increase in firing rate near the interictal discharge peak. Surprisingly, a subset of 30% of units showed significant firing rate variations several hundred of milliseconds before the ID. In a few patients, we observed that this neuronal firing pattern was related with elevated LFP oscillations at 40–120 Hz. Finally, based on two statistical methods that identify spike synchronization, we found that only a very small subset of cells showed significant coincident firing before or during ID.

Our observations of neuronal firing during the interictal discharge are consistent with the paroxysmal depolarizing shift (PDS) mechanism—a large depolarization phase followed by a long hyperpolarization—that have been studied in animal models of epilepsy (Matsumoto and Marsan, 1964; Prince, 1968). The first part of the depolarization phase is believed to be generated by intrinsic membrane conductance (de Curtis et al., 1999), and the later from feedback recurrent synaptic excitation mediated by AMPA and NMDA receptor subtypes, and glutamate receptorcoupled calcium conductances (Traub et al., 1993). Thus, PDS has been shown to be the result of giant excitatory postsynaptic potentials. The PDS is usually followed by a hyperpolarization, which represents GABA-mediated recurrent inhibition, as well as Ca2+-dependent K+ currents. Interestingly, consistent with *in vitro* studies on hippocampal slices from human patients with temporal lobe epilepsy (Cohen et al., 2002; Wittner et al., 2009), the presence of a similar suppression of unit activities in our *in vivo* data suggests that IDs can occur in cortical regions maintaining substantial inhibitory function.

However, in contrast to simple models of PDS and in line with other observations in human epileptic neocortex (Keller et al., 2010), we found that ID, rather than requiring a large synchronization of neurons, can occur with relatively sparse single neuron participation (estimated at about 30% of the cells). Furthermore, a small subset of the units significantly increased or decreased their firing well before ID. Concomitant with changes in firing rate for certain neurons, at least in some patients, high-frequency oscillations at 40–120 Hz can be seen in the LFP preceding the ID and correlate closely with the changes in neuronal firing. Because interneurons are involved in the generation of high frequency oscillations through mechanisms of post-inhibition resetting of neuronal firing (Cobb et al., 1995; Ylinen et al., 1995; Le Van Quyen et al., 2008; Le Van Quyen, 2012), it is here tempting to speculate that GABA-mediated events may contribute to enhance synchronization of local epileptic networks before ID. Interestingly, emerging evidence indicates that GABA promotes epileptiform synchronization (Pavlov et al., 2013). For instance, GABA receptor-mediated inhibition can facilitate thalamocortical processes leading to the occurrence of generalized spike and wave discharges that occur during absence seizures (Danober et al., 1998). Following a similar mechanism, ID may be caused by a rebound synchronization of cells that may start firing

synchronizations (red circles) were defined as coincidences between the two units (green and blue points) occurring over a 5-ms interval. Note the significant increase in coincidences during ID **(A)** and the pre-ID period **(B)**, over the statistical threshold defined by a random jitter of the original data.

synchronously shortly after inhibition ceases and permit the fast component of the ID. Moreover, intense synaptic activation of GABAA receptors in the hippocampus can lead to a shift in GABAergic neurotransmission from inhibitory to excitatory, contributing to epileptic discharges (Kohling et al., 2000; Cohen et al., 2002). Interestingly, pre-event changes have also been seen in advance of seizures in an animal model of temporal lobe epilepsy (Bower and Buckmaster, 2008) and around seizure onset in human epilepsy (Babb and Crandall, 1976; Truccolo et al., 2011), suggesting a possible similar mechanism before seizures.

Taken together, our data suggest that ID in patients with temporal lobe epilepsy is not a simple paroxysm of hypersynchronous

## **REFERENCES**

properties of human neocortical neurons maintained *in vitro*. *Prog. Neurobiol*. 48, 519–554. doi: 10.1016/0301-0082(95)00050-X

excitatory activity, but rather represents a heterogeneous synchronization process possibly initiated by GABAergic responses in small subsets of cells and emerging over hundreds of milliseconds before the paroxysmal discharges.

## **ACKNOWLEDGMENTS**

Catalina Alvarado-Rojas was supported by the Administrative Department for Science, Technology and Innovation (COLCIENCIAS), Colombia. Vincent Navarro was supported by a Contrat Inferface INSERM. This work was supported by funding from the program "Investissements d'avenir" ANR-10-IAIHU-06 and from the ICM and OCIRP.

by time-frequency analysis. *Clin. Neurophysiol*. 122, 32–42. doi: 10.1016/j.clinph.2010.05.033

*Applications, and Related Fields*, eds E. Niedermeyer and F. Lopes da Silva (Philadelphia, PA: Lippincott Williams and Wilkins).

brainweb: phase synchronization and large-scale integration. *Nat. Rev. Neurosci*. 2, 229–239. doi: 10.1038/35067550

**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: 27 September 2013; published online: 18 October 2013.*

*Citation: Alvarado-Rojas C, Lehongre K, Bagdasaryan J, Bragin A, Staba R, Engel Jr J, Navarro V and Le Van Quyen M (2013) Single-unit activities during epileptic discharges in the human hippocampal formation. Front. Comput. Neurosci. 7:140. doi: 10.3389/ fncom.2013.00140*

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

*Copyright © 2013 Alvarado-Rojas, Lehongre, Bagdasaryan, Bragin, Staba, Engel, Navarro and Le Van Quyen. 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.*

## A unified view on weakly correlated recurrent networks

#### *Dmytro Grytskyy1 \*, Tom Tetzlaff 1, Markus Diesmann1,2 and Moritz Helias 1*

*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 Medical Faculty, RWTH Aachen University, Germany*

#### *Edited by:*

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

#### *Reviewed by:*

*Brent Doiron, University of Pittsburgh, USA Shigeru Shinomoto, Kyoto University, Japan*

#### *\*Correspondence:*

*Dmytro Grytskyy, Jülich Research Centre and JARA, Institute of Neuroscience and Medicine (INM-6) and Institute for Advanced Simulation (IAS-6), Building 15.22, 52425 Jülich, Germany e-mail<EMAIL_ADDRESS>
The diversity of neuron models used in contemporary theoretical neuroscience to investigate specific properties of covariances in the spiking activity raises the question how these models relate to each other. In particular it is hard to distinguish between generic properties of covariances and peculiarities due to the abstracted model. Here we present a unified view on pairwise covariances in recurrent networks in the irregular regime. We consider the binary neuron model, the leaky integrate-and-fire (LIF) model, and the Hawkes process. We show that linear approximation maps each of these models to either of two classes of linear rate models (LRM), including the Ornstein–Uhlenbeck process (OUP) as a special case. The distinction between both classes is the location of additive noise in the rate dynamics, which is located on the output side for spiking models and on the input side for the binary model. Both classes allow closed form solutions for the covariance. For output noise it separates into an echo term and a term due to correlated input. The unified framework enables us to transfer results between models. For example, we generalize the binary model and the Hawkes process to the situation with synaptic conduction delays and simplify derivations for established results. Our approach is applicable to general network structures and suitable for the calculation of population averages. The derived averages are exact for fixed out-degree network architectures and approximate for fixed in-degree. We demonstrate how taking into account fluctuations in the linearization procedure increases the accuracy of the effective theory and we explain the class dependent differences between covariances in the time and the frequency domain. Finally we show that the oscillatory instability emerging in networks of LIF models with delayed inhibitory feedback is a model-invariant feature: the same structure of poles in the complex frequency plane determines the population power spectra.

**Keywords: correlations, linear response, Hawkes process, leaky integrate-and-fire model, binary neuron, linear rate model, Ornstein–Uhlenbeck process**

## **1. INTRODUCTION**

The meaning of correlated neural activity for the processing and representation of information in cortical networks is still not understood, but evidence for a pivotal role of correlations increases (recently reviewed in Cohen and Kohn, 2011). Different studies have shown that correlations can either decrease (Zohary et al., 1994) or increase (Sompolinsky et al., 2001) the signal to noise ratio of population signals, depending on the readout mechanism. The architecture of cortical networks is dominated by convergent and divergent connections among the neurons (Braitenberg and Schüz, 1991) causing correlated neuronal activity by common input from shared afferent neurons in addition to direct connections between pairs of neurons and common external signals. It has been shown that correlated activity can faithfully propagate through convergent-divergent feed forward structures, such as synfire chains (Abeles, 1991; Diesmann et al., 1999), a potential mechanism to convey signals in the brain. Correlated firing was also proposed as a key to the solution of the binding problem (von der Malsburg, 1981; Bienenstock, 1995; Singer, 1999), an idea that has been discussed controversially (Shadlen and Movshon, 1999). Independent of a direct functional role of correlations in cortical processing, the covariance function between the spiking activity of a pair of neurons contains the information about time intervals between spikes. Changes of synaptic coupling, mediated by spike-timing dependent synaptic plasticity (STDP, Markram et al., 1997; Bi and Poo, 1999), are hence sensitive to correlations. Understanding covariances in spiking networks is thus a prerequisite to investigate the evolution of synapses in plastic networks (Burkitt et al., 2007; Gilson et al., 2009, 2010).

On the other side, there is ubiquitous experimental evidence of correlated spike events in biological neural networks, going back to early reports on multi-unit recordings in cat auditory cortex (Perkel et al., 1967; Gerstein and Perkel, 1969), the observation of closely time-locked spikes appearing at behaviorally relevant points in time (Kilavik et al., 2009; Ito et al., 2011) and collective oscillations in cortex [recently reviewed in Buzsáki and Wang (2012)].

The existing theories explaining correlated activity use a multitude of different neuron models. Hawkes (1971) developed the theory of covariances for linear spiking Poisson neurons (Hawkes processes). Ginzburg and Sompolinsky (1994) presented the approach of linearization to treat fluctuations around the point of stationary activity and to obtain the covariances for networks of non-linear binary neurons. The formal concept of linearization allowed Brunel and Hakim (1999) and Brunel (2000) to explain fast collective gamma oscillations in networks

of spiking leaky integrate-and-fire (LIF) neurons. Correlations in feed-forward networks of LIF models are studied in Moreno-Bote and Parga (2006), exact analytical solutions for such network architectures are given in Rosenbaum and Josic (2011) for the case of stochastic random walk models, and threshold crossing neuron models are considered in Tchumatchenko et al. (2010) and Burak et al. (2009). Covariances in structured networks are investigated for Hawkes processes (Pernice et al., 2011), and in linear approximation for LIF (Pernice et al., 2012) and exponential integrate-and-fire neurons (Trousdale et al., 2012). The latter three works employ an expansion of the propagator (time evolution operator) in terms of the order of interaction. Finally Buice et al. (2009) investigate higher order cumulants of the joint activity in networks of binary model neurons.

Analytical insight into a neuroscientific phenomenon based on correlated neuronal activity often requires a careful choice of the neuron model to arrive at a solvable problem. Hence a diversity of models has been proposed and is in use. This raises the question which features of covariances are generic properties of recurrent networks and which are specific to a certain model. Only if this question can be answered one can be sure that a particular result is not an artifact of oversimplified neuronal dynamics. Currently it is unclear how different neuron models relate to each other and whether and how results obtained with one model carry over to another. In this work we present a unified theoretical view on pairwise correlations in recurrent networks in the asynchronous and collective-oscillatory regime, approximating the response of different models to linear order. The joint treatment allows us to answer the question of genericness and moreover naturally leads to a classification of the considered models into only two categories, as illustrated in **Figure 1**. The classification in addition enables us to extend existing theoretical results to biologically relevant parameters, such as synaptic delays and the presence of inhibition, and to derive explicit expressions for the timedependent covariance functions, in quantitative agreement with direct simulations, which can serve as a starting point for further work.

The remainder of this article is organized as follows. In the first part of our results in "Covariance structure of noisy rate models" we investigate the activity and the structure of covariance functions for two versions of linear rate models (LRM); one with input the other with output noise. If the activity relaxes exponentially after application of a short perturbation, both models coincide with the OUP. We mainly consider the latter case, although most results hold for arbitrary kernel functions. We extend the analytical solutions for the covariances in networks of OUP (Risken, 1996) to the neuroscientifically important case of synaptic conduction delays. Solutions are derived first for general forms of connectivity in "Solution of the convolution equation with input noise" for input noise and in "Solution of convolution equation with output noise" for output noise. After analyzing the spectral properties of the dynamics in the frequency domain in "Spectrum of the dynamics," identifying poles of the propagators

**linear rate models (LRM).** The arrows indicate analytical methods which enable a mapping from the original spiking (LIF model, Hawkes model) or binary neuron dynamics to the analytically more tractable linear rate models. Depending on the original dynamics (spiking or binary) the resulting LRM contains an additive noise component *x* either on the output side (left) or on the input side (right).

and their relation to collective oscillations in neuronal networks, we show in "Population-averaged covariances"' how to obtain pairwise averaged covariances in homogeneous Erdös-Rényi random networks. We explain in detail the use of the residue theorem to perform the Fourier back-transformation of covariance functions to the time domain in "Fourier back transformation" for general connectivity and in "Explicit expression for the population averaged cross covariance in the time domain" for averaged covariance functions in random networks, which allows us to obtain explicit results and to discuss class dependent features of covariance functions.

In the second part of our results in "Binary neurons," "Hawkes processes," and "Leaky integrate-and-fire neurons" we consider the mapping of different neuronal dynamics on either of the two flavors of the linear rate models discussed in the first part. The mapping procedure is qualitatively the same for all dynamics as illustrated in **Figure 1**: Starting from the dynamic equations of the respective model, we first determine the working point described in terms of the mean activity in the network. For unstructured homogeneous random networks this amounts to a mean-field description in terms of the population averaged activity (i.e., firing rate in spiking models). In the next step, a linearization of the dynamical equations is performed around this working point. We explain how fluctuations can be considered in the linearization procedure to improve its accuracy and we show how the effective linear dynamics maps to the LRM. We illustrate the results throughout by a quantitative comparison of the analytical results to direct numerical simulations of the original non-linear dynamics. The appendices "Implementation of noisy rate models," "Implementation of binary neurons in a spiking simulator code," and "Implementation of Hawkes neurons in a spiking simulator code." describe the model implementations and are modules of our long-term collaborative project to provide the technology for neural systems simulations (Gewaltig and Diesmann, 2007).

#### **2. COVARIANCE STRUCTURE OF NOISY RATE MODELS**

#### **2.1. DEFINITION OF MODELS**

Let us consider a network of linear model neurons, each characterized by a continuous fluctuating rate *r* and connections from neuron *j* to neuron *i* given by the element *wij* of the connectivity matrix **w**. We assume that the response of neuron *i* to input can be described by a linear kernel *h* so that the activity in the network fulfills

$$\mathbf{r}(t) = h(\circ) \* [\mathbf{w}\mathbf{r}(\circ - d) + b\mathbf{x}(\circ)](t),\tag{1}$$

where *f*(◦ − *d*) denotes the function *f* shifted by the delay *d*, **x** is an uncorrelated noise with

$$
\langle \mathbf{x}\_i(t) \rangle = 0, \qquad \langle \mathbf{x}\_i(s) \mathbf{x}\_j(t) \rangle = \delta\_{ij} \delta(s - t) \rho^2, \tag{2}
$$

e.g., a Gaussian white noise and (*f* ∗ *g*)(*t*) = *t* −∞ *f*(*t* − *t* ) *g*(*t* ) *dt* is the convolution. With the particular choice *b* = **w**δ(◦ − *d*)∗ we obtain

$$\mathbf{r}(t) = [h(\circ) \ast \mathbf{w}(\mathbf{r}(\circ - d) + \mathbf{x}(\circ - d))](t). \tag{3}$$

We call the dynamics (3) the linear noisy rate model (LRM) with noise applied to output, as the sum *r* + *x* appears on the right hand side. Alternatively, choosing *b* = **1** we define the model with input noise as

$$\mathbf{r}(t) = h(\circ) \* [\mathbf{w}\mathbf{r}(\circ - d) + \mathbf{x}(\circ)](t). \tag{4}$$

Hence, Equations (3) and (4) are special cases of (1). In the following we consider the particular case of an exponential kernel

$$h(\mathbf{s}) = \frac{1}{\mathfrak{r}} \theta(\mathbf{s}) e^{-\mathbf{s}/\mathfrak{r}},\tag{5}$$

where θ denotes the Heaviside function, θ(*t*) = 1 for *t* > 0, 0 else. Applying to (1) the operator *O* = τ *d ds* + 1 which has *h* as a Green's function (i.e., *Oh* = δ) we get

$$\mathbf{r}\frac{d}{dt}\mathbf{r}(t) + \mathbf{r}(t) = \mathbf{w}\mathbf{r}(t-d) + \mathbf{b}\mathbf{x}(t),\tag{6}$$

which is the equation describing a set of delay coupled Ornstein-Uhlenbeck-processes (OUP) with input or output noise for *b* = **1** or *b* = **w**δ(◦ − *d*)∗, respectively. We use this representation in "Binary neurons" to show the correspondence to networks of binary neurons.

#### **2.2. SOLUTION OF THE CONVOLUTION EQUATION WITH INPUT NOISE**

The solution for the system with input noise obtained from the definition (4) after Fourier transformation is

$$\mathbf{R} = H\_d \mathbf{w} \mathbf{R} + H \mathbf{X},\tag{7}$$

where the delay is consumed by the kernel function *hd*(*s*) = 1 τ θ(*s* − *d*)*e*−(*s*−*d*)/τ. We use capital letters throughout the text to denote objects in the Fourier domain and lower case letters for objects in the time domain. Solved for **R** = (1 − *Hd***w**)−1*H***X** the covariance function of **r** in the Fourier domain is found with the Wiener–Khinchin theorem (Gardiner, 2004) as **R**(ω)**R***T*(−ω), also called the cross spectrum

$$\mathbf{C}(\boldsymbol{\omega}) = \langle \mathbf{R}(\boldsymbol{\omega}) \mathbf{R}^{\top}(-\boldsymbol{\omega}) \rangle \tag{8}$$

$$= (1 - H\_d(\boldsymbol{\omega}) \mathbf{w})^{-1} H(\boldsymbol{\omega}) \langle \mathbf{X}(\boldsymbol{\omega}) \mathbf{X}^{T}(-\boldsymbol{\omega}) \rangle$$

$$H(-\boldsymbol{\omega}) (1 - H\_d(-\boldsymbol{\omega}) \mathbf{w}^{T})^{-1}$$

$$= (H\_d(\boldsymbol{\omega})^{-1} - \mathbf{w})^{-1} \mathbf{D} (H\_d(-\boldsymbol{\omega})^{-1} - \mathbf{w}^{T})^{-1},$$

where we introduced the matrix **D** = **X**(ω)**X***T*(−ω). From the second to the third line we used the fact that the non-delayed kernels *H*(ω) can be replaced by delayed kernels *Hd*(ω) and that the corresponding phase factors *ei*ω*d* and *e*−*i*ω*d* cancel each other. If **x** is a vector of pairwise uncorrelated noise, **D** is a diagonal matrix and needs to be chosen accordingly in order for the cross spectrum (8) to coincide (neglecting non-linear effects) with the cross spectrum of a network of binary neurons, as described in "Equivalence of binary neurons and Ornstein–Uhlenbeck processes".

#### **2.3. SOLUTION OF CONVOLUTION EQUATION WITH OUTPUT NOISE**

For the system with output noise we consider the quantity *yi* = *ri* + *xi* as the dynamic variable representing the activity of neuron *i* and aim to determine pairwise correlations. It is easy to get from (3) after Fourier transformation

$$\mathbf{R} = H\_d \mathbf{w}(\mathbf{R} + \mathbf{X}),\tag{9}$$

which can be solved for **R** = (1 − *Hd***w**)−1*Hd***wX** in order to determine the Fourier transform of **Y** as

$$\mathbf{Y} = \mathbf{R} + \mathbf{X} = (1 - H\_d \mathbf{w})^{-1} \mathbf{X}.\tag{10}$$

The cross spectrum hence follows as

$$\mathbf{C}(\boldsymbol{\omega}) = \langle \mathbf{Y}(\boldsymbol{\omega}) \mathbf{Y}^T(-\boldsymbol{\omega}) \rangle \tag{11}$$

$$= (1 - H\_d(\boldsymbol{\omega}) \mathbf{w})^{-1} \langle \mathbf{X}(\boldsymbol{\omega}) \mathbf{X}^T(-\boldsymbol{\omega}) \rangle (1 - H\_d(-\boldsymbol{\omega}) \mathbf{w}^T)^{-1}$$

$$= (1 - H\_d(\boldsymbol{\omega}) \mathbf{w})^{-1} \mathbf{D} (1 - H\_d(-\boldsymbol{\omega}) \mathbf{w}^T)^{-1},$$

with **D** = **X**(ω)**X***T*(−ω). **D** is a diagonal matrix with the *i*-th diagonal entry ρ2 *i* . For the correspondence to spiking models **D** must be chosen appropriately, as discussed in "Hawkes processes" and "Leaky integrate-and-fire neurons" for Hawkes processes and LIF neurons, respectively.

#### **2.4. SPECTRUM OF THE DYNAMICS**

For both linear rate dynamics, with output and with input noise, the cross spectrum **C**(ω) has poles at certain frequencies ω in the complex plane. These poles are defined by the zeros of det(*Hd*(ω)−1 − **w**) and the corresponding term with the opposite sign of ω. The zeros of det(*Hd*(ω)−1 − **w**) are solutions of the equation

$$\left(H\_d(\omega)\right)^{-1} = (1 + i\omega\mathfrak{r})e^{i\omega d} = L\_j$$

where *Lj* is the *j*-th eigenvalue of **w**. The same set of poles arises from (1) when solving for **R**. For *d* > 0 and the exponential kernel (5), the poles can be expressed as

$$z\_k(L\_j) = \frac{i}{\mathfrak{r}} - \frac{i}{d} \mathcal{W}\_k \left( L\_j \frac{d}{\mathfrak{r}} e^{\frac{d}{\mathfrak{r}}} \right), \tag{12}$$

where *Wk* is the *k*-th of the infinitely many branches of the Lambert-W function (Corless et al., 1996). For vanishing synaptic delay *d* = 0 there is obviously only one solution for every *Lj* given by *z* = −*i* τ (*Lj* − 1).

Given the same parameters *d*, **w**, τ, the pole structures of the cross spectra of both systems (8) and (11) are identical, since the former can be obtained from the latter by multiplication with (*Hd*(ω)*Hd*(−ω))−1 = (*H*(ω)*H*(−ω))−1, which has no poles. The only exception causing a different pole structure for the two models is the existence of an eigenvalue *Lj* = 0 of the connectivity matrix **w**, corresponding to a pole *z*(0) = *i* τ . However, this pole corresponds to an exponential decay of the covariance for input noise in the time domain and hence does not contribute to oscillations. For output noise, the multiplication with the term (*H*(ω)*H*(−ω))−1, vanishing at ω = *i* τ , cancels this pole in the covariance. Consequently both dynamics exhibit similar oscillations. A typical spectrum of poles for a negative eigenvalue *Lj* < 0 is shown in **Figures 2B,D**.

#### **2.5. POPULATION-AVERAGED COVARIANCES**

Often it is desirable to consider not the whole covariance matrix but averages over subpopulations of pairs of neurons. For instance the average over the whole network would result in a single scalar value. Separately averaging pairs, distinguishing excitatory and inhibitory neuron populations, yields a 2 by 2 matrix of covariances. For these simpler objects closed form solutions can be obtained, which already preserve some useful information and show important features of the network. Averaged covariances are also useful for comparison with simulations and experimental results.

In the following we consider a recurrent random network of *Ne* = *N* excitatory and *Ni* = γ*N* inhibitory neurons with synaptic weight *w* for excitatory and −*gw* for inhibitory synapses. The probability *p* determines the existence of a connection between

**FIGURE 2 | Pole structure determines dynamics.** Autocovariance of the population activity **(A,C)** measured in ρ2/τ and its Fourier transform called power spectrum **(B,D)** of the rate models with output noise (dots) (3) and input noise (diagonal crosses) (4) for delays *d* = 3 ms **(A,B)**, and *d* = 1 ms **(C,D)**. Black symbols show averages over the excitatory population activity and gray symbols over the inhibitory activity obtained by direct simulation. Light gray curves show theoretical predictions for the spectrum (20) and the covariance (21) for output noise and the spectrum (17) and the covariance (18) for input noise. Black crosses (12) in **(B,D)** denote the locations of the poles of the cross spectra - with the real parts corresponding to the damping (vertical axis), and the imaginary parts to oscillation frequencies (horizontal axis). The detailed parameters for this and following figures are given in "Parameters of simulations".

two randomly chosen neurons. We study the dynamics averaged over the two subpopulations by introducing the quantities *ra* = 1 *Na j* ∈ *a rj* and noise terms *xa* = 1 *Na j* ∈ *a xj* for *a* ∈ {*E*, *I*}; indices *I* and *E* stand for inhibitory and excitatory neurons and corresponding quantities. Calculating the average local input *N*−1 *a j* ∈ *a wjkrk* to a neuron of type *a*, we obtain

$$N\_a^{-1} \sum\_{j \in a} \sum\_k \boldsymbol{w}\_{jk} \boldsymbol{r}\_k = N\_a^{-1} \left( \sum\_{j \in a} \sum\_{k \in \mathcal{E}} \boldsymbol{w}\_{jk} \boldsymbol{r}\_k + \sum\_{j \in a} \sum\_{k \in \mathcal{I}} \boldsymbol{w}\_{jk} \boldsymbol{r}\_k \right) (13)$$

$$= N\_a^{-1} \left( p N\_a \boldsymbol{w} \sum\_{k \in \mathcal{E}} \boldsymbol{r}\_k - p N\_a \boldsymbol{g} \boldsymbol{w} \sum\_{k \in \mathcal{I}} \boldsymbol{r}\_k \right)$$

$$= p \boldsymbol{w} N(\boldsymbol{r}\_\mathcal{E} - \boldsymbol{\chi} \boldsymbol{g} \boldsymbol{r} \boldsymbol{\chi}),$$

where, from the second to the third line we used the fact that in expectation a given neuron *k* has *pNa* targets in the population *a*. The reduction to the averaged system in (13) is exact if in every column *k* in **w***jk* there are exactly *K* non-zero elements for *j* ∈ *E* and γ*K* for *j* ∈ *I*, which is the case for networks with fixed out-degree (number of outgoing connections of a neuron to the neurons of a particular type is kept constant), as noted earlier (Tetzlaff et al., 2012). For fixed in-degree (number of connections to a neuron coming in from the neurons of a particular type is kept constant) the substitution of *rj* ∈ *a* by *ra* is an additional approximation, which could be considered as an average over possible realizations of the random connectivity. In both cases the effective population-averaged connectivity matrix **M** turns out to be

$$\mathbf{M} = K\nu \begin{pmatrix} 1 \ -\gamma \mathbf{g} \\ 1 \ -\gamma \mathbf{g} \end{pmatrix},\tag{14}$$

with *K* = *pN*. So the averaged activities fulfill the same Equations (3) and (4) with the non-averaged quantities **r**, **x**, and **w** replaced by their averaged counterparts **r**¯ = (*rE* ,*rI*)*T*, **x**¯ = (*xE* , *xI*)*T*, and **M**. The population averaged activities *ra* are directly related to the block-wise averaged covariance matrix **c**¯ = *cEE cEI cIE cII* , with *cab* = *N*−1 *a N*−1 *b i* ∈ *a j* ∈ *b cij*. With

$$\bar{D}\_{ab} = N\_a^{-1} N\_b^{-1} \left\langle \sum\_{i \in a} \mathbf{x}\_i \sum\_{j \in b} \mathbf{x}\_j \right\rangle \tag{15}$$

$$= N\_a^{-1} N\_b^{-1} \sum\_{i \in a} \sum\_{j \in b} D\_{ij}$$

$$= \delta\_{ab} N\_a / N\_a^2 \rho^2 = \delta\_{ab} N\_a^{-1} \rho^2$$

we replace **D** by **D**¯ = ρ2 *N*−1 0 0 (γ*N*)−1 and **c** by **c**¯ so that the same Equations (11) and (8) and their general solutions also hold for the block-wise averaged covariance matrices.

The covariance matrices separately averaged over pairs of excitatory, inhibitory or mixed pairs are shown in **Figure 2** for both linear rate dynamics (3) and (4). (Parameters for all simulations presented in this article are collected in "Parameters of simulations," the implementation of LRM is described in "Implementation of noisy rate models"). The poles of both models shown in **Figure 2B** are given by (12) and coincide with the peaks in the cross spectra (8) and (11) for output and input noise, respectively. The results of direct simulation and the theoretical prediction are shown for two different delays, with the longer delay leading to stronger oscillations.

**Figure 3C** shows the distribution of eigenvalues in the complex plane for two random connectivity matrices with different synaptic amplitudes *w*. The model exhibits a bifurcation, if at least one eigenvalue assumes a zero real part. For fixed out-degree the averaging procedure (13) is exact, reflected by the precise agreement of theory and simulation in **Figure 3D**. For fixed in-degree, the averaging procedure (13) is an approximation, which is good only for parameters far from the bifurcation. Even in this regime still small deviations of the theory from the simulation results are visible in **Figure 3B**. On the stable side close to a bifurcation, the appearance of long living modes causes large fluctuations. These weakly damped modes appearing in one particular realization of the connectivity matrix are not represented after the replacement of the full matrix **w** by the average **M** over matrix realizations. The eigenvalue spectrum of the connectivity matrix provides an alternative way to understand the deviations. By the averaging the set of *N* eigenvalues of the connectivity matrix is replaced withby the two eigenvalues of the reduced matrix **M**, one of which is zero due to identical rows of **M**. The eigenvalue spectrum of the full matrix

**FIGURE 3 | Limits of the theory for fixed in-degree and fixed out-degree.** Autocovariance **(A)** and covariance **(B)** in random networks with fixed in-degree (dots) and fixed out-degree (crosses). Simulation results for *cEE* , *cEI* , and *cII* are shown in dark gray, black and light gray, respectively for synaptic weight *w* = 0.011 far from bifurcation. For larger synaptic weight *w* = 0.018 close to bifurcation (see text at the end of "Population-averaged covariances"), *cEE* is also shown in **(D)** for fixed in-degree (dark gray dots) and for fixed out-degree (black dots). Corresponding theoretical predictions for the autocovariance (34) **(A)** and the covariance (18) **(B,D)** are plotted as light gray curves throughout. The set of eigenvalues is shown as black dots in panel **(C)** for the smaller weight. The gray circle denotes the spectral radius *w Np*(1 − *p*)(1 + γ*g*2) (Rajan and Abbott, 2006; Kriener et al., 2008) confining the set of eigenvalues for the larger weight. The small filled gray circle and the triangle show the effective eigenvalues *L* of the averaged systems for small and large weight, respectively.

is illustrated in **Figure 3C**. Even if the eigenvalue(s) *L***M** of **M** are far in the stable region (corresponding to (*z*(*L***M**)) > 0) some eigenvalues *L***w** of the full connectivity matrix in the vicinity of the bifurcation region may still have an imaginary part becoming negative and the system can feel their influence, shown in **Figure 3D**.

#### **2.6. FOURIER BACK TRANSFORMATION**

Although the cross spectral matrices (8) and (11) for both dynamics look similar in the Fourier domain, the procedures for back transformation differ in detail. In both cases, the Fourier integral along the real ω-axis can be extended to a closed integration contour by a semi-circle with infinite radius centered at 0 in the appropriately chosen half-plane. The half-plane needs to be selected such that the contribution of the integration along the semi-circle vanishes. By employing the residue theorem (Bronstein et al., 1999) the integral can be replaced by a sum over residua of the poles encircled by the contour. For a general covariance matrix we only need to calculate **c**(*t*) for *t* ≥ 0, as for *t* < 0 the solution can be found by symmetry **c**(*t*) = **c***T*(−*t*).

For input noise it is possible to close the contour in the upper half-plane where the integrand **C**(ω)*ei*ω*t* vanishes for |ω| → ∞ for all *t* > 0, as |*Cij*(ω)| decays as |ω| −2. This can be seen from (8), because the highest order of *H*−1 *d* ∝ ω appearing in det(*H*−1 *d* − **w**) is equal to the dimensionality *N* of **w** (*N* = 2 for **M**), and in det(adjugate matrix *ij* of *H*−1 *d* − **w**) it is *N* − 1 (*i* = *j*) or *N* − 2 (*i* = *j*). So |(*H*−1 *d* − **w**)−1| is proportional to |ω| −1|*e*−*i*ω*d*| and |**C**(ω)|∝|ω| −2 for large |ω|.

For the case of output noise (11) **C**(ω) can be obtained from the **C**(ω) for input noise (8) multiplied with (*Hd*(ω)*Hd*(−ω))−1 ∼ |ω| 2 for large |ω|. The multiplication with this factor changes the asymptotic behavior of the integrand, which therefore contains terms converging to a constant value and terms decaying like |ω| −1 for |ω| → ∞. These terms result in non-vanishing integrals over the semicircle in the upper half-plane and have to be considered separately. To this end we rewrite (11) as

$$\begin{split} \mathbf{C}(\boldsymbol{\omega}) &= ((1 - H\_d(\boldsymbol{\omega})\mathbf{w})^{-1} H\_d(\boldsymbol{\omega}) \mathbf{w} + 1) \\ &\quad \mathbf{D}(\mathbf{w}^T H\_d(-\boldsymbol{\omega})(1 - H\_d(-\boldsymbol{\omega})\mathbf{w}^T)^{-1} + 1) \\ &= (1 - H\_d(\boldsymbol{\omega})\mathbf{w})^{-1} H\_d(\boldsymbol{\omega}) \mathbf{w} \mathbf{D} \mathbf{w}^T H\_d(-\boldsymbol{\omega})(1 - H\_d(-\boldsymbol{\omega})\mathbf{w}^T)^{-1} \\ &\quad + (1 - H\_d(\boldsymbol{\omega})\mathbf{w})^{-1} H\_d(\boldsymbol{\omega}) \mathbf{w} \mathbf{D} \\ &\quad + \mathbf{D} \mathbf{w}^T H\_d(-\boldsymbol{\omega})(1 - H\_d(-\boldsymbol{\omega})\mathbf{w}^T)^{-1} \\ &\quad + \mathbf{D}, \end{split}$$

and find the constant term **D** which turns into a δ-function in the time domain. The first term in the second line of (16) decays like |ω| −2 and can be transformed just as **C**(ω) for input noise closing the contour in the upper half-plane. The second and third term are the transposed complex conjugates of each other, because of the dependence of *H* on −ω instead of ω, and require a special consideration. Multiplied by *ei*ω*t* under the Fourier integral, the first term is proportional to *Hdei*ω*t* ∼ ω−1*ei*ω(*t*−*d*) and vanishes faster than |ω| −1 for large |ω| in the upper half-plane for *t* > *d* and in the lower half plane for *t* < *d*. For the second term the half planes are interchanged. The application of the residue theorem requires closing the integration contour in the half-plane where the integral over the semi-circle vanishes faster than |ω| −1. For **w** = **M** and in the general case of a stable dynamics all poles of the first term are in the upper half-plane (*zk*(*Lj*)) > 0, and have no contribution to **c**(*t*) for *t* < *d*. For the second term the same is true for *t* > −*d*; these terms correspond to the jumps of **c**(*t*) after one delay, caused by the effect of the sending neuron arriving at the other neurons in the system after one synaptic delay. These terms correspond to the response of the system to the impulse of the sending neuron – hence we call them "echo terms" in the following (Helias et al., 2013). The presence of such discontinuous jumps at time points *d* and −*d* in the case of output noise is reflected in the convolution of *h***w** with **D** in the time domain in (37). For input noise the absence of discontinuities can be inferred from the absence of such terms in (33), where the derivative of the correlation function is equal to the sum of finite terms. The first summand in (16) corresponds to the covariance evoked by fluctuations propagating through the system originating from the same neuron and we call it "correlated input term". In the system with input noise a similar separation into effective echo and correlated input terms can be performed. We obtain the correlated input term as the covariance in an auxiliary population without outgoing connections and echo terms as the difference between the full covariance between neurons within the network and the correlated input term.

### **2.7. EXPLICIT EXPRESSION FOR THE POPULATION AVERAGED CROSS COVARIANCE IN THE TIME DOMAIN**

We obtain the population averaged cross spectrum in a recurrent random network of Ornstein–Uhlenbeck processes (OUP) with input noise by inserting the averaged connectivity matrix **w** = **M** (14) into (8). The explicit expression for the covariance function follows by taking into account all (both) eigenvalues of **M** with values 0 and *L* = *Kw*(1 − γ*g*). The detailed derivation of the results presented in this section are documented in "Calculation of the Population Averaged Cross Covariance in Time Domain". The expression for the cross spectrum (8) takes the form

$$\mathbf{C}(\boldsymbol{\omega}) = f(\boldsymbol{\omega})f(-\boldsymbol{\omega}) \left( 1 + K\nu \begin{pmatrix} \gamma \text{g } -\gamma \text{g} \\ 1 & -1 \end{pmatrix} H\_d(\boldsymbol{\omega}) \right)$$

$$\mathbf{D} \begin{pmatrix} 1 + K\nu \begin{pmatrix} \gamma \text{g} & 1 \\ -\gamma \text{g} & -1 \end{pmatrix} H\_d(-\boldsymbol{\omega}) \end{pmatrix}, \tag{17}$$

where we introduced *f*(ω) = (*Hd*(ω)−1 − *L*)−1 as a short hand. Sorting the terms by their dependence on ω, introducing the functions -1(ω), . . . , -4(ω) for this dependence, and ϕ1(*t*), . . . , ϕ4(*t*) for the corresponding functions in the time domain, the covariance in the time domain **c**(*t*) = 1 2π - +∞ −∞ **C**(ω)*ei*ω*t d*ω takes the form

$$\begin{split} \mathbf{c}(t) &= \mathbf{D}\boldsymbol{\varphi}\_{1}(t) \\ &+ K\boldsymbol{\nu} \left( \begin{pmatrix} \boldsymbol{\chi}\mathbf{g} & -\boldsymbol{\chi}\mathbf{g} \\ 1 & -1 \end{pmatrix} \mathbf{D}\boldsymbol{\varphi}\_{2}(t) + \mathbf{D} \begin{pmatrix} \boldsymbol{\chi}\mathbf{g} & 1 \\ -\boldsymbol{\chi}\mathbf{g} & -1 \end{pmatrix} \boldsymbol{\varphi}\_{3}(t) \right) \\ &+ K^{2}\boldsymbol{\nu}^{2} \begin{pmatrix} \boldsymbol{\chi}\mathbf{g} & -\boldsymbol{\chi}\mathbf{g} \\ 1 & -1 \end{pmatrix} \mathbf{D} \begin{pmatrix} \boldsymbol{\chi}\mathbf{g} & 1 \\ -\boldsymbol{\chi}\mathbf{g} & -1 \end{pmatrix} \boldsymbol{\varphi}\_{4}(t) .\end{split}$$

The previous expression is valid for arbitrary **D**. In simulations presented in this article we consider identical marginal input statistics for all neurons. In this case the averaged activities for excitatory and inhibitory neurons are the same, so we can insert the special form of **D** given in (15), which results in

$$
\mathbf{c}(t) = \frac{\rho^2}{N} \begin{pmatrix} 1 & 0 \\ 0 \ \mathbf{y}^{-1} \end{pmatrix} \boldsymbol{\varphi}\_1(t) \tag{18}
$$

$$
+ \frac{\rho^2}{N} \text{Kw} \begin{pmatrix} \text{\"g} & -\text{g} \\ 1 & -\text{\"} \end{pmatrix} \boldsymbol{\varphi}\_2(t) + \frac{\rho^2}{N} \text{Kw} \begin{pmatrix} \text{\"g} & 1 \\ -\text{g} & -\text{\"} \end{pmatrix} \boldsymbol{\varphi}\_3(t)
$$

$$
+ \frac{\rho^2}{N} (\text{\"\/} + 1) \boldsymbol{K}^2 \boldsymbol{w}^2 \begin{pmatrix} \text{\"\"g}^2 & \text{\"\"g} \\ \text{\"g} & \text{\"\"} \end{pmatrix} \boldsymbol{\varphi}\_4(t).
$$

The time-dependent functions ϕ1,..., ϕ4 are the same in both cases. Using the residue theorem ϕ*i*(*t*) = 1 2π - +∞ −∞ *i*(ω)*ei*ω*t d*ω = *i z*∈poles of *i* Res(*i*,*z*)*eizt* for *t* - 0 they can be expressed as a sum over the poles *zk*(*L*) given by (12) and the pole *z* = *i* τ of *Hd*(ω). At ω = *zk*(*L*) the residue of *f*(ω) is Res(*f* , ω = *zk*(*L*)) = *idL* + *i*τ*ei*ω*d* −1 , the residue of *Hd*(ω) at *z* = *i* τ is − *i* τ *ed*/τ, so that the explicit forms of ϕ1,..., ϕ4 follow as

$$\begin{split} \varphi\_{1}(t) &= \sum\_{\boldsymbol{\alpha} = \boldsymbol{z}\_{k}(L)} i \text{Res}(f, \boldsymbol{\alpha}) f(-\boldsymbol{\alpha}) e^{i \boldsymbol{\alpha} t} \\ \varphi\_{2}(t) &= \sum\_{\boldsymbol{\alpha} = \boldsymbol{z}\_{k}(L)} i \text{Res}(f, \boldsymbol{\alpha}) f(-\boldsymbol{\alpha}) H\_{d}(\boldsymbol{\alpha}) e^{i \boldsymbol{\alpha} t} \\ &+ \frac{e^{(d-t)/\tau}}{\tau} f\left(\frac{i}{\tau}\right) f\left(-\frac{i}{\tau}\right) \end{split}$$

$$\begin{split} \varphi\_{3}(t) &= \sum\_{\boldsymbol{\alpha} = \boldsymbol{z}\_{k}(L)} i \text{Res}(f, \boldsymbol{\alpha}) f(-\boldsymbol{\alpha}) H\_{d}(-\boldsymbol{\alpha}) e^{i \boldsymbol{\alpha} t} \\ \varphi\_{4}(t) &= \sum\_{\boldsymbol{\alpha} = \boldsymbol{z}\_{k}(L)} i \text{Res}(f, \boldsymbol{\alpha}) f(-\boldsymbol{\alpha}) H\_{d}(\boldsymbol{\alpha}) H\_{d}(-\boldsymbol{\alpha}) e^{i \boldsymbol{\alpha} t} \\ &+ \frac{e^{-t/\tau}}{2\tau} f\left(\frac{i}{\tau}\right) f\left(-\frac{i}{\tau}\right). \end{split}$$

The corresponding expression for **C**(ω) for output noise is obtained by multiplying (17) with *H*−1 *d* (ω)*H*−1 *d* (−ω) = (1 + ω2τ2)

$$\mathbf{C}(\omega) = H\_d^{-1}(\omega) H\_d^{-1}(-\omega) f(\omega) f(-\omega) \tag{20}$$

$$\times (\mathbf{1} + \mathbf{K}\mathbf{w} \begin{pmatrix} \gamma \mathbf{g} & -\gamma \mathbf{g} \\ 1 & -1 \end{pmatrix} H\_d(\omega)) \mathbf{D}(\mathbf{1} + \mathbf{K}\mathbf{w} \begin{pmatrix} \gamma \mathbf{g} & 1 \\ -\gamma \mathbf{g} & -1 \end{pmatrix} H\_d(-\omega)),$$

which, after Fourier transform, provides the expression for **c**(*t*) in the time domain for *t* -0

$$\begin{split} \mathbf{c}(t) &= \mathbf{M} \mathbf{D} \mathbf{M}^T \boldsymbol{\varphi}\_1(t) + \mathbf{M} \mathbf{D} \boldsymbol{\varphi}\_0(t) + \mathbf{D} \boldsymbol{\delta}(t) \\ &= K^2 \boldsymbol{\nu}^2 \frac{\boldsymbol{\rho}^2}{N} (1 + \boldsymbol{\gamma} \mathbf{g}^2) \begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix} \boldsymbol{\varphi}\_1(t) + K \boldsymbol{\nu} \frac{\boldsymbol{\rho}^2}{N} \begin{pmatrix} 1 & - \mathbf{g} \\ 1 & - \mathbf{g} \end{pmatrix} \boldsymbol{\varphi}\_0(t) \\ &+ \frac{\boldsymbol{\rho}^2}{N} \begin{pmatrix} 1 & 0 \\ 0 & \boldsymbol{\gamma}^{-1} \end{pmatrix} \boldsymbol{\delta}(t). \end{split} \tag{21}$$

As in (18), the first line holds for arbitrary **D**, and the second for **D** given by (15), valid if the firing rates are homogeneous. ϕ1 is defined as before, and

$$\varphi\_0(t) = \theta(t - d) \sum\_{\omega = z\_k(L)} \left( dL + \tau e^{i\alpha d} \right)^{-1} e^{i\alpha t} \tag{22}$$

vanishes for *t* < *d*. All matrix elements of the first term in (21) are identical. Therefore all elements of **c**(*t*) are equal for 0 < |*t*| < *d*. Both rows of the matrix in front of ϕ0 are identical, so for *t* > 0 the off diagonal term *cIE* coincides with *cEE* and *cEI* with *cII* and vice versa for *t* < 0.

As an illustration we show the functions ϕ0,..., ϕ4 for one set of parameters in **Figure 4**. The left panels **(A,C)** correspond to contributions to the covariance caused by common input to a pair of neurons, the right panels **(B,D)** to terms due to the effect of one of the neurons' activities on the remaining network (echo terms). The upper panels **(A,B)** belong to the model with input noise, the lower panel **(C,D)** to the one with output noise.

**FIGURE 4 | Functions ϕ0 (D), ϕ1 (C), ϕ2, ϕ3 (B), ϕ4 (A) introduced in (19) and (22) for decomposition of covariance c***(t)***.** In **(B)** ϕ3(−*t*) is shown in gray and ϕ2(*t*) in black. The two functions are continuations of each other, joint at *t* = 0. Both functions appear in the echo term for input noise. The function ϕ0 in **(D)** describing the corresponding echo term in the case of output noise is shifted to be aligned with the function in **(B)** to facilitate the comparison of **(B,D)**. Parameters in all panels are *d* = 3 ms, τ = 10 ms, *L* = −1.72.

For the rate dynamics with output noise, the term with ϕ1 in (21) (shown in **Figure 4C**) is symmetric and describes the common input covariance and the term with ϕ0 (shown in **Figure 4D**) is the echo part of the covariance. For the rate dynamics with input noise (18) the term containing ϕ4 (shown in **Figure 4A**) is caused by common input and is hence also symmetric, the terms with ϕ2 and ϕ3 (shown in **Figure 4B**) correspond to the echo part and have hence their peak outside the origin. The second echo term in (18) is equal to the first one transposed and with opposite sign of the time argument, so we show ϕ2(*t*) and ϕ3(−*t*) together in one panel in **Figure 4B**. Note that for input noise, the term with ϕ1 describes the autocovariance, which corresponds to the term with the δ-function in case of output noise.

The solution (18) is visualized in **Figure 6**, the solution (21) in **Figure 7**, and the decomposition into common input and echo parts is also shown and compared to direct simulations in **Figure 8**.

#### **3. BINARY NEURONS**

In the following sections we study, in turn, the binary neuron model, the Hawkes model and the LIF model and show how they can be mapped to one of the two OUPs; either the one with input or the one with output noise, so that the explicit solutions (18) and (21) for the covariances derived in the previous section can be applied. In the present section, we start with the binary neuron model (Ginzburg and Sompolinsky, 1994; Buice et al., 2009).

Following Ginzburg and Sompolinsky (1994) the state of the network of *N* binary model neurons is described by a binary vector **n** ∈ {0, 1}*N* and each neuron is updated at independently drawn time points with exponentially distributed intervals of mean duration τ. This stochastic update constitutes a source Grytskyy et al. Unified correlations

of noise in the system. Given the *i*-th neuron is updated, the probability to end in the up-state (*ni* = 1) is determined by the gain function *Fi*(**n**) which depends on the activity **n** of all other neurons. The probability to end in the down state (*ni* = 0) is 1 − *Fi*(**n**). Here we implemented the binary model in the NEST simulator (Gewaltig and Diesmann, 2007) as described in "Implementation of Binary Neurons in a Spiking Simulator Code". Such systems have been considered earlier (Ginzburg and Sompolinsky, 1994; Buice et al., 2009), and here we follow the notation employed in the latter work. In the following we collect results that have been derived in these works and refer the reader to these publications for the details of the derivations. The zero-time lag covariance function is defined as *cij*(*t*) = *ni*(*t*)*nj*(*t*) − *ai*(*t*)*aj*(*t*), with the expectation value taken over different realizations of the stochastic dynamics. Here **a**(*t*) = (*a*1(*t*), . . . , *aN*(*t*))*T* is the vector of mean activities *ai*(*t*) = *ni*(*t*). *cij*(*t*) fulfills the differential equation

$$\mathbf{u}\frac{d}{dt}\mathbf{c}\_{\vec{\eta}}(t) = -2\mathbf{c}\_{\vec{\eta}}(t) + \langle (\mathbf{u}\_{\vec{\jmath}}(t) - a\_{\vec{\jmath}}(t))F\_{\vec{\imath}}(\mathbf{n}) \rangle$$

$$+ \langle (\mathbf{n}\_{\vec{\imath}}(t) - a\_{\vec{\imath}}(t))F\_{\vec{\jmath}}(\mathbf{n}) \rangle.$$

In the stationary state, the covariance therefore fulfills

$$c\_{\vec{\eta}} = \frac{1}{2} \langle (\eta\_{\vec{\eta}} - a\_{\vec{\eta}}) F\_{\vec{\imath}}(\mathbf{n}) \rangle + \frac{1}{2} \langle (\eta\_{\vec{\imath}} - a\_{\vec{\imath}}) F\_{\vec{\jmath}}(\mathbf{n}) \rangle. \tag{23}$$

The time lagged covariance *cij*(*t*, *s*) = *ni*(*t*)*nj*(*s*) − *ai*(*t*)*aj*(*s*) fulfills for *t* > *s* the differential equation

$$\text{tr}\frac{d}{dt}\boldsymbol{c}\_{\vec{\eta}}(t,\mathbf{s}) = -\boldsymbol{c}\_{\vec{\eta}}(t,\mathbf{s}) + \langle F\_i(\mathbf{n},\mathbf{t})(\eta\_j(\mathbf{s}) - a\_j(\mathbf{s})) \rangle. \tag{24}$$

This equation is also true for *i* = *j*, the autocovariance. The term *Fi*(**n**, *t*)(*nj*(*s*) − *aj*(*s*)) has a simple interpretation: it measures the influence of a fluctuation of neuron *j* at time *s* around its mean value on the gain of neuron *i* at time *t* (Ginzburg and Sompolinsky, 1994). We now assume a particular form for the coupling between neurons

$$F\_l(\mathbf{n}, t) = \phi(\mathbf{J}\_l \mathbf{n}(t - d)) = \phi\left(\sum\_{k=1}^{N} J\_{ik} \eta\_k(t - d)\right), \quad (25)$$

where **J***i* is the vector of incoming synaptic weights into neuron *i* and φ is a non-linear gain function. Assuming that the fluctuations of the total input **J***i***n** into the *i*-th neuron are sufficiently small to allow a linearization of the gain function φ, we obtain the Taylor expansion

$$F\_i(\mathbf{n}, t) = F\_i(\mathbf{a}) + \phi'(\mathbf{J}\_i \mathbf{a}) \mathbf{J}\_i(\mathbf{n}(t - d) - \mathbf{a}(t - d)),$$

where

$$\phi'(\mathbf{J}\_i \mathbf{a})\tag{26}$$

is the slope of the gain function at the point of mean input.

Up to this point the treatment of the system is identical to the work of Ginzburg and Sompolinsky (1994). Now we present an alternative approach for the linearization which takes into account the effect of fluctuations in the input. For sufficiently asynchronous network states, the fluctuations in the input **J***i***n**(*t* − *d*) to neuron *i* can be approximated by a Gaussian distribution *N* (μ, σ). In the following we consider a homogeneous random network with fixed in-degree as described in "Populationaveraged covariances". As each neuron receives the same number *K* of excitatory and γ*K* inhibitory synapses, the marginal statistics of the summed input to each neuron is identical. The mean input to a neuron then is μ = *KJ*(1 − γ*g*)*a*, where *a* is the mean activity of a neuron in the network. If correlations are small, the variance of this input signal distribution can be approximated as the sum of the variances of the individual contributions from the incoming signals, resulting in σ2 = *KJ*2(1 + γ*g*2) *a*(1 − *a*), where we used the fact that the variance of a binary variable with mean *a* is *a*(1 − *a*). This results from a direct calculation: since *n* ∈ {0, 1}, *n*2 = *n*, so that the variance is *n*2−*n*2 = *n*−*n*2 = *a*(1 − *a*). Averaging the slope φ of the gain function over the distribution of the input variable results in the averaged slope

$$
\langle \phi' \rangle \cong \int\_{-\infty}^{\infty} \mathcal{N}(\mu, \sigma, \mathbf{x}) \, \phi'(\mathbf{x}) \, d\mathbf{x} \tag{27}
$$

$$
\mathbf{1} \qquad \Big( \quad (\mathbf{x} - \mu)^2 \Big)
$$

$$\text{with}\,\mathcal{N}(\mu,\sigma,\mathbf{x}) = \frac{1}{\sqrt{2\pi}\sigma} \exp\left(-\frac{(\mathbf{x}-\mu)^2}{2\sigma^2}\right).$$

The two alternative methods of linearization of φ are illustrated in **Figure 5**. In the given example, the linearization procedure taking into account the fluctuations of the input signal results in a smaller effective slope φ than taking the slope φ (*a*) at the mean activity *a* near its maximum. Averaging the slope φ over this distribution fits simulation results better than φ (*a*) calculated at the mean of *a*, as shown in **Figure 6**.

The finite slope of the non-linear gain function can be understood as resulting from the combination of a hard threshold with an intrinsic local source of noise. The inverse strength of this noise

**FIGURE 5 | Alternative linearizations of the binary neuron model.** The black curve represents the non-linear gain function φ(*x*) = 1 2 + 1 2 tanh(β*x*). The dashed gray line is its tangent at the mean input value (denoted by the diagonal cross). The solid curve is the slope φ averaged over the distribution of the fluctuating input (27). This distribution estimated from direct simulation is presented by black dots, the corresponding theoretical prediction of a normal distribution *N*(μ, σ) (27) is shown as the light gray curve.

determines the slope parameter β (Ginzburg and Sompolinsky, 1994). In this sense, the network model contains two sources of noise, the explicit local noise, quantified by β and the fluctuating synaptic input interpreted as self-generated noise on the network level, quantified by σ. Even in the absence of local noise (β → ∞), the above mentioned linearization is applicable and yields a finite effective slope φ (27). In the latter case the resulting effective synaptic weight is independent of the original synapse strength (Grytskyy et al., 2013).

We now extend the classical treatment of covariances in binary networks (Ginzburg and Sompolinsky, 1994) by synaptic conduction delays. In (25) *Fi*(**n**, *t*) must therefore be understood as a functional acting on the function **n**(*t* ) for *t* ∈ [−∞, *t*], so that also synaptic connections with time delay *d* can be realized. We define an effective weight vector to absorb the gain factor as **w***i* = β*i***J***i*, with either β*i* = φ (μ) or β*i* = φ depending on the linearization procedure, and expand the right hand side of (24) to obtain

$$
\langle F\_{\vec{l}}(\mathbf{n},t)(n\_{\vec{j}}(s) - a\_{\vec{j}}(s))\rangle = \sum\_{k=1}^{N} \,\omega\_{ik}c\_{kj}(t-d,s).
$$

Thus the cross-covariance fulfills the matrix delay differential equation

$$\mathfrak{tr}\frac{d}{dt}\mathbf{c}(t,s) + \mathbf{c}(t,s) = \mathbf{w}\mathbf{c}(t-d,s). \tag{28}$$

This differential equation is valid for *t* > *s*. For the stationary solution, the differential equation only depends on the relative timing *u* = *t* − *s*

$$\mathfrak{tr}\frac{d}{du}\mathfrak{c}(u) + \mathfrak{c}(u) = \mathfrak{wc}(u-d). \tag{29}$$

The same linearization applied to (23) results in the boundary condition for the solution of the previous equation

$$\mathbf{2c}(0) = \mathbf{wc}(-d) + \left(\mathbf{wc}(-d)\right)^T \tag{30}$$

or, if we split **c** into its diagonal and its off-diagonal parts **c***a* and **c**=

$$\mathbf{2c}\_{\neq}(0) = \mathbf{w}\mathbf{c}\_{\neq}(-d) + \left(\mathbf{w}\mathbf{c}\_{\neq}(-d)\right)^{T} + \mathbf{O} \tag{31}$$

$$\text{with } \mathbf{O} = \mathbf{w}\mathbf{c}\_{a}(-d) + \left(\mathbf{w}\mathbf{c}\_{a}(-d)\right)^{T}.$$

In the following section we use this representation to demonstrate the equivalence of the covariance structure of binary networks to the solution for OUP with input noise.

#### **3.1. EQUIVALENCE OF BINARY NEURONS AND ORNSTEIN–UHLENBECK PROCESSES**

In the following subsection we show that the same Equations (29) and (31) for binary neurons also hold for the Ornstein-Uhlenbeck process (OUP) with input noise. In doing so here we also extend the existing framework of OUP (Risken, 1996) to synaptic conduction delays *d*. A network of such processes is described by

$$\mathbf{r}\frac{d}{dt}\mathbf{r}(t) + \mathbf{r}(t) = \mathbf{w}\mathbf{r}(t-d) + \mathbf{x}(t),\tag{32}$$

where **x** is a vector of pairwise uncorrelated white noise with **x**(*t*)*x* = 0 and *xi*(*t*)*xj*(*t* + *t* )*x* = δ*ij*δ(*t* )ρ2. With the help of the Green's function *G* satisfying (τ *d dt* + 1)*G*(*t*) = δ(*t*), namely *G*(*t*) = 1 τ θ(*t*)*e*−*t*/τ, we obtain the solution of Equation (32) as

$$\mathbf{r}(t) = \mathbf{r}G(t)\mathbf{r}(0) + \int\_0^t \mathbf{G}(t - t')(\mathbf{w}\mathbf{r}(t' - d) + \mathbf{x}(t'))\,dt'.$$

The equation for the fluctuations δ**r**(*t*) = **r**(*t*) − **r**(*t*)*x* around the expectation value

$$\delta \mathbf{r}(t) = \int\_0^t G(t - t') (\mathbf{w} \delta \mathbf{r}(t' - d) + \mathbf{x}(t')) \, dt'$$

coincides with the noisy rate model with input noise (4) with delay *d* and convolution kernel *h* = *G*. In the next step we investigate the covariance matrix *cij*(*t*, *s*) = δ*ri*(*t* + *s*)δ*rj*(*t*)*x* to show for which choice of parameters the covariance matrices for the binary model and the OUP with input noise coincide. To this end we derive the differential equation with respect to the time lag *s* for positive lags *s* > 0

$$\mathbf{r}\frac{d}{ds}\mathbf{c}(t,s) = \left\langle \mathbf{r}\frac{d}{ds}\delta\mathbf{r}(t+s)\delta\mathbf{r}^T(t) \right\rangle\_{\mathbf{x}} \tag{33}$$

$$= \left\langle (\mathbf{w}\delta\mathbf{r}(t+s-d) - \delta\mathbf{r}(t+s) + \mathbf{x}(t+s))\delta\mathbf{r}^T(t) \right\rangle\_{\mathbf{x}}$$

$$= \mathbf{w}\mathbf{c}(t,s-d) - \mathbf{c}(t,s),$$

where we used **x**(*t* + *s*))δ**r**(*t*)*x* = 0, because the noise is realized independently for each time step and the system is causal. Equation (33) is identical to the differential equation satisfied by the covariance matrix (28) for binary neurons (Ginzburg and Sompolinsky, 1994). To determine the initial condition of (33) we need to take the limit **c**(*t*, 0) = lim*s*→+0 **c**(*t*,*s*). This initial condition can be obtained as the stationary solution of the following differential equation

$$\begin{split} \pi \frac{d}{dt} \mathbf{c}(t, \mathbf{0}) &= \lim\_{s \to +0} \left( \left\langle \mathbf{r} \frac{d}{dt} \delta \mathbf{r}(t+s) \delta \mathbf{r}^T(t) \right\rangle\_{\mathbf{x}} + \left\langle \delta \mathbf{r}(t+s) \mathbf{r} \frac{d}{dt} \delta \mathbf{r}^T(t) \right\rangle\_{\mathbf{x}} \right) \\ &= \lim\_{s \to +0} \left( \langle (\mathbf{w} \delta \mathbf{r}(t+s-d) - \delta \mathbf{r}(t+s) + \mathbf{x}(t+s)) \delta \mathbf{r}^T(t) \rangle\_{\mathbf{x}} \right) \\ &\quad + \langle \delta \mathbf{r}(t+s) (\delta \mathbf{r}^T(t-d) \mathbf{w}^T - \delta \mathbf{r}^T(t) + \mathbf{x}^T(t)) \rangle\_{\mathbf{x}} \right) \\ &= -2\mathbf{c}(t, 0) + \mathbf{w} \mathbf{c}(t, -d) + \mathbf{c}(t-d, d) \mathbf{w}^T + \mathbf{D}. \end{split}$$

Here we used that **x**(*t* + *s*)δ**r***T*(*t*) vanishes due to independent noise realizations and causality and

$$\begin{split} \mathbf{D} &= \lim\_{s \to +0} \langle \delta \mathbf{r}(t+s) \mathbf{x}^T(t) \rangle\_{\mathbf{x}} \\ &= \lim\_{s \to +0, \, s \prec d} \int\_0^{t+s} G(t+s-t') (\mathbf{w} \underbrace{\langle \delta \mathbf{r}(t'-d) \mathbf{x}^T(t) \rangle\_{\mathbf{x}}}\_{=0 \text{ causality}} + \underbrace{\langle \mathbf{x}(t') \mathbf{x}^T(t) \rangle\_{\mathbf{x}}}\_{=1 \delta(t-t')\rho^2} dt' \\ &= \lim\_{s \to +0, \, s \prec d} \int\_0^{t+s} G(t+s-t') 1 \delta(t-t') \rho^2 dt' \\ &= \lim\_{s \to +0, \, s \prec d} G(s) \mathbf{1} \boldsymbol{\rho}^2 = \frac{1}{\pi} \mathbf{1} \boldsymbol{\rho}^2. \end{split}$$

In the stationary state, **c** only depends on the time lag *s* and is independent of the first time argument *t*, which, with the symmetry **c**(−*d*)*T* = **c**(*d*) yields the additional condition for the solution of (33)

$$\mathbf{2c}(0) = \mathbf{wc}(-d) + (\mathbf{wc}(-d))^T + \mathbf{D}$$

or, if **c** is split in diagonal and off-diagonal parts **c***a* and **c**=, respectively,

$$2\mathbf{c}\_{\neq}(0) = \mathbf{w}\mathbf{c}\_{\neq}(-d) + \left(\mathbf{w}\mathbf{c}\_{\neq}(-d)\right)^{T} + \mathbf{O}$$

$$2\mathbf{c}\_{a}(0) = \mathbf{w}\mathbf{c}\_{\neq}(-d) + \left(\mathbf{w}\mathbf{c}\_{\neq}(-d)\right)^{T} + \mathbf{D}$$