Maximum likelihood estimation of spatially dependent interactions in large populations of cortical neurons

Understanding how functional connectivity between cortical neurons varies with spatial distance is crucial for characterizing large-scale neural dynamics. However, inferring these spatial patterns is challenging when spike trains are collected from large populations of neurons. Here, we present a maximum likelihood estimation (MLE) framework to quantify distance-dependent functional interactions directly from observed spiking activity. We validate this method using both synthetic spike trains generated from a linear Poisson model and biologically realistic simulations performed with Izhikevich neurons. We then apply the approach to large-scale electrophysiological recordings from V1 cortical neurons. Our results show that the proposed MLE approach robustly captures spatial decay in functional connectivity, providing insights into the spatial structure of population-level neural interactions.

Introduction

Understanding the organization of brain networks is a fundamental challenge in neuroscience (Honey et al., 2010). A key principle of neuronal wiring is spatial embedding, where the probability of synaptic contact between neurons decreases as a function of physical distance (Ercsey-Ravasz et al., 2013; Horvát et al., 2016; Perin et al., 2011; Song et al., 2005; Stetter et al., 2012; Zingg et al., 2014). This principle holds across multiple spatial scales, from local neuronal networks to large-scale brain regions (Markov et al., 2013), and has been leveraged in computational models to explain various properties of cortical connectivity, neuronal communication (Ercsey-Ravasz et al., 2013), and functions including learning and working memory (Keller et al., 2024).

Advances in experimental neuroscience have enabled large-scale recordings of neural activity, capturing the dynamics of broad neuronal populations (Stosiek et al., 2003; Stringer et al., 2019a,b). These advances provide unprecedented insight into the complex interactions that govern brain function (Chen et al., 2022; Honey et al., 2010; Seth, 2008; Suárez et al., 2020). However, characterizing distance-dependent interactions within these populations remains challenging, as the underlying synaptic connectivity is typically unknown in vivo (Tian et al., 2024). Addressing this challenge would provide a rich understanding of neuronal organization and enable comparisons across brain regions that play distinct roles in behavioral and mental states.

To examine distance-dependent interactions, most studies have relied on functional connectivity to measure statistical dependencies or mutual information between patterns of neuronal activity (Friston et al., 2013; Stetter et al., 2012). While measuring structural connectivity requires access to the physical connections between neurons, based on invasive or post-mortem techniques, functional connectivity can be derived from recorded neural activity. This makes functional connectivity a pivotal approach for investigating the dynamic relationship between interconnected neural elements. Thus, estimating functional connectivity is critical for understanding information flow in the brain and the organization of neuronal circuits.

A variety of measures have been proposed to infer functional connectivity, including pairwise correlations, transfer entropy, and Granger causality (Tian et al., 2024). Distance-based interactions estimated from pairwise correlations (Smith and Kohn, 2008; Yatsenko et al., 2015) and spike-based functional connectivity (English et al., 2017) have been measured in cortex. These interactions show specificity across layers (Senzai et al., 2019), cell types (Khoury et al., 2025), brain regions (Horvát et al., 2016), as well as a dependence on sensory features such as stimulus orientation (Smith and Kohn, 2008; Yatsenko et al., 2015).

However, measures of functional connectivity have notable limitations. First, estimating functional connectivity is complicated by the inherently noisy and irregular nature of neuronal signals (Aertsen et al., 1989; Ito et al., 2011; Izzi et al., 2024; Liao et al., 2010; Vicente et al., 2011). Second, functional connectivity does not directly measure synaptic efficacy, making interpretation difficult (Ursino et al., 2020). Third, conventional functional connectivity measures do not explicitly account for spatial interactions, often overlooking the impact of physical distance on neuronal communication (Holmgren et al., 2003; Jiang et al., 2015; Perin et al., 2011; Vegué et al., 2017). As a result, these methods fail to capture the underlying spatial structure of neuronal networks, limiting our ability to understand their functional organization. Given the critical role of spatially structured interactions in neural processing (Ohiorhenuan et al., 2010), new approaches are needed to extract meaningful spatial information from large-scale recordings.

Here, we introduce a novel maximum likelihood estimation (MLE) approach to characterize distance-dependent functional connectivity. Instead of estimating individual interaction strengths between all pairs of neurons, this method estimates a global parameter that governs the spatial decay of connectivity across a neuronal population. The proposed MLE assumes an exponential relationship between connectivity and physical distance (Ercsey-Ravasz et al., 2013; Horvát et al., 2016) and infers this decay parameter from observed spike trains.

We validate this approach using synthetic networks of Poisson and Izhikevich neurons (Izhikevich, 2003) and apply it to large-scale recordings from mouse V1 neurons obtained from two-photon calcium imaging (Stringer et al., 2019a). Our results demonstrate that MLE successfully extracts distance-dependent interactions, revealing spatially structured neuronal communication.

Results

Distance-based interactions in a linear Poisson model

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,

$p\left(S_i\left(t\right)\right)=r_0+\sum _j\left(J_{ij}\cdot \mathrm{\alpha }\cdot S_j\left(t-1\right)\right)$

FIGURE 1

Figure 1. Spatially embedded Poisson neurons. (A) Schema of a stationary Poisson neuron receiving spiking activity S_j(t−1) from j^th presynaptic neurons at time t-1 via weighted connections α. Presynaptic activity is summed and a stochastic process generates spikes p(S_i(t)). (B) Illustration of spiking activity over time-steps t for a population of neurons randomly positioned on a two-dimensional grid. (C) The probability of connection between pairs of neurons is determined by an exponential function with slope parameter λ_decay.

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),

$\mathbb{P}(J_{ij}=1)=exp(-{\lambda }_{\text{decay}}{d}_{ij})(1)$

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.

FIGURE 2

Figure 2. Estimating spatially dependent interactions in a network of Poisson neurons. (A) Spike raster generated by 1,000 interconnected neurons. (B) The log likelihood function was maximal near the true value of λ_decay. (C) The maximum likelihood estimate of λ_decay across different runs of the Poisson model. Dashed line: unity line showing the correspondence between the true and estimated values of λ_decay. Five approaches were employed to estimate λ_decay, from left to right: MLE, transfer entropy, pairwise correlations, Granger causality, and spatial GLM. (D) Estimated value of exponential decay. Filled black circles are individual data points obtained from connection probability (MLE and spatial GLM), transfer entropy, spike-count correlations, Granger causality.

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.

