As a starting point, we designed a simple linear model where spike trains from individual neurons are generated by stationary Poisson processes, implying that the number of spikes in non-overlapping intervals is conditionally independent and that the spike rate remains constant over time (Figure 1A) (see section “Materials and methods”). The spike probability p(S_i(t)) of a neuron i at time t depends on the baseline rate (r₀) and the spiking activity of surrounding neurons,

where S_j(t−1) is the spiking activity of neuron j at the previous time-step, α is a fixed coupling strength, and J_ij = 1 if neuron j is connected to neuron i, and 0 otherwise, where the connection probability follows an exponential function with a decay parameter (λ_decay) (Figure 1C),

where d_ij denotes the Euclidean distance between neurons i and j, expressed in units of millimeters. In this model, neurons are recurrently connected and embedded in two-dimensional space, meaning that each neuron is assigned a unique position along a finite two-dimensional grid where boundaries do not wrap around. Hence, spiking activity in this model is described by both spatial and temporal coordinates (Figure 1B).

An example of a spike raster generated with this model using λ_decay = 5 mm^–1 is shown in Figure 2A (1 time-step set as equivalent to 1 ms). The central question addressed here is whether we can reliably estimate λ_decay based on this spike raster. This was performed by maximum likelihood estimation (see section “Materials and methods”). Because λ_decay is estimated from spike trains and not synaptic connections, it captures the functional interactions between spatially embedded neurons, and is not a proxy for physical connectivity, though the two may be linked, as exemplified below.

A numerical example of the log likelihood obtained from 1,000 Poisson neurons simulated for 100 s of activity is shown in Figure 2B. The log likelihood was maximal near the true value of λ_decay = 5 mm^–1. A similar result held across a range of λ_decay values (Figure 2C). Hence, MLE yielded a close approximation of the λ_decay parameter, providing a characterization of distance-dependent interactions based on the neural data. MLE also approximated distance-dependent interactions characterized by half-Gaussian, linear, and lognormal functions, albeit poorer results were obtained with an inverse square (Supplementary Figure 1).

To evaluate the extent to which spike data from the Poisson model could be used to recover underlying distance-dependent functions, we applied a spatial Generalized Linear Model (GLM, see section “Materials and methods”) (Pillow et al., 2008). Five different generating functions were employed to produce spike trains (exponential, half-Gaussian, linear, inverse square, and lognormal). The maximum likelihood of each generating function was computed to identify the most likely candidate function. Results show that the exponential, half-Gaussian, linear, and inverse square functions could be identified based on their maximum log likelihood (Table 1). The lognormal function, however, was confounded with the linear function and was therefore excluded from further analyses. Thus, maximum likelihood estimation based on a spatial GLM enabled the identification of several distance-dependent functions that determined spatial interactions in the Poisson model.

Simulation results obtained with MLE were compared to three related approaches. First, functional connectivity between all pairs of simulated neurons was computed using transfer entropy (Kobayashi and Shinomoto, 2024; Vicente et al., 2011). A least-squares approximation was employed to find the best-fitting value of λ_decay relating transfer entropy to spatial distances between pairs of neurons. This approach offered a poor fit to the true value of λ_decay from which spike trains were generated (Figure 2C).

In a second approach, spike-count correlations were computed between all pairs of neurons using non-overlapping time bins (10 ms), then least-squares was applied to find the λ_decay that best describes the relation between correlations and spatial distances. This approach offered a marginally better fit than transfer entropy, but its accuracy remained lower than MLE. Third, Granger causality (Tian et al., 2024) was computed followed by least-squares fitting, yielding a poor fit between neural interactions and spatial distances. A fourth approach employed the spatial GLM to estimate values of λ_decay, again yielding a poor fit compared to the MLE. Examples of fit obtained from MLE, transfer entropy, Granger causality, spatial GLM and spike-count correlations are shown in Figure 2D. In summary, the fit obtained with MLE was not attained with the alternative approaches.

Next, we examined the computational run time required for MLE, transfer entropy, Granger causality, spatial GLM, and spike-count correlations. Run times were computed on an Intel i9 processor at 2500 Mhz. For small networks consisting of 100–500 neurons, run times were similar across the five approaches (Figure 3A). However, as the number of neurons increased, the run time of transfer entropy and spatial GLM increased rapidly, whereas the runtimes of MLE, Granger causality, and spike-count correlations remained comparatively low. Overall, MLE required only marginally higher run time than spike-count correlations.

