Synchronization in spiking neural networks with short and long connections and time delays

Brain activity is characterized by rich spatiotemporal dynamics across scales. Among these, synchronization¹ between rhythms of neuronal populations² has been shown to be relevant for brain function³ and dysfunction.⁴ Different spatial patterns of synchronization have been hallmarked both during tasks⁵ and rest^6,7 and through various neuroimaging modalities such as functional Magnetic Resonance Imaging (fMRI), Magneto-Electro-Encephalography (MEG/EEG), and Electrocorticography (ECoG). Similarly, neurons are known to synchronize locally through short-range connections, as seen in Local Field Potentials (LFPs), which also bear relevance for the brain function.⁸ Hence, even if it is established that the macroscale rhythms emerge from the microscopic spiking neuronal activity,^9,10 the interplay between the two scales is often poorly understood.

Generative computational models have been commonly used^11–14 to understand how the brain’s anatomical structure, as defined by its connectome,¹⁵ impacts structured macroscopic patterns of brain activity. Here, we explore how these emergent dynamics depend not only on the weights of the connectome but also on the time it takes to propagate through the fibers.^16–18 Together, they define a space–time structure that is determinant for the network behavior, including resonance, synchronization, and emergence of brain rhythms.¹⁹ 

The space-invariant short-range cortical connections in the brain are characterized by negligible delays.¹³ Alongside these are heterogeneous links from the connectome, distinguished by their time delays and weights.²⁰ To investigate the interplay between these two types of connections, this study simplifies the scenario by focusing on the impact of a single long-range link. Through systematic alteration of its length, and consequently its delay by keeping a constant transmission velocity, the study aims to elucidate the effects of these long-range connections on neural interactions.

The effects of network heterogeneity, including that of time delays, were studied on the macroscopic level,^21–23 as well as for the combined macro- and mesoscopic representations^24–28 using a neural field approach, including both global and local connectivities.^13,29,30 However, neural fields in those cases were either modeled by mean-field representations of the neuronal masses or by phenomenological neuronal models that correspond to reduced representations of Hodgkin–Huxley neurons. In the case of integrate-and-fire neuron models, there is a lack of a precise and self-consistent mathematical framework to link microscopic and macroscopic dynamics within neural field theory.³¹ In particular, the validation of macroscopic phenomena emerging from spiking neuronal networks is missing.

In this study, we focus on the emergence of synchronization in neural fields coupled through long-range connections (LRCs) with finite transmission velocity.^{24,26,27,32–35} We use a methodology similar to Jirsa and Stefanescu³² but with a more detailed and physiologically grounded model of neurons including adaptation phenomena. This last aspect is critical since the previously used simpler neuron model embodied strong slow attractive manifolds known to be the key determinant in network synchronization. As a consequence, the synchronization results reported in the literature may not always translate into real-world brain networks.

To address this, we introduce a detailed description of the spiking network followed by a description of a reduced analytical model to capture the observed synchronization. In the results, we first present the effects of the LRC parameters and noise intensity on synchronization. Then, we explore the joint effects of neuron parameters (e.g., adaptation) and LRC delay on the modulation of synchronization. We then compare the spiking network and the analytical model, and finally, we describe the parameters (E/I balance, network size) that tend to break synchronization, always in relation to the LRC delay and the emergent dynamics.

The LRC is defined by a source cell and a target cell. In the source cell, on top of the short-range connections, 40% of the excitatory neurons project to the distant target following a two-dimensional Gaussian probability density (⁠ $σ l o n g=2.0$⁠) [see Figs. 1(c) and 1(d)] centered in the target cell. The weight of the LRC is independent of the local ones (by default: 15.0 nS), and their delay is uniform between neurons and independent of the distance.

In addition, we used this theoretical model to check for co-existence of multiple stable solutions for $Ω$⁠, as it is common in delay-coupled phase oscillators.^23,56–58 To do so, we solved geometrically the system of transcendental Eqs. (5a) and (5b) for the recovered $K h$ and $K$⁠, where the frequency of synchronization in our system is confined in the interval $[ω− k h,ω+ k h]$⁠, as a consequence of (5a). Next, we calculated the linear stability of the observed solutions by determining the largest real eigenvalue of the Jacobian of the system.

We seek to recover and understand the conditions in which the network of spiking neurons is synchronized and to characterize its spatiotemporal dynamics, particularly with respect to the long-range connection. To this end, we conducted analyses of parameter spaces of synchronization as a function of global and local network parameters, i.e., weight and delay of unidirectional LRC for the former, and excitatory weights and variance of white noise for the latter.

Overall, time series analysis reveals that different levels of synchronization are associated with specific spatiotemporal patterns of activity. When the network is asynchronous [Fig. 3(a), corresponding to mark A in Fig. 4(b)], the activity corresponds to a traveling wave along the longest axis of the grid. As expected, the histogram is flat, and the spike trains form a diagonal of spikes due to the ordering of the neurons. In contrast, when the network is highly synchronized [Fig. 3(c) corresponding to mark C in Fig. 4(b)], its activity alternates between bursts of activation and quiet periods. In this case, activity bursts are very short (less than 20 ms) and composed of multiple bumps (locally synchronized firing of neurons) mainly starting from the target cell, but also from other random locations due to the noise. The eight snapshots indicate that when two waves collide, they collapse. Between these two extremes, spatiotemporal patterns present a gradient where synchronization is linked to shorter activity bursts (on the scale of the whole network) and a higher number of bumps from random locations [see Fig. 3(b) corresponding to mark B in Fig. 4(b)].

When the noise is turned off, the network is generally highly synchronized, as shown by the coefficient of variation being systematically higher than 1.0 [Fig. 4(a)]. Time series (Fig. S1 in the supplementary material) are characterized by periodic bumps starting from the target cell, forming a circular wave, and eventually colliding with itself. We find an almost proportional relationship between the level of synchronization and the delay. This relation is preserved when the weight of the LRC, the number of connections, and the position of the target cell are varied. However, the network is silent when the LRC weight is less than 5 nS, when the delay is less than 70 ms, or when the standard deviation of the Gaussian kernel is lower than 0.5.

When the LRC are removed, the parameter space [Fig. 4(c)] shows that the network is synchronized only when the noise is strong enough (variance greater than 1500 pA) and if the excitatory synaptic weight is sufficiently large (higher than 8 nS). When these conditions are met, synchronization decreases as noise increases with little effect of excitatory synaptic weights. A related phenomenon to this is the appearance of synchronized bumps displaying chaotic-like behavior⁵⁹ (see Figs. S2A and S2B in the supplementary material). When the excitatory synaptic weight is too low, below 8 nS, the system remains in an asynchronous regime (Fig. S2 in the supplementary material).

Under the combination of noise and the presence of LRC, the results show that the level of synchronization depends on the delay along the LRC [Fig. 4(b)], with a regular pattern of higher synchronization occurring at 0, 110, 220, or 330 ms. This indicates inherent oscillations with a period of 110 ms. The observed impact of the delay on synchronization is similar to the one found in populations of oscillators.^56,60 The same effect was also observed for continuous spiking neurons.³² 

If the weight of the LRC or the number of connections is too low, this effect is lost, and the network remains synchronized regardless of the delay. Like before, we found that the position of the target cell does not have a substantial impact on synchronization. In other words, if the contribution of the LRC is too small, then the noise is the dominant cause for the emergence of synchronized activity, Figs. 4(b) and 4(c). Otherwise, LRC modulates synchronization through the ratio between the delay and mean firing rate, determining the increased/decreased synchronization ranges.

When the synchronization is facilitated, e.g., point C in Fig. 4(b), the wave is homotopic and propagates faster [Fig. 3(b) compared to Fig. 3(c)], compared to the case when the synchronization is decreased, e.g., point B in Fig. 4(b). In certain conditions, point A in Fig. 4(b), regular traveling waves appear. In this case, the source and the target cells are anti-phase (close to half of the period), and the bump does not collapse but continues to travel across the network.

The adaptive current is a mechanism that affects the frequency of the spikes. When adaptation parameters increase, new spikes are harder to trigger so neurons receive less input from their neighbors. This impacts the mean-firing rate, which serves as a “natural frequency” of self-sustained network oscillations in the case of synchronization, ultimately leading to a similar effect in synchronized oscillators as reducing the time delays (Fig. 5).

The adaptation has three main parameters: the time constant $τ w$⁠, the spike-triggered adaptation $b$⁠, and the sub-threshold adaptation $a$⁠. As seen in Fig. 5, synchronization can be found for all ranges of these parameters. The periodic effect of the LRC delay on synchronization found previously is preserved when adaptation is varied but tends to fade out as the adaptation increases. The spike-triggered adaptation and the time constant are associated with a gradient of $C V F R$⁠, not found with subthreshold adaptation. Low values of $τ w$ and $b$ induce an increase of the mean firing rate. $a$ also impacts firing rate but to a lesser degree (details in Fig. S3 in the supplementary material).

We applied heterogeneous external currents, $I e$⁠, distributed across neurons, to create heterogeneity of excitability within the population (see Sec. II). The voltage membrane excitability has opposite effects on adaptation as it facilitates the generation of spikes. Therefore, as seen in Fig. 6(a), increasing heterogeneity tends to reduce the $C V F R$ and the periodic effect of the delay also fades out for low values, mirroring the effect of adaptation. Although not presented here, if heterogeneity is too extreme, it affects the balance of the network and pushes it into an asynchronous irregular state. The relationship between voltage reset and firing rate can be found in Fig. S7 in the supplementary material.

Synchronization is mainly present for a voltage reset between $−$63 and $−$45 mV [see Fig. 6(b)], otherwise the network is mostly asynchronous. In the bottom part of the panel, there is almost no activity in the network (close to zero mean firing rate and CV, Figs. S11E and S11F in the supplementary material) except for some high-frequency bursts for specific values of the LRC delay, yielding narrow stripes of synchronized activity. This separation between synchronous and asynchronous regimes may be explained by a strong non-linear relationship between the mean firing rate and the voltage reset (Fig. S4D in the supplementary material).

When the refractory time is varied [Fig. 6(c)], the parameter space is separable in two subspaces: under and over 1.7 ms, approximately. For low refractory times, more than 95% of the neurons have irregular burst activities (Fig. S7 in the supplementary material). In this case, reducing refractory time increases the variation of the $C V F R$ proportionally to the delay. However, large values of $C V F R$ (above 3) capture the presence of large waves (see Fig. S5A in the supplementary material). When the refractory period exceeds 1.7 ms, neurons regularly spike with few bursts with little to no effect of the delay along the LRC. In this scenario, very low values of CV (below 0.5) capture traveling wave activity (shown in Fig. S5C in the supplementary material).

In this section, we compare synchronization results obtained in the spiking network and the analytical model derived for two uni-directionally delay-coupled phase oscillators. We now characterize the synchronization for the spiking neurons using the Kuramoto order parameter (R) instead of the $CV F R$⁠, while for the pairwise coherence of the phase oscillators, we use the absolute value of the phase locking value (PLV). Specifically, we select [Figs. 7(a)–7(c)] the LRC weight, sub-threshold adaptation, and heterogeneity between neurons for a range of the LRC delay as the parameters to be compared with $K,ω, K h$⁠, and $τ$ of the analytical model [Figs. 7(d)–7(f)]. We find that the model of phase oscillators consistently reproduces the activity of the spiking neuronal network. The weight of the LRC has the same effect on synchronization in the network as $K$ (long-range coupling) between the oscillators, despite the changes in the frequency of synchronization. The pattern of synchronization found when varying the sub-threshold adaptation parameter is also retrieved when decreasing the natural frequency (⁠ $ω$⁠) of the oscillators. Lastly, decreasing the heterogeneity of the neurons in the network tends to facilitate synchronization by increasing the probability of simultaneously having a greater number of neurons close to the firing threshold. Hence, this has the same effect as increasing the local coupling $K h$ that also brings neurons closer to synchrony and we find a similar pattern of synchronization for different values of $τ$⁠.

Additionally, solving geometrically the system of transcendental Eqs. (5a) and (5b) for the recovered values of $K h$ and $K$⁠, allowed us to calculate the linear stability of the observed solutions. Even though we observed a small range in the values of $τ$ where two stable fixed points coexist (see Fig. S6A and S6B in the supplementary material), we always observed only one of them in the numerical simulations of the phase oscillators. It is, thus, likely that this co-existence is a numerical artifact of the method for geometrically solving the transcendental equations since the range of co-existence is within the margins of the numerical error for the solver. However, another possibility is that the other solution was not observed because it is less stable than the one that was observed (see Fig. S6C in the supplementary material).

As seen in previous studies,⁶¹ synchronizations occurs in a network of spiking neurons for a specific balance between excitation and inhibition that allows for wave propagation. Other factors, such as the number of neurons and the synaptic weights also affect network synchronization.³² Plots in Fig. 8 show the effect of the relative strength between excitatory and inhibitory synaptic weights on synchronization for different noise intensities (0, 1800, and 2600 pA). Without noise, only the LRC delay has a significant effect on synchronization, similarly to the previous results of Fig. 4 with the same characteristic proportional relationship. The main difference is that for a low E/I ratio (approximately under 0.06) asynchronous activity is seen [but no activity at all for low $k$ in the plot (a)]. This regime is also independent of noise intensity and simulations always have a firing rate over 50 Hz.

When noise is introduced, there appears a specific range of E/I balance that generates synchronization and the parameter space is divided into two regimes: one asynchronous regime, either dominated by excitation (top) or by inhibition (bottom), and a synchronized “balanced” regime in the middle.^62,63 The increase in noise intensity narrows this range of synchronization by pushing the upper limit down to a lower E/I balance.

To explore the effect of the number of neurons on synchronization, we varied the percentage of active neurons in the network. As seen from left to right in Fig. 9, the sparsening of neurons is associated with an increase in the lower bound of excitatory synaptic weights needed to reach synchronization. However, when the number of neurons is too low [plot (d)] increasing excitation can no longer compensate for the loss of activity. This balance between the number of neurons and excitatory synaptic weights is also reflected in the firing rate as the increase of synaptic weights can compensate for the reduction in the number of neurons to keep it constant⁶⁴ (Fig. S8 in the supplementary material).

In addition, we explored different grid sizes without changing the density of the neurons and found invariant results except for very small grids (less than five cells) (see Fig. S14 in the supplementary material).

In the previously shown results with unidirectional LRC [see Figs. S10A and 4(b)] , varying the weights had very little or no effect on the synchronization for any fixed delay. However, when the LRC becomes bidirectional, synchronization is lost for LRC weights higher than 4 nS (Fig. S10B in the supplementary material).

In the absence of synchronization, activity resembles that of Fig. 3(a) with an unstable traveling wave across the network. When tuned near synchronization, we again find bumps of activity, but this time originating from both sides of the LRC, Figs. S12B and S12C in the supplementary material. If the weight of the LRC is increased above 4 nS (from point C to D in Fig. S10C in the supplementary material), bumps of activation appear with origins alternating between both ends of the LRC, with roughly one-third of them never propagating through the network. The principal cause of the disappearance of the periodic dependence on the delay in the case of bidirectional LRC is that regardless of the noise level, the network is not synchronized. Still, a traveling wave occurs (see Fig. S9A in the supplementary material), while the same condition for a unidirectional LRC would have caused synchronization instead (see Fig. S9B in the supplementary material).

The present study demonstrates the conditions necessary for synchronization in a two-dimensional surface comprised of adaptive exponential integrate-and-fire neurons,³⁸ incorporating long-range connections (LRC) with transmission delays. In particular, we have identified how the global synchronization is affected by the delay along the LRC. This relationship typically follows a regular pattern, contingent upon other parameters that we have examined. These results align with earlier works, confirming this phenomenon in a neural field derived from spiking neural networks,³² as well as with broader literature illustrating Arnold’s tongues for synchronization.⁶⁵ This relationship is contingent upon the ratio between natural frequencies and time delays.^22,56 Based on the latter, we also introduced a reduced analytical model. Besides its simplicity, it omits all the higher-order terms in the Fourier series of the phase-interaction function, and it still qualitatively captures and describes the observed synchronization phenomena, thus strengthening the validity of the numerically observed results.

We identified four key conditions that influence synchronization modulation by the LRC in the network. First, the network with LRC naturally tends to synchronize at a frequency that is dependent on the delay along the LRC. Adding an external input (here in the form of random noise) is also crucial for achieving synchronization and revealing a more intricate relationship between synchronization and delay. Second, the LRC must be sufficiently strong (measured by the weight and width of the kernel) to affect the network. However, synchronization no longer occurs if the LRC are bidirectional and their influence is too strong. The third condition is that the network must maintain a balance between inhibition and excitation^61–63 for the LRC to modulate the synchronization effectively. Finally, the number of neurons must be large enough³⁴ or sufficiently interconnected for this phenomenon to take place.

Through investigating the parameters space, we have done thorough parametrizations that correspond to various neuron sub-types. We started with values corresponding to adapting neurons in sub-threshold activities.⁴² It corresponds to type II neuron models (as shown in Fig. S4C in the supplementary material) exhibiting resonator regimes³⁹ and type I Phase Response Curve (PRC).⁴¹ When modifying neuronal properties, we changed firing rate dynamics, neuron type (I/II), voltage membrane regime, and PRC (as shown in Figs. S4 and S17 in the supplementary material). For instance, when the sub-threshold is shifted to low values, the adapting neuron in sub-threshold activities becomes a type I neuron with integrator regimes. Interestingly, we found that the type of neuron is not a significant condition for the effect of the LRC on synchronization. Moreover, when setting the spike-triggered adaptation to 0 pA, the neuron model becomes equivalent to an exponential integrate-and-fire neuron without adaptation. We observed that the results remain consistent for this neuron type as well. This indicates that mean field models derived from various neuron types would still conform to our findings. Additionally, considering that for every delay value along the LRC, there exists at least one set of neuron parameters (i.e., a neuron type) enabling network synchronization, it would be possible to deduce neuron properties from knowledge about the LRC, under certain conditions about the spiking activity.

In terms of generative mechanisms, we identified that the bumps are the main cause of synchronization captured by the coefficient of variation of the firing rate $C V F R$⁠, and that the bumps correspond to a burst of activity where many neurons are collectively triggered in a short period of time. These bursts are generated periodically and their generative process resembles the one of action potentials by following a sequence of voltage accumulation, threshold crossing, spike triggering, and a refractory period. When the noise drive is high enough to trigger a spike at the single neuron level, the synaptic currents will accumulate locally through the synaptic weights (mediated by the local E/I balance). If the local voltage becomes large enough, it generates circular waves that propagate through the network. In the case of a homogeneous network, the probability of these events forming anywhere on the grid is uniform, potentially generating multiple waves from different locations at the same time. The collapsing of waves is followed by a silent period due to the refractory time ensuing collective firing of neurons. Indeed, the local increase of adaptation current prevents neurons from generating new spikes. The adaptation current then decays until new events occur. As a result, the frequency of bursts depends on a complex interaction between the adaptive current, synaptic current, voltage reset, refractory time, E/I balance, and noise.

Finally, the characteristics of both, neurons and the network, influence how emergent phenomena like synchronization and waves develop as spatiotemporal patterns in the population. This study establishes connections between these emergent dynamics and the properties of long-range fiber-like connections. Understanding and describing this self-organization is crucial for understanding complex systems like the brain. In particular, both the delay and the weights of the fibers play crucial roles in network synchronization, alongside neuronal properties. Specifically, variations in adaptive current and voltage reset can affect the synchronization of certain neurons, impacting the emergent rhythm of the neuronal population.⁶⁶ Furthermore, wave propagation, as demonstrated in this study, is crucial for synchronization as an emergent dynamic. Our results corroborate previous findings by showing that local modifications such as a change in the E/I balance or a loss of neurons can disrupt higher-order phenomena (e.g., wave propagation and synchronization) and, thus, could potentially be responsible for altered brain activity observed in pathology.

One major limitation of our work lies in the choice of the synchronization metrics. Here, synchronization is quantified using the coefficient of variation of firing rate (⁠ $C V F R$⁠), which is less precise in the case of bursts (the cause of variation can be the high activities of a sub-population). Nevertheless, we were able to reproduce some of the main results using the time-averaged Kuramoto order-parameter,⁵³ when we compared the analytical model to the network activity. The time-averaged Kuramoto order-parameter has the particularity to be robust to the bursting effect, as shown by Borges et al.⁵⁴ Finally, for the analytical model, a better fit would have been obtained by a more detailed model containing higher-order terms in the Fourier series of the coupling function.⁶⁷ Together with the higher-order network interactions⁶⁸ these play an important role in brain activity,^69,70 but the current parsimonious parametrization strengthens the universality of the observed phenomena.

Among other limitations, we can note that each network configuration was only tested once, thus bringing nuances to the conclusions that can be drawn about multiple regimes potentially coexisting for the same set of parameters. Indeed, some simulations showed transitions between regimes (see Fig. S13 in the supplementary material), making it difficult to determine if the simulation represents a complete stationary regime or only one regime in a more complex landscape. However, simulations were run long enough to reach a steady regime whenever it existed. Moreover, the analysis and especially the simulations of the reduced model did not give a strong indication for the coexistence of different solutions, and hence we did not pursue a more thorough exploration of dynamical regimes, since it would have been tedious and beyond the focus of this study.

Lastly, the analysis was not fully exhaustive; synaptic parameters such as Gaussian connectivity profiles,⁷¹ time-constants, or reversal potentials, for example, were not explored because we chose to focus on neuronal properties important for mean field derivations.

In summary, the synchronization of a network of biological neurons with long-range connections and finite transmission velocities is modulated by the delay of these connections. The main conditions for this effect are the pre-existence of spontaneous activity driven by an external input (or noise), a significant contribution of the LRC, a balanced network (excitation/inhibition), and a large number of neurons. Given that our objective was to provide evidence for the validity of neural field models, future work could include testing the relaxation of neural field model assumptions (e.g., plasticity, delays, and weight distributions) and using multiple heterogeneous neuron models. Additionally, further analysis with bidirectional LRC is needed to understand spatiotemporal patterns in this case. This analysis should consider larger grids or modified connectivity, as recent findings suggest that traveling waves are based on a fraction of neurons in a population.⁶⁶ Validation of this phenomenon with more recent neural field models (e.g., Zerlaut et al.⁷² or Augustin et al.⁷³) is also recommended.

See the supplementary material for additional information about other network metrics that we provide in open-source databases and details of the network simulation for the reproducibility of the results. The two supplementary notes describe the databases and the additional metrics applied to the network. The supplementary table is based on the proposition of a table by Nordlie et al.74 for the reproducibility of the neural network model. The figures show the network topology, selection of different time series, and parameter explorations with selected additional metrics.

We acknowledge the use of Fenix Infrastructure resources, which are partially funded by the European Union’s Horizon 2020 Research and Innovation Programme through the ICEI project under Grant Agreement No. 800858. This research has received funding from the European Union’s Horizon 2020 Framework Programme for Research and Innovation under the Specific Grant Agreement Nos. 945539 (Human Brain Project SGA3), 101147319 (EBRAINS 2.0 Project), and 101137289 (Virtual Brain Twin Project), by the French government grant managed by the Agence Nationale de la Recherch (ANR) Reference No. ANR-22-PESN-0012 (France 2030 program), and from the ANR under Grant No. ANR-17-CE37-0001-CONNECTOME. The authors gratefully acknowledge the support team of the supercomputer Piz Daint at CSCS—Swiss National Supercomputing Centre.

The authors have no conflicts to disclose.

Spase Petkoski and Viktor K. Jirsa contributed equally to this work as the last author.

Lionel Kusch: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Martin Breyton: Writing – original draft (equal); Writing – review & editing (equal). Damien Depannemaecker: Investigation (supporting); Methodology (supporting); Writing – review & editing (equal). Spase Petkoski: Conceptualization (equal); Formal analysis (lead); Investigation (lead); Methodology (lead); Project administration (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Viktor K. Jirsa: Conceptualization (equal); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (lead); Resources (equal); Writing – original draft (equal); Writing – review & editing (equal).

The code for the simulation and analysis of this paper is freely available under v2 Apache license at https://bitbucket.org/lionelkusch/torus_long_range_connection/src/paper/. Each exploration’s analysis results are registered in databases available on Zenodo (http://dx.doi.org/10.5281/zenodo.7629156) with the associate code. Supplementary material Note 1 describes the contents of them.