We first analyze data from motor cortex of macaques during rest, recorded with $4\times 4{\mathrm{mm}}^2$, 100-electrode Utah arrays with 400 µm inter-electrode distance (Figure 1A). The resting condition of motor cortex in monkeys is ideal to assess intrinsic coordination between neurons during ongoing activity. In particular, our analyses focus on true resting state data, devoid of movement-related transients in neuronal firing (see Materials and methods). Parallel single-unit spiking activity of $\approx 130$ neurons per recording session, sorted into putative excitatory and inhibitory cells, shows strong spike-count correlations across the entire Utah array, well beyond the typical scale of the underlying short-range connectivity profiles (Figure 1B and D). Positive and negative correlations form patterns in space that are furthermore seemingly unrelated to the neuron types. All populations show a large dispersion of both positive and negative correlation values (Figure 1C).

Figure 1

Download asset Open asset

Figure 1—source data 1

Code and data.

https://cdn.elifesciences.org/articles/68422/elife-68422-fig1-data1-v1.zip

Download elife-68422-fig1-data1-v1.zip

The classical view on pairwise correlations in balanced networks (Ginzburg and Sompolinsky, 1994; Renart et al., 2010; Pernice et al., 2011; Pernice et al., 2012; Tetzlaff et al., 2012; Helias et al., 2014) focuses on averages across many pairs of cells: average correlations are small if the network dynamics is stabilized by an excess of inhibitory feedback; dynamics known as the ‘balanced state’ arise (van Vreeswijk and Sompolinsky, 1996; Amit and Brunel, 1997; van Vreeswijk and Sompolinsky, 1998): Negative feedback counteracts any coherent increase or decrease of the population-averaged activity, preventing the neurons from fluctuating in unison (Tetzlaff et al., 2012). Breaking this balance in different ways leads to large correlations (Rosenbaum and Doiron, 2014; Darshan et al., 2018; Baker et al., 2019). Can the observation of significant correlations between individual cells across large distances be reconciled with the balanced state? In the following, we provide a mechanistic explanation.

Connections mediate interactions between neurons. Many studies therefore directly relate connectivity and correlations (Pernice et al., 2011; Pernice et al., 2012; Trousdale et al., 2012; Brinkman et al., 2018; Kobayashi et al., 2019). From direct connectivity, one would expect positive correlations between excitatory neurons and negative correlations between inhibitory neurons and a mix of negative and positive correlations only for excitatory-inhibitory pairs. Likewise, a shared input from inside or outside the network only imposes positive correlations between any two neurons (Figure 2A). The observations that excitatory neurons may have negative correlations (Figure 1D), as well as the broad distribution of correlations covering both positive and negative values (Figure 1C), are not compatible with this view. In fact, the sign of correlations appears to be independent of the neuron types. So how do negative correlations between excitatory neurons arise?

Figure 2

Download asset Open asset

The view that equates connectivity with correlation implicitly assumes that the effect of a single synapse on the receiving neuron is weak. This view, however, regards each synapse in isolation. Could there be states in the network where, collectively, many weak synapses cooperate, as perhaps required to form low-dimensional neuronal manifolds? In such a state, interactions may not only be mediated via direct connections but also via indirect paths through the network (Figure 2B). Such effective multi-synaptic connections may explain our observation that far apart neurons that are basically unconnected display considerable correlation of arbitrary sign.

Let us here illustrate the ideas first and corroborate them in subsequent sections. Direct connections yield correlations of a predefined sign, leading to correlation distributions with multiple peaks, for example a positive peak for excitatory neurons that are connected and a peak at zero for neurons that are not connected. Multi-synaptic paths, however, involve both excitatory and inhibitory intermediate neurons, which contribute to the interaction with different signs (Figure 2B). Hence, a single indirect path can contribute to the total interaction with arbitrary sign (Pernice et al., 2011). If indirect paths dominate the interaction between two neurons, the sign of the resulting correlation becomes independent of their type. Given that the connecting paths in the network are different for any two neurons, the resulting correlations can fall in a wide range of both positive and negative values, giving rise to the broad distributions for all combinations of neuron types in Figure 1C. This provides a hypothesis why there may be no qualitative difference between the distribution of correlations for excitatory and inhibitory neurons. In fact, their widths are similar and their mean is close to zero (see Materials and methods for exact values); the latter being the hallmark of the negative feedback that characterizes the balanced state. The subsequent model-based analysis will substantiate this idea and show that it also holds for networks with spatially organized heterogeneous connectivity.

To play this hypothesis further, an important consequence of the dominance of multi-synaptic connections could be that correlations are not restricted to the spatial range of direct connectivity. Through interactions via indirect paths the reach of a single neuron could effectively be increased. But the details of the spatial profile of the correlations in principle could be highly complex as it depends on the interplay of two antagonistic effects: On the one hand, signal propagation becomes weaker with distance, as the signal has to pass several synaptic connections. Along these paths mean firing rates of neurons are typically diverse, and so are their signal transmission properties (de la Rocha et al., 2007). On the other hand, the number of contributing indirect paths between any pair of neurons proliferates with their distance. With single neurons typically projecting to thousands of other neurons in cortex, this leads to involved combinatorics; intuition here ceases to provide a sensible hypothesis on what is the effective spatial profile and range of coordination between neurons. Also it is unclear which parameters these coordination patterns depend on. The model-driven and analytical approach of the next section will provide such a hypothesis.

We first note that the large magnitude and dispersion of individual correlations in the data and their spatial structure primarily stem from features in the underlying covariances between neuron pairs (Appendix 1—figure 1). Given the close relationship between correlations and covariances (Appendix 1—figure 1D and E), in the following we analyze covariances, as these are less dependent on single neuron properties and thus analytically simpler to treat. To gain an understanding of the spatial features of intrinsically generated covariances in balanced critical networks, we investigate a network of excitatory and inhibitory neurons on a two-dimensional sheet, where each neuron receives external Gaussian white noise input (Figure 3A). We investigate the covariance statistics in this model by help of linear-response theory, which has been shown to approximate spiking neuron models well (Pernice et al., 2012; Trousdale et al., 2012; Tetzlaff et al., 2012; Helias et al., 2013; Grytskyy et al., 2013; Dahmen et al., 2019). To allow for multapses, the connections between two neurons are drawn from a binomial distribution, and the connection probability decays with inter-neuronal distance on a characteristic length scale $d$ (for more details see Materials and methods). Previous studies have used linear-response theory in combination with methods from statistical physics and field theory to gain analytic insights into both mean covariances (Ginzburg and Sompolinsky, 1994; Lindner et al., 2005; Pernice et al., 2011; Tetzlaff et al., 2012) and the width of the distribution of covariances (Dahmen et al., 2019). Field-theoretic approaches, however, were so far restricted to purely random networks devoid of any network structure and thus not suitable to study spatial features of covariances. To analytically quantify the relation between the spatial ranges of covariances and connections, we therefore here develop a theory for spatially organized random networks with multiple populations. The randomness in our model is based on the sparseness of connections, which is one of the main sources of heterogeneity in cortical networks in that it contributes strongly to the variance of connections (see Appendix 1 Section 15).