In follow-up analyses, we compared the various approaches in terms of estimation error (computed as the absolute difference between the true and estimated values of λ_decay). We fixed λ_decay = 5 mm^–1 and systematically increased the number of simulated neurons. MLE attained markedly lower error compared to related approaches (Figure 3B). While the spatial GLM achieved similar performance to MLE with >1,500 neurons, this came at the cost of a high CPU time. Thus, MLE provided an accurate estimation of spatial interactions based on spiking activity and required conservative run times compared to alternative approaches.

Additional analyses revealed that while several factors impacted MLE accuracy, including firing rates, subsampling, and synaptic delays, error remained low under typical experimental conditions (Supplementary Figures 2–4).

Going beyond a linear Poisson model, we examined a network of Izhikevich neurons with distance-dependent connectivity (Figure 4A). The model generated sparse, uncorrelated activity (Figure 4B). Importantly, because the MLE approach assumes that spike trains are close to Poisson statistics, we computed the Fano factor (FF) – the ratio of spike count variance to its mean – for all neurons and ensured a distribution with mean FF ∼ 1 (Figure 4C). In addition, the mean spike-count correlation across all pairs of neurons was close to zero (mean: 0.0017) (Figure 4D).

To test the proposed MLE approach, a series of simulations were performed where a population of 1,000 neurons generated spike trains for 100 s with varying values of λ_decay. Examples of fit between the true and estimated λ_decay are shown in Figure 5A. Close correspondence between these values was found over a range of parameters (Figure 5B).

Next, the Izhikevich model was employed to investigate the impact of distance-dependent interactions on neural dynamics. Broadening interactions by increasing λ_decay from 1 mm^–1 to 5 mm^–1 (Figure 6A) resulted in a decrease in synchronized activity (Figure 6B). This makes intuitive sense given that longer-range connectivity promotes communication between a broader network of neurons. In turn, broadening spatial interactions by lowering λ_decay in a range between 0.5 and 3.0 resulted in a high average spike-count correlations (Figure 7A).

Interestingly, an abrupt transition was observed between a state of high and low correlations, occurring around λ_decay = 3.0. To formally examine this transition, we considered an approximation for the correlation C_ij between pairs of neurons receiving a shared input under the assumption of weak coupling (see section “Materials and Methods”),

where J is the mean synaptic strenght and σ² is the noise variance of the shared input. Plotting Equation 2 across values of λ_decay shows a non-linear relation between spatially dependent connectivity and pairwise correlations (Figure 7B). Thus, particularly for stronger coupling strength J, correlations remained high for low values of spatial decay, with a sharp decline as decay increased further.

We further examined this effect in a scenario with randomized connectivity by shuffling the distance-dependent connections of the Izhikevich model. A similar transition in correlation was observed in the random model (Figure 7C). However, in the distance-dependent model, shorter-range correlations (< 0.1 mm) retained higher values than longer-range connections (> 1 mm) following the transition. This effect did not arise in the random model due to the lack of spatial dependencies (Figure 7D). Thus, while an abrupt shift in correlation may be explained by synaptic density rather than spatial dependency (Roxin et al., 2004), the phase transition exhibited a distance-dependent effect that was not present in the randomized model.

We applied MLE to estimate spatial interactions in two-photon calcium imaging data (Stringer et al., 2019a). These data are particularly challenging given that they are composed of approximately 10,000 neurons, making approaches based on pairwise computations prohibitive. Instead, the approach taken by MLE is to estimate the best-fitting λ_decay on the basis of the spike raster and the physical distance amongst neurons, without needing to compute the strength of interactions between all pairs of neurons. Thus, this approach differs from attempts at reconstructing connectivity based on calcium data (Orlandi et al., 2014; Stetter et al., 2012). In our case, the goal is strictly focused on estimating the decay parameter of a distance-based function, without reconstructing the full connectivity.

A total of 14 datasets were analyzed, including 9 rasters of spontaneous activity and 5 rasters of evoked activity with drifting gratings (see section “Materials and methods” for details). The mean number of neurons was 11,602. A sample of neurons arrayed according to their three-dimensional coordinates is shown in Figure 8A, showing a spatially homogeneous distribution of cell locations. A raster of spiking activity for this sample shows a weak degree of coordination across neurons (Figure 8B) with a Fano factor near FF ∼ 1, indicative of approximate Poisson statistics (Figure 8C). Mean firing rates were similar for spontaneous (2.76 Hz) and evoked (2.87 Hz) activity (Figure 8D). Mean spike-count correlations (100 ms non-overlapping time bins) were low for both spontaneous (C = 0.0072) and evoked (C = 0.0121) activity.

An example of a log-likelihood function obtained from MLE applied to a single raster is shown in Figure 9A. The function shows a defined peak at λ_decay = 3.53 mm^–1. Hence, the characteristic length scale of spatial interactions is 3.53^–1 = 0.28 mm, which means that the probability of pairwise interactions decays significantly over approximately 0.28 mm (280 μm). To illustrate this result, we considered a two-dimensional neuronal space and assumed a putative neuron located in the center of this space. We then plotted an exponentially-distributed probability of pairwise interactions according to the function p(d) = ce^−3.53d (Figure 9B). This function shows that most of the interactions are contained within a delimited spatial area surrounding each neuron, in line with work showing that a majority of cortical interactions occur within a 100–300 μm radius (Ercsey-Ravasz et al., 2013; Holmgren et al., 2003; Horvát et al., 2016; Levy and Reyes, 2012; Markov et al., 2013). This result is consistent with distance-dependent interactions observed in related measures including spike-count correlations, transfer entropy, and Granger causality for spontaneous and evoked activity (Supplementary Figure 5).

Next, we examined whether the number of neurons and duration of recordings impacted the estimation of λ_decay. Despite the large number of neurons available, a subset of ∼500 neurons was sufficient to obtain a reliable estimate of λ_decay (Figure 9C, left). This is encouraging in terms of applying MLE to smaller datasets. Similarly, stable estimates were obtained with less than 5,000 ms of data (Figure 9C, right). This result is likely impacted by the stationarity of spike trains, with more data points required in datasets that exhibit complex fluctuations such as metastable states, a phenomenon that is not explored here.

We found that estimates of λ_decay overlapped across datasets of spontaneous and evoked activity (Figure 9D). This result is consistent with an exponential decay in spatial interactions that reflects structural properties of neuronal circuits and not fluctuations across different dynamical states. Simulations of the Izhikevich model showed that restricted spatial interactions (λ_decay = 3 mm^–1) result in low spike-count correlations (Figure 7A). This prediction agrees with cortical recordings, where the mean λ_decay was high for both spontaneous (λ_decay = 3.49 mm^–1) and evoked (λ_decay = 4.12 mm^–1) activity. Further, across all datasets, a negative correlation was found between values of λ_decay and mean spike-count correlations, showing that spatially restricted interactions yielded lower spike-count correlations (N = 14, R = −0.5477, P = 0.0527). Interestingly, simulated activity with values of λ_decay between 3–4 mm^–1 puts the Izhikevich model at the edge of a phase transition between high and low spike-count correlations. In this state, subtle changes in distance-based interactions can have a drastic effect of neuronal dynamics.

Finally, we compared the fit of an exponential decay to a half-Gaussian and an inverse square decay (Figure 10A). Log likelihood estimates for these functions were derived across all 14 datasets using the spatial GLM (Materials and Methods) (Figure 10B). Overall, the exponential decay attained a higher log likelihood, indicative of a superior fit (Figure 10C). Of the 14 datasets analyzed, 5 yielded a higher log-likelihood for the half-Gaussian function (three datasets of spontaneous activity, and two evoked). When comparing the average best-fitting exponential and half-Gaussian functions, the key difference was in the steepness of the decay occurring within distances <1 mm, where the half-Gaussian function suggested a shallower decay (Figure 10D). The best fitting half-Gaussian function had a mean standard deviation across datasets of σ = 0.37, suggesting that 95% of functional interactions resided within a distance of 2σ = 0.74 mm. This value is consistent with estimates of synaptic connectivity in rodent visual cortex, where the slope of the half-Gaussian function varied between 0.49 (layer 2) and 0.3 (layer 3) (Hellwig, 2000).

Overall, these results suggest that within local circuits of the visual cortex, distance-dependent interactions are best characterized by an exponential decay. To be sure, this finding may depend on several factors not explored here. For broader spatial scales, including inter-areal interactions, the decay function would be broader to reflect longer-range connections (Campagnola et al., 2022; Ercsey-Ravasz et al., 2013; Horvát et al., 2016). Further, we cannot rule out the possibility that distance-based interactions reflect unique connectivity patterns across cortical layers and may be altered under different regimes of neuronal activity that reflect cognitive states (Huang et al., 2019).

To estimate λ_decay based on the spiking activity of a linear Poisson model with baseline firing rate r₀ and coupling strength α, a maximum likelihood estimator (MLE) was employed as follows. Assuming Bernoulli spike data S_i(t) ∈ {0,1}, the probability of spiking is given by

where

assuming a large population of neurons, and σ(x)=1/(1 + exp⁡(−x)). The likelihood is

Taking the log-likelihood,

Defining P_i(t)=σ(h_i(t)) and applying the chain rule, the derivative of Equation 4 is obtained as

with the full derivative given by

Practically, finding the best-fitting λ_decay, denoted λdecay*, involves numerically minimizing the negative log-likelihood:

using gradient descent. We assume that the log-likelihood function (Equation 4) is maximized at the true parameter λ_decay employed to generate spike trains, which holds for long recordings (T→∞) given the law of large numbers. We further assume that the parameter λ_decay is identifiable, meaning that distinct values of λ_decay lead to distinct distributions of the observed data. This holds true given that λ_decay uniquely determines the distance-dependent interactions, and thus spike probabilities. Under reasonable conditions of smoothness for ℒ_MLE(λ_decay) and finite Fisher information, the MLE is asymptotically normal and converges in probability to the true parameter λ_decay.

The MLE approach employed to estimate λ_decay from spike trains was based on the same method as described above for Poisson neurons but adapted to a scenario where only the spike raster is known, and not the baseline rate or the interaction strength between neurons. The log likelihood ℒ_MLE is the joint probability of observing the spiking activity given λ_decay (Equation 4). Given that the baseline firing rate and coupling strength for the Izhikevich model and calcium imaging data are assumed to be unknown, we substituted the spike probability of Equation 3 for

where N is the number of neurons. The interaction term I_i(t) was computed as

where

and S_j(t) is the spike activity of the presynaptic neuron j at time t. Finding the best-fitting λ_decay involves numerically minimizing the negative log-likelihood.

In the case of a half-Gaussian function, W(i,j) (Equation 5) is replaced with

where we aim to estimate σ (the width of the half-Gaussian function). The log likelihood function is the same as before. We can also substitute in a log-normal function:

with dimensionless parameters σ and μ. To focus on estimating σ, we fix μ as follows:

For a linear distance-dependent function,

and for an inverse square function,

For the spatial GLM, the expression for the spike probability P(S_i(t)=1) follows Equation 3. The log⁡ℒ_GLM(λ_decay,α,r₀) follows Equation 4, but is estimated over three parameters instead of only λ_decay. Gradients for these parameters are obtained as follows,

where

The spatial GLM captures the spatiotemporal dynamics of the system, including time-dependence and directional interactions (i.e., from neuron j to i). By comparison, MLE is not a full generative model — instead, summary statistics are extracted from the spike raster as a whole, then a spatially-dependent function, such as an exponential decay, is fitted to capture how functional connectivity declines with spatial distance.

Transfer entropy is a powerful information-theoretic measure used to quantify directed interactions between time series (Izzi et al., 2024; Kajiwara et al., 2021; Lizier et al., 2008; Lungarella et al., 2007; Novelli et al., 2019; Schreiber, 2000; Ursino et al., 2020; Verdes, 2005; Vicente et al., 2011; Voges et al., 2024). It quantifies causal influences by assessing the degree to which the past activity of one neuron (potentially pre-synaptic) improves the predictability of the future activity of another neuron (potentially post-synaptic) (Kobayashi and Shinomoto, 2024; Vicente et al., 2011). Given two neurons X (source) and Y (target), transfer entropy is calculated as:

where y_t+1 is the state of the target neuron at the next time step, ytk are the past states of the target neuron Y up to k time steps, xt(l-△⁢t) are the past states of the source neuron up to l−△t time steps, and △t is a time delay (set to 1 ms).

To compute Granger causality for a pair of neurons, we assumed two spike trains X(t) and Y(t) over a time window [0,T] and binned them using bins of size △t = 10 ms. We then chose a history window ℋ = 100 ms and for each time t, create covariates using the past spiking activity of both neurons:

Then, two models were fitted to represent the instantaneous probability of neuron X. A reduced model based only on the history of X is provided by

and a full model combining the history of X and Y is given by

where μ_X is the mean firing rate of neuron X, and w_XX and wY⁢XT are weights to be estimated. Both models were fitted numerically using maximum likelihood estimation. The likelihood of the reduced model is given by

where Δ is the bin width. Similarly, for the full model,

The strength of causal interaction is given by

which yields a unitless ratio of how much better the model performs at predicting the activity of neuron X by including the history of neuron Y.

A network of Izhikevich neurons with N = 1,000 excitatory neurons was simulated (Figure 3A). The dynamics of each neuron’s membrane potential (v_i) and recovery variable (u_i) are described by

Spikes are emitted when v_i=v_thresh, where v_thresh = 30 mV. Then, the membrane potential and recovery variable are reset to v_i←c_i and u_i←u_i + d_i, where a_i, b_i, c_i, and d_i are parameters set as follows: a_i = 0.02, b_i = 0.2, c_i = −65, and d_i = 8. The input current is

where J_ij is the synaptic weight from neuron j to i, and S_j = 1 if v_i=v_threshold, and S_j = 0 otherwise. The external input is Iie⁢x⁢t=𝒩⁢(μ,σ2) and the total input is Ii=Iis⁢y⁢n+Iie⁢x⁢t. Each neuron i was assigned a position x_i=(x_i,y_i) in two-dimensional space, and the probability of pairwise connections decayed exponentially with distance (Equation 1), with c = 30. Sparse connectivity was generated by connecting any given pair of neurons with a 10% probability, while the remainder of weights were set to zero. While inhibitory neurons are beyond our scope, they could be incorporated in further refinements of the model.

Using mean field theory, we define the mean firing rate of the network as

where ⟨⋅⟩ denotes averaging over time and neurons. We then replace the exact spike trains with their statistical averages,

where J is the mean synaptic strength and N is the number of presynaptic neurons. Synaptic fluctuations are approximated by a stationary Gaussian process with mean JNr and variance J²r. Thus, the effective input current follows

where η_i(t) is a Gaussian noise term with variance

Let us consider two neurons i and j in the network. These neurons receive inputs from a subset of common presynaptic neurons. Let S_ij denote the fraction of presynaptic neurons shared by both i and j. The total input currents are:

where pre(i) denotes the set of presynaptic neurons projecting to i and j. The fraction of shared presynaptic inputs between neurons i and j is given by

The shared presynaptic inputs introduce correlations in the inputs I_i and I_j. The covariance of inputs between neurons i and j is

with variance

The correlation coefficient is given by

Thus, the correlation simplifies to

For weak coupling, we can approximate the correlation as

Assuming a homogeneous spatial network where the shared input fraction decays exponentially with distance,

we can rewrite Equation 6 as

Given that neurons receive synaptic inputs from a common pool of presynaptic neurons, the fraction of shared inputs between two neurons separated by distance d is approximately

We can then rewrite the approximate correlations as

As λ_decay→0, connectivity becomes all-to-all, hence C_ij approaches its upper bound,

Conversely, when λ_decay is large, shared inputs become negligible, hence C_ij→0.

Activity from mouse V1 cortex was obtained using two-photon calcium imaging. Details are found in related work (Stringer et al., 2019a). Each recording site measured 12 × 12 μm², and the vertical spacing between adjacent sites was 20 μm. The recording area was 0.9 mm by 0.935 mm, with plane depth of 0.35 mm. Spike were extracted using OASIS (Friedrich et al., 2017; Pachitariu et al., 2018).