TABLE 1

Table 1. Identification of various distance-dependent functions.

Comparison to related approaches

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.

FIGURE 3

Figure 3. Comparison of run time (A) and mean error (B) of four approaches to estimate exponential spatial decay. GLM, generalized linear model; TE, transfer entropy; GC, Granger causality; corr, 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).

Izhikevich model

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).

FIGURE 4

Figure 4. Recurrent Izhikevich model with distance-dependent connectivity. (A) Schema of the recurrent model. The connections (black lines) of three neurons (red circles) are highlighted to show spatially-delimited interactions with neighboring units. (B) Spike raster obtained from a simulation of 1,000 neurons. Distance-dependent connectivity followed an exponential distribution with λ_decay = 5 mm^{– 1}. (C) Distribution of Fano factors obtained from all individual neurons using 10 ms non-overlapping bins. (D) Distribution of spike-count correlations.

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).

FIGURE 5

Figure 5. Estimation of distance-dependent connectivity using MLE. (A) Examples of exponential functions estimated from different runs of the Izhikevich model that varied the true value of λ_decay. (B) Estimation of λ_decay across various runs of the model. Dashed line: best-fitting regression.

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).

FIGURE 6

Figure 6. Impact of distance-dependent connectivity on population dynamics in the Izhikevich model. (A) Two exponential decay functions tested on the model. (B) Spike rasters (left) and distribution of synchronous spikes using 10 ms bins (right) obtained with two different decay functions.

FIGURE 7

Figure 7. Spike-count correlations are modulated by distance-dependent connectivity. (A) Exponentially decaying connectivity modulates correlations. Above: schema of connectivity patterns. Short-range: mean of correlations between neurons with distance <0.1 mm. Long-range: neurons with distance >1 mm. Average: mean over all neurons. (B) Analytical expression relating correlations to spatial decay and coupling strength (Equation 2). (C) Randomized network obtained after shuffling weights from panel “A.” (D) Zoom of panels “A” and “C” for higher values of λ_decay.

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”),

$C_{ij}\approx \frac{J^2}{J^2+{\mathrm{\sigma }}^2}{e}^{-2{\lambda }_{\text{decay}}{d}_{ij}}(2)$

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.

Application of MLE to calcium imaging data

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.

FIGURE 8

Figure 8. Two-photon calcium imaging data from rodent visual cortex. (A) Three-dimensional position of individual neurons from a subset of the population. (B) Raster obtained from calcium imaging data (subset of neurons from panel “A”). (C) Distribution of Fano factors, with mean indicated by the arrow. (D) Mean firing rate (left) and spike-count correlations (right) across 9 datasets of spontaneous activity and 5 datasets of evoked 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).

FIGURE 9

Figure 9. Maximum likelihood estimation (MLE) applied to two-photon calcium data. (A) Log-likelihood function obtained from a single raster of spontaneous activity (10,164 neurons, 15,761 ms). Vertical dashed line shows the optimal value of λ_decay obtained with MLE. (B) Best-fitting exponential function based on the log-likelihood from panel “A” shown on a two-dimensional spatial grid. (C) Fluctuations in the estimated λ_decay as the number of neurons and time increased. Dashed line: λ_decay value obtained with all neurons and time-steps included. (D) Estimates of λ_decay across different recording sessions of spontaneous and stimulus-evoked activity. Filled circles are individual sessions; horizontal lines are within-group averages.

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).

FIGURE 10

Figure 10. Comparison of MLE across distance-dependent functions. (A) Examples of half-Gaussian and lognormal functions (μ = 1.5) characterizing distance-based interactions. (B) Examples of log likelihood estimates obtained on single calcium imaging datasets. (C) Comparison of log likelihood across three candidate functions and 14 experimental datasets. Values shown were obtained by subtracting the minimum log likelihood. (D) Average best-fitting exponential (λ_decay = 3.8 mm^{– 1}) and half-Gaussian (σ = 0.37 mm) functions across all experimental datasets.

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).

Discussion

This study introduced a statistical method based on maximum likelihood estimation (MLE) to investigate distance-dependent interactions from observed neuronal population activity. A key advantage of this approach is its ability to estimate spatial decay directly from spike data without requiring a full reconstruction of all pairwise interactions, leading to considerable savings in computational time when scaled to large networks, hence enabling the characterization of spatial interactions across large populations of neurons. Combined with simulations of spatially-embedded networks, this approach draws links between the structure and dynamics of neuronal circuits by relating spatial interactions and synchronous activity.

When applied to neuronal activity obtained from calcium imaging (Stringer et al., 2019a), MLE revealed spatially-delimited interactions that were consistent across both spontaneous and evoked states, suggesting an underlying network architecture that is strongly biased toward stronger interactions amongst nearby neurons. This result is consistent with a growing understanding that neural circuits are fundamentally shaped by their spatial embedding (Wang and Kennedy, 2016). This embedding emerges across a wide range of spatial scales, ranging from local microcircuits (Campagnola et al., 2022) to large inter-areal projections (Markov et al., 2013), and carries significant consequences for both the cytoarchitectural organization of the brain as well as the dynamical patterns of activity that emerge from this architecture.

One consequence of spatially constrained connectivity is the relative rarity of “giant hubs” characterized by neurons (or whole areas) that receive dense afferences from both nearby and distant projection sites (Deco et al., 2021). While disproportionally rare, these hubs may fulfill important functional roles in relaying information across broad networks, with implications for information processing (Ohiorhenuan et al., 2010). As such, a full understanding of cortical organization relies on several requisites, including not only distance-dependent interactions, but also hubs and locally clustered connectivity (Litwin-Kumar and Doiron, 2012).

The prevalence of an exponential decay in spatially embedded brain networks may be attributed to a trade-off: on one side, the need for efficient communication, and on the other, the energetic and spatial costs associated with axonal extension, myelination, and maintenance. Theoretical work shows that an exponential decay optimizes this trade-off in a scale invariant fashion, making spatially-embedded models a widely applicable framework for describing brain networks at various levels of organization (Ercsey-Ravasz et al., 2013; Horvát et al., 2016).