Figure 3

Download asset Open asset

Figure 3—source data 1

Code and data.

https://cdn.elifesciences.org/articles/68422/elife-68422-fig3-data1-v1.zip

Download elife-68422-fig3-data1-v1.zip

A distance-resolved histogram of the covariances in the spatially organized E-I network shows that the mean covariance is close to zero but the width or variance of the covariance distribution stays large, even for large distances (Figure 3C). Analytically, we derive that, despite the complexity of the various indirect interactions, both the mean and the variance of covariances follow simple exponential laws in the long-distance limit (see Appendix 1 Section 4 - Section 12). These laws are universal in that they do not depend on details of the spatial profile of connections. Our theory shows that the associated length scales are strikingly different for means and variances of covariances. They each depend on the reach of direct connections and on specific eigenvalues of the effective connectivity matrix. These eigenvalues summarize various aspects of network connectivity and signal transmission into a single number: Each eigenvalue belongs to a ‘mode’, a combination of neurons that act collaboratively, rather than independently, coordinating neuronal activity within a one-dimensional subspace. To start with, there are as many such subspaces as there are neurons. But if the spectral bound in Figure 3B is close to one, only a relatively small fraction of them, namely those close to the spectral bound, dominate the dynamics; the dynamics is then effectively low-dimensional. Additionally, the eigenvalue quantifies how fast a mode decays when transmitted through a network. The eigenvalues of the dominating modes are close to one, which implies a long lifetime. The corresponding fluctuations thus still contribute significantly to the overall signal, even if they passed by many synaptic connections. Therefore, indirect multi-synaptic connections contribute significantly to covariances if the spectral bound is close to one, and in that case we expect to see long-range covariances.

To quantify this idea, for the mean covariance $\overline{c}$ we find that the dominant behavior is an exponential decay $\overline{c}\sim \exp (-x/\overline{d})$ on a length scale $\overline{d}$. This length scale is determined by a particular eigenvalue, the population eigenvalue, corresponding to the mode in which all neurons are excited simultaneously. Its position solely depends on the ratio between excitation and inhibition in the network and becomes more negative in more strongly inhibition-dominated networks (Figure 3B). We show in Appendix 1 Section 9.4 that this leads to a steep decay of mean covariances with distance. The variance of covariances, however, predominantly decays exponentially on a length scale d_eff that is determined by the spectral bound $R$, the largest real part among all eigenvalues (Figure 3B and D). In inhibition-dominated networks, $R$ is determined by the heterogeneity of connections. For $R\lesssim 1$ we obtain the effective length scale

What this means is that precisely at the point where $R$ is close to one, when neural activity occupies a low-dimensional manifold, the length scale d_eff on which covariances decay exceeds the reach of direct connections by a large factor (Figure 3D). As the network approaches instability, which corresponds to the spectral bound $R$ going to one, the effective decay constant diverges (Figure 3D inset) and so does the range of covariances.

Our population-resolved theoretical analysis, furthermore, shows that the larger the spectral bound the more similar the decay constants between different populations, with only marginal differences for $R\lesssim 1$ (Figure 3E). This holds strictly if connection weights only depend on the type of the presynaptic neuron but not on the type of the postsynaptic neuron. Moreover, we find a relation between the squared effective decay constants and the squared anatomical decay constants of the form

This relation is independent of the eigenvalues of the effective connectivity matrix, as the constant of order ${\displaystyle \mathcal{O}(1)}$ does only depend on the choice of the connectivity profile. For $R\simeq 1$, this means that even though the absolute value of both effective length scales on the left hand side is large, their relative difference is small because it equals the small difference of anatomical length scales on the right hand side.

To check if these predictions are confirmed by the data from macaque motor cortex, we first observe that, indeed, covariances in the resting state show a large dispersion over almost all distances on the Utah array (Figure 4). Moreover, the variance of covariances agrees well with the predicted exponential law: Performing an exponential fit reveals length constants above 1 mm. These large length constants have to be compared to the spatial reach of direct connections, which is about an order of magnitude shorter, in the range of 100-400 μm (Schnepel et al., 2015), so below the 400 μm inter-electrode distance of the Utah array. The shallow decay of the variance of covariances is, next to the broad distribution of covariances, a second indication that the network is in the dynamically balanced critical regime, in line with the prediction by Equation (1).

Figure 4

Download asset Open asset

Figure 4—source data 1

Code and data.

https://cdn.elifesciences.org/articles/68422/elife-68422-fig4-data1-v1.zip

Download elife-68422-fig4-data1-v1.zip

The population-resolved fits to the data show a larger length constant for excitatory covariances than for inhibitory ones (Figure 4A). This is qualitatively in line with the prediction of Equation (2) given the – by tendency – longer reach of excitatory connections compared to inhibitory ones, as derived from morphological constraints (Reimann et al., 2017, Fig. S2). In the dynamically balanced critical regime, however, the predicted difference in slope for all three fits is practically negligible. Therefore, we performed a second fit where the slope of the three exponentials is constrained to be identical (Figure 4B). The error of this fit is only marginally larger than the ones of fitting individual slopes (Figure 4C). This shows that differences in slopes are hardly detectable given the empirical evidence, thus confirming the predictions of the theory given by Equation (1) and Equation (2).

Since covariances measure the coordination of temporal fluctuations around the individual neurons’ mean firing rates, they are determined by how strong a neuron transmits such fluctuations from input to output (Abeles, 1991). To leading order this is explained by linear-response theory (Ginzburg and Sompolinsky, 1994; Lindner et al., 2005; Pernice et al., 2011; Tetzlaff et al., 2012): How strongly a neuron reacts to a small change in its input depends on its dynamical state, foremost the mean and variance of its total input, called ‘working point’ in the following. If a neuron receives almost no input, a small perturbation in the input will not be able to make the neuron fire. If the neuron receives a large input, a small perturbation will not change the firing rate either, as the neuron is already saturated. Only in the intermediate regime the neuron is susceptible to small deviations of the input. Mathematically, this behavior is described by the gain of the neuron, which is the derivative of the input-output relation (Abeles, 1991). Due to the non-linearity of the input-output relation, the gain is vanishing for very small and very large inputs and non-zero in the intermediate regime. How strongly a perturbation in the input to one neuron affects one of the subsequent neurons therefore not only depends on the synaptic weight ${\displaystyle \mathit{J}}$ but also on the gain ${\displaystyle \mathit{S}}$ and thereby the working point. This relation is captured by the effective connectivity ${\displaystyle \mathit{W}=\mathit{S}\cdot \mathit{J}}$. What is the consequence of the dynamical interaction among neurons depending on the working point? Can it be used to reshape the low-dimensional manifold, the collective coordination between neurons?

The first part of this study finds that long-range coordination can be achieved in a network with short-range random connections if effective connections are sufficiently strong. Alteration of the working point, for example by a different external input level, can affect the covariance structure: The pattern of coordination between individual neurons can change, even though the anatomical connectivity remains the same. In this way, routing of information through the network can be adapted dynamically on a mesoscopic scale. This is a crucial difference of such coordination as opposed to coordination imprinted by complex but static connection patterns.

Here, we first illustrate this concept by simulations of a network of 2000 sparsely connected threshold-linear (ReLU) rate neuron models that receive Gaussian white noise inputs centered around neuron-specific non-zero mean values (see Materials and methods and Appendix 1 Section 14 for more details). The ReLU activation function thereby acts as a simple model for the vanishing gain for neurons with too low input levels. Note that in cortical-like scenarios with low firing rates, neuronal working points are far away from the high-input saturation discussed above, which is therefore neglected by the choice of the ReLU activation function. For independent and stationary external inputs covariances between neurons are solely generated inside the network via the sparse and random recurrent connectivity. External inputs only have an indirect impact on the covariance structure by setting the working point of the neurons.

We simulate two networks with identical structural connectivity and identical external input fluctuations, but small differences in mean external inputs between corresponding neurons in the two simulations (Figure 5A). These small differences in mean external inputs create different gains and firing rates and thereby differences in effective connectivity and covariances. Since mean external inputs are drawn from the same distribution in both simulations (Figure 5B), the overall distributions of firing rates and covariances across all neurons are very similar (Figure 5E1, F1). But individual neurons’ firing rates do differ (Figure 5E2). For the simple ReLU activation used here, we in particular observe neurons that switch between non-zero and zero firing rate between the two simulations. This resulting change of working points substantially affects the covariance patterns (Figure 5F2): Differences in firing rates and covariances between the two simulations are significantly larger than the differences across two different epochs of the same simulation (Figure 5C). The larger the spectral bound, the more sensitive are the intrinsically generated covariances to the changes in firing rates (Figure 5D). Thus, a small offset of individual firing rates is an effective parameter to control network-wide coordination among neurons. As the input to the local network can be changed momentarily, we predict that in the dynamically balanced critical regime coordination patterns should be highly dynamic.

Figure 5

Download asset Open asset

Figure 5—source data 1

Code and data.

https://cdn.elifesciences.org/articles/68422/elife-68422-fig5-data1-v1.zip

Download elife-68422-fig5-data1-v1.zip

In order to test the theoretical prediction in experimental data, we analyze parallel spiking activity from macaque motor cortex, recorded during a reach-to-grasp experiment (Riehle et al., 2013; Brochier et al., 2018). In contrast to the resting state, where the animal was in an idling state, here the animal is involved in a complex task with periods of different cognitive and behavioral conditions (Figure 6A). We compare two epochs in which the animal is requested to wait and is sitting still but which differ in cognitive conditions. The first epoch is a starting period (S), where the monkey has self-initiated the behavioral trial and is attentive because it is expecting a cue. The second epoch is a preparatory period (P), where the animal has just received partial information about the upcoming trial and is waiting for the missing information and the GO signal to initiate the movement.

Figure 6

Download asset Open asset

Figure 6—source data 1

Code and data.

https://cdn.elifesciences.org/articles/68422/elife-68422-fig6-data1-v1.zip

Download elife-68422-fig6-data1-v1.zip

Within each epoch, S or P, the neuronal firing rates are mostly stationary, likely due to the absence of arm movements which create relatively large transient activities in later epochs of the task, which are not analyzed here (see Appendix 1 Section 3). The overall distributions of the firing rates are comparable for epochs S and P, but the firing rates are distributed differently across the individual neurons: Figure 6C shows one example session of monkey N, where the changes in firing rates between the two epochs are visible in the spread of markers around the diagonal line in panel C2. To assess the extent of these changes, we split each epoch, S and P, into two disjoint sub-periods, S1/S2 and P1/P2 (Figure 6A). We compute the correlation coefficient between the firing rate vectors of two sub-periods of different epochs (‘between’ markers in Figure 6E) and compare it to the correlation coefficient between the firing rate vectors of two sub-periods of the same epoch (‘within’ markers): Firing rate vectors in S1 are almost perfectly correlated with firing rate vectors in S2 ($\rho \approx 1$ for all of the five/eight different recording sessions from different recording days for monkey E/N, similarly for P1 and P2), confirming stationarity investigated in Appendix 1 Section 3. Firing rate vectors in S1 or S2, however, show significantly lower correlation to firing rate vectors in P1 and P2, confirming a significant change in network state between epochs S and P (Figure 6E).

The mechanistic model in the previous section shows a qualitatively similar scenario (Figure 5C and E). By construction it produces different firing rate patterns in the two simulations. While the model is simplistic and in particular not adapted to quantitatively reproduce the experimentally observed activity statistics, its simulations and our underlying theory make a general prediction: Differences in firing rates impact the effective connectivity between neurons and thereby evoke even larger differences in their coordination if the network is operating in the dynamically balanced critical regime (Figure 5D). To check this prediction, we repeat the correlation analysis between the two epochs, which we described above for the firing rates, but this time for the covariance patterns. Despite similar overall distributions of covariances in S and P (Figure 6D1), covariances between individual neuron pairs are clearly different between S and P: Figure 6B shows the covariance pattern for one representative reference neuron in one example recording session of monkey N. In both epochs, this covariance pattern has a salt-and-pepper structure as for the resting state data in Figure 1D. Yet, neurons change their individual coordination: a large number of neuron pairs even changes from positive covariance values to negative ones and vice versa. These neurons fire cooperatively in one epoch of the task while they show antagonistic firing in the other epoch. The covariances of all neuron pairs of that particular recording session are shown in Figure 6D2. Markers in the upper left and lower right quadrant show neuron pairs that switch the sign of their coordination (45 % of all neuron pairs). The extent of covariance changes between epochs is again quantified by correlation coefficients between the covariance patterns of two sub-periods (Figure 6F). As for the firing rates, we find rather large correlations between covariance patterns in S1 and S2 as well as between covariance patterns in P1 and P2. Note, however, that correlation coefficients are around 0.8 rather than 1, presumably since covariance estimates from 200 ms periods are noisier than firing rate estimates. The covariance patterns in S1 or S2 are, however, significantly more distinct from covariance patterns in P1 and P2, with correlation coefficients around 0.5 (Figure 6F). This more pronounced change of covariances compared to firing rates is predicted by a network whose effective connectivity has a large spectral bound, in the dynamically balanced critical state. In particular, the theory provides a mechanistic explanation for the different coordination patterns between neurons on the mesoscopic scale (range of a Utah array), which are observed in the two states S and P (Figure 6B). The coordination between neurons is thus considerably reshaped by the behavioral condition.

Two adult macaque monkeys (monkey E - female, and monkey N - male) are recorded in behavioral experiments of two types: resting state and reach-to-grasp. The recordings of neuronal activity in motor and pre-motor cortex (hand/arm region) are performed with a chronically implanted $4\mathrm{x}4{\mathrm{mm}}^2$ Utah array (Blackrock Microsystems). Details on surgery, recordings, spike sorting and classification of behavioral states can be found in Riehle et al., 2013; Riehle et al., 2018; Brochier et al., 2018; Dąbrowska et al., 2020. All animal procedures were approved by the local ethical committee (C2EA 71; authorization A1/10/12) and conformed to the European and French government regulations.

The mean and variance of covariances are calculated for a two-dimensional network consisting of one excitatory and one inhibitory population of neurons. The connectivity profile $p(\mathit{x})$, describing the probability of a neuron having a connection to another neuron at distance $\mathit{x},$ decays with distance. We assume periodic boundary conditions and place the neurons on a regular grid (Figure 3A), which imposes translation and permutation symmetries that enable the derivation of closed-form solutions for the distance-dependent mean and variance of the covariance distribution. These simplifying assumptions are common practice and simulations show that they do not alter the results qualitatively.

Our aim is to find an expression for the mean and variance of covariances as functions of distance between two neurons. While the theory in Dahmen et al., 2019 is restricted to homogeneous connections, understanding the spatial structure of covariances here requires us to take into account the spatial structure of connectivity. Field-theoretic methods, combined with linear-response theory, allow us to obtain expressions for the mean covariance ${\displaystyle \overline{\mathit{c}}}$ and variance of covariance ${\displaystyle \overline{\delta {\mathit{c}}^2}}$

with identity matrix 1, mean ${\displaystyle \mathit{M}}$ and variance $\mathit{S}$ of connectivity matrix $\mathit{W}$, input noise strength $D$, and spectral bound $R$. Since $\mathit{M}$ and $\mathit{S}$ have a similar structure, the mean and variance can be derived in the same way, which is why we only consider variances in the following.

To simplify Equation (4), we need to find a basis in which $\mathit{S},$ and therefore also $\mathit{A}=\mathbf{1}-\mathit{S}$, is diagonal. Due to invariance under translation, the translation operators $\mathit{T}$ and the matrix $\mathit{S}$ have common eigenvectors, which can be derived using that translation operators satisfy ${\mathit{T}}^N=1$, where $N$ is the number of lattice sites in $x$- or $y$-direction (see Appendix 1). Projecting onto a basis of these eigenvectors shows that the eigenvalues $s_{\mathit{K}}$ of $\mathit{S}$ are given by a discrete two-dimensional Fourier transform of the connectivity profile$s_{\mathit{K}}\propto \sum _{x}p(x){\mathrm{e}}^{-\mathrm{i}\mathbf{KX}}\mathit{\hspace{1em}}.$

Expressing ${\mathit{A}}^{-1}$ in the eigenvector basis yields ${\mathit{A}}^{-1}(\mathit{x})=\mathbf{1}+\mathit{B}(\mathit{x})$, where $\mathit{B}(\mathit{x})$ is a discrete inverse Fourier transform of the kernel $s_{\mathit{K}}/(1-s_{\mathit{K}})$. Assuming a large network with respect to the connectivity profiles allows us to take the continuum limit

As we are only interested in the long-range behavior, which corresponds to $|\mathit{x}|\to \mathrm{\infty }$, or $|\mathit{K}|\to 0$, respectively, we can approximate the Fourier kernel around $|\mathit{K}|\approx 0$ by a rational function, quadratic in the denominator, using a Padé approximation. This allows us to calculate the integral which yields

where $K_0(x)$ denotes the modified Bessel function of second kind and zeroth order (Olver et al., 2010), and the effective decay constant d_eff is given by Equation (1). In the long-range limit, the modified Bessel function behaves like

Writing Equation (4) in terms of $B(\mathit{x})$ gives

with the double asterisk denoting a two-dimensional convolution. $(B\ast \ast B)(\mathit{x})$ is a function proportional to the modified Bessel function of second kind and first order (Olver et al., 2010), which has the long-range limit

Hence, the effective decay constant of the variances is given by d_eff. Note that further details of the above derivation can be found in the Appendix 1 Section 4 - Section 12.

The explanation of the network state dependence of covariance patterns presented in the main text is based on linear-response theory, which has been shown to yield results quantitatively in line with non-linear network models, in particular networks of spiking leaky integrate-and-fire neuron models (Tetzlaff et al., 2012; Trousdale et al., 2012; Pernice et al., 2012; Grytskyy et al., 2013; Helias et al., 2013; Dahmen et al., 2019). The derived mechanism is thus largely model independent. We here chose to illustrate it with a particularly simple non-linear input-output model, the rectified linear unit (ReLU). In this model, a shift of the network’s working point can turn some neurons completely off, while activating others, thereby leading to changes in the effective connectivity of the network. In the following, we describe the details of the network model simulation.

We performed a simulation with the neural simulation tool NEST (Jordan, 2019) using the parameters listed in Appendix 1—table 4. We simulated a network of $N$ inhibitory neurons (threshold_lin_rate_ipn, Hahne, 2017), which follow the dynamical equation

where z_i is the input to neuron i, $\nu$ the output firing rate with (threshold linear activation function)

time constant $\tau$, connectivity matrix $\mathit{J}$, a constant external input ${\mu }_{\mathrm{ext},i}$, and uncorrelated Gaussian white noise $⟨{\xi }_i(t)⟩=0$, $⟨{\xi }_i(s){\xi }_j(t)⟩={\delta }_{ij}\delta (s-t)$, with noise strength $\sqrt{\tau }{\sigma }_{\mathrm{noise},i}$. The neurons were connected using the fixed_indegree connection rule, with connection probability $p$, indegree $K=p\cdot N$, and delta-synapses (rate_connection_instantaneous) of weight $w$.

The constant external input ${\mu }_{\mathrm{ext},i}$ to each neuron was normally distributed, with mean ${\mu }_{\mathrm{ext}}$, and standard deviation ${\sigma }_{\mathrm{ext}}$. It was used to set the firing rates of neurons, which, via the effective connectivity, influence the intrinsically generated covariances in the network. The two parameters ${\mu }_{\mathrm{ext}}$ and ${\sigma }_{\mathrm{ext}}$ were chosen such that, in the stationary state, half of the neurons were expected to be above threshold. Which neurons are active depends on the realization of ${\mu }_{\mathrm{ext},i}$ and is therefore different for different networks.

To assess the distribution of firing rates, we first considered the static variability of the network and studied the stationary solution of the noise-averaged input ${⟨z⟩}_{\mathrm{noise}}$, which follows from Equation (5) as

Note that ${⟨{\nu }_j⟩}_{\mathrm{noise}}={⟨\varphi (z_j)⟩}_{\mathrm{noise}}$, through the nonlinearity $\varphi$, in principle depends on fluctuations of the system. This dependence is, however, small for the chosen threshold linear $\varphi$, which is only nonlinear in the point $z=0$.

The derivation of ${\mu }_{\mathrm{ext}}$ is based on the following mean-field considerations: according to Equation (6) the mean input to a neuron in the network is given by the sum of external input and recurrent input

The variance of the input is given by

The mean firing rate can be calculated using the diffusion approximation (Tuckwell, 2009; Amit and Tsodyks, 2009), which is assuming a normal distribution of inputs due to the central-limit theorem, and the fact that a linear threshold neuron only fires if its input is positive

where $P$ denotes the probability density of the firing rate $\nu$. The variance of the firing rates is given by

The number of active neurons is the number of neurons with a positive input, which we set to be equal to $N/2$

which is only fulfilled for $\mu =0$. Inserting this condition simplifies the equations above and leads to

For the purpose of relating synaptic weight $w$ and spectral bound $R$, we can view the nonlinear network as an effective linear network with half the population size (only the active neurons). In the latter case, we obtain

For a given spectral bound $R$, this relation allows us to derive the value

that, for a arbitrarily fixed ${\sigma }_{\mathrm{ext}}$ (here ${\sigma }_{\mathrm{ext}}=1$), makes half of the population being active. We were aiming for an effective connectivity with only weak fluctuations in the stationary state. Therefore, we fixed the noise strength for all neurons to the small value ${\sigma }_{\mathrm{noise}}=0.1\ll {\sigma }_{\mathrm{ext}}$ compared to the external input, such that the noise fluctuations did not have a large influence on the calculation above that determines which neurons were active.

To show the effect of a change in the effective connectivity on the covariances, we simulated two networks with identical connectivity, but supplied them with slightly different external inputs. This was realized by choosing

with

$ϵ\ll 1$, and $\alpha \in \{1,2\}$ indexing the two networks. The main component ${\mu }_{\mathrm{ext},i}$ of the external input was the same for both networks. But, the small component ${\mu }_{\mathrm{ext},i}^{(\alpha )}$ was drawn independently for the two networks. This choice ensures that the two networks have a similar external input distribution (Figure 5B1), but with the external inputs distributed differently across the single neurons (Figure 5B2). How similar the external inputs are distributed across the single neurons is determined by $ϵ$.

The two networks have a very similar firing rate distribution (Figure 5E1), but, akin to the external inputs, the way the firing rates are distributed across the single neurons differs between the two networks (Figure 5E2). As the effective connectivity depends on the firing rates

this leads to a difference in the effective connectivities of the two networks and therefore to different covariance patterns, as discussed in Figure 5.

We performed the simulation for spectral bounds ranging from 0.1 to 0.9 in increments of 0.1. We calculated the correlation coefficient of firing rates and the correlation coefficient of time-lag integrated covariances between $N_{\mathrm{sample}}$ neurons in the two networks (Figure 5D) and studied the dependence on the spectral bound.

To check whether the simulation was long enough to yield a reliable estimate of the rates and covariances, we split each simulation into two halves, and calculate the correlation coefficient between the rates and covariances from the first half of the simulation with the rates and covariances from the second half. They were almost perfectly correlated (Figure 5C). Then, we calculated the correlation coefficients comparing all halves of the first simulation with all halves of the second simulation, showing that the covariance patterns changed much more than the rate patterns (Figure 5C).

A typical measure for the strength of neuronal coordination is the Pearson correlation coefficient, here applied to spike counts in $1\mathrm{s}$ bins. Correlation coefficients, however, comprise features of both auto- and cross-covariances. From a theoretical point of view, it is simpler to study cross-covariances separately. Indeed, linear-response theory has been shown to faithfully predict cross-covariances in spiking leaky integrate-and-fire networks (Tetzlaff et al., 2012; Pernice et al., 2012; Trousdale et al., 2012; Helias et al., 2013; Dahmen et al., 2019; Grytskyy et al., 2013). Appendix 1—figure 1 justifies the investigation of cross-covariances instead of correlation coefficients for the purpose of this study. It shows that the spatial organization of correlations closely matches the spatial organization of cross-covariances.

Appendix 1—figure 1

Download asset Open asset

The analysis of the experimental data involves a number of preprocessing steps, which may affect the resulting statistics. In our study one such critical step is the separation of putative excitatory and inhibitory units, which is partially based on setting thresholds on the widths of spike waveform, as described in the Methods section. We tested the robustness of our conclusions with respect to these thresholds.

As mentioned in the Methods, two thresholds for the width of a spike waveform are chosen, based on all SU average waveforms: A width larger than the ‘‘broadness’’ threshold indicates a putative excitatory neuron. A width lower than the ‘‘narrowness’’ threshold indicates a putative inhibitory neuron. Units with intermediate widths are unclassified. Additionally, to increase the reliability of the classification, we perform it in two steps: first on the SU’s average waveform, and second on all its single waveforms. We calculate the percentage of single waveforms classified as either type. Finally, only SUs showing a high enough percentage of single waveforms classified the same as the average waveform are sorted as the respective type. The minimal percentage required, referred to as consistency $c$, is initially set to the lowest value which ensures no contradictions between average- and single-waveform thresholding results. While the ‘‘broadness’’ and ‘‘narrowness’’ thresholds are chosen based on all available data for a given monkey, the required consistency is determined separately for each recording session. For monkey N $c$ is set to 0.6 in all but one sessions: In resting state session N1 it is increased to 0.62. For monkey E the values of $c$ equals 0.6 in the resting state recordings and take the following values in five analyzed reach-to-grasp sessions: 0.6, 0.89, 0.65, 0.61, 0.64.

The only step of our analysis for which the separation of putative excitatory and inhibitory neurons is crucial is the fitting of exponentials to the distance-resolved covariances. This step only involves resting state data. To test the robustness of our conclusions, we manipulate the required consistency value for sessions E1, E2, N1, and N2 by setting it to 0.75. Appendix 1—figure 2 and Appendix 1—table 1 summarize the resulting fits.

Appendix 1—figure 2

Download asset Open asset

It turns out that increasing $c$ to 0.75, which implies disregarding about 20-25 percent of all data, does not have a strong effect on the fitting results. The obtained decay constants are smaller than for a lower $c$ value, but they stay in a range about an order of magnitude larger than the anatomical connectivity. We furthermore see that fitting individual slopes to different populations in some sessions leads to unreliable results (cf. yellow lines in Appendix 1—figure 2A, I and blue lines in Appendix 1—figure 2C,D,K,L). Therefore, the data is not sufficient to detect differences in decay constants for different neuronal populations. Fitting instead a single decay constant yields trustworthy results (cf. yellow lines in Appendix 1—figure 2E,M and blue lines in Appendix 1—figure 2G,H,O,P). Our data thus clearly expose that decay constants of covariances are in the millimeter range.

The linear-response theory, with the aid of which we develop our predictions about the covariance structure in the network, assumes that the processes under examination are stationary in time. However, this assumption is not necessarily met in experimental data, especially in motor cortex during active behavioral tasks. For this reason we analyzed the stationarity of average single unit firing rate and pairwise zero time-lag covariance throughout a reach-to-grasp trial, similarly to Dahmen et al., 2019. Although the spiking activity becomes highly non-stationary during the movement, those epochs that are chosen for the analysis in our study (S and P) show only moderate variability in time (Appendix 1—figure 3). An analysis on the level of single-unit resolved activity also shows that the majority of neurons has stationary activity statistics within the relevant epochs S and P, especially when comparing to their whole dynamic range that is explored during movement transients towards the end of the task (Appendix 1—figure 5). Appendix 1—figure 6 shows that there are, however, a few exceptions (e.g. units 11, 84 in this session) that show moderate transients also within an epoch. Nevertheless, these transients are small compared to changes between the two epochs S and P.

Appendix 1—figure 3

Download asset Open asset

Appendix 1—figure 4

Download asset Open asset

Appendix 1—figure 5

Download asset Open asset

Appendix 1—figure 6

Download asset Open asset

Thus both the population-level and single-unit level analyses are in line with the second test for stationarity that we show in Figure 6. There we compare the firing rate and covariance changes between two 200 ms segments of the same epoch to the firing rate and covariance changes between two 200 ms segments of different epochs. If the neural activity was not stationary within an epoch then we would not obtain correlation coefficients of almost one between firing rates in Figure 6E and correlation coefficients up to 0.9 between covariance patterns within one epoch in Figure 6F. In summary, the analyses together make us confident that assuming stationarity within an epoch is a good approximation to show that there are significant behaviorally related changes in covariances across epochs of the reach-to-grasp experiment.

We are considering neuronal network models with isotropic and distance-dependent connection profiles. Ultimately, we are interested in describing cortical networks with two-dimensional sheet-like structure. But, for developing the theory, we first consider the simpler case of a one-dimensional ring and subsequently develop the theory on a two-dimensional torus, ensuring periodic boundary conditions in both cases. $N$ equidistantly distributed neurons form a grid on these manifolds. The position of neuron $i\in \{1,\mathrm{\dots },N\}$ is described by the vector ${\displaystyle {\mathit{r}}_i\in {\mathbb{R}}^D}$, with ${\displaystyle D\in \{1,2\}}$. The connections $W_{ij}$ from neuron $j$ to neuron i are drawn randomly with a connection probability that decays with distance between neurons $\left|{\mathit{R}}_i-{\mathit{R}}_j\right|$, described by the normalized connectivity profile $p(\mathit{R})$, $\int p(\mathit{R}){\mathrm{d}}^{D}r=1$, which we assume to obey radial symmetry. The connection probability decays on a characteristic length scale $d$. As we are working on discrete lattices, we introduce the probability of two neurons being connected $p_{ij}$, which is defined by the relation $p({\mathit{R}}_i-{\mathit{R}}_j)={lim}_{a\to 0}{p}_{ij}/a$, with lattice spacing $a$. We set the synaptic weights for connections of a single type to a fixed value $w$, but allow for multiple connections between neurons, that is $W_{ij}\in \{0,w,2w,\mathrm{\dots }\}=n_{ij}\cdot w$ for all sending neurons $j$ of a given type, where $n_{ij}$ is binomially distributed. Such multapses are required to simultaneously meet biological constraints on neuronal indegrees, neuron densities, and spatial ranges of connections. If instead one assumed Bernoulli connectivity, an analysis analogous to Eq. 7 of Senk et al., 2018 would yield a connection probability exceeding unity.

We introduce two populations of neurons, excitatory (E) and inhibitory (I) neurons. The number of neurons of a given population $a\in \{\mathrm{E},\mathrm{I}\}$ is $N_a$, and their ratio is $q=N_{\mathrm{E}}/N_{\mathrm{I}}$, which, for convenience, we assume to be an even number (see permutation symmetry below). The connection from population $b$ to population $a$ has the synaptic weight $w_{ab}$ and characteristic decay length of the connectivity profile $d_{ab}$. The average number of inputs drawn per neuron is fixed to $K_{ab}$. In order to preserve translation symmetry, $q$ excitatory neurons and one inhibitory neuron are put onto the same lattice point, as shown in Figure 3A in the main text.

Linear-response theory has been shown to faithfully capture the statistics of fluctuations in asynchronous irregular network states (Lindner et al., 2005). Here we follow Grytskyy et al., 2013, who show that different types of neuronal network models can be mapped to an Ornstein-Uhlenbeck process and that the low-frequency limit of this simple rate model describes spike count covariances of spiking models well (Tetzlaff et al., 2012). In particular, Dahmen et al., 2019 showed quantitative agreement of linear-response predictions for the statistics of spike-count covariances in leaky integrate-and-fire networks for the full range of spectral bounds $R\in [0,1)$. Therefore, we consider a network of linear rate neurons, whose activity $\mathit{x}\in {ℝ}^N$ is described by

with uncorrelated Gaussian white noise $\mathit{\xi }$, ${\displaystyle ⟨{\xi }_i(t)⟩=0}$, ${\displaystyle ⟨{\xi }_i(s){\xi }_j(t)⟩=D_i{\delta }_{ij}\delta (s-t)}$. The solution to this differential equation can be found by multiplying the whole equation with the left eigenvectors ${\mathit{u}}_{\alpha }$ of $\mathit{W}$

where $y_{\alpha }={\mathit{U}}_{\alpha }\cdot \mathit{x}$, ${\xi }_{\alpha }={\mathit{u}}_{\alpha }\cdot \mathit{\xi }$, and ${\lambda }_{\alpha }$ is denoting the corresponding eigenvalue of $\mathit{W}$. Neglecting the noise term, the solutions are given by

with Heaviside function $\mathrm{\Theta }(t)$. These are the eigenmodes of the linear system and they are linear combinations of the individual neuronal rates

Note that the weights ${\left(u_{\alpha }\right)}_i$ of these linear combinations depend on the details of the effective connectivity matrix $\mathit{W}$. The stability of an eigenmode is determined by the corresponding eigenvalue ${\lambda }_{\alpha }$. If $\mathrm{Re}\left({\lambda }_{\alpha }\right)\S lt;1$, the eigenmode is stable and decays exponentially. If $\mathrm{Re}\left({\lambda }_{\alpha }\right)\S gt;1$, the eigenmode is unstable and grows exponentially. If $\mathrm{Im}\left({\lambda }_{\alpha }\right)\ne 0,$ the eigenmode is oscillatory with an exponential envelope. $\mathrm{Re}\left({\lambda }_{\alpha }\right)=1$ is here referred to as the critical point. This type of stability is also called linear stability to stress that these considerations are only valid in the linear approximation. Realistic neurons have a saturation at high rates, which prevents activity from diverging indefinitely. A network is called linearly stable if all modes are stable. This is determined by the real part of the largest eigenvalue of $\mathit{W}$, called spectral bound $R.$ In inhibition-dominated networks, the spectral bound is determined by the heterogeneity in connections and $R⪅1$ defines the dynamically balanced critical state (Dahmen et al., 2019).

The different noise components ${\xi }_{\alpha }$ excite the corresponding eigenmodes of the system and act as a driving force. A noise vector $\mathit{\xi }$ that is not parallel to a single eigenvector ${\mathit{U}}_{\alpha }$ excites several eigenmodes, each with the corresponding strength ${\xi }_{\alpha }$.

Note that the different eigenmodes do not interact, which is why the total activity $\mathit{x}$ is given by a linear combination, or superposition, of the eigenmodes

where ${\mathit{V}}_{\alpha }$ denotes the α-th right eigenvector of the connectivity matrix $\mathit{W}$.

Time-lag integrated covariances $c_{ij}=\int D\tau ⟨x_i(t)x_j(t+\tau )⟩-⟨x_i(t)⟩⟨x_j(t+\tau )⟩$ can be computed analytically for the linear dynamics (Gallego et al., 2020). They follow from the connectivity $\mathit{W}$ and the noise strength $D$ as (Pernice et al., 2011; Trousdale et al., 2012; Grytskyy et al., 2013; Lindner et al., 2005)

with identity matrix 1. These covariances are equivalent to covariances of spike counts in large time windows, given by the zero-frequency component of the Fourier transform of ${\displaystyle \mathit{x}}$ (sometimes referred to as Wiener-Khinchin theorem Gardiner, 1985; even though the theorem proper applies in cases where the Fourier transforms of the signals ${\displaystyle \mathit{x}}$ do not exist). Spike count covariances (Figure 1B in the main text) can be computed from trial-resolved spiking data (Dahmen et al., 2019). This equivalence allow us to directly relate theoretical predictions for covariances to the experimentally observed ones.

While Equation (10) provides the full information on covariances between any two neurons in the network, this information is not available in the experimental data. Only a small subset of neuronal activities can be recorded such that inference of connectivity parameters from Equation 10 is unfeasible. We recently proposed in Dahmen et al., 2019 to instead consider the statistics of covariances as the basis for comparison between models and data. Using Equation 8 and Equation 10 as a starting point, field theoretical techniques allow the derivation of equations for the mean ${\displaystyle \overline{\mathit{c}}}$ and variance ${\displaystyle \overline{\delta {\mathit{c}}^2}}$ of cross-covariances in relation to the mean $\mathit{M}$ and variance $\mathit{S}$ of the connectivity matrix $\mathit{W}$ (Dahmen et al., 2019):

$\mathit{M}$and $\mathit{S}$ are defined in the subsequent section. The renormalized input noise strength is given by

with input noise covariance $\mathit{D}$, and the all-ones vector $\mathit{I}={(1,\mathrm{\dots },1)}^{\mathrm{T}}\in {ℝ}^N$. Note that Equation (12) only holds for cross-covariances ($i\ne j$). The diagonal terms ${\displaystyle {\left[\overline{\delta {\mathit{c}}^2}\right]}_{ii}}$, that is the variance of auto-covariances, do get a second contribution, which is negligible for the cross-covariances considered here.

For calculating the mean and variance of the covariances of the network activity ((11) and (12)) we need mean $\mathit{M}$ and variance $\mathit{S}$ of connectivity $\mathit{W}$. In the following, we derive the cumulant generating function (Gardiner, 1985) of $W_{ij}$.

The number of connections $n$ from neuron $j$ to neuron i is a binomial random variable with $K$ trials with the probability of success given by $p_{ij}$ (in the following, for brevity, we ignore the index i, $p_{ij}\equiv p_j$)

The average number of connections from neuron $j$ to neuron i is $K_j=p_{j}K$, which assures the correct average total indegree

The moment generating function of a connectivity matrix element $W_j\equiv W_{ij}\in \{0,w,2w,\mathrm{\dots }\}$ is given by

In a realistic network, $K$ is very large. In the limit $K\to \mathrm{\infty }$, while keeping $Kp=\mathrm{const}.$, the binomial distribution converges to a Poisson distribution and we can write

Taking the logarithm leads to the cumulant generating function

and the first two cumulants

Note that $\mathit{M}$ and $\mathit{S}$ have an identical structure determined by the connectivity profile and the structure of the covariance equation is identical for the mean Equation (11) and variance Equation (12) as well. This is why in the following we only derive the results for the mean of covariances. The results for the variance of covariances is obtained by substituting $w$ by $w^2$ and ${\mathit{D}}_{\mathrm{r}}$ by ${\mathit{D}}_{\mathrm{r}}^2$. As we show, divergences in expressions related to the mean covariances arise if the population eigenvalue ${\lambda }_0$ of the effective connectivity matrix approaches one. In expressions related to the variance of covariances, the divergences are caused by the squared spectral bound $R^2$ being close to one. In general expressions, we sometimes write $\zeta$ in order to denote either the population eigenvalue or the spectral bound, corresponding to the context of mean or variance of covariances.

For real neuronal networks, the anatomical connectivity is never known completely, let alone the effective connectivity. This is why we are considering disorder-averaged systems. They are described by the mean $\mathit{M}$ and variance $\mathit{S}$ of the connectivity. The latter inherit the underlying symmetries of the network, like for example the same radially symmetric connectivity profile for all neurons of one type. As neuronal networks are high dimensional systems, calculating covariances from Equation (11) and Equation (12) first seems like a daunting task. But, leveraging the aforementioned symmetries similarly as in Kriener et al., 2013 allows for an effective reduction of the dimensionality of the system, thereby rendering the problem manageable.

As a demonstrative example of how this is done, consider a random network of $N$ neurons on a one-dimensional ring, in which a neuron can form a connection with weight $w$ to any other neuron with probability p₀. In that case, $\mathit{M}$ is a homogeneous matrix, with all entries given by the same average connectivity weight

This corresponds to an all-to-all connected ring network. Due to the symmetry of the system, moving all neurons by one lattice constant does not change the system. The translation operator $\mathit{T}$, representing this operation mathematically, is defined via its effect on the vector of neuron activity $\mathit{x}$

Applying $\mathit{T}N$ ${\displaystyle N}$-times yields the identity operation

Hence, its eigenvalues are given by complex roots of one

with $L=Na$ denoting the circumference of the ring. This shows that $\mathit{T}$ has $N$ one-dimensional eigenspaces. Since the system is invariant under translation, $\mathit{M}$ is invariant under the transformation ${\displaystyle {\displaystyle \mathit{T}\mathit{M}{\mathit{T}}^{-1}=\mathit{M},}}$ and thus $\mathit{M}$ and $\mathit{T}$ commute. As $\mathit{M}$ leaves eigenspaces of $\mathit{T}$ invariant (if $\mathit{V}$ is an eigenvector of $\mathit{T}$, $\mathit{M}\mathit{V}$ is an eigenvector with the same eigenvalue, so they need to be multiples of each other), all eigenvectors of $\mathit{T}$ must be eigenvectors of $\mathit{M}$. Accordingly, knowing the eigenvectors of $\mathit{T}$ allows diagonalizing $\mathit{M}$. The normalized (left and right) eigenvectors of $\mathit{T}$ are given by

We get the eigenvalues of $\mathit{M}$ by multiplying it with the eigenvectors of$\mathit{T}$

which is always zero, except for $l=0$, which corresponds to the population eigenvalue ${\lambda }_0:=m_{k_0}=N{p}_{0}w$ of $\mathit{W}$ (Figure 3C in the main text). Now, we can simply write down the diagonalized form of $\mathit{M}$

and we effectively reduced the $N$-dimensional to a one dimensional problem. Inverting $\mathit{A}:=\mathbf{1}-\mathit{M}$ in Equation (11) is straightforward now, since it is diagonal in the new basis. Its eigenvalues can be written as ${\displaystyle {\displaystyle a_k=1-m_k,}}$ where we suppressed the index $l$. Therefore its inverse is given by

The renormalized noise can be evaluated using that the all-ones vector occurring in equation Equation (13) is the eigenvector ${\mathit{V}}_0$ of $\mathit{S}$. After identifying the eigenvalue s₀ with the squared spectral bound $R^2$, we find

which allows us to express the mean cross-covariances $\overline{c}$ (see Equation (11)) and the variance of cross-covariances ${\displaystyle \overline{\delta c^2}}$ (see Equation (12)) in terms of the eigenvectors of $\mathit{M}$ and $\mathit{S}$ respectively

The simplest network with spatial connectivity is a one-dimensional ring of neurons with one population of neurons. Following section Section 6, the mean connectivity matrix has the form

As $p_{ij}$ only depends on the distance of two neurons, the rows in $\mathit{M}$ are identical, but shifted by one index.