Spatial and temporal correlations in neural networks with structured connectivity

Correlated fluctuations in the activity of neural populations reflect the network's dynamics and connectivity. The temporal and spatial dimensions of neural correlations are interdependent. However, prior theoretical work mainly analyzed correlations in either spatial or temporal domains, oblivious to their interplay. We show that the network dynamics and connectivity jointly define the spatiotemporal profile of neural correlations. We derive analytical expressions for pairwise correlations in networks of binary units with spatially arranged connectivity in one and two dimensions. We find that spatial interactions among units generate multiple timescales in auto- and cross-correlations. Each timescale is associated with fluctuations at a particular spatial frequency, making a hierarchical contribution to the correlations. External inputs can modulate the correlation timescales when spatial interactions are nonlinear, and the modulation effect depends on the operating regime of network dynamics. These theoretical results open new ways to relate connectivity and dynamics in cortical networks via measurements of spatiotemporal neural correlations.


I. INTRODUCTION
Neocortical activity fluctuates endogenously on multiple spatial and temporal scales. These intrinsic fluctuations are usually quantified by correlations in neural activity. The spatial scale of correlations is measured by equal-time cross-correlations between spike counts in pairs of neurons [1]. The spatial correlations decrease with lateral distance between neurons in the cortex [2][3][4][5][6][7]. The temporal scale of correlations is measured by the decay rate of time-delayed auto-correlation of activity in single neurons and time-delayed cross-correlations between pairs of neurons. Timescales of spontaneous neural activity range widely from tens of milliseconds [8] up to several seconds [9] and increase from sensory to association and prefrontal cortical areas [8,10,11]. Spatial and temporal correlations of neural activity can be modulated during changes in behavioral states, such as selective attention [1,5,[12][13][14][15][16][17][18] or working memory maintenance [10], and relate to computations across different cognitive tasks [19][20][21]. Hence, understanding how neural correlations arise from the network connectivity and dynamics will help to identify mechanisms of neural computations in the brain.
Theoretical models suggest that spatial and temporal correlations in neural activity originate from the connectivity structure of biological circuits. In mammalian neocortex, the wiring of neural circuits is highly structured in space. Neurons in primate cortex are organized in minicolumns which consist of ∼80 − 100 vertically connected neurons spanning all cortical layers [22,23]. Minicolumns form local spatial clusters through short-range horizontal connections tiling the lateral dimension of the cortex [22]. The spatial organization of local intracortical connectivity is consistent with the dependence of cross-correlations on distance [4][5][6][24][25][26]. Similarly, network models suggest that differences in timescales across cortical areas may be directly related to areal differences in recurrent connectivity strength in the primate cortex [27].
These prior theoretical studies considered either spatial or temporal dimensions of neural correlations separately. In the spatial domain, network models with spatially arranged connectivity can produce spatial patterns of neural correlations with realistic distance dependence. This mechanism have been demonstrated in different types of network models including networks of spiking model neurons [4,24,25], binary units [18,26], and rate units [5,6]. In the temporal domain, theoretical studies of neural correlations focused primarily on randomly connected networks. Recurrent interactions in these networks generate slow timescales in autocorrelation, which can be significantly longer than the membrane time constant of individual neurons [28][29][30][31][32]. In these models, slow timescales can arise from operating in the transition to chaos regime [28][29][30] or from metastable transitions between finite randomly-connected clusters [32]. Moreover, a heterogeneous distribution of self-coupling strengths can generate heterogeneous timescales across network units [33,34].
Temporal and spatial correlations arise from the same spatiotemporal dynamics in the network and are therefore intertwined. However, prior theoretical work did not explore the relationship between correlations in these two domains, especially in networks with spatially arranged connectivity. Theoretical understanding of how the interplay between temporal and spatial correlations arises from the network's dynamics and connectivity will provide tight constraints on models of cortical dynamics.
We show that spatial and temporal correlations are tightly interdependent in networks with stochastic dy-namics and spatially arranged connectivity. We study analytically and in numerical simulations the spatiotemporal correlations in networks of binary units with connectivity arranged in one-and two-dimensional space. These networks generate rich spatiotemporal patterns of activity varying across multiple temporal and spatial scales. Using Fourier transformation, we show that each spatial frequency mode of fluctuations is related to a specific timescale, contributing hierarchically to the overall patterns of correlations. We find that both spatial and temporal scales of correlations depend on the spatial connectivity range, and timescales associated with nonzero spatial-frequency components are heterogeneously reduced with broader spatial connectivity. Moreover, the external input in networks with nonlinear interactions can modulate the fluctuations of mean activity and timescales of correlations with effects depending on the operating regime of network dynamics.
The organization of the paper is as follows. In Sec. II, we define the network models and derive general forms of dynamical equations for correlations. We use these equations to compute the spatiotemporal structure of correlations in one-(Sec. III) and two-dimensional (Sec. IV) models with spatially arranged connectivity. In Section V, we investigate how the input current modulates timescales of correlations in different operating regimes of network dynamics.

A. The network architecture
We consider networks of binary interacting units [35][36][37][38][39] with spatially structured connectivity. We study oneand two-dimensional networks with different ranges of spatial connectivity. In the one-dimensional model, N units are evenly spaced on a ring with the periodic boundary condition (Fig. 1a). In the two-dimensional model, N 2 units are evenly placed on the nodes of N × N square lattice with periodic boundary conditions (Fig. 1b). In both models, units receive directed connections from their neighbors within a ball of the radius R in Chebyshev distances (L ∞ norm) in one or two dimensions. In models with R = 1, each unit only receives inputs from its nearest neighbors, which we refer to as nearest neighbor connectivity. We refer to models with R > 1 as models with long-range connectivity.
In one-dimensional models, the perimeter of the ring is L, so the distance between neighboring nodes on the ring is a = L/N , where a denotes the lattice constant. Thus, the spatial position of unit i is x i = a · i, i = 0, ..., N − 1. For connectivity radius R, each target unit i receives directed connections from 2R nearby units ranging from i−R, i−R+1, ..., i+R (R = 1, 2, ..., N/2). The strength of connectivity is uniform across these 2R units and scaled by 1/R. When R = 1, each unit i only receives inputs from two nearest neighbors i − 1 and i + 1. The network architecture. (a) One-dimensional network with nearest neighbor (R = 1, left) and long-range connectivity (R > 1, right). (b) Two-dimensional network with nearest neighbor (R = 1, left) and long-range connectivity (R > 1, right).

B. Dynamics of binary units
In the network model, each unit i can be in one of two states S i ∈ {0, 1}. These states could represent the presence (S i = 1) or absence (S i = 0) of a spike in a time-bin for a single neuron or the high or low activity state in a local group of neurons such as a cortical minicolumn [40]. For simplicity, we call these states active (S i = 1) and inactive (S i = 0).
The state of each unit S i ∈ {0, 1} is updated based on transition rates between active-to-inactive and inactiveto-active states, given by ω(1 → 0) and ω(0 → 1), respectively. We parametrize the transition rates as: Here α 1 and α 2 are the intrinsic transition rates of one unit in the absence of network interactions. The second term in these expressions represents modulation of transition rates due to interactions with connected units in the network. In the interaction terms, j S j represents the sum of activity of units directly connected to unit i. We assume that all connected units uniformly contribute to the transition rate of the target unit (i.e. have uniform connectivity strength to the target unit). A nonlinear activation function F is a monotonically increasing function of x that satisfies conditions F(0) = 0, F(∞) = 1.
In previous models of binary-unit networks [35][36][37][38][39], F(x) is usually approximated by a Heaviside function with a fixed threshold. Here we consider F of the form: where θ is a positive constant that controls the gain of recurrent inputs, and n is the number of connected neighbors to each target unit. The parameters β 1 and β 2 control the interaction strength. To satisfy the condition of transition rate being positive, we require α 2 − β 2 0. In our models, we assume the connectivity is excitatory, hence β 1 > 0 and β 2 > 0. Thus, inputs from active neighbors will increase the transition rate from inactive to active state ω(0 → 1) and suppress the transition rate from active to inactive state ω(1 → 0). Since the connectivity is spatially organized, nearby units are more likely to become active simultaneously. Therefore, the recurrent interaction tends to enhance the spatial clustering of high activity states. Combining ω(0 → 1) and ω(1 → 0), the general expression of transition rate is given by

Linearized approximation
We study correlated patterns of activity fluctuations at the steady state. Outside the steady state regime, the activity is unstable and we cannot define activity fluctuations as perturbations around a fixed point. Hence, we focus on a region of model parameters in which after a sufficiently long time, the global network activity reaches an equilibrium state at a fixed point. In this case, linearization of the dynamical equations around the fixed point provides a good approximation and simplifies derivations of correlation functions. In our main analyses, we use the linear approximation of the interaction terms where we defined the effective interaction strengths Here we assumed that the mean global activity is close to zero, so F (S) ≈ F (0) and F(S) ≈ 0. With these conditions, the linearized transition rates become We discuss the case when mean activityS is nonnegligible in Sect. V.

Simulations of network dynamics
We verify our analytical derivations using numerical simulations of the network models (Fig. 2). We simulate the networks in discrete time using transition probabilities instead of transition rates. Specifically, the state of each unit is updated at each time step t k (k indexes time steps with t k − t k−1 = ∆t) based on the transition probabilities:

FIG. 2.
Simulations of a two-dimensional network model (N × N = 10, 000 units) with nearest-neighbor connectivity. (a) Snapshots of population activity. (b) Time-series of activity states for six example units sampled from two local neighborhoods indicated with circles in a.
Here p s is the self-excitation probability, p ext is the probability of external excitation, and p r is the probability of recurrent excitation from active neighbors. j S j denotes the number of active neighbors for the target unit. This transition probability scheme is a discrete-time approximation of continuous-time dynamics with the transition rates Eq. 7, linked via the parameter transformation (Appendix A): C. Dynamical equations for the mean activity and correlations We use the master equation to derive the dynamical equation for the mean activity and the general forms of time-evolution equations for the correlation functions (Appendix B) [35][36][37][38][39]. We assume β 1 = β 2 in all calculations and model simulations unless stated otherwise.
The mean activity S i (t) of unit i at time t (where · denotes averaging over the distribution of all possible configurations, Eq. B2) obeys the equation: with τ 0 is given by Eq. 12 shows that in the absence of network interactions (β 1 = 0), the activity of each unit drifts toward the same mean value α 1 /(α 1 + α 2 ) with the intrinsic timescale τ 0 . For finite interaction strength, we find the steady-state solution for the mean global activityS by averaging over all unitsS = i lim t→∞ S i (t) /N d , which yields where n r is the number of units connected to a target unit. With the connectivity radius R, n r = 2R for onedimensional models, and n r = (2R + 1) 2 − 1 for twodimensional models. Thus, the mean activity is scaled by a factor of 1/[1 − n r β 1 /(α 1 + α 2 )], which describes the effect of network interactions.
The value ofS sets the upper bound on the interaction strength, sinceS is a non-negative number, which implies n r β 1 /(α 1 + α 2 ) < 1. When the interaction strength exceeds this bound, the network activity becomes unstable and the mean-field approximation fails. We focus on the strong interaction regime which is close to the threshold of instability, i.e. n r β 1 /(α 1 + α 2 ) ≈ 1. For one-and twodimensional models with nearest-neighbor connectivity, the strong interaction limit is 2β 1 /(α 1 + α 2 ) ≈ 1 and 8β 1 /(α 1 + α 2 ) ≈ 1, respectively. In Sec. III C and IV C, we show that in this regime, spatial recurrent interactions generate slow timescales in auto-and cross-correlations that are much longer than intrinsic timescale τ 0 .
To compute neural correlations, we analyze the dynamics of fluctuations around the fixed point of the mean global activityS. We define the activity fluctuation of unit i as δS i = S i −S. The equal-time crosscorrelation function is then defined as δS i (t)δS j (t) (i = j), the time-delayed cross-correlation function is δS i (t)δS j (t + τ ) (i = j), and the auto-correlation function is δS i (t)δS i (t + τ ) .
The mean global activityS also determines the average variance of activity, which is the average auto-correlation at zero time lag: Using the property of binary units S i = S 2 i , we can express A(0) viaS: To obtain the analytical expressions for correlations, we used the general form time-evolution equations for correlation functions (Appendix B) derived based on the master equation formalism [26,[35][36][37][38][39]. We then applied Fourier expansion of these time-evolution equations to solve for the average equal-time and time-delayed crosscorrelations and autocorrelations. Fourier expansion was used in previous work but only to study equal-time crosscorrelations in one-dimensional binary network models [26] and firing-rate networks [6] with spatial connectivity. Next, we obtained the steady-state solution based on the Fourier transformation of time-evolution of equal-time cross-correlation function. Finally, we solved the timeevolution equation of time-delayed cross-correlations and autocorrelations, where initial conditions are given by the steady state of equal-time cross-correlations. In Sec. III and IV, we discuss the analytical solutions and numerical simulations of these correlations in different network configurations.

III. ONE-DIMENSIONAL MODELS
In this section we study spatiotemporal patterns of correlations in the one-dimensional models (Fig. 1a). We use Fourier transformation to derive the analytical form of the auto-and cross-correlations and study their dependence on the spatial network structure.

A. Correlation functions in Fourier space
The one-dimensional model contains N units. Each unit is located in position x, x = j · a, j = 0, ..., N − 1, with periodic boundary condition: x + N · a = x (Fig.  1a). We can expand the state S(x) of the unit at position x in Fourier space (here we omit the time index of S(x) for notation clarity): whereS(k) denotes the state variable in Fourier space with the wave number k. The periodic boundary condition requires the state variable to be invariant under translation S(x + N a) = S(x), which restricts the allowed values of the wave number in Fourier space: hence, e ikN a = 1, kN a = 2πm, m = 0, ±1, ±2, . . . . (18) Without loss of generality, we can define the Fourier mode spectrum to be N discrete values: k = 2πn/(N a) = 2πn/L, where n = 0, 1, ..., N − 1, which is analogous to the first Brillouin zone in solid state physics [41]. Similarly, we can expand the equal-time pairwise correlation function between the units located at x 1 and x 2 as (omitting the time index) e ik1x1 e ik2x2 δS(k 1 )δS(k 2 ) . (19) We are interested in the average correlation function C(x 1 , x 2 ) with a fixed difference between x 1 and x 2 : This correlation function can be expanded in Fourier space: whereC(k) is the amplitude of k-th Fourier mode of correlation function,C(k) = δS(k)δS(−k) = x C(x)e −ikx /N . Here, we focus on the case when the correlation function is symmetric C(x) = C(−x), in which case the correlation function can be expressed as a function of distance ∆ = |x 1 − x 2 |. The distance ∆ takes N/2 discrete values: ∆ = na, n = 1, 2, ..., N/2. In this case, the Fourier modes of correlation function are restricted to take N/2 values: k = 0, ..., (N/2 − 1): Without loss of generality, here and below we assume N to be an even integer. By using the identity C(∆) = C(−∆), we can rewrite the Eq. 22 as whereC(k) is the inverse Fourier transformation of C(∆): For the time-delayed cross-correlation function, we can also define the average correlation: and expand it in Fourier space as a function of distance using time-dependent Fourier amplitudesC(k, t): C(∆, t) has the initial condition C(∆, t = 0) ≡ C(∆), which givesC(k, t = 0) ≡C(k).
The average autocorrelation is defined as Since the average autocorrelation does not have spatial dependence, we do not directly apply Fourier transformation.

B. Time-evolution equations for the average correlation functions
We can simplify the general equations for correlations (Appendix B) to obtain the time-evolution equations for the average correlation functions C(x), C(x, t), A(t) defined in Sec. III A. By solving these equations in Fourier space, we obtain the Fourier amplitudesC(k) and then compute C(∆), C(∆, t) using Eqs. 23, 26. In the case of nearest-neighbor connectivity (R = 1), the steady-state equation for equal-time cross-correlation function reads The time evolution equation for the time-delayed crosscorrelation function is The time evolution equation for the time-delayed autocorrelation function is Here, we can see that the cross-correlations contribute to the auto-correlation function.
For the one-dimensional model with connectivity radius R > 1 (i.e. long-range connectivity), we denote the average equal-time cross-correlation with fixed position difference x as C(x; R), the average time-delayed crosscorrelation C(x, t; R), and the auto-correlation A(t; R). Similar to the case of nearest-neighbor connectivity, we obtain the steady-state equation for C(x; R): and the time-evolution equation for A(t; R): C(ma, t; R) . (33) C. Spatiotemporal structure of correlation functions

Nearest-neighbor connectivity
Here we study the spatiotemporal structure of correlation functions in the case of nearest-neighbor connectivity (Fig. 1a). Eq. 30 describes the time evolution of autocorrelation function. The first term on the right-hand side of the equation represents the decay of autocorrelation with the rate given by the intrinsic timescale τ 0 . In the limit of weak interactions β 1 → 0, we can neglect the contribution of cross-correlation to the auto-correlation and obtain the solution for autocorrelation as For finite interaction strength β 1 > 0, the crosscorrelation C(a, t) acts as an external source term that brings additional temporal structures into A(t). Therefore, A(t) contains two types of timescales: the intrinsic timescale τ 0 that is independent of network interactions, and the interaction timescales that are shared with crosscorrelation C(a, t).
To get the analytical form of C(a, t), we first solve Eq. 28 and get C(∆), which provides the initial condition for C(∆, t) at t = 0. Then, we can solve Eq. 29 to find C(∆, t). Eq. 28 and Eq. 29 are coupled equations for C(x) and C(x ± a), but they can be decoupled in Fourier Analytics Simulation space. Using Eq. 28, for each Fourier mode k, we obtaiñ

(35) Then,C(k) is given bỹ
In this expression, the factor 1/N comes from the normalization of the discrete Fourier transformation. The inverse Fourier transformation ofC(k) leads to C(∆): Relative weight, 2 Spatial frequency where we defined the correlation length L c : Our analytical calculation of C(∆) agrees well with the equal-time cross-correlation function computed from the model simulations (Fig. 3a). In the limit of strong interactions, 2β 1 /(α 1 + α 2 ) → 1, the correlation length can be approximated as Next, we compute the time-delayed cross-correlation. Eq. 29 includes both auto-and cross-correlations. The autocorrelation A(t) contains the intrinsic timescale τ 0 , which, as we will show, is faster than dominant timescales in the cross-correlation. Therefore, here we neglect A(t) on the right-hand side of Eq. 29. Under this approximation, the Fourier transformation of Eq. 29 is given by Here τ (k) is the interaction timescale for mode k (Fig.  4a) defined as Eq. 40 shows that each spatial Fourier modeC(k) fluctuates independently with the timescale τ (k): Thus, the time-dependence of C(∆, t) is described by a superposition of N/2 Fourier modes where each mode has a characteristic timescale τ (k) with the weight C(k) cos(k∆): Our analytical calculation of C(∆, t) agrees well with the results from numerical simulations (Fig. 3b). We find that C(∆, t) decay in time much slower than the intrinsic timescale τ 0 , indicating that interaction timescales are much longer than τ 0 . At short time lags, the decay rate of C(∆, t) decreases with increasing distance ∆ (seen as flattening profile of C(∆, t) at short time lags).
To understand the structure of weights for interaction timescales τ (k) in C(∆, t), we define the average timescale τ (∆) of the correlation function This expression shows that the average timescale of crosscorrelation τ (∆) is a weighted sum of N/2 timescales, where the normalized relative weight for mode k is given by [2C(k) cos(k∆)/C(∆)]. These weights define the relative contribution of different Fourier modes to the crosscorrelation.
These components lead to a slow decay of correlations at short time lags, flattening the temporal profile of correlations (Fig. 3b).
Using the analytical approximation of C(∆, t), we can solve Eq. 30 to obtain the analytical form of autocorrelation: This equation shows that A(t) contains N/2 + 1 timescales: the intrinsic timescale τ 0 (Eq. 13) and N/2 interaction timescales τ (k) (Eq. 41) inherited from the cross-correlation. The mixture of these timescales defines the temporal profile of autocorrelation. At short time lags, the decay of autocorrelation is dominated by the intrinsic timescale (Fig. 5a). At intermediate time lags, the autocorrelation decays with a characteristic timescale similar to τ (∆) which is between τ 0 and τ global (Fig. 5b).
In the limit of long time lags, the timescale of decay approaches the global timescale τ global (Fig. 5c). In addition, at time lags much larger than τ 0 , the autocorrelation and cross-correlations decay at a similar rate (Fig. 5d), confirming the effects of shared Fourier amplitudesC(k) in auto-and cross-correlations (Eq. 46).

Long-range connectivity
Here we study correlations in one-dimensional models with long-range connectivity (R > 1, Fig. 1a). We investigate how the connectivity radius R affects the spatiotemporal patterns of correlations. Same as for the nearest-neighbor connectivity, we solve the steady-state equation for cross-correlation Eq. 31 in Fourier space. The Fourier amplitudes of equal-time cross-correlatioñ C(k; R) are given bỹ (47) This equation shows that k = 0 mode is independent of R. For all other modes, the magnitude ofC(k; R) decreases with increasing connectivity radius R, especially for high-k modes (short-range correlations, Fig. 7). Thus, increasing R leads to more spatially homogeneous correlations (i.e. reduces the distance-dependence of correlations). This effect is evident in the position space, where equal-time cross-correlation C(∆; R) is given by with the correlation length L R To compute C(∆; R), here we used the approximation sin(Rka/2)/[sin(ka/2)R] ≈ 1 when R 1 to simplifỹ C(k; R). L c is the correlation length for the model with nearest-neighbor interactions (R = 1, Eq. 38). Eq. 49 shows that the correlation length L R is proportional to connectivity radius R. With increasing R, the network activity is more homogeneous, which is reflected in an increase of the correlation length. C(∆, R) estimated from the model simulations exhibits an increase in the correlation length (measured as the slope of C(∆, R) in the logarithmic-linear coordinates) with increasing R that is in agreement with the analytical prediction (Fig. 6).
To understand how the connectivity radius R affects the temporal structure of correlations, we solve the equation for the time-delayed cross-correlation (Eq. 32) in Fourier space, under the approximation of neglecting A(t; R) (similar to the case R = 1, Eq. 40). The timescale of each modeC(k, t; R) is determined by the equation: where the interaction timescales are . (51) Eq. 51 shows that τ (k; R) depends on the interaction radius R in heterogeneous manner, depending on the value of k (Fig. 7). Specifically, for k = 0, the global timescale is invariant to the change of R:  which means the timescale of global activity fluctuations is the same in all networks with different R. For finite-k modes, the associated timescales decrease with increasing R , and approach the intrinsic timescale τ 0 . In the limit of large R, the cross-correlation has only two nondegenerate timescales: τ 0 and τ global .
Next, we compute the time-delayed cross-correlation and auto-correlation functions for networks with longrange connectivity. Summing over all Fourier modes, the time-delayed cross-correlation in the position space is given by: .
(53) Combining Eq. 53 and Eq. 33, we obtain an analytical form of autocorrelation for networks with long-range connectivity: .
Similar to the case of nearest-neighbor interactions, A(t; R) contains N/2 + 1 timescales: the intrinsic time scale τ 0 (Eq. 13) and N/2 interaction timescales τ (k; R) (Eq. 51) inherited from cross-correlation C(∆, t; R). In A(t; R), the amplitude of τ (k; R) isC(k; R). Changing the connectivity radius R affects both the amplitudes and the corresponding timescales τ (k; R) ( Fig. 7), leading to R-dependent changes in autocorrelation. With increasing R, the relative weight of τ 0 is enhanced due to reduction ofC(k; R) for finite-k models. Accordingly, in autocorrelation A(t; R), the crossover from τ 0 to the average interaction timescale (mixture of τ (k; R)) occurs at a larger time lag t (Fig. 8a). Since in the large time-lag limit, the autocorrelation is dominated by the largest interaction timescale τ global , the autocorrelations decay at the same rate for all values of R in this region (Fig. 8b).
In summary, the connectivity radius R affects both the spatial and temporal structure of correlation functions. Increasing R diminishes the non-zero spatial-frequency components in equal-time cross-correlation and suppresses the amplitude of interaction timescales (except for the global timescale) in the autocorrelation. In the large R limit, cross-correlations become spatially homogeneous, and autocorrelations contain only two residual timescales, the intrinsic timescale and the global timescale.

IV. TWO-DIMENSIONAL MODELS
In this section, we generalize the analytical methods used for the one-dimensional models to study the spatiotemporal correlations in the two-dimensional models (Fig. 1b). Similar to the one-dimensional model, we can expand the state of each unit in Fourier space. We denote the location of units on the lattice as (x 1 , x 2 ), where x 1 = n 1 a, x 2 = n 2 a, n 1,2 = 0, ..., N − 1. The periodic boundary conditions are x 1 + N a = x 1 , x 2 + N a = x 2 . Similar to Sect. III A, the periodic boundary conditions lead to discrete modes in Fourier space: k 1 = 2πm 1 /(N a) = 2πm 1 /L, k 2 = 2πm 2 /(N a) = 2πm 2 /L, where m 1,2 = 0, ..., N − 1. Then, the activity state of the unit at (x 1 , x 2 ) is
The average autocorrelation in two-dimensional models is defined as x,y δS(x, y, t 0 )δS(x, y, t 0 + t) /N 2 . (63) B. Spatial structure of correlations Next, we study the dependence of cross-correlation on the spatial distance. In the case of nearest-neighbor interactions, we solve the steady state equation for equal-time cross-correlation in Fourier space (Appendix C) and find the amplitude of each spatial Fourier mode (k 1 , k 2 ) as An inverse Fourier transformation ofC 2 (k 1 , k 2 ) gives rise to C 2 (∆ 1 , ∆ 2 ): Here L c,2 is defined as the correlation length in two dimensions: In the limit of strong interaction, 8β 1 /(α 1 + α 2 ) → 1, the correlation length can be approximated as Equal-time cross-correlations decay exponentially with increasing distance. Fig. 10a shows C 2 (∆ 1 , ∆ 2 ) as a function of distance for ∆ 1 = ∆ 2 , where L c,2 is given by Eq. 66. The spatial profile of the average correlation C 2 (∆) can be obtained by combining Eqs. 60 and 64 (Fig.  10b).

C. Timescales of correlations
Here we explore timescales of auto-and crosscorrelations in the two-dimensional model with the nearest-neighbor connectivity.
To quantify how the average temporal profile of crosscorrelations depends on distance, we compute the average interaction timescale for cross-correlation C(∆ 1 , ∆ 2 , t): As we can see from numerical values of τ (∆ 1 , ∆ 2 ) (Fig. 13c), when ∆ increases from the minimal distance ∆ = a, the average interaction timescale increases and reaches a peak at ∆ = 7a. When ∆ increases further, the average timescale decreases and approaches a value close to τ global,2D , because at large distances ∆, C 2 (∆, ∆, t) is dominated by the homogeneous (distance-independent) component, which has a global timescale τ global,2D .
The average time-delayed correlation C(∆, t) can also be written as a summation ofC 2 (k 1 , k 2 , t): Here we defined M 2 (∆; k 1 , k 2 ) to be the weight of each mode (k 1 , k 2 ) in C 2 (∆, t). The patterns of M 2 in k-space are shown in Fig. 14 for different distances ∆. Qualitatively, M 2 has a similar behavior as M 1 . When ∆ = a, M 2 is always positive, creating a super-linear correlation C 2 (∆, t) in the logarithmic scale. When ∆ > a, M 2 oscillates between negative and positive values in the small-k region, generating a slow time decay (plateau) of C 2 (∆, t) at short time lags (Fig. 13b).

Autocorrelation
We obtain the analytical form of autocorrelation A 2 (t) by solving the time-evolution equation (Appendix Eq. C3). In the limit of weak interactions β 1 /(α 1 +α 2 ) → 0, A 2 (t) is an exponential decay function with time constant τ 0 : A 2 (t) = A(0) exp(−t/τ 0 ). For finite interaction strength, time dependence of A 2 (t) is influenced by the cross-correlation terms C 2 (a, a, t), C 2 (a, 0, t) and C 2 (0, a, t). Therefore, A 2 (t) inherits N 2 /4 interaction timescales τ (k 1 , k 2 ) (Eq. 69) from the time-delayed crosscorrelation. Altogether, A 2 (t) contains N 2 /4 interaction timescales and an intrinsic timescale τ 0 , with the following analytical expression: The temporal decay pattern of A 2 (t) is dominated by different timescales at different ranges of time lags. At short time-lags, A 2 (t) decays with the intrinsic timescale τ 0 (Fig. 15a). At intermediate time lags, A 2 (t) decays with an intermediate timescale that is in between τ 0 and τ global,2D (Fig. 15a). This intermediate timescale comes from a superposition of all interaction timescales τ (k 1 , k 2 ) and is similar to the timescales of cross-correlations C 2 (∆ 1 , ∆ 2 , t) (Fig. 15c), which reflects the link between auto-and cross-correlations. In the limit of large time lags, the time decay of A 2 (t) is dominated by the largest timescale τ global,2D and contributions of all other timescales are negligible (Fig. 15b).
To quantify the average temporal profile of the interaction part of auto-correlation and its relation to the average timescales of cross-correlations, we define the average interaction timescale τ AC as a time integral of autocorrelation after subtracting the component associated with the intrinsic timescale: Here C 2 (0, 0) = 4 2π(N/2−1) L k1,k2=0C 2 (k 1 , k 2 ). As we can see from numerical values of τ AC and τ (∆ 1 , ∆ 2 ) (Fig. 15d), the average interaction timescale of autocorrelation is similar to the average interaction timescale of cross-correlation C 2 (∆ 1 = a, ∆ 2 = a, t), indicating the link between auto-and cross-correlations.
The maximal value f (k 1 = 0, k 2 = 0; R) ≡ 1 does not depend on R, whereas the band width scales with 1/R. Hence increasing R acts to reduce the number of k modes that contribute to f (k 1 , k 2 ; R). The dependence of f (k 1 , k 2 ; R) on R is reflected inC(k 1 , k 2 ; R). The amplitude of zero-k modeC(0, 0; R) does not depend on R, and non-negligible values ofC(k 1 , k 2 ; R) are restricted to the region k 1,2 ∈ [0, π/(Ra)] (Fig. 16).
As spatial scale of correlation C 2 (∆ 1 , ∆ 2 ; R) scales approximately with the inverse of Fourier wave number k, the correlation length should scale approximately with R. Indeed, we find that f (k 1 , k 2 ; R) ≈ cos((R + 1)k 1 a/2) cos((R + 1)k 2 a/2) in the limit of R 1, and the equal-time cross-correlation is then where the correlation length is proportional to R: Numerical values of C 2 (∆ 1 , ∆ 2 ; R) for intermediate R confirm that the decay rate of correlations with distance (inverse of the correlation length) decreases with increasing R (Fig. 17), which indicates the diminishing amplitudes of high wave-number modes. When R reaches the maximal value R = N/2, only the zero-k mode has a non-zero amplitude, hence C 2 (∆ 1 , ∆ 2 ; R = N/2) becomes a homogeneous (distance independent) function.
In summary, the increase of interaction-radius R smooths the spatial profile of the equal-time cross-correlation by reducing the amplitudes of all non-zero wave number modes.
To understand how temporal patterns of correlations depend on the connectivity range, we solve the timeevolution equation for the time-delayed cross-correlation C 2 (∆ 1 , ∆ 2 , t; R). We solve this equation in Fourier space with the approximation of neglecting A 2 (t; R) terms (Appendix C). We find that each modeC 2 (k 1 , k 2 , t; R) is an exponential decay function of time-lag t with an interaction timescale Then, the time-delayed cross-correlation can be written as a weighted sum of N 2 /4 modes, where each mode carries an interaction timescale τ (k 1 , k 2 ; R): Eq. 83 shows that the magnitude of timescales τ (k 1 , k 2 ; R) depend on R. In particular, the global timescale associated with k 1 = k 2 = 0 mode does not depend on R: τ (k 1 = 0, k 2 = 0; R)=τ glabal,2D . All other interaction timescales decrease with the increasing R and are pushed towards the value of the intrinsic timescale τ 0 (Fig. 18). Autocorrelation A 2 (t; R) is given by the combination of a component with the intrinsic timescale τ 0 and N 2 /4 components inherited from the cross-correlation modes with interaction timescales τ (k 1 , k 2 ; R): The temporal profile of autocorrelation is influenced by the relative weights of intrinsic timescale τ 0 and the interaction timescales. At short time lags t ≈ τ 0 , A 2 (t; R) decays with the timescale τ 0 . At intermediate time lags, A 2 (t; R) decays with an intermediate timescale which reflects the cumulative effect of all interaction timescales. In between these two regions, the autocorrelation slope (in the logarithmic-linear coordinates) changes abruptly indicating a crossover from decay rate dominated by τ 0 to the intermediate timescales.
Since the amplitudes C 2 (k 1 , k 2 ; R) of interaction timescales decrease with increasing R (except the zero-k mode), the time lag where the crossover occurs increases monotonically with R (Fig. 19a). At large time lags, the overall decay rate is governed by τ global,2D . Since the amplitude of zero-k mode with the global timescale τ global,2D is independent of R (Fig. 16), autocorrelations of models with different R exhibit the same slope at large time lags, with different intercepts reflecting the R-dependence of components associated with other interaction timescales (Fig. 19b).

V. OPERATING REGIME OF NETWORK DYNAMICS AND TIMESCALES OF CORRELATIONS
In previous sections, we focused on the case where the mean global activityS was very close to zero. Here we discuss how the mean global activity affects correlations. We show that increasing mean global activity can increase or decrease the intrinsic timescale of correlations, depending on the sign of β 1 − β 2 . Previous studies of binary neuron models [35][36][37][38][39] analyzed only the special case β 1 = β 2 , in which the intrinsic timescale does not depend on the mean global activity. The mean global activity also affects effective interaction strengths β 1,2 and therefore influences the interaction timescales of correlations. We also show how the external input affects the magnitude of the mean global activity, its stability, and timescales of correlations in different operating regimes of network dynamics.

A.
The mean global activity and the intrinsic timescale of correlation The mean global activity modulates the transition rates and therefore can affect the intrinsic timescale. In the derivation of transition rates ω(0 → 1) and ω(1 → 0) (Eq. 1), we can perform Taylor expansion around the mean global activityS. The expansion for the interaction terms are given by where F 0 is defined as Here F denotes the derivative of F. Since the activation function is the explicit forms of F 0 and dF 0 /dS are The interaction strengths in the linearized approximation are When the mean global activityS 1, we can neglect F 0 term and replace F (nS) in the expressions for the interaction strengths by F (0). However, when the mean global activityS is of order one, we have to include the contribution from F 0 as well as modulations of interaction strengths due to F (nS). In this case we can rewrite the transition rates as Thus, the effective intrinsic transition rates are activity dependent: With these effective intrinsic transition rates, the intrinsic timescale also becomes activity dependent. Based on Eq. 13, the equation for the intrinsic timescale can be rewritten as According to this equation, increasing mean global activ-ityS leads to a decrease of τ 0 when (β 1 − β 2 ) > 0 and to an increase of τ 0 when (β 1 − β 2 ) < 0. The changes of the intrinsic timescale result from a non-linear activation function and large values of the mean global activitȳ S. In the linear networks, F 0 is zero and the intrinsic timescale is constant.

B. Influence of external input on the operating regime of network dynamics and interaction timescales
In models with non-linear interactions, the linear response of the system (derivative of the activation function) depends on the mean global activity and inputs. We show that for fixed interaction strength, the external input current changes the operating regime of network dynamics, which affects the magnitude of the mean global activity, its stability, and the intrinsic and interaction timescales.
The activation function with the external input is defined as: where I represents a constant global input current (here we only consider the case I 0). In the steady-state, the mean global activityS follows the equation: Since 0 < F 1 and F(∞) = 1, in the large input limit I → ∞, we haveS(I → ∞) = α1+β 1 α1+α2+(β 1 −β 2 ) . Eq. 99 can have different solutions forS depending on the sign of (β 1 − β 2 ). To find these solutions, we define the function g(S; I): g(S; I) = α 1 + β 1 F(nS + I) α 1 + α 2 + (β 1 − β 2 )F(nS + I) . (100) The solutions of Eq. 99 are the intersections between the curve x =S, y = g(S; I) and the straight line x =S, y =S in the (x, y) plane (Fig. 20). The number and locations of the intersections depend on the first and second derivatives of g(S; I). The first derivative of g(S; I) is and the second derivative is To determine the stability of solutions forS, we consider a small deviation δS around the solutionS = g(S; I).
The magnitude of fluctuation of the mean global activity is equal to g (nS + I)δS.
Using the solution forS, we can determine the effect of external input on the intrinsic (Eq. 97) and interaction timescales. For simplicity, we consider only a representative interaction timescale, the global timescale, which is the largest interaction timescales. For a given mean global activityS and external input I, the global timescale is To study the dependence of τ global onS, we take the derivative of the denominator with respect toS: Since both timescales τ 0 and τ global depend onS, the external input can affect the timescales by changingS. If F is a linear function ofS, then both F 0 and F are independent ofS and hence the input does not influence the timescales [42]. Depending on the parameters α 1,2 , β 1,2 , θ, there are two classes of solutions forS depending on the sign of (β 1 − β 2 ). In the following, we discuss these possible solutions and how they affect the intrinsic and global timescales.
In this configuration, there is only one solution forS, in the range from α 1 /(α 1 + α 2 ) toS(I → ∞) (Fig. 20a). With increasing current I, the global activityS increases, leading to a reduction in both the intrinsic timescale τ 0 and the global interaction timescale τ global . b) One solution forS in the sub-linear or supra-linear region (Fig. 20c,d). When conditions in case (a) are not satisfied, there is always one solution forS in sub-linear or supra-linear region.

VI. DISCUSSION
We studied the spatial and temporal scales of neural correlations in binary-unit networks with connectivity arranged in one-and two-dimensional space. We used the time-evolution equations for correlation functions derived from the master equation. We solved these equations using the discrete Fourier transform and translational symmetry of the model and obtained analytical solutions for spatiotemporal correlations. We found that the spatial and temporal scales of correlations are related to each other and shaped by the spatial profile of the recurrent connectivity. Finally, we showed that external inputs can control the operating regime of the global network activity and thus influence the timescales of correlations. To confirm our theoretical results, we performed numerical simulations and found a good agreement between analytical solutions and simulation results.
One of our key findings is that spatial recurrent interactions generate multiple timescales in network dynamics. The spatial interactions we considered are similar to the spatial connectivity structure in the primate cortex. The distance-dependent connectivity perseveres the translational symmetry, hence in Fourier space, each spatial Fourier mode of correlations is approximately decoupled and evolves with a unique characteristic interaction timescale. In the strong interaction limit, the interaction timescales can be significantly larger than the intrinsic timescale. The overall temporal profile of correlations arises from a superposition of all Fourier modes with distinct timescales. These interaction timescales depend on the spatial range of connectivity in heterogeneous manner. In particular, local spatial connectivity tends to enhance a broad spectrum of interaction timescales, while homogeneous all-to-all connectivity eliminates all interaction timescales except for the global timescale associated with the spatially homogeneous component of correlations. Therefore, in our network models, multiple timescales are inherently coupled to the spatial connectivity, which is different from other models where heterogeneous timescales are generated by single-cell proprieties such as self-couplings [33]. The relation between timescales and connectivity has been analyzed in a deterministic linear network model [34], where timescales are defined by the eigenvalues of the connectivity matrix. Here we study the relation between structural connectivity and timescales of neural correlations in stochastic networks of binary-units.
Another major contribution of our work is to establish the link between spatial and temporal scales of correlations. Our theory predicts that slow interaction timescales in autocorrelations of networks with spatial connectivity are generated by correlations between the activity of units at different distances. In these networks, correlations at different distance have distinct amplitude spectra of their spatial Fourier modes. Since each Fourier mode carries a unique interaction timescale, the overall temporal structure of correlations depends on the the spatial distance. In particular, the average interaction timescale tends to be larger for correlations between pairs of neurons with a larger distance. This feature is supported by a recent analysis of spiking activity in primate visual cortex [18].
We showed that when the interaction between the network units is non-linear, the external input current changes the operating regime of network dynamics and modifies the intrinsic and interaction timescales of correlations. This mechanism of modulating timescales through external input may have implications in biological circuits. For example, in neocortex top-down inputs from higher cortical areas can regulate dynamics of cortical states in sensory areas [12], which may modulate timescales of fluctuations [5,18,40].
In this paper, we considered models with spatial connectivity patterns where each unit connects to all its neighbors within a radius R. In the future, it would be interesting to extend the current framework to study network models with random spatial connectivity. In this case, the connectivity patterns can be described by random band matrices. According to theories of random band matrices [43], the spatial correlation undergoes a transition from localization to delocalization phases, when the range of spatial connectivity exceeds certain thresholds. In addition, we focused here on the regime with a stable activity. Another future extension is to explore the spatiotemporal correlations in the dynamical regime where rate dynamics are chaotic [28,44]. we discuss how these two different representations of dynamics are related to each other.
In discrete-time dynamics, the state of a binary unit S i ∈ {0, 1} is updated at each time step ∆t based on the transition probabilities, which depend on the sum of states of its directly connected neighbours (denoted by j S j ): Here, we define the interaction term as a linear function Generally, F should satisfy the condition F(0) = 0, F(∞) = 1, and F(x) is a monotonically increasing function of x. When the mean global activity of the network is much smaller than 1, the above linear definition serves as a good approximation. Thus, for each unit we define a transition matrix between binary states: .
Using equations A1-A4, we can write the transition matrix as: In our analytical calculations of binary-unit dynamics, we use the instantaneous transition rates α 1,2 and β 1,2 to describe the changes in the probability density of the states. To link transition rate parameters to transition probabilities (p ext , p s , p r ), we use the fact that the transition matrix P (∆t) can be approximated by the matrix exponential of transition rate matrix e Q∆t , where the transition rate matrix Q is given by Here, ω(0 → 1) describes the transition rate from state 0 to 1, and ω(1 → 0) describes the transition rate from state 1 to 0. Then, the matrix exponential of transition rate matrix can be written as Solving e Q∆t = P (∆t) and using Eq. A1-A5, we have where, transition rates α 1 , β 1 are given by: and where transition rates α 2 , β 2 are given by We see that for the transition rates corresponding to the parameters of discrete model, β 1 = β 2 .
We denote the probability of the network to be in a certain configuration {S} = {S 1 , S 2 , ..., S N } at time t by P ({S}, t). The master equation describes the time evolution of P ({S}, t), which is given by [35][36][37][38][39]: where {S} i * = {S 1 , S 2 , ..., 1−S i , ..., S N } and w(S i ) is the transition rate from state S i to 1 − S i . Using the master equation, one can write down the equation for the time evolution of arbitrary moments. For example, the average activity of a unit i is defined as where we sum over all configurations of variables {S} at a given time t. The time evolution of the average activity is given by Similarly, the rate of change of a second moment for each pair of units is The time evolution of time-delayed second moment can be computed as [35][36][37][38][39]: where P ({σ}, t + τ |{S}, t) is conditional probability of finding the system in configuration {σ} at time t + τ , given that it was in configuration {S} at time t. Since the conditional probability obeys the same master equation, we have d dτ S i (t)S j (t + τ ) = S i (t)(1 − 2S j (t + τ ))w(S j (t + τ )) .
[2] M. A. Smith and A. Kohn, Spatial and Temporal Scales