Spatially embedded interactions carry profound implications for the dynamical state of neuronal circuits, shaping the scope of activity that can be expressed at the population level. Various forms of spatiotemporal dynamics, including traveling waves (Huang et al., 2019) and spiral waves (Paraskevov and Zendrikov, 2017) are supported by spatially ordered connectivity and form a hallmark of information propagation across brain areas (Keller et al., 2024). Further, theoretical results presented here demonstrate a direct impact of distance-based interactions on spike-count correlations (Figure 7), showing a key role of spatially constrained connectivity in shaping population activity. This prediction of the model will need to be validated in experiments where the strength of distance-based interactions can be modulated, for instance by altering sensory input (Ohiorhenuan et al., 2010).

While our analyses revealed that an exponential decay function prevailed over alternatives including a half-Gaussian and a lognormal function, this result may depend upon several factors including the dynamical state of the network, the spatial scale of analysis, the laminar origin of the data, and the brain area considered (Horvát et al., 2016; Paraskevov and Zendrikov, 2017; Pastorelli et al., 2018). In particular, understanding how connectivity varies across cortical layers—with superficial layers showing more lateral, long-range connections and deep layers exhibiting more localized connectivity—could provide important insights into the organization of neural circuits (Holmgren et al., 2003; Levy and Reyes, 2012). As a caveat, candidate functions including an exponential and half-Gaussian decay may yield similar predictions at shorter spatial scales (< 1 mm) (Figure 10D).

While an MLE approach provided accurate estimates of spatially-dependent interactions in both a linear Poisson model (Figure 2C) and an Izhikevich network (Figure 5B), it is worth re-iterating sources of error that may compromise accuracy, including the number of neurons considered, duration of activity, firing rate, subsampling of neurons, spatial edge effects, domain size, and neuronal delays. Further, because MLE evaluates functional interactions and not structural connectivity, it reflects a combination of mono- and poly-synaptic interactions, and hence may overestimate the spatial extent of direct synaptic connections. In recent work, the connection probability of mouse neocortex as a function of lateral intersomatic distance was estimated by fitting a Gaussian decay with σ values ranging from between 0.1 mm and 0.13 mm (Campagnola et al., 2022). These values are lower than our findings (mean of σ = 0.37 mm), but comparable to visual cortical reconstruction in other work (Hellwig, 2000).

In general, estimates of functional connectivity inferred from spike data may serve as a proxy for synaptic connectivity when several conditions are met, including: (1) spike trains are stationary over the time window of estimation; (2) functional connectivity predominantly targets monosynaptic interactions; (3) the common input between pairs of neurons is low; (4) synaptic connectivity is sparse; (5) excitation and inhibition is balanced, to avoid high synchronization or inactivity; (6) neural dynamics are approximately linear (Cai et al., 2017; Chen et al., 2011; Endo et al., 2021; Graber et al., 2024; Kobayashi and Shinomoto, 2024; Masud and Borisyuk, 2011; Okatan et al., 2005). While these conditions are not always met in practice, they increase the reliability of synaptic connectivity inferred from functional interactions.

One limitation of the proposed approach is that while the spatial GLM successfully identified distance-based interactions generated from exponential, half-Gaussian, linear, and inverse square functions, it confounded linear and lognormal functions. Therefore, we did not include these functions in the analysis of calcium imaging data. As a result, we could not assess the fit of these functions and cannot rule out the possibility that other functions not considered here may offer a more accurate characterization of cortical interactions. Ultimately, however, experimentalists may be more interested in how these interactions are altered across layers (Senzai et al., 2019), cell types (Khoury et al., 2025), brain regions (Horvát et al., 2016), and stimuli (Ohiorhenuan et al., 2010; Smith and Kohn, 2008; Yatsenko et al., 2015), particularly given that several distance-dependent functions have a similar shape at short distances.

Another consideration is that while the proposed MLE approach provides a useful framework for investigating spatially dependent interactions, it describes a global decay function across an entire network and cannot resolve individual interactions between pairs of neurons. This limitation turns out to be an advantage in large-scale networks where we aim to characterize a spatial decay function without having to estimate all pairwise connections. In addition, MLE is not designed to distinguish between spatial interactions driven by local recurrent connections versus those driven by external inputs to a neuronal population.

A future direction is to apply MLE to distinguish between the spatial scales of different neuronal subtypes, including pyramidal neurons as well as subclasses of inhibitory interneurons (parvalbumin, somatostatin, and vasopressin) (Fino and Yuste, 2011; Jiang et al., 2015; Tremblay et al., 2016). Characterizing subtype-specific spatial patterns would deepen our understanding of how local and longer-range interactions contribute to network function.

More broadly, the proposed framework may offer insights into neurodevelopmental and neurodegenerative disorders marked by alterations in spatial organization. For instance, Alzheimer′s disease is associated with synaptic loss and reduced connectivity density (DeKosky and Scheff, 1990), while conditions such as epilepsy and autism spectrum disorder involve distinct disruptions in spatial circuit structure (Stoner et al., 2014). Understanding these differences at a subtype-resolved spatial scale could help clarify the circuit-level mechanisms underlying such pathologies.

To conclude, the proposed MLE approach offers an efficient and accessible tool for investigating spatially dependent interactions in neuronal circuits. This analysis provides a promising avenue to deepen our understanding of how structural constraints shape neural dynamics. Further, it opens the door to characterizing the relation between spatial organization and dynamics in large-scale data obtained from healthy and pathological brain circuits.

Materials and methods

Maximum likelihood estimation in a linear Poisson model

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

$P(S_i\left(t\right)=1)=\mathrm{\sigma }\left(h_i(t)\right)(3)$

where

$h_i\left(t\right)=r_0+\sum _{j\ne i}\mathrm{\alpha }\cdot S_j\left(t-1\right)\cdot \exp \left(-{\lambda }_{\text{decay}}\cdot d_{ij}\right)$

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

${ℒ}_{MLE}\left({\lambda }_{\text{decay}}\right)=\prod _{t=2}^T\prod _{i=1}^N\mathrm{\sigma }{\left(h_i\left(t\right)\right)}^{S_i(t)}\cdot {\left(1-\mathrm{\sigma }\left(h_i\left(t\right)\right)\right)}^{1-S_i(t)}$

Taking the log-likelihood,

$\log {ℒ}_{MLE}\left({\lambda }_{\text{decay}}\right)=\sum _{t=2}^T\sum _{i=1}^N[S_i\left(t\right)\log \left(\mathrm{\sigma }\left(h_i\left(t\right)\right)\right)+(4)$

$(1-S_i(t))log(1-\mathrm{\sigma }\left(h_i\left(t\right)\right))]$

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

$\frac{d}{d{\lambda }_{\text{decay}}}\log {ℒ}_{MLE}({\lambda }_{\text{decay}})=$

$\sum _{t=2}^T\sum _{i=1}^N\left[\frac{S_i\left(t\right)-P_i\left(t\right)}{P_i\left(t\right)\left(1-P_i\left(t\right)\right)}\cdot \frac{d{P}_i\left(t\right)}{d{\lambda }_{\text{decay}}}\right]$

with the full derivative given by

$\frac{d}{d{\lambda }_{\text{decay}}}\log {ℒ}_{MLE}({\lambda }_{\text{decay}})=-\mathrm{\alpha }\sum _{t=2}^T\sum _{i=1}^N(S_i\left(t\right)-P_i\left(t\right))\cdot$

$\sum _{j\ne i}{S}_j\left(t-1\right)d_{ij}\cdot e^{-{\lambda }_{\text{decay}}{d}_{ij}}$

Practically, finding the best-fitting λ_decay, denoted ${\lambda }_{\text{decay}}^{*}$, involves numerically minimizing the negative log-likelihood:

${\lambda }_{\text{decay}}^{*}=arg\underset{{\lambda }_{\text{decay}}}{min}-\log {ℒ}_{MLE}({\lambda }_{\text{decay}})\mathit{\hspace{1em}}$

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.

Maximum likelihood estimation: exponential function

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

$P(S_i\left(t\right)=1)\frac{I_i\left(t\right)}{N}$

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

$I_i\left(t\right)=\sum _j{S}_j\left(t\right)W(i,j)$

where

$W(i,j)=e^{-{\lambda }_{\text{decay}}{d}_{ij}}(5)$

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.

Alternative functions for maximum likelihood estimation

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

$W(i,j)=exp\left(-\frac{d_{ij}^2}{2{\mathrm{\sigma }}^2}\right)$

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:

$W(i,j)=\frac{1}{d_{ij}\mathrm{\sigma }\sqrt{2\mathrm{\pi }}}exp\left(-\frac{{\left(\log d_{ij}^2-\mathrm{\mu }\right)}^2}{2{\mathrm{\sigma }}^2}\right)$

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

$\mathrm{\mu }=\frac{1}{N\left(N-1\right)}\sum _{i\ne j}log\left(d_{i\ne j}\right)$

For a linear distance-dependent function,

$W(i,j)=1-\mathrm{\sigma }{d}_{ij}$

and for an inverse square function,

$W(i,j)=\frac{1}{d_{ij}^2+\mathrm{\sigma }}$

Spatial generalized linear model

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,

$\frac{\partial {ℒ}_{GLM}}{\partial r_0}=\sum _{t=2}^T\sum _{i=1}^N\left(S_i\left(t\right)-P_i(t)\right)$

$\frac{\partial {ℒ}_{GLM}}{\partial \mathrm{\alpha }}=\sum _{t=2}^T\sum _{i=1}^N\left(S_i\left(t\right)-P_i(t)\right)\cdot h_i(t)$

$\frac{\partial {ℒ}_{GLM}}{\partial {\lambda }_{\text{decay}}}=\sum _{t=2}^T\sum _{i=1}^N\left(S_i\left(t\right)-P_i(t)\right)\cdot g_i(t)$

where

$g_i\left(t\right)=\sum _{j\ne i}{S}_j\left(t-1\right)\cdot [-\mathrm{\alpha }{d}_{ij}\exp (-{\lambda }_{\text{decay}}{d}_{ij})]$

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

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:

$T{E}_{X\to Y}(\mathrm{△}t)=$

$\sum _{y_{t+1,y_t^{(k)},x_t^{(l)}}}p(y_{t+1},y_t^{(k)},x_t^{(l-\mathrm{△}t)})log\frac{p(y_{t+1}|y_t^{(k)},x_t^{(l-\mathrm{△}t)})}{p(y_{t+1}|y_t^{(k)})}$

where y_t+1 is the state of the target neuron at the next time step, $y_t^k$ are the past states of the target neuron Y up to k time steps, $x_t^{(l-\mathrm{△}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).

Granger causality

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:

${\mathbf{\text{h}}}_X\left(t\right)=[x_{t-1},x_{t-2},,x_{t-\mathscr{H}}]$

${\mathbf{\text{h}}}_Y\left(t\right)=[y_{t-1},y_{t-2},,y_{t-\mathscr{H}}]$

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

$\log {\lambda }_X^{reduced}(t)={\mathrm{\mu }}_X+{\mathbf{\text{w}}}_{XX}{\mathbf{\text{h}}}_X\left(t\right)$

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

$\log {\lambda }_X^{full}(t)={\mathrm{\mu }}_X+{\mathbf{\text{w}}}_{XX}{\mathbf{\text{h}}}_X\left(t\right)+{\mathbf{\text{w}}}_{YX}^T{\mathbf{\text{h}}}_Y\left(t\right)$

where μ_X is the mean firing rate of neuron X, and w_XX and ${\mathbf{\text{w}}}_{YX}^T$ are weights to be estimated. Both models were fitted numerically using maximum likelihood estimation. The likelihood of the reduced model is given by

${ℒ}_{reduced}=\sum _t\left[x(t)\cdot \log \left({\lambda }_X^{reduced}(t)\cdot \mathrm{\Delta }\right)-{\lambda }_X^{reduced}(t)\cdot \mathrm{\Delta }\right]$

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

${ℒ}_{full}=\sum _t\left[x(t)\cdot \log \left({\lambda }_X^{full}(t)\cdot \mathrm{\Delta }\right)-{\lambda }_X^{full}(t)\cdot \mathrm{\Delta }\right]$

The strength of causal interaction is given by

$G{C}_{Y\to X}=\frac{{ℒ}_{full}-{ℒ}_{reduced}}{\left|{ℒ}_{reduced}\right|}$

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.

Izhikevich model

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

$\frac{d{v}_i}{dt}=0.04v_i^2+5v_i+140-u_i+I_i$

$\frac{d{u}_i}{dt}=a_i\left(b_i{v}_i-u_i\right)$

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

$I_i^{syn}=\sum _{j=1}^N{J}_{ij}{S}_j$

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 $I_i^{ext}=𝒩(\mathrm{\mu },{\mathrm{\sigma }}^2)$ and the total input is $I_i=I_i^{syn}+I_i^{ext}$. 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

$r=⟨S_j(t)⟩$

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

$I_i^{syn}\left(t\right)\approx J\sum _j{S}_j\left(t\right)\approx JNr$

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

$I_i\left(t\right)\approx JNr+{\mathrm{\eta }}_i\left(t\right)$

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

$Var\left({\mathrm{\eta }}_i\right)=J^{2}r+{\mathrm{\sigma }}^2$

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:

$I_i=\sum _{k\in pre\left(i\right)}J{S}_k$

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

$K_{ij}=\frac{\left|pre(i)pre(j)\right|}{\left|pre(i)\right|}$

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

$Cov(I_i,I_j)\approx J^{2}r{K}_{ij}$

with variance

$Var\left(I_i\right)=Var\left(I_i\right)=J^{2}r+{\mathrm{\sigma }}^2$

The correlation coefficient is given by

$C_{ij}=\frac{Cov(I_i,I_j)}{\sqrt{Var\left(I_i\right)Var\left(I_j\right)}}$

Thus, the correlation simplifies to

$C_{ij}=\frac{J^{2}r{K}_{ij}}{J^{2}r+{\mathrm{\sigma }}^2}$

For weak coupling, we can approximate the correlation as

$C_{ij}\approx \frac{J^2}{J^2+{\mathrm{\sigma }}^2}{K}_{ij}(6)$

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

$K_{ij}\approx e^{-2{\lambda }_{\text{decay}}{d}_{ij}}$

we can rewrite Equation 6 as

$C_{ij}\approx \frac{J^2}{J^2+{\mathrm{\sigma }}^2}{e}^{-2{\lambda }_{\text{decay}}{d}_{ij}}$

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

$K_{ij}\approx {\int }_{0}P\left(d\right)e^{-{\lambda }_{\text{decay}}d}\mathit{d}d=e^{-2{\lambda }_{\text{decay}}d}$

We can then rewrite the approximate correlations as

$C_{ij}\approx \frac{J^2{e}^{-2{\lambda }_{\text{decay}}d}}{J^2{e}^{-2{\lambda }_{\text{decay}}d}+{\mathrm{\sigma }}^2}$

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

$C_{ij,{\lambda }_{\text{decay}}\to 0}\to \frac{J^2}{J^2+{\mathrm{\sigma }}^2}$

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

Calcium data

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).

Data availability statement

Publicly available datasets were analyzed in this study. These data can be found here: https://janelia.figshare.com/articles/dataset/Recordings_of_ten_thousand_neurons_in_visual_cortex_during_spontaneous_behaviors/6163622.

Author contributions

CG: Data curation, Formal analysis, Investigation, Methodology, Software, Writing – review & editing. JT: Data curation, Formal analysis, Investigation, Methodology, Software, Writing – review & editing, Conceptualization, Funding acquisition, Project administration, Resources, Supervision, Validation, Visualization, Writing – original draft.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This work was funded by a grant to J.P.T. from the Natural Sciences and Engineering Research Council of Canada (NSERC Grant No. 210977).

Conflict of interest

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

Generative AI statement

The authors declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncom.2025.1639829/full#supplementary-material

References

Aertsen, A. M., Gerstein, G. L., Habib, M. K., and Palm, G. (1989). Dynamics of neuronal firing correlation: Modulation of “effective connectivity.” J. Neurophysiol. 61, 900–917. doi: 10.1152/jn.1989.61.5.900

PubMed Abstract | Crossref Full Text | Google Scholar

Cai, Z., Neveu, C. L., Baxter, D. A., Byrne, J. H., and Aazhang, B. (2017). Inferring neuronal network functional connectivity with directed information. J. Neurophysiol. 118, 1055–1069. doi: 10.1152/jn.00086.2017

PubMed Abstract | Crossref Full Text | Google Scholar

Campagnola, L., Seeman, S. C., Chartrand, T., Kim, L., Hoggarth, A., Gamlin, C., et al. (2022). Local connectivity and synaptic dynamics in mouse and human neocortex. Science 375:eabj5861. doi: 10.1126/science.abj5861

PubMed Abstract | Crossref Full Text | Google Scholar

Chen, Y., Rosen, B. Q., and Sejnowski, T. J. (2022). Dynamical differential covariance recovers directional network structure in multiscale neural systems. Proc. Natl. Acad. Sci. U.S.A. 119:e2117234119. doi: 10.1073/pnas.2117234119

PubMed Abstract | Crossref Full Text | Google Scholar

Chen, Z., Putrino, D. F., Ghosh, S., Barbieri, R., and Brown, E. N. (2011). Statistical inference for assessing functional connectivity of neuronal ensembles with sparse spiking data. IEEE Trans. Neural Syst. Rehabil. Eng. Publ. IEEE Eng. Med. Biol. Soc. 19, 121–135. doi: 10.1109/TNSRE.2010.2086079

PubMed Abstract | Crossref Full Text | Google Scholar

Deco, G., Sanz Perl, Y., Vuust, P., Tagliazucchi, E., Kennedy, H., and Kringelbach, M. L. (2021). Rare long-range cortical connections enhance human information processing. Curr. Biol. 31, 4436–4448.e5. doi: 10.1016/j.cub.2021.07.064.

PubMed Abstract | Crossref Full Text | Google Scholar

DeKosky, S. T., and Scheff, S. W. (1990). Synapse loss in frontal cortex biopsies in Alzheimer’s disease: Correlation with cognitive severity. Ann. Neurol. 27, 457–464. doi: 10.1002/ana.410270502

PubMed Abstract | Crossref Full Text | Google Scholar

Endo, D., Kobayashi, R., Bartolo, R., Averbeck, B. B., Sugase-Miyamoto, Y., Hayashi, K., et al. (2021). A convolutional neural network for estimating synaptic connectivity from spike trains. Sci. Rep. 11:12087. doi: 10.1038/s41598-021-91244-w

PubMed Abstract | Crossref Full Text | Google Scholar

English, D. F., McKenzie, S., Evans, T., Kim, K., Yoon, E., and Buzsáki, G. (2017). Pyramidal cell-interneuron circuit architecture and dynamics in hippocampal networks. Neuron 96, 505–520.e7. doi: 10.1016/j.neuron.2017.09.033.

PubMed Abstract | Crossref Full Text | Google Scholar

Ercsey-Ravasz, M., Markov, N. T., Lamy, C., Van Essen, D. C., Knoblauch, K., Toroczkai, Z., et al. (2013). A predictive network model of cerebral cortical connectivity based on a distance rule. Neuron 80, 184–197. doi: 10.1016/j.neuron.2013.07.036

PubMed Abstract | Crossref Full Text | Google Scholar

Fino, E., and Yuste, R. (2011). Dense inhibitory connectivity in neocortex. Neuron 69, 1188–1203. doi: 10.1016/j.neuron.2011.02.025

PubMed Abstract | Crossref Full Text | Google Scholar

Friedrich, J., Zhou, P., and Paninski, L. (2017). Fast online deconvolution of calcium imaging data. PLoS Comput. Biol. 13:e1005423. doi: 10.1371/journal.pcbi.1005423

PubMed Abstract | Crossref Full Text | Google Scholar

Friston, K., Moran, R., and Seth, A. K. (2013). Analysing connectivity with Granger causality and dynamic causal modelling. Curr. Opin. Neurobiol. 23, 172–178. doi: 10.1016/j.conb.2012.11.010

PubMed Abstract | Crossref Full Text | Google Scholar

Graber, C., Vlasov, Y., and Schwing, A. (2024). Dynamic organization of visual cortical networks revealed by machine learning applied to massive spiking datasets. eLife 13, 1–48. doi: 10.7554/eLife.95449.1

Crossref Full Text | Google Scholar

Hellwig, B. (2000). A quantitative analysis of the local connectivity between pyramidal neurons in layers 2/3 of the rat visual cortex. Biol. Cybern. 82, 111–121. doi: 10.1007/PL00007964

PubMed Abstract | Crossref Full Text | Google Scholar

Holmgren, C., Harkany, T., Svennenfors, B., and Zilberter, Y. (2003). Pyramidal cell communication within local networks in layer 2/3 of rat neocortex. J. Physiol. 551, 139–153. doi: 10.1113/jphysiol.2003.044784

PubMed Abstract | Crossref Full Text | Google Scholar

Honey, C. J., Thivierge, J.-P., and Sporns, O. (2010). Can structure predict function in the human brain? NeuroImage 52, 766–776. doi: 10.1016/j.neuroimage.2010.01.071

PubMed Abstract | Crossref Full Text | Google Scholar

Horvát, S., Gămănuţ, R., Ercsey-Ravasz, M., Magrou, L., Gămănuţ, B., Essen, D. C. V., et al. (2016). Spatial embedding and wiring cost constrain the functional layout of the cortical network of rodents and primates. PLoS Biol. 14:e1002512. doi: 10.1371/journal.pbio.1002512

PubMed Abstract | Crossref Full Text | Google Scholar

Huang, C., Ruff, D. A., Pyle, R., Rosenbaum, R., Cohen, M. R., and Doiron, B. (2019). Circuit models of low-dimensional shared variability in cortical networks. Neuron 101, 337–348.e4. doi: 10.1016/j.neuron.2018.11.034.

PubMed Abstract | Crossref Full Text | Google Scholar

Ito, S., Hansen, M. E., Heiland, R., Lumsdaine, A., Litke, A. M., and Beggs, J. M. (2011). Extending transfer entropy improves identification of effective connectivity in a spiking cortical network model. PLoS One 6:e27431. doi: 10.1371/journal.pone.0027431

PubMed Abstract | Crossref Full Text | Google Scholar

Izhikevich, E. M. (2003). Simple model of spiking neurons. IEEE Trans. Neural Netw. 14, 1569–1572. doi: 10.1109/TNN.2003.820440

PubMed Abstract | Crossref Full Text | Google Scholar

Izzi, J. V. R., Ferreira, R. F., Girardi, V. A., and Pena, R. F. O. (2024). Identifying effective connectivity between stochastic neurons with variable-length memory using a transfer entropy rate estimator. Brain Sci. 14:442. doi: 10.3390/brainsci14050442

PubMed Abstract | Crossref Full Text | Google Scholar

Jiang, X., Shen, S., Cadwell, C. R., Berens, P., Sinz, F., Ecker, A. S., et al. (2015). Principles of connectivity among morphologically defined cell types in adult neocortex. Science 350:aac9462. doi: 10.1126/science.aac9462

PubMed Abstract | Crossref Full Text | Google Scholar

Kajiwara, M., Nomura, R., Goetze, F., Kawabata, M., Isomura, Y., Akutsu, T., et al. (2021). Inhibitory neurons exhibit high controlling ability in the cortical microconnectome. PLoS Comput. Biol. 17:e1008846. doi: 10.1371/journal.pcbi.1008846

PubMed Abstract | Crossref Full Text | Google Scholar

Keller, T. A., Muller, L., Sejnowski, T. J., and Welling, M. (2024). A spacetime perspective on dynamical computation in neural information processing systems. arXiv [Preprint] doi: 10.48550/arXiv.2409.13669

Crossref Full Text | Google Scholar

Khoury, C. F., Ferrone, M., and Runyan, C. A. (2025). Local differences in network organization in the auditory and parietal cortex, revealed with single neuron activation. J. Neurosci. 45:e1385242025. doi: 10.1523/JNEUROSCI.1385-24.2025

PubMed Abstract | Crossref Full Text | Google Scholar

Kobayashi, R., and Shinomoto, S. (2024). Inference of monosynaptic connections from parallel spike trains: A review. Neurosci. Res. 215, 37–46. doi: 10.1016/j.neures.2024.07.006

PubMed Abstract | Crossref Full Text | Google Scholar

Levy, R. B., and Reyes, A. D. (2012). Spatial profile of excitatory and inhibitory synaptic connectivity in mouse primary auditory cortex. J. Neurosci. 32, 5609–5619. doi: 10.1523/JNEUROSCI.5158-11.2012

PubMed Abstract | Crossref Full Text | Google Scholar

Liao, W., Mantini, D., Zhang, Z., Pan, Z., Ding, J., Gong, Q., et al. (2010). Evaluating the effective connectivity of resting state networks using conditional granger causality. Biol. Cybern. 102, 57–69. doi: 10.1007/s00422-009-0350-5

PubMed Abstract | Crossref Full Text | Google Scholar

Litwin-Kumar, A., and Doiron, B. (2012). Slow dynamics and high variability in balanced cortical networks with clustered connections. Nat. Neurosci. 15, 1498–1505. doi: 10.1038/nn.3220

PubMed Abstract | Crossref Full Text | Google Scholar

Lizier, J. T., Prokopenko, M., and Zomaya, A. Y. (2008). Local information transfer as a spatiotemporal filter for complex systems. Phys. Rev. E 77:026110. doi: 10.1103/PhysRevE.77.026110

PubMed Abstract | Crossref Full Text | Google Scholar

Lungarella, M., Pitti, A., and Kuniyoshi, Y. (2007). Information transfer at multiple scales. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 76:056117. doi: 10.1103/PhysRevE.76.056117

PubMed Abstract | Crossref Full Text | Google Scholar

Markov, N. T., Ercsey-Ravasz, M., Lamy, C., Ribeiro Gomes, A. R., Magrou, L., Misery, P., et al. (2013). The role of long-range connections on the specificity of the macaque interareal cortical network. Proc. Natl. Acad. Sci. U. S. A. 110, 5187–5192. doi: 10.1073/pnas.1218972110

PubMed Abstract | Crossref Full Text | Google Scholar

Masud, M. S., and Borisyuk, R. (2011). Statistical technique for analysing functional connectivity of multiple spike trains. J. Neurosci. Methods 196, 201–219. doi: 10.1016/j.jneumeth.2011.01.003

PubMed Abstract | Crossref Full Text | Google Scholar

Novelli, L., Wollstadt, P., Mediano, P., Wibral, M., and Lizier, J. T. (2019). Large-scale directed network inference with multivariate transfer entropy and hierarchical statistical testing. Netw. Neurosci. Camb. Mass 3, 827–847. doi: 10.1162/netn_a_00092

PubMed Abstract | Crossref Full Text | Google Scholar

Ohiorhenuan, I. E., Mechler, F., Purpura, K. P., Schmid, A. M., Hu, Q., and Victor, J. D. (2010). Sparse coding and high-order correlations in fine-scale cortical networks. Nature 466, 617–621. doi: 10.1038/nature09178

PubMed Abstract | Crossref Full Text | Google Scholar

Okatan, M., Wilson, M. A., and Brown, E. N. (2005). Analyzing functional connectivity using a network likelihood model of ensemble neural spiking activity. Neural Comput. 17, 1927–1961. doi: 10.1162/0899766054322973

PubMed Abstract | Crossref Full Text | Google Scholar

Orlandi, J. G., Stetter, O., Soriano, J., Geisel, T., and Battaglia, D. (2014). Transfer entropy reconstruction and labeling of neuronal connections from simulated calcium imaging. PloS One 9:e98842. doi: 10.1371/journal.pone.0098842

PubMed Abstract | Crossref Full Text | Google Scholar

Pachitariu, M., Stringer, C., and Harris, K. D. (2018). Robustness of spike deconvolution for neuronal calcium imaging. J. Neurosci. 38, 7976–7985. doi: 10.1523/JNEUROSCI.3339-17.2018

PubMed Abstract | Crossref Full Text | Google Scholar

Paraskevov, A. V., and Zendrikov, D. K. (2017). A spatially resolved network spike in model neuronal cultures reveals nucleation centers, circular traveling waves and drifting spiral waves. Phys. Biol. 14:026003. doi: 10.1088/1478-3975/aa5fc3

PubMed Abstract | Crossref Full Text | Google Scholar

Pastorelli, E., Paolucci, P. S., Simula, F., Biagioni, A., Capuani, F., Cretaro, P., et al. (2018). “Gaussian and exponential lateral connectivity on distributed spiking neural network simulation,” in Proceedings of the 2018 26th Euromicro International Conference on Parallel, Distributed and Network-Based Processing (PDP), (Piscataway, NJ: IEEE), 658–665.

Google Scholar

Perin, R., Berger, T. K., and Markram, H. (2011). A synaptic organizing principle for cortical neuronal groups. Proc. Natl. Acad. Sci. U. S. A. 108, 5419–5424. doi: 10.1073/pnas.1016051108

PubMed Abstract | Crossref Full Text | Google Scholar

Pillow, J. W., Shlens, J., Paninski, L., Sher, A., Litke, A. M., Chichilnisky, E. J., et al. (2008). Spatio-temporal correlations and visual signalling in a complete neuronal population. Nature 454, 995–999. doi: 10.1038/nature07140

PubMed Abstract | Crossref Full Text | Google Scholar

Roxin, A., Riecke, H., and Solla, S. A. (2004). Self-Sustained Activity in a Small-World Network of Excitable Neurons. Phys. Rev. Lett. 92:198101. doi: 10.1103/PhysRevLett.92.198101

PubMed Abstract | Crossref Full Text | Google Scholar

Schreiber, T. (2000). Measuring information transfer. Phys. Rev. Lett. 85, 461–464. doi: 10.1103/PhysRevLett.85.461

PubMed Abstract | Crossref Full Text | Google Scholar

Senzai, Y., Fernandez-Ruiz, A., and Buzsáki, G. (2019). Layer-specific physiological features and interlaminar interactions in the primary visual cortex of the mouse. Neuron 101, 500–513.e5. doi: 10.1016/j.neuron.2018.12.009.

PubMed Abstract | Crossref Full Text | Google Scholar

Seth, A. K. (2008). Causal networks in simulated neural systems. Cogn. Neurodyn. 2, 49–64. doi: 10.1007/s11571-007-9031-z

PubMed Abstract | Crossref Full Text | Google Scholar

Smith, M. A., and Kohn, A. (2008). Spatial and temporal scales of neuronal correlation in primary visual cortex. J. Neurosci. 28, 12591–12603. doi: 10.1523/JNEUROSCI.2929-08.2008

PubMed Abstract | Crossref Full Text | Google Scholar

Song, S., Sjöström, P. J., Reigl, M., Nelson, S., and Chklovskii, D. B. (2005). Highly nonrandom features of synaptic connectivity in local cortical circuits. PLoS Biol. 3:e68. doi: 10.1371/journal.pbio.0030068

PubMed Abstract | Crossref Full Text | Google Scholar

Stetter, O., Battaglia, D., Soriano, J., and Geisel, T. (2012). Model-free reconstruction of excitatory neuronal connectivity from calcium imaging signals. PLoS Comput. Biol. 8:e1002653. doi: 10.1371/journal.pcbi.1002653

PubMed Abstract | Crossref Full Text | Google Scholar

Stoner, R., Chow, M. L., Boyle, M. P., Sunkin, S. M., Mouton, P. R., Roy, S., et al. (2014). Patches of disorganization in the neocortex of children with autism. N. Engl. J. Med. 370, 1209–1219. doi: 10.1056/NEJMoa1307491

PubMed Abstract | Crossref Full Text | Google Scholar

Stosiek, C., Garaschuk, O., Holthoff, K., and Konnerth, A. (2003). In vivo two-photon calcium imaging of neuronal networks. Proc. Natl. Acad. Sci. U.S.A. 100, 7319–7324. doi: 10.1073/pnas.1232232100

PubMed Abstract | Crossref Full Text | Google Scholar

Stringer, C., Pachitariu, M., Steinmetz, N., Carandini, M., and Harris, K. D. (2019a). High-dimensional geometry of population responses in visual cortex. Nature 571, 361–365. doi: 10.1038/s41586-019-1346-5

PubMed Abstract | Crossref Full Text | Google Scholar

Stringer, C., Pachitariu, M., Steinmetz, N., Reddy, C. B., Carandini, M., and Harris, K. D. (2019b). Spontaneous behaviors drive multidimensional, brainwide activity. Science 364: 255. doi: 10.1126/science.aav7893

PubMed Abstract | Crossref Full Text | Google Scholar

Suárez, L. E., Markello, R. D., Betzel, R. F., and Misic, B. (2020). Linking structure and function in macroscale brain networks. Trends Cogn. Sci. 24, 302–315. doi: 10.1016/j.tics.2020.01.008

PubMed Abstract | Crossref Full Text | Google Scholar

Tian, Z.-Q. K., Chen, K., Li, S., McLaughlin, D. W., and Zhou, D. (2024). Causal connectivity measures for pulse-output network reconstruction: Analysis and applications. Proc. Natl. Acad. Sci. U. S. A. 121:e2305297121. doi: 10.1073/pnas.2305297121

PubMed Abstract | Crossref Full Text | Google Scholar

Tremblay, R., Lee, S., and Rudy, B. (2016). GABAergic interneurons in the neocortex: From cellular properties to circuits. Neuron 91, 260–292. doi: 10.1016/j.neuron.2016.06.033

PubMed Abstract | Crossref Full Text | Google Scholar

Ursino, M., Ricci, G., and Magosso, E. (2020). Transfer entropy as a measure of brain connectivity: A critical analysis with the help of neural mass models. Front. Comput. Neurosci. 14:45. doi: 10.3389/fncom.2020.00045

PubMed Abstract | Crossref Full Text | Google Scholar

Vegué, M., Perin, R., and Roxin, A. (2017). On the structure of cortical microcircuits inferred from small sample sizes. J. Neurosci. 37, 8498–8510. doi: 10.1523/JNEUROSCI.0984-17.2017

PubMed Abstract | Crossref Full Text | Google Scholar

Verdes, P. F. (2005). Assessing causality from multivariate time series. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 72:026222. doi: 10.1103/PhysRevE.72.026222

PubMed Abstract | Crossref Full Text | Google Scholar

Vicente, R., Wibral, M., Lindner, M., and Pipa, G. (2011). Transfer entropy–a model-free measure of effective connectivity for the neurosciences. J. Comput. Neurosci. 30, 45–67. doi: 10.1007/s10827-010-0262-3

PubMed Abstract | Crossref Full Text | Google Scholar

Voges, N., Lima, V., Hausmann, J., Brovelli, A., and Battaglia, D. (2024). Decomposing Neural Circuit Function into Information Processing Primitives. J. Neurosci. 44:e0157232023. doi: 10.1523/JNEUROSCI.0157-23.2023

PubMed Abstract | Crossref Full Text | Google Scholar

Wang, X.-J., and Kennedy, H. (2016). Brain structure and dynamics across scales: In search of rules. Curr. Opin. Neurobiol. 37, 92–98. doi: 10.1016/j.conb.2015.12.010

PubMed Abstract | Crossref Full Text | Google Scholar

Yatsenko, D., Josić, K., Ecker, A. S., Froudarakis, E., Cotton, R. J., and Tolias, A. S. (2015). Improved estimation and interpretation of correlations in neural circuits. PLoS Comput. Biol. 11:e1004083. doi: 10.1371/journal.pcbi.1004083

PubMed Abstract | Crossref Full Text | Google Scholar

Zingg, B., Hintiryan, H., Gou, L., Song, M. Y., Bay, M., Bienkowski, M. S., et al. (2014). Neural networks of the mouse neocortex. Cell 156, 1096–1111. doi: 10.1016/j.cell.2014.02.023

PubMed Abstract | Crossref Full Text | Google Scholar