How synaptic function controls critical transitions in spiking neuron networks: insight from a Kuramoto model reduction

The dynamics of synaptic interactions within spiking neuron networks play a fundamental role in shaping emergent collective behavior. This paper studies a finite-size network of quadratic integrate-and-fire neurons interconnected via a general synaptic function that accounts for synaptic dynamics and time delays. Through asymptotic analysis, we transform this integrate-and-fire network into the Kuramoto-Sakaguchi model, whose parameters are explicitly expressed via synaptic function characteristics. This reduction yields analytical conditions on synaptic activation rates and time delays determining whether the synaptic coupling is attractive or repulsive. Our analysis reveals alternating stability regions for synchronous and partially synchronous firing, dependent on slow synaptic activation and time delay. We also demonstrate that the reduced microscopic model predicts the emergence of synchronization, weakly stable cyclops states, and non-stationary regimes remarkably well in the original integrate-and-fire network and its theta neuron counterpart. Our reduction approach promises to open the door to rigorous analysis of rhythmogenesis in networks with synaptic adaptation and plasticity.

1 Introduction

Cooperative rhythms play a pivotal role in brain functioning. Fully or partially synchronized oscillations, observed across various frequency bands, underlie fundamental processes such as perception, cognition, and motor control Churchland and Sejnowski, 1992; Mizuseki and Buzsaki, 2014; Kopell et al., 2000. Extensive research has focused on the emergence of cooperative rhythms in networks of spiking and bursting neurons, encompassing synchronization Kopell et al., 2000; Brunel, 2000; Börgers and Kopell, 2003; Somers and Kopell, 1993; Izhikevich, 2007; Belykh et al., 2005; Ermentrout and Terman, 2010, partial and cluster synchronization Achuthan and Canavier 2009; Shilnikov et al., 2008; Belykh and Hasler, 2011; Schöll, 2016; Berner et al., 2021a, neural bumps Laing and Chow, 2001; Gutkin et al., 2001, and chimera states Olmi et al., 2011; Omelchenko et al., 2013.

Networks of spiking neurons with fast synaptic connections are often modeled via pulsatile on-off coupling, which sharply activates upon the arrival of a spike from a pre-synaptic cell. Such interactions are conveniently represented by networks of quadratic integrate-and-fire (QIF) models particularly suitable for large-scale simulations and analysis of cooperative dynamics Gerstner and Kistler (2002). The macroscopic dynamics of QIF networks have received extensive attention through the reduction to low-dimensional model descriptions, especially in the thermodynamic limit of infinite-dimensional networks Montbrió et al., 2015; Pazó and Montbrió, 2016; Devalle et al., 2017; Esnaola-Acebes et al., 2017; Devalle et al., 2018; Schmidt et al., 2018; Pietras et al., 2019; Montbrió and Pazó 2020; Lin et al., 2020; Gast et al., 2020; Taher et al., 2020; Byrne et al., 2022; Clusella et al., 2022; Clusella and Montbrió 2024; Ratas and Pyragas, 2018; Pyragas and Pyragas 2022, 2023; Coombes 2023; Ferrara et al., 2023. Notably, Montbrió et al. (2015) derived exact macroscopic equations for QIF networks, uncovering an effective coupling between firing rate and mean membrane potential governing network dynamics. Pietras et al. (2023) offered an analytical description of QIF network macroscopic dynamics, extending beyond the Ott-Antosen ansatz Ott and Antonsen (2008) and exploring various fast synaptic pulse profile choices. The impact of synaptic time delay on the collective dynamics of integrate-and-fire networks with sharply activated synaptic coupling, modeled by the Dirac delta function, was also extensively explored (Ernst et al., 1995; Devalle et al., 2018; Ratas and Pyragas, 2018; Pyragas and Pyragas 2022, 2023). In particular, Devalle et al. (2018) reduced a QIF model with synaptic delay to a set of firing rate equations to analyze the existence and stability of partially synchronous states. Ratas and Pyragas (2018) employed a Lorenzian ansatz to characterize macroscopic oscillations of a QIF network with heterogeneous time-delayed delta function synapses. However, there is a lack of analytical studies on the role of slower synaptic activation, potentially in the presence of time delays, in controlling critical phase transitions in QIF networks. Nevertheless, since the seminal paper by Van Vreeswijk et al. (1994), it has been recognized that slow inhibitory and excitatory synapses can reverse their roles, with slow inhibition favoring synchronization Golomb and Rinzel 1993; Terman et al., 1998; Elson et al., 2002. While predicting the exact rates of synaptic activation inducing such critical transitions in conductance-based spiking models may be challenging, analytically tractable QIF networks offer promising avenues for such exploration.

Toward this goal, this paper investigates a finite-size network of QIF neurons globally connected via a general kernel function that governs both synaptic activation and synaptic time delay. We analytically illustrate how the shape of the kernel function impacts neuron interaction, significantly altering the microscopic and macroscopic behavior of QIF networks representing Type I neuron populations. This is achieved by reducing QIF networks and their phase analog, theta neuron networks, to the Kuramoto-Sakaguchi (KS) model. Here, oscillator frequencies, coupling strength, and the Sakaguchi phase lag parameter are determined by the pulse profile’s first and second terms in the Fourier expansion. We conduct this reduction under the weak coupling assumption, utilizing the intermediate step of representing the QIF network as a generalized Winfree model, subsequently reduced to the KS model.

In our recent study Munyayev et al. (2023), we elucidated the qualitative connection between the dynamics of QIF networks incorporating synaptic dynamics and neuronal refractoriness, and the second-order Kuramoto model with high-order mode coupling. Here, we use multi-scale analysis to derive exact relationships between the QIF network with an arbitrary synaptic activation function and the KS model. Specifically, we establish explicit conditions on the parameters of the general kernel function that lead to critical transitions, determining whether the coupling is attractive or repulsive. Consequently, these conditions dictate the emergence of stable synchronization or nonstationary generalized splay states Berner et al. (2021b) and cyclops states Munyayev et al. (2023). Our analysis reveals alternating stability regions for network synchronization, dependent on both the (slow) synaptic activation and time delay. With some important caveats, this finding can be interpreted as an analogous stability criterion for synchronization in time-delayed phase oscillator networks Earl and Strogatz (2003).

Our approach serves as a connecting link between two alternative methodologies for describing macroscopic dynamics: QIF networks and theta neurons, and Winfree-type models Pazó and Montbrió, 2014; Gallego et al., 2017; Montbrió and Pazó, 2018; Pazó et al., 2019; Pazó and Gallego 2020; Manoranjani et al., 2021; Bick et al., 2020. Our KS model reduction of the generalized Winfree model with a general synaptic activation function can be seen as an extension of the work Montbrió and Pazó (2018), where a two-population Kuramoto model was derived from a network of Winfree oscillators featuring a feedback loop between fast excitation and slow inhibition.

The structure of this paper is outlined as follows. Section 2 presents the QIF network model, its theta neuron equivalent, and the general synaptic activation function. Section 3 details transforming the theta neuron model into the generalized Winfree model. We expand the pulse profile as a Fourier series and further simplify the model to the KS model using weak coupling-enabled averaging techniques. Section 4 focuses on a specific example of synaptic activation, presenting a class of kernel functions. We establish exact conditions determining whether the synaptic coupling is attractive, promoting synchronization, or repulsive, favoring splay and cyclops states. Section 5 offers numerical validation of the derived conditions and presents a comparison between the dynamics of the QIF network, the theta neuron model, and the reduced KS model. We demonstrate that the KS model accurately predicts firing rates and times, capturing the emergence of synchronization, weakly stable cyclops states, and non-stationary regimes. Section 6 contains concluding remarks and discussions.

2 The general QIF network and its theta neuron representation

Physiologically, excitable neurons are commonly categorized into two types. We focus here on Type I neurons, a group encompassing cortical excitatory pyramidal neurons. When subjected to a sufficiently large input stimulus, these neurons exhibit action potentials at an arbitrarily low rate, signaling the disappearance of a resting state through a saddle-node bifurcation. The canonical model used to describe Type I neurons is the QIF neuron model, which characterizes neurons’ dynamics near the spiking threshold Izhikevich (2007).

This study investigates a globally coupled network of $N$ QIF neurons interacting through chemical synapses. Each neuron’s microscopic state is characterized by its individual membrane potential $v_n$, governed by the following ordinary differential equation Izhikevich (2007):

$\begin{array}{ccc}\hfill {\dot{v}}_n& =v_n^2+{\eta }_n+\varkappa S\left(t\right)\hfill & \hfill \text{if\hspace{0.17em}}{v}_n<v_{th},\\ \hfill v_n& =v_r\hfill & \hfill \text{if\hspace{0.17em}}{v}_n\ge v_{th}.\end{array}(1)$

Here, ${\eta }_n$, $n=\mathrm{1,2},\dots ,N$ represents external constant currents applied to neurons, $\varkappa$ is a common synaptic weight controlling the total strength of synaptic inputs, and $S\left(t\right)$ is a time-varying input drive. When the membrane potential $v_n$ of the $n$th neuron reaches the threshold value $v_{th}$, the neuron generates a spike, and its voltage resets to $v_r$. In the absence of the input drive $(S\left(t\right)=0)$, the intrinsic applied current ${\eta }_n=0$ places the corresponding $n$th neuron at a saddle-node bifurcation, marking the onset of periodic firing. Thus, if ${\eta }_n<0$, the neuron is in the excitable regime, while if ${\eta }_n>0$, it is in the oscillatory regime.

The last term on the right-hand side of (Eq. 1) represents synaptic interactions characterized by the coupling strength $\varkappa$ of the global synaptic drive. In general, the mean population synaptic activity $S\left(t\right)$ can be expressed by the following recurrent input equation:

$S\left(t\right)=\frac{1}{N}{\displaystyle \sum _{n^{\prime }=1}^N}{\displaystyle \sum _{m\backslash t_{n^{\prime }}^m<t}}\underset{-\infty }{\overset{t}{\int }}d\tilde{t}G\left(\left(t-\tilde{t}\right)/\tau \right)\delta \left(\tilde{t}-t_{n^{\prime }}^m\right).(2)$

This equation accounts for relaxation processes and describes a specific type of neuron activation and its sensitivity to stimuli from other cells, including signal duration and post-spike latency. Here, $t_{n^{\prime }}^m$ denotes the time of the $m$th spike of the $n^{\prime }$th neuron, $\delta \left(t\right)$ represents the Dirac delta function, and $G\left(t/\tau \right)$ is the normalized synaptic activation caused by a single presynaptic spike with a time scale $\tau$. Notably, the integral transformation with the kernel $G\left(t/\tau \right)$ acts as a low-pass filter.

The QIF-neuron model (Eq. 1) describes the membrane potential $v_n$ and operates as a hybrid dynamical system, incorporating instantaneous resets to a base value $v_r$ upon spike emission. While this formulation provides a direct physical interpretation, discontinuities can pose challenges for certain applications. Fortunately, a smooth change of coordinates exists, transforming the QIF-neuron dynamics into a space where the membrane potential $v_n$ is represented by a phase variable ${\theta }_n$ on the unit circle. This representation captures nonlinear spike-generating mechanisms of Type I neurons, ensuring smooth solutions within a compact domain. In the limit $v_{th}=-v_r\to \infty$, the transformation $v_n\left(t\right)=\tan \left({\theta }_n\left(t\right)/2\right)$ Ermentrout and Kopell (1986) converts the membrane potential description of the QIF-neuron model (Eq. 1) into a canonical theta neuron model of a population of Type I neurons coupled by an excitatory or inhibitory synaptic drive. In this case, each neuron’s dynamics is governed by the following equation:

$\frac{d{\theta }_n}{dt}=\left(1-\cos {\theta }_n\right)+{\eta }_n\left(1+\cos {\theta }_n\right)+\varkappa \left(1+\cos {\theta }_n\right)S\left(t\right),(3)$

where $n=1,\dots ,N$ represents the index of the $n$th neuron, and its state is characterized by the phase angle ${\theta }_n$. We assume that the constant excitability parameters ${\eta }_n$, akin to fixed input currents, have slightly different values across the network elements. Moreover, we consider ${\eta }_n>0$, placing each neuron in an oscillatory regime indicative of periodic spiking. A neuron is deemed to spike or produce an action potential when ${\theta }_n$ crosses $\pi$ while increasing. Consequently, in addition to a common external input ${\eta }_n$, each cell receives stimulation from other cells. Thus, the neurons are recurrently coupled via synaptic current $S\left(t\right)$.

The last term on the right-hand side of (Eq. 3) accounts for chemical interactions among neurons. The coupling strength $\varkappa$ is assumed to be uniform across all neurons. We model the synaptic activity $S\left(t\right)$ acting on a neuron with the following expression:

$S\left(t\right)=\underset{-\infty }{\overset{t}{\int }}d\tilde{t}G\left(\left(t-\tilde{t}\right)/\tau \right)\left[\frac{1}{N}{\displaystyle \sum _{n=1}^N}{P}_{\nu }\left({\theta }_n\left(\tilde{t}\right)\right)\right],(4)$

where the function

$P_{\nu }\left({\theta }_n\right)={\mathit{\text{p}}}_{\nu }{\left(1-\cos {\theta }_n\right)}^{\nu },\nu \in \mathbb{N}(5)$

determines the shape of the pulsatile chemical synapse. The positive integer parameter $\nu$ controls the sharpness of $P_{\nu }({\theta }_n)$, with higher values yielding sharper peaks. Note that as $\nu \to \infty$, the smooth profile $P_{\nu }({\theta }_n)$ converges to $P_{\infty }({\theta }_n)$, effectively representing $\delta$-pulses so that

$P_{\nu }\left({\theta }_n\right)\stackrel{\nu \to \infty }{\to }{P}_{\infty }\left({\theta }_n\right)=2{\text{C}}_{2\pi }\left({\theta }_n-\pi \right)=2{\displaystyle \sum _{k=-\infty }^{+\infty }}\delta \left({\theta }_n-\pi -2\pi k\right),(6)$

where ${\text{C}}_T(\cdot )$ denotes the Dirac comb with period $T$. In this limit and under the assumption $v_{th}=-v_r\to \infty$, the theta neuron model (Eqs 3, 4) fully aligns with the QIF-neuron model (Eqs 1, 2). It is noteworthy that these models exhibit an unconventional characteristic: in contrast to conductance-based models, inhibitory $\delta$-pulse coupling $(\varkappa <0)$ promotes synchronization, thus being considered attractive, while excitatory $\delta$-pulse coupling $(\varkappa >0)$ is repulsive Izhikevich (2007).

Note that the limiting $\delta$-pulse case with $P_{\infty }({\theta }_n)$ offers a convenient framework in which function $G(t)$ controls synaptic activation and deactivation after the instantaneous occurrence of a presynaptic action potential. The use of $P_{\infty }({\theta }_n)$ simplifies the analytical expressions in the KS model to be more explicit and manageable (see next section). However, the general Eqs 4, 5 for synaptic input $S(t)$ provide a more comprehensive and accurate description that captures the nuances of synaptic transmission processes. The release of neurotransmitters from the presynaptic neuron that induces a chemically activated synaptic current is non-instantaneous. Thus, the general function $P_{\nu }({\theta }_n)$ more adequately describes this gradual neurotransmitter release that takes place as the presynaptic neuron state approaches a spike activation threshold.

In this work, we primarily focus on the pulse shape defined by (Eq. 5), originally proposed in Ariaratnam and Strogatz (2001) and widely adopted in recent studies of pulse-coupled phase oscillators O’Keeffe and Strogatz 2016; Pazó and Montbrió, 2014; Gallego et al., 2017; Pazó et al., 2019; Bick et al., 2020 and populations of theta neurons Luke et al., 2013; So et al., 2014; Laing 2014, Laing 2015; Montbrió et al., 2015; Pazó and Montbrió, 2016; Chandra et al., 2017; Goel and Ermentrout 2002; Bick et al., 2020; Pietras et al., 2023. However, our approach is directly applicable to alternative pulse shapes satisfying common properties such as unimodality, normalization, symmetry, and localization around ${\theta }_n=\pi$ and considered in previous studies Gallego et al., 2017; Pietras et al., 2023.

In the following, we explore how the shape of the kernel function $G\left(t/\tau \right)$, governing synaptic activation, impacts the network’s collective behavior. To tackle this analytically, we will demonstrate that the model (Eqs 3, 4) consideration can be reduced to a more analytically-tractable KS model. The process involves several key steps: firstly, leveraging the assumption ${\eta }_n>0$, we introduce an alternative phase representation for (Eqs 3, 4). Subsequently, we will derive the KS model of phase oscillators by employing multiple time scale analyses in weak chemical synaptic coupling scenarios.

3 Deriving the KS model from the theta-neuron model: an asymptotic analysis

We begin by assuming that each neuron operates within an oscillatory regime in the absence of interaction, i.e., ${\eta }_n>0$. Hence, we can use a dynamical variable transformation:

$2\tan \left({\theta }_n/2\right)=\mathrm{\Omega }\tan \left({\phi }_n/2\right),(7)$

where $\mathrm{\Omega }$ is an unknown parameter to be determined subsequently. It is notable that $\mathrm{\Omega }$ is close to $2\sqrt{⟨{\eta }_n⟩}$, where $⟨{\eta }_n⟩$ denotes the mean value of the distributed external constant currents ${\eta }_n$. However, a more precise determination of $\mathrm{\Omega }$ is feasible for finite-size networks, as will be shown later in this section.

The transformation (Eq. 7) transitions the model (Eqs 3, 4) to an alternative phase representation:

$\frac{d{\phi }_n}{dt}=\mathrm{\Omega }+\frac{2\varkappa }{\mathrm{\Omega }}\left(1+\cos {\phi }_n\right)\left(\frac{{\eta }_n-{\mathrm{\Omega }}^2/4}{\varkappa }+S\left(t\right)\right),(8)$

$S\left(t\right)=\underset{-\infty }{\overset{t}{\int }}d\tilde{t}G\left(\left(t-\tilde{t}\right)/\tau \right)\left[\frac{1}{N}{\displaystyle \sum _{n=1}^N}{Q}_{\nu }\left({\phi }_n\left(\tilde{t}\right)\right)\right],(9)$

where

$Q_{\nu }\left({\phi }_n\right)=P_{\nu }\left(2\arctan \left(\mathrm{\Omega }\tan \left[{\phi }_n/2\right]/2\right)\right).(10)$

This representation remains consistent with the original description; notably, ${\phi }_n\in [-\pi ,\pi ]$, and the occurrence of spikes for the $n$th neuron is still defined by ${\phi }_n$ crossing $\pi$. However, given that in the model (Eqs 3, 4), all cells receive constant external inputs ${\eta }_n>0$, and all units operate within the oscillatory regime, each phase ${\phi }_n$ uniformly rotates in the absence of network interaction. Hence, the representation (Eqs 8, 9) proves more convenient for subsequent analysis. Thus, we derive the governing equation for the introduced dynamical variables ${\phi }_n$ that may be viewed as the generalized Winfree model, which, in turn, can be reduced to the KS model. Further elaboration on this approach’s derivation and technical intricacies are presented below.

For further analysis, it is convenient to expand the symmetric pulse $Q_{\nu }\left({\phi }_n\right)$ into a Fourier series with respect to ${\phi }_n$. This expansion takes the following form:

$Q_{\nu }\left({\phi }_n\right)={\displaystyle \sum _{\ell =-\infty }^{+\infty }}{{\mathcal{Q}}_{\nu }}_{\ell }{e}^{i\ell {\phi }_n},{{\mathcal{Q}}_{\nu }}_{\ell }=\frac{1}{2\pi }\underset{-\pi }{\overset{+\pi }{\int }}d\phi e^{-i\ell \phi }{Q}_{\nu }\left(\phi \right),{\mathcal{Q}}_{\nu ,-\ell }={\mathcal{Q}}_{\nu \ell }\in \mathbb{R}.(11)$

The coefficients ${{\mathcal{Q}}_{\nu }}_{\ell }$ in this series can be expressed analytically as follows:

$\begin{array}{cc}\hfill {{\mathcal{Q}}_{\nu }}_{\ell }& =\frac{{\mathit{\text{p}}}_{\nu }{2}^{\nu }}{\pi }{\displaystyle \sum _{m=0}^{\ell }}{\left(-1\right)}^m{C}_{2\ell }^{2m}{\left(\frac{\left|\mathrm{\Omega }\right|}{2}\right)}^{2\left(\ell -m\right)+1}\frac{\mathrm{\Gamma }\left(\nu +m+\frac{1}{2}\right)\mathrm{\Gamma }\left(\ell -m+\frac{1}{2}\right)}{\mathrm{\Gamma }\left(\nu +\ell +1\right)}\hfill \\ \hfill & {\times }_2{F}_1\left(\ell +1,\ell -m+\frac{1}{2};\nu +\ell +1;1-\frac{{\mathrm{\Omega }}^2}{4}\right),\hfill \end{array}$

where $C_{2\ell }^{2m}$ represents combinations, $\mathrm{\Gamma }(\cdot )$ denotes the gamma function, ${\mathrm{a}\mathrm{n}\mathrm{d}}_2{F}_1(\cdot ,\cdot ;\cdot ;\cdot )$ is the hypergeometric function. Specifically, for the first two Fourier series coefficients ${{\mathcal{Q}}_{\infty }}_0$ and ${{\mathcal{Q}}_{\infty }}_1$ of the symmetric pulse $Q_{\infty }({\phi }_n)$ in (Eq. 10), we obtain:

${{\mathcal{Q}}_{\nu }}_0={\frac{\left|\mathrm{\Omega }\right|}{2\pi }}_2{F}_1\left(1,\frac{1}{2};\nu +1;1-\frac{{\mathrm{\Omega }}^2}{4}\right),(12)$

${{\mathcal{Q}}_{\nu }}_1=\frac{\left|\mathrm{\Omega }\right|}{4\pi \left(\nu +1\right)}\left[{\frac{{\mathrm{\Omega }}^2}{4}}_2{F}_1\left(2,\frac{3}{2};\nu +2;1-\frac{{\mathrm{\Omega }}^2}{4}\right)-{\left(2\nu +1\right)}_2{F}_1\left(2,\frac{1}{2};\nu +2;1-\frac{{\mathrm{\Omega }}^2}{4}\right)\right].(13)$

Noteworthy, in the limit $\nu \to \infty$, i.e., for $Q_{\infty }\left({\phi }_n\right)$, all coefficients in the Fourier series (Eq. 11) converge to the same value ${{\mathcal{Q}}_{\infty }}_{\ell }={\left(-1\right)}^{\ell }\left|\mathrm{\Omega }\right|/2\pi$.

We proceed by assuming that the synaptic coupling is weak, allowing us to express it as $2\varkappa /\mathrm{\Omega }=\epsilon \kappa$, where $\epsilon \ll 1$ is a small parameter. Similarly, we assume that the deviations of external inputs ${\eta }_n$ from the value ${\mathrm{\Omega }}^2/4$ are small, i.e., ${\eta }_n-{\mathrm{\Omega }}^2/4=\epsilon \mathrm{\Omega }{\sigma }_n/2$.

These assumptions enable a multiple-time scale analysis. To facilitate this analysis, we introduce a separation of time scales:

$t_k={\epsilon }^{k}t,k=\mathrm{0,1},\dots ,\infty (14)$

and represent each phase variable, ${\phi }_n(t)$, as an asymptotic series with respect to the small parameter $\epsilon$:

${\phi }_n\left(t\right)={\displaystyle \sum _{k=0}^{\infty }}{\epsilon }^k{\phi }_n^{\left(k\right)}\left(t_0,t_1,t_2,\dots \right).(15)$

Substituting the series (Eq. 15) with times (Eq. 14) and (Eq. 11) into (Eqs 8, 9), and considering the zeroth order in $\epsilon$, we obtain for ${\phi }_n^{(0)}(t_0,t_1,t_2,\dots )$:

${\phi }_n^{\left(0\right)}\left(t_0,t_1,t_2,\dots \right)=\mathrm{\Omega }{t}_0+{\varphi }_n\left(t_1,t_2,\dots \right)(16)$

and, taking into account (Eq. 16), for $S(t_0,t_1,t_2,\dots )$ we arrive at

$\begin{array}{cc}\hfill S\left(t_0,t_1,t_2,\dots \right)& =\tau \underset{0}{\overset{+\infty }{\int }}d\xi G\left(\xi \right)\left[\frac{1}{N}{\displaystyle \sum _{n^{\prime }=1}^N}{Q}_{\nu }\left(\mathrm{\Omega }{t}_0-\mathrm{\Omega }\tau \xi +{\varphi }_{n^{\prime }}\left(t_1-\epsilon \tau \xi ,t_2-{\epsilon }^2\tau \xi ,\dots \right)\right.\right.\hfill \\ \hfill & +\left.\left.\epsilon {\phi }_{n^{\prime }}^{\left(1\right)}\left(t_0-\tau \xi ,t_1-\epsilon \tau \xi ,\dots \right)+\cdots \right)\right]\approx \frac{1}{N}{\displaystyle \sum _{n^{\prime }=1}^N}{\displaystyle \sum _{\ell =-\infty }^{+\infty }}{\mathcal{G}}_{\ell }{{\mathcal{Q}}_{\nu }}_{\ell }{e}^{i\ell \mathrm{\Omega }{t}_0+i\ell {\varphi }_{n^{\prime }}\left(t_1,t_2,\dots \right)},\hfill \end{array}(17)$

where each corresponding complex coefficient ${\mathcal{G}}_{\ell }$ is determined as follows:

${\mathcal{G}}_{\ell }=\tau \underset{0}{\overset{+\infty }{\int }}d\xi G\left(\xi \right)e^{-i\ell \mathrm{\Omega }\tau \xi }.(18)$

In (Eq. 16), the first term $\mathrm{\Omega }{t}_0$ describes the fast, free-running rotation of period $2\pi /\mathrm{\Omega }$, while the slow phase drifts induced by synaptic interaction are characterized by the set of slow variables ${\varphi }_n\left(t_1,t_2,\dots \right)$. Hence, ${\varphi }_n\left(t_1,t_2,\dots \right)$ can be considered constant over time scales comparable to the period of the corresponding fast rotation. Consequently, the standard averaging method can be applied to derive the KS model corresponding to (Eqs 3, 4). Notably, in this case, the Sakaguchi phase shift emerges naturally due to the complex nature of the coefficient ${\mathcal{G}}_{\ell }$ determined in (Eq. 18).

In accordance with the averaging procedure after substituting expression (Eq. 17) into (Eq. 8), the next step of our asymptotic approach involves considering all terms that are $o(\epsilon )$ and, in the first order in $\epsilon$, we obtain a set of equations for ${\phi }_n^{(1)}(t_0,t_1,\dots )$.

To eliminate the secular terms that grow without bounds as $t_0\to \infty$, we impose the conditions

$\frac{\partial {\varphi }_n}{\partial t_1}={\sigma }_n+\kappa {{\mathcal{Q}}_{\nu }}_0+\frac{\kappa |{\mathcal{G}}_1|{{\mathcal{Q}}_{\nu }}_1}{N}{\displaystyle \sum _{n^{\prime }=1}^N}\cos \left({\varphi }_{n^{\prime }}-{\varphi }_n+\arg \left({\mathcal{G}}_1\right)\right).(19)$

This yields a solution for ${\phi }_n^{(1)}(t_0,t_1,\dots )$ without the secular terms. Note that (Eq. 19) measures the rate of change of ${\varphi }_n(t_1,t_2,\dots )$ with respect to the slow time scale $t_1$. Finally, taking into account the relation $d{\varphi }_n/dt\approx \epsilon \partial {\varphi }_n/\partial t_1$, we find that the dynamics of the slow phases ${\varphi }_n\left(t\right)$ is approximately described by the KS model:

$\frac{d{\varphi }_n}{dt}={\omega }_n+\frac{K}{N}{\displaystyle \sum _{n^{\prime }=1}^N}\sin \left({\varphi }_{n^{\prime }}-{\varphi }_n-\alpha \right),(20)$

where

${\omega }_n=\epsilon {\sigma }_n+\epsilon \kappa {{\mathcal{Q}}_{\nu }}_0=2\left({\eta }_n-{\mathrm{\Omega }}^2/4+\varkappa {{\mathcal{Q}}_{\nu }}_0\right)/\mathrm{\Omega },(\mathrm{21a})$

$K=\epsilon \kappa |{\mathcal{G}}_1|{{\mathcal{Q}}_{\nu }}_1=2\varkappa |{\mathcal{G}}_1|{{\mathcal{Q}}_{\nu }}_1/\mathrm{\Omega },(\mathrm{21b})$

$\alpha =-\arg \left({\mathcal{G}}_1\right)-\pi /2.(\mathrm{21c})$

To determine the unknown parameter $\mathrm{\Omega }$, we set the mean value of the intrinsic frequencies ${\omega }_n$ in the KS model (Eq. 20) to zero. Hence, the value of $\mathrm{\Omega }$ can be found by solving the nonlinear algebraic equation:

$⟨{\eta }_n⟩-{\mathrm{\Omega }}^2/4+\varkappa {\mathcal{Q}}_{\nu 0}\left(\mathrm{\Omega }\right)=0.(22)$

This choice of the optimal value of parameter $\mathrm{\Omega }$ yields a better quantitative match between the numerical simulation results of the theta neuron model (Eqs 3, 4) and the KS model (Eq. 20), compared to the conventional choice of $2\sqrt{⟨{\eta }_n⟩}$. In the limiting case $\nu \to +\infty$, where the shape of the pulsatile chemical synapse $P_{\infty }({\theta }_n)$ is determined by Eq. 6, the first two Fourier series coefficients ${{\mathcal{Q}}_{\infty }}_0$ and ${{\mathcal{Q}}_{\infty }}_1$ in (Eqs 12, 13) of the symmetric pulse $Q_{\infty }({\phi }_n)$ in Eqs 8, 9 converge to ${{\mathcal{Q}}_{\infty }}_0=-{{\mathcal{Q}}_{\infty }}_1=|\mathrm{\Omega }|/2\pi$ and, hence, Eq. 22 for $\mathrm{\Omega }$ can be solved analytically. This yields the simpler expressions for the parameters ${\omega }_n$, $K$ and $\alpha$ of the KS model (Eq. 20) with coefficients (Eqs 21a, 21b, 21c):

${\omega }_n=\frac{2\pi \left({\eta }_n-\langle {\eta }_n\rangle \right)}{\varkappa +\sqrt{4{\pi }^2\langle {\eta }_n\rangle +{\varkappa }^2}},(\mathrm{23a})$

$K=-\varkappa \left|{\mathcal{G}}_1\right|/\pi ,(\mathrm{23b})$

$\alpha =-\arg \left({\mathcal{G}}_1\right)-\pi /2,(\mathrm{23c})$

where the complex coefficient ${\mathcal{G}}_1$ is determined as follows:

${\mathcal{G}}_1=\tau \underset{0}{\overset{+\infty }{\int }}d\xi G\left(\xi \right)e^{-i\left(\varkappa +\sqrt{4{\pi }^2⟨{\eta }_n⟩+{\varkappa }^2}\right)\tau \xi /\pi }.(24)$

Note that oscillator frequencies ${\omega }_n$, coupling strength $K$, and the Sakaguchi phase lag parameter $\alpha$ in the KS model (Eq. 20) are explicitly defined through by the pulse profile’s first and second Fourier series terms ${\mathcal{Q}}_{\nu 0}$ and ${\mathcal{Q}}_{\nu 1}$. In particular, the expressions (Eq. 23) clearly indicate that for the $\delta$-pulses, the sign of the coupling strength $K$ in the KS model is solely determined by and is opposite to the sign of the coupling $\kappa$ in the original QIF model. The following section will show that this property carries over to the general finite-width pulses defined by (Eq. 5). For a specific class of the activation function $G(t/\tau )$, we will also demonstrate that increasing the pulse width decreases the coupling strength $K$ and practically does not affect the Sakachichi phase lag parameter $\alpha$. The direct dependence on the properties of synaptic activation enables a straightforward assessment of the role of synaptic interactions in critical phase transitions and dynamics. Specifically, the coupling strength $K$ and Sakachichi parameter $\alpha$ directly reflect the impact of synaptic activation on critical phase transitions and the synchronization behavior of the neuron population. In particular, determining whether the coupling in the theta neuron model (Eqs 3, 4) is attractive or repulsive presents a challenge due to the complexity of the system, whereas the KS model (Eq. 20) makes this process straightforward. In the subsequent section, we delve into this process, focusing on a particular synaptic activation profile.

4 The role of synaptic profile: a combined effect of activation, deactivation, and time delays

To demonstrate the important effects arising from the specific selection of the shape of a “low-pass filter” in synaptic activation, we examine the following class of kernel functions:

$G\left(t/\tau \right)=\frac{{\left(t/\tau \right)}^q{e}^{-t/\tau }}{q!\tau }H\left(t/\tau \right),(25)$

where $H(t/\tau )$ is the Heaviside step function and $q$ represents an integer parameter pivotal to the network dynamics, as will become evident subsequently. The chosen kernel corresponds to the Green’s function for a non-homogeneous linear differential equation of order $(q+1)$. In this case, the equation for the mean synaptic activity $S(t)$ is generated by the $(q+1)$th order linear differential operator $\widehat{\mathcal{L}}={\left(\tau d/dt+1\right)}^{q+1}$ and contains an external signal that represents an average profile of spike pulses:

${\left(\tau \frac{d}{dt}+1\right)}^{q+1}S\left(t\right)=\left\{\begin{array}{c}\frac{1}{N}{\displaystyle \sum _{n^{\prime }=1}^N}{\displaystyle \sum _{m\backslash t_{n^{\prime }}^m<t}}\delta \left(t-t_{n^{\prime }}^m\right),\hfill \\ \frac{1}{N}{\displaystyle \sum _{n^{\prime }=1}^N}{P}_{\nu }\left({\theta }_{n^{\prime }}\left(t\right)\right),\hfill \\ \frac{1}{N}{\displaystyle \sum _{n^{\prime }=1}^N}{Q}_{\nu }\left({\phi }_{n^{\prime }}\left(t\right)\right).\hfill \end{array}\right.(26)$

The source term on the right-hand side of (Eq. 26) is presented in three interchangeable forms corresponding to the QIF model (Eqs 1, 2), theta neuron (Eqs 3, 4), and their averaged representation through the KS model (Eq. 20). This term can be interpreted as the population firing rate, which induces a post-synaptic current in response to the arrival of spikes. In the limit $\tau \to 0$, one might assume that the interaction between neurons becomes instantaneous. However, if the characteristic time scale $\tau$ of a post-synaptic response is not negligibly small, it becomes imperative to take into account the individual adaptation dynamics of the synaptic variable $S(t)$, encompassing its activation, deactivation, and time delay.

To describe the adaption dynamics of $S(t)$, it is common to assume that the synaptic variable $S(t)$ follows the first-order Devalle et al., 2017; Bick et al., 2020; Afifurrahman et al., 2020; Afifurrahman et al., 2021; Pietras et al., 2023 or second-order ordinary differential equation Mohanty and Politi 2006; Zillmer et al., 2007; Bolotov et al., 2016; Chen et al., 2017, corresponding to $q=0$ and $q=1$ in (Eq. 26), respectively.

In the case $q=0$, the mean population synaptic activity $S(t)$ is governed by the standard relaxation rule. Here, when the $n^{\prime }$th neuron fires at time $t_{n^{\prime }}^m$, generating the $m$th spike in the form of the Dirac delta function, the mean activity $S(t)$ instantaneously changes and subsequently decays exponentially in the absence of further firings. The parameter $\tau$ acts as a synaptic time constant.

For $q=1$, when the $n^{\prime }$th neuron fires at time $t_{n^{\prime }}^m$ and the $m$th Dirac delta pulse is generated, the variable $S(t)$ is augmented by the function $G\left((t-t_{n^{\prime }}^m)/\tau \right)/N$ defined by (Eq. 25) with $q=1$, coinciding with the so-called alpha-function pulse Mohanty and Politi 2006; Zillmer et al., 2007; Bolotov et al., 2016; Chen et al., 2017. For this alpha-pulse created by a spiking neuron, $\tau$ determines both the signal’s width and the time at which it attains its maximum value (Figure 1).

Figure 1

Figure 1. Synaptic dynamics $S\left(t\right)$, defined by (Eq. 4), induced by presynaptic spikes $p\left(t\right)$, taking the form of (A,B) $P_{\infty }\left(t\right)$ as defined by (Eq. 6) and (C,D) $P_{\nu }\left(t\right)$ as defined by (Eq. 5) with $\nu =40$ and linearly increasing phase ${\theta }_1(t)=t$ $(N=1)$. The parameter $\tau =0.18$ is used for all cases.

We extend the argument to arbitrary $q$ and $\tau$ that make the model of synaptic adaptation (Eq. 26) encompass a broad spectrum of realistic biophysical scenarios, ranging from fast non-delayed to slow delayed activation. These scenarios include neurotransmitter release in the synaptic cleft and the opening/closing of postsynaptic ion channels, characterized by distinct time scales such as latency (time delay), rise, and decay times, as reflected in the kernel function (Eq. 25).

Our choice of the kernel function $G\left(t/\tau \right)$ for synaptic activation aligns with the gamma distribution Zwillinger (1992), a continuous probability distribution characterized by two parameters: $\tau >0$, the scale parameter, and $q>0$, the shape parameter. This distribution reaches its maximum value at $t_{\max }=q\tau$, with mean value $t_{\text{mean}}=\left(q+1\right)\tau$, variance ${\sigma }_t^2=\left(q+1\right){\tau }^2$, and skewness ${\gamma }_t=2/\sqrt{q+1}$. The skewness, reflecting the symmetry of the distribution about its mean, is maximal for the exponential case $q=0$ and diminishes for larger values of $q$, indicating increased symmetry. While both $q$ and $\tau$ influence synaptic dynamics, $q$ has a more significant impact on synaptic time delay than $\tau$, whereas $\tau$ predominantly shapes the synaptic activation profile. Figures 1, 2 demonstrate how the parameters $\tau$ and $q$ determine the time profile of the post-synaptic response and its characteristics.

Figure 2

Figure 2. Characteristics of the synaptic dynamic profile for $S\left(t\right)$ (Eq. 4) as a function of $q$. (A) Peak latency, $t_m$, (B) maximum value, $S_m$, and (C) full width at half maximum $(\text{FWHM})$, induced by spikes $P_{\nu }\left(t\right)$ with $\tau =0.1$ and different $\nu$ (cyan markers – $\nu =10$, red markers – $\nu =100$, orange markers – $\nu =1000$).

Towards our objective of deriving explicit conditions for the attractiveness or repulsiveness of synaptic coupling governed by (Eq. 24) with the kernel function (Eq. 25), we calculate the complex coefficient ${\mathcal{G}}_1$ and its modulus and argument from (Eq. 24) as follows:

${\mathcal{G}}_1=\frac{1}{{\left(1+i\mathrm{\Omega }\tau \right)}^{q+1}},\left|{\mathcal{G}}_1\right|=\frac{1}{{\left(1+{\left(\mathrm{\Omega }\tau \right)}^2\right)}^{\left(q+1\right)/2}},\arg \left({\mathcal{G}}_1\right)=-\left(q+1\right)\arctan \left(\mathrm{\Omega }\tau \right),(27)$

where $\mathrm{\Omega }$ is determined from (Eq. 22). While the hypergeometric function ${{\mathcal{Q}}_{\nu }}_1$, a factor defining the coupling strength $K$ in the KS model (Eq. 20) with coefficients (Eq. 27), cannot be expressed via elementary functions (except for the liming $\delta$-pulse case with $\nu \to \infty$), it is evident that ${\mathcal{Q}}_{\nu 1}\le 0$ for $\nu \ge 0$. Consequently, the coupling strength $K=2\varkappa |{\mathcal{G}}_1|{{\mathcal{Q}}_{\nu }}_1/\mathrm{\Omega }$ in the KS model (Eq. 20) for $\mathrm{\Omega }>0$ and the coupling strength $\kappa$ in the QIF model (Eq. 1) have opposite signs. Therefore, for $\varkappa <0$, a positive coupling strength $K$ corresponds to attractive coupling, provided that the Sakaguchi phase lag parameter $\alpha <\pi /2$. According to (Eq. 21), $\cos \alpha >0$ if

$\sin \left[\left(q+1\right)\arctan \left(\tau \mathrm{\Omega }\right)\right]>0.(28)$

Solving the inequality (Eq. 28), we obtain the following $⌊q/2⌋+1$ intervals of parameters that correspond to the attractive coupling in the KS model (Eq. 20) and therefore in the QIF (Eqs 1, 2) and theta-neuron models (Eqs 3, 4) with $\kappa <0$:

$\begin{array}{ccc}\hfill & \tau \mathrm{\Omega }\in D_0\cup \left(\tan \left(\frac{⌊q/2⌋}{q+1}\pi \right),+\infty \right),\hfill & \hfill \text{for\hspace{0.17em}even\hspace{0.17em}}\lfloor q/2\rfloor ,\\ \hfill & \tau \mathrm{\Omega }\in \left(-\infty ,-\tan \left(\frac{⌊q/2⌋}{q+1}\pi \right)\right)\cup D_0,\hfill & \hfill \text{for\hspace{0.17em}odd\hspace{0.17em}}\lfloor q/2\rfloor ,\end{array}(29)$

where $D_0=\left[\bigcup _{n=-⌊q/4⌋}^{⌊\left(q-2\right)/4⌋}\left(\tan \left(\frac{2n}{q+1}\pi \right),\tan \left(\frac{2n+1}{q+1}\pi \right)\right)\right]$.

Figure 3 displays the regions, as defined by (Eq. 29), where the coupling in the KS model is attractive (blue) or repulsive (red). These regions exhibit an alternating pattern as functions of the synaptic time constant $\tau$ and the common external input $\theta$ for a fixed $q$. Increasing the parameter $q$, which primarily controls the time delay, makes the synaptic coupling more sensitive and results in more alternating zones. This alternating pattern of attractive and repulsive coupling resembles the stability criterion for synchronization in time-delayed phase oscillator networks, where the time delay controls the sign of the derivative of the periodic coupling function Earl and Strogatz (2003).

Figure 3

Figure 3. Regions of attractive (blue) and repulsive (red) coupling in the theta neuron model (Eqs 3, 4), corresponding to the KS model regions defined by (Eq. 29). The colors represent the coupling strength $K$ sign as a function of the synaptic time constant $\tau$ and the common external input $\theta$ for a fixed $q$. (A) $q=2$, (B) $q=4$, and (C) $q=12$. Other parameters: $\varkappa =-0.2\pi$, $\nu =20$, ${\eta }_1={\eta }_2=\cdots ={\eta }_N=\eta$. The yellow points $A$, $B$, $C$, and $D$ indicate the parameter values used for numerical simulations of Figures 6–9.

Figure 4 provides further insight into how the synchronization properties of the KS model, critically controlled by the coupling strength $K$ and Sakaguchi parameter $\alpha$, depend on the original parameters of the QIF model. Notably, both $K$ and $\alpha$ exhibit a weak dependence on $\nu$, suggesting that the pulse width has minimal impact on synchronization dynamics, except for the case of biologically irrelevant wide pulses where $\nu$ approaches 0, making the coupling strength $K$ significantly weaker. In contrast, both $K$ and $\alpha$ are highly sensitive to variations in $q$, with the latter producing the alternating regions of attractive and repulsive coupling, as illustrated in Figure 3.

Figure 4

Figure 4. Coupling strength $K$ (A) and Sakaguchi parameter $\alpha$ (B) as functions of synaptic adaptation/delay (parameter $q$) and the finite pulse width (parameter $\nu$). The dependencies are calculated analytically via (Eq. 21) and (Eq. 27) for $\tau =0.15$, $\varkappa =-0.2\pi$, $\bar{\eta }=2$. The coupling strength $K$ increases as the pulse width decreases (via increasing $\nu$) and decreases as synaptic adaptation slows down and experiences larger time delays (via increasing $q$). The Sakaguchi parameter $\alpha$ is highly sensitive to $q$ and only weakly dependent on $\nu$.

In the following, we offer additional evidence supporting the predictive power of the derived KS model. We show numerically that it effectively captures the emergence of robust dynamical regimes like synchronization and more intricate partially synchronized dynamics such as weakly stable cyclops states and non-stationary generalized splay states in both the QIF and theta neuron models.

5 Dynamical equivalence of the models: numerical validation

We conduct numerical computations using a widely accepted fifth-order Runge–Kutta method with a fixed time step of 0.01, providing additional validation for our analytical findings and predictions.

To characterize the dynamical regimes, we utilize both microscopic measures (pairwise phase differences and firing times) and macroscopic indicators such as the first- and second-order complex Kuramoto parameters Daido 1992; Skardal et al., 2011:

$R_{\ell }\left(t\right)=\frac{1}{N}{\displaystyle \sum _{n^{\prime }=1}^N}{e}^{i\ell {\phi }_{n^{\prime }}}=r_{\ell }{e}^{i{\psi }_{\ell }},(30)$

where $r_{\ell }$ and ${\psi }_{\ell }$, $l=\mathrm{1,2}$ define the magnitude and the phase of the $\ell$th moment Kuramoto order parameter $R_{\ell }(t)$, respectively. The first-order scalar parameter $r_1=|R_1|$ characterizes the degree of phase synchrony with $r_1=1$ corresponding to full phase synchrony. The second-order scalar parameter $r_2=|R_2|$ determines the degree of cluster synchrony, where $r_2=0$ corresponds to generalized splay states Berner et al. (2021b) and their particular case of cyclops states Munyayev et al. (2023). We will also use the firing rate, which measures the average rate at which neurons emit spikes, as another important macroscopic observable to characterize dynamical regimes. We will use the following formula Pietras et al. (2023) in simulations later in the section:

$\rho \left(t\right)=\frac{1}{N}{\displaystyle \sum _{n^{\prime }=1}^N}{\displaystyle \sum _{m\backslash t_{n^{\prime }}^m<t}}\frac{1}{\mathrm{\Delta }t}\underset{t-\mathrm{\Delta }t}{\overset{t}{\int }}d\tilde{t}\delta \left(\tilde{t}-t_{n^{\prime }}^m\right).(31)$

To identify the time steps corresponding to neuron spike events in both the QIF-neuron model and the theta neuron model, we monitored the sign changes of $v_n(t)-v_{th}$ and $\left({\theta }_n(t)\mathrm{mod}2\pi \right)-\pi$, along with their time derivative signs. Subsequently, we determined the spike moment $t_n^m$ using linear interpolation within each time step.

Figures 5, 6 illustrate the perfect correspondence between the emergence of full synchronization and non-stationary generalized splay states in the theta neuron model (Eqs 3, 4) and the KS model (Eq. 20) within the range of attractive coupling (point A in Figure 3) and repulsive coupling (point B in Figure 3), respectively. In the case of full synchronization (Figure 5), the first-order and second-order scalar parameters (Eq. 30), $|R_1|$ and $|R_2|$, converge to unity but cannot reach 1 due to intrinsic parameter mismatch in ${\eta }_n$. Likewise, the first-order scalar parameter, $|R_1|$, associated with the non-stationary generalized splay state oscillates closely around 0 (Figure 6). Figure 7 demonstrates that the KS model maintains its excellent predictive power for the emergence of full synchronization and non-stationary generalized splay states in large networks of 100 and 1,000 theta neurons.

Figure 5

Figure 5. Dynamical equivalence between the theta neuron (Eqs 3, 4) and KS models (Eq. 20), demonstrated via the onset of full synchronization. (A) The evolution of the first $|R_1|$ (solid curves) and second $|R_2|$ (dashed curves) order parameters for the theta neuron (green curves) and KS model (red curves), including the transient period. Initial phases ${\theta }_n$, $n=1,\dots ,N=21$ are uniformly distributed over the interval $[-\pi ;\pi ]$. (B) The colors depict the phase differences ${\theta }_n-{\theta }_{21}$ in the theta neuron model converging to imperfect full synchronization, subject to mismatched parameters ${\eta }_n$ that are uniformly distributed on the segment $\left[\bar{\eta }-\delta \eta /2;\bar{\eta }+\delta \eta /2\right]$, $\bar{\eta }=2.0$, $\delta \eta =6\times 10^{-3}$. (C) Comparison between the dynamics of the phase differences ${\theta }_n-{\theta }_{21}$ for $n=n_1=17$ (thick red curve) and $n=n_2=20$ (thick blue curve) in the theta neuron model and ${\theta }_n-{\theta }_{21}$ recalculated from phases ${\phi }_n$ using the relation (Eq. 7) for $n=n_1=17$ (thin cyan curve) and $n=n_2=20$ (thin red curve) in the KS model. Note the perfect alignment of the phase-difference dynamics in the two models. Parameters: $q=2$, $\varkappa =-0.2\pi$, $\tau =0.5$, $\nu =20$, $\bar{\eta }=2.0$ correspond to point $A$ on Figure 3.

Figure 6

Figure 6. Dynamical equivalence between the theta neuron (Eqs 3, 4) and KS models (Eq. 20), demonstrated via the onset of non-stationary generalized splay state with an oscillating $|R_1|\approx 0$. Notations are as in Figure 5. (A) The evolution of the first $|R_1|$ (solid curves) and second $|R_2|$ (dashed curves) order parameters for the theta neuron (green curves) and KS model (red curves), including the transient period. (B) The colors depict the phase differences θ_n–θ₂₁ in the theta neuron model. (C) Comparison between the dynamics of the phase differences θ_n–θ₂₁ for n = n₁ = 17 (thick red curve) and n = n₁ = 20 (thick blue curve) in the theta neuron model and θ_n–θ₂₁ for n = n₁ = 17 (thin cyan curve) and n = n₁ = 20 (thin red curve) in the KS model. Mismatch parameters ${\eta }_n$ are chosen from a uniform distribution $\left[\bar{\eta }-\delta \eta /2;\bar{\eta }+\delta \eta /2\right]$, $\bar{\eta }=2.0$, $\delta \eta =10^{-3}$. Other parameters: $N=21$, $q=2$, $\varkappa =-0.2\pi$, $\tau =0.8$, and $\nu =20$ correspond to point $B$ on Figure 3 and yield the frequency parameter $\mathrm{\Omega }\approx 2.639$, calculated from (Eq. 22) [not shown].

Figure 7

Figure 7. Large-size networks. Dynamical equivalence between the theta neuron (Eqs 3, 4) and KS models (Eq. 20), demonstrated via the onset of full synchronization (A,B) and non-stationary generalized splay state with an oscillating $|R_1|\approx 0$ (C,D) for $N=100$ (A,C) and $N=\mathrm{1,000}$ (B,D). Notations are as in Figures 5, 6. Parameters: $q=2$, $\varkappa =-0.2\pi$, $\tau =0.5$, $\nu =20$, $\bar{\eta }=2.0$ (A,B); $q=2$, $\varkappa =-0.2\pi$, $\tau =0.8$, and $\nu =20$ (C,D).

Figures 8–10 illustrate the remarkable agreement between cooperative dynamics in the QIF and KS models. Specifically, Figure 8 depicts the onset of full synchronization, as evidenced by synchronized firing rates and times determined via (Eq. 31). The slight discrepancy in the firing times between the QIF and KS models may stem from various sources, such as accumulated numerical errors and the approximate calculation of the frequency parameter $\mathrm{\Omega }$ derived from (Eq. 22) for selecting the KS model parameters. Figure 9 provides evidence for the capability of the KS model to perfectly predict even non-stationary, asynchronous firing in the QIF model. Figure 10 illustrates the predictive power of the derived KS model in discerning stable complex cluster patterns like cyclops states in the QIF model. Introduced in Munyayev et al. (2023), cyclops states are formed by two distinct, coherent clusters, and a solitary oscillator reminiscent of the Cyclops’ eye. While detecting stable cyclops states can be challenging in the QIF model, the KS model provides a more convenient and constructive approach.

Figure 8

Figure 8. Onset of full synchronization in the QIF (Eq. 1) and KS models (Eq. 20). (A) Firing rate and (B) firing times of QIF neurons (cyan curves and round markers) and oscillators of the KS model (black curves and cross markers) with the firing times recorded at ${\theta }_n(t_f)=\pi$. Each row in (B) represents the firing times of a neuron/oscillator. Inset (C) zooms-in on the firing time pattern from (B). Initial conditions are as in Figure 5. Parameters: $N=21$, $q=4$, $\tau =0.15$, $\varkappa =-0.2\pi$, $v_{th}=-v_r=10^5$, $\bar{\eta }=2$, $\delta \eta =6\times 10^{-3}$, $\nu =10^5$ correspond to point $C$ in Figure 3. The firing rate was calculated within a sliding time window of $\delta t=5\times 10^{-2}$.

Figure 9

Figure 9. Diagram similar to Figure 8 showing a nearly perfect match for asynchronous firing rate (A) and firing times (B) in the QIF network (cyan curves and round markers) and the KS model (black curves and cross markers). Parameters: $N=21$, $q=4$, $\tau =0.41$, $\varkappa =-0.2\pi$, $v_{th}=-v_r=10^5$, $\bar{\eta }=2$, $\delta \eta =6\times 10^{-3}$, $\nu =10^5$ correspond to point $D$ in Figure 3. Other notations and settings are as in Figure 8.

Figure 10

Figure 10. Firing rate (A) and firing times (B) of a three-cluster cyclops state in the QIF network (cyan curves and round markers) and the KS model (black curves and cross markers). (C) Snapshots of the cyclops state phase distributions ${\phi }_n$ in the KS model at two time instants. The oscillators’ coloring represents their phase. The cyclops states is formed by a solitary oscillator (red) and two coherent clusters, each composed of 10 oscillators (orange and blue). The initial phases are chosen near a cyclops state. Parameters: $N=21$, $q=2$, $\tau =0.8$, $\varkappa =-0.2\pi$, $v_{th}=-v_r=10^5$, $\bar{\eta }=2$, $\delta \eta =0$, $\nu =10^5$. Other notations and settings are as in Figure 8.

In the numerical calculations of relatively small-size networks of QIF neurons presented in Figures 8–10, we employed the fourth–order Runge–Kutta method using the procedure for identifying the spike events described above. To ensure better consistency between the QIF-neuron model and the theta neuron model, we used sufficiently large values of $v_{th}$ and $v_r$, specifically $v_{th}=-v_r=10^5$. We did not account for the time interval required for the membrane potential to reach $\pm \infty$. However, this approach becomes computationally expensive for simulating large networks of QIF neurons. Therefore, in the numerical calculations presented in Figures 11, 12, we employed the algorithm described in Pazó and Montbrió (2016). This algorithm, based on the Euler method, accounts for the time it takes for the dynamic variable associated with the membrane potential to pass through the singularity $\pm \infty$ after exceeding $v_{th}$ and reach the final value $v_r$. A detailed description of the algorithm and its parameters is available in the Supplemental Material for Pazó and Montbrió (2016). Figures 11, 12 confirm that the main analytical predictions for the dynamics of the QIF model remain valid for large network sizes and regardless of the calculation scheme used. Specifically, Figure 11 validates the prediction of coupling attractiveness at point $C$ in Figure 3 and demonstrates the onset of full synchronization in the QIF network of 1,000 neurons. Similarly, Figure 12 confirms the coupling repulsiveness at point $D$ in Figure 3 and illustrates the emergence of asynchronous dynamics. Notably, the asynchronous firing rate in the 21-node QIF network shown in Figure 9, corresponding to point $D$, may resemble weak coherence patterns. In contrast, its 1,000-node counterpart in Figure 12, with parameters also corresponding to point $D$, exhibits nearly perfect asynchronous dynamics. Our numerous additional simulations of 1,000-node QIF networks for a large set of parameters, including those near the boundary of the alternating regions of attractive and repulsive coupling in Figure 3, further validated the consistency of the predictions for cooperative dynamics and the critical transitions in the QIF networks [not shown].

Figure 11

Figure 11. Firing rate (A) and firing times (B) in the QIF network demonstrating the transition to full synchronization starting from random initial conditions. Parameters $N=\mathrm{1,000}$, $q=4$, $\tau =0.15$, $\varkappa =-0.2\pi$, $v_{th}=-v_r=100$, $\bar{\eta }=2$, and $\delta \eta =6\times 10^{-3}$correspond to point $C$ in Figure 3. The simulations use the algorithm from Pazó and Montbrió (2016), which accounts for the neurons’ refractory time. The mean synaptic activation was calculated within a sliding time window of $\delta t=5\times 10^{-2}$.

Figure 12

Figure 12. Firing rate (A) and firing times (B) in the QIF network accompanying the transition from synchronous to asynchronous dynamics. Parameters $N=\mathrm{1,000}$, $q=4$, $\tau =0.41$, $\varkappa =-0.2\pi$, $v_{th}=-v_r=100$, $\bar{\eta }=2$, and $\delta \eta =6\times 10^{-3}$correspond to point $D$ in Figure 3. Other notations and settings are as in Figure 11.

6 Conclusions

Understanding the influence of synaptic dynamics, including activation rates, deactivation processes, and latency, on collective dynamics in neuronal networks is of significant importance. Considerable advancements have been made in analyzing the role of fast or time-delayed synapses in integrate-and-fire neuron networks. However, there remains a scarcity of analytical studies exploring the influence of slower synaptic dynamics, potentially in the presence of time delays, on controlling critical phase transitions in neuronal networks.

In this paper, we have made substantial contributions to advancing analytical methods in this field. We studied a finite-size network of QIF neurons globally interconnected via a generalized kernel function governing both synaptic activation and time delay. Our analytical exploration demonstrated how the shape of the kernel function profoundly affects neuron interaction, thereby significantly modifying the microscopic and macroscopic behavior of QIF networks. To achieve this, we reduced the QIF and theta neuron network models to the Kuramoto–Sakaguchi model. In this model, oscillator frequencies, coupling strength, and the Sakaguchi phase lag parameter are determined by the Fourier terms of the pulse profile series expansion.

We established exact conditions determining whether synaptic coupling is attractive, fostering synchronization, or repulsive, promoting splay and cyclops states. Furthermore, we demonstrated a remarkable correspondence between the dynamics of the derived KS model and the original QIF and theta neuron models. Specifically, the KS model accurately predicted firing rates and times, capturing the emergence of synchronization, weakly stable cyclops states, and non-stationary regimes in the QIF model. Our reduction approach complements the work by Ratas and Pyragas (2018), which assumed a Lorentzian distribution of the input currents and employed the thermodynamic limit to reduce the QIF model with time-delayed coupling to macroscopic equations that characterize the mean membrane potential, the spiking rate, and the mean synaptic current. The bifurcation analysis of these macroscopic equations, performed in Ratas and Pyragas (2018), revealed alternating parameter regions where the QIF network exhibits macroscopic self-oscillations as a function of input current heterogeneity and time-delayed coupling. While sharing some similarities and goals, our approach is fundamentally different. It reduces finite-size QIF networks with an arbitrary current distribution and arbitrary synaptic activation function $G(t/\tau )$ to the finite-size microscopic KS model. This allows for characterizing fine dynamical patterns, including cluster and cyclops states, which might be out of reach for a macroscopic description. In light of this, our reduction approach to an analytically tractable Kuramoto model holds promise in facilitating constructive analysis of rhythmogenesis in QIF networks. By utilizing the reduced KS model, a variety of methods and analytical machinery can be applied to the analysis of collective dynamics in QIF models, including the constructive selection of complex patterns, such as chimeras, whose existence and emergence might be easier to deduce from the reduced Kuramoto model description. Furthermore, the reduction approach holds the potential for extensions to incorporate synaptic adaptation, Hebbian learning, and a complex network structure by employing node-degree block approximation.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

LAS: Conceptualization, Formal Analysis, Investigation, Methodology, Visualization, Writing–review and editing. VOM: Formal Analysis, Investigation, Methodology, Visualization, Writing–review and editing. MIB: Formal Analysis, Investigation, Methodology, Visualization, Writing–review and editing. GVO: Conceptualization, Investigation, Methodology, Supervision, Writing–review and editing. IB: Conceptualization, Investigation, Methodology, Supervision, Writing–original draft.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This work was supported by the Ministry of Science and Higher Education of the Russian Federation under project No. 0729-2020-0036 (to GO and MB.), the Russian Science Foundation under project No. 22-12-00348 (to VM and LS), the Georgia State University Brains and Behavior Program, the National Science Foundation (United States) under Grants No. CMMI-2009329 and CMMI-1953135 (to IB).

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.

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.

References

Achuthan, S., and Canavier, C. C. (2009). Phase-resetting curves determine synchronization, phase locking, and clustering in networks of neural oscillators. J. Neurosci. 29, 5218–5233. doi:10.1523/JNEUROSCI.0426-09.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Afifurrahman, A., Ullner, E., and Politi, A. (2021). Collective dynamics in the presence of finite-width pulses. Chaos An Interdiscip. J. Nonlinear Sci. 31, 043135. doi:10.1063/5.0046691

CrossRef Full Text | Google Scholar

Afifurrahman, , Ullner, E., and Politi, A. (2020). Stability of synchronous states in sparse neuronal networks. Nonlinear Dyn. 102, 733–743. doi:10.1007/s11071-020-05880-4

CrossRef Full Text | Google Scholar

Ariaratnam, J. T., and Strogatz, S. H. (2001). Phase diagram for the Winfree model of coupled nonlinear oscillators. Phys. Rev. Lett. 86, 4278–4281. doi:10.1103/PhysRevLett.86.4278

PubMed Abstract | CrossRef Full Text | Google Scholar

Belykh, I., De Lange, E., and Hasler, M. (2005). Synchronization of bursting neurons: what matters in the network topology. Phys. Rev. Lett. 94, 188101. doi:10.1103/PhysRevLett.94.188101

PubMed Abstract | CrossRef Full Text | Google Scholar

Belykh, I., and Hasler, M. (2011). Mesoscale and clusters of synchrony in networks of bursting neurons. Chaos An Interdiscip. J. Nonlinear Sci. 21, 016106. doi:10.1063/1.3563581

PubMed Abstract | CrossRef Full Text | Google Scholar

Berner, R., Vock, S., Schöll, E., and Yanchuk, S. (2021a). Desynchronization transitions in adaptive networks. Phys. Rev. Lett. 126, 028301. doi:10.1103/PhysRevLett.126.028301

PubMed Abstract | CrossRef Full Text | Google Scholar

Berner, R., Yanchuk, S., Maistrenko, Y., and Scholl, E. (2021b). Generalized splay states in phase oscillator networks. Chaos An Interdiscip. J. Nonlinear Sci. 31, 073128. doi:10.1063/5.0056664

CrossRef Full Text | Google Scholar

Bick, C., Goodfellow, M., Laing, C. R., and Martens, E. A. (2020). Understanding the dynamics of biological and neural oscillator networks through exact mean-field reductions: a review. J. Math. Neurosci. 10 (9), 9. doi:10.1186/s13408-020-00086-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Bolotov, M., Osipov, G., and Pikovsky, A. (2016). Marginal chimera state at cross-frequency locking of pulse-coupled neural networks. Phys. Rev. E 93, 032202. doi:10.1103/PhysRevE.93.032202

PubMed Abstract | CrossRef Full Text | Google Scholar

Börgers, C., and Kopell, N. (2003). Synchronization in networks of excitatory and inhibitory neurons with sparse, random connectivity. Neural comput. 15, 509–538. doi:10.1162/089976603321192059

PubMed Abstract | CrossRef Full Text | Google Scholar

Brunel, N. (2000). Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons. J. Comput. Neurosci. 8, 183–208. doi:10.1023/a:1008925309027

PubMed Abstract | CrossRef Full Text | Google Scholar

Byrne, Á., Ross, J., Nicks, R., and Coombes, S. (2022). Mean-field models for EEG/MEG: from oscillations to waves. Brain Topogr. 35, 36–53. doi:10.1007/s10548-021-00842-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Chandra, S., Hathcock, D., Crain, K., Antonsen, T. M., Girvan, M., and Ott, E. (2017). Modeling the network dynamics of pulse-coupled neurons. Chaos An Interdiscip. J. Nonlinear Sci. 27, 033102. doi:10.1063/1.4977514

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, B., Engelbrecht, J. R., and Mirollo, R. (2017). Cluster synchronization in networks of identical oscillators with α-function pulse coupling. Phys. Rev. E 95, 022207. doi:10.1103/PhysRevE.95.022207

PubMed Abstract | CrossRef Full Text | Google Scholar

Churchland, P. S., and Sejnowski, T. J. (1992). The computational brain. MIT press.

Google Scholar

Clusella, P., and Montbrió, E. (2024). Exact low-dimensional description for fast neural oscillations with low firing rates. Phys. Rev. E 109, 014229. doi:10.1103/PhysRevE.109.014229

PubMed Abstract | CrossRef Full Text | Google Scholar

Clusella, P., Pietras, B., and Montbrió, E. (2022). Kuramoto model for populations of quadratic integrate-and-fire neurons with chemical and electrical coupling. Chaos An Interdiscip. J. Nonlinear Sci. 32, 013105. doi:10.1063/5.0075285

CrossRef Full Text | Google Scholar

Coombes, S. (2023). Next generation neural population models. Front. Appl. Math. Statistics 9, 1128224. doi:10.3389/fams.2023.1128224

CrossRef Full Text | Google Scholar

Daido, H. (1992). Order function and macroscopic mutual entrainment in uniformly coupled limit-cycle oscillators. Prog. Theor. Phys. 88, 1213–1218. doi:10.1143/ptp.88.1213

CrossRef Full Text | Google Scholar

Devalle, F., Montbrió, E., and Pazó, D. (2018). Dynamics of a large system of spiking neurons with synaptic delay. Phys. Rev. E 98, 042214. doi:10.1103/physreve.98.042214

CrossRef Full Text | Google Scholar

Devalle, F., Roxin, A., and Montbrió, E. (2017). Firing rate equations require a spike synchrony mechanism to correctly describe fast oscillations in inhibitory networks. PLoS Comput. Biol. 13, e1005881. doi:10.1371/journal.pcbi.1005881

PubMed Abstract | CrossRef Full Text | Google Scholar

Earl, M. G., and Strogatz, S. H. (2003). Synchronization in oscillator networks with delayed coupling: a stability criterion. Phys. Rev. E 67, 036204. doi:10.1103/PhysRevE.67.036204

PubMed Abstract | CrossRef Full Text | Google Scholar

Elson, R. C., Selverston, A. I., Abarbanel, H. D., and Rabinovich, M. I. (2002). Inhibitory synchronization of bursting in biological neurons: dependence on synaptic time constant. J. Neurophysiology 88, 1166–1176. doi:10.1152/jn.2002.88.3.1166

PubMed Abstract | CrossRef Full Text | Google Scholar

Ermentrout, B., and Kopell, N. (1986). Parabolic bursting in an excitable system coupled with a slow oscillation. SIAM J. Appl. Math. 46, 233–253. doi:10.1137/0146017

CrossRef Full Text | Google Scholar

Ermentrout, G. B., and Terman, D. H. (2010). Mathematical foundations of neuroscience, 35. Springer Science and Business Media.

Google Scholar

Ernst, U., Pawelzik, K., and Geisel, T. (1995). Synchronization induced by temporal delays in pulse-coupled oscillators. Phys. Rev. Lett. 74, 1570–1573. doi:10.1103/PhysRevLett.74.1570

PubMed Abstract | CrossRef Full Text | Google Scholar

Esnaola-Acebes, J. M., Roxin, A., Avitabile, D., and Montbrió, E. (2017). Synchrony-induced modes of oscillation of a neural field model. Phys. Rev. E 96, 052407. doi:10.1103/PhysRevE.96.052407

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferrara, A., Angulo-Garcia, D., Torcini, A., and Olmi, S. (2023). Population spiking and bursting in next-generation neural masses with spike-frequency adaptation. Phys. Rev. E 107, 024311. doi:10.1103/PhysRevE.107.024311

PubMed Abstract | CrossRef Full Text | Google Scholar

Gallego, R., Montbrió, E., and Pazó, D. (2017). Synchronization scenarios in the Winfree model of coupled oscillators. Phys. Rev. E 96, 042208. doi:10.1103/PhysRevE.96.042208

PubMed Abstract | CrossRef Full Text | Google Scholar

Gast, R., Schmidt, H., and Knösche, T. R. (2020). A mean-field description of bursting dynamics in spiking neural networks with short-term adaptation. Neural Comput. 32, 1615–1634. doi:10.1162/neco_a_01300

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerstner, W., and Kistler, W. M. (2002). Spiking neuron models: single neurons, populations, plasticity. Cambridge University Press.

Google Scholar

Goel, P., and Ermentrout, B. (2002). Synchrony, stability, and firing patterns in pulse-coupled oscillators. Phys. D. Nonlinear Phenom. 163, 191–216. doi:10.1016/s0167-2789(01)00374-8

CrossRef Full Text | Google Scholar

Golomb, D., and Rinzel, J. (1993). Dynamics of globally coupled inhibitory neurons with heterogeneity. Phys. Rev. E 48, 4810–4814. doi:10.1103/physreve.48.4810

PubMed Abstract | CrossRef Full Text | Google Scholar

Gutkin, B. S., Laing, C. R., Colby, C. L., Chow, C. C., and Ermentrout, G. B. (2001). Turning on and off with excitation: the role of spike-timing asynchrony and synchrony in sustained neural activity. J. Comput. Neurosci. 11, 121–134. doi:10.1023/a:1012837415096

PubMed Abstract | CrossRef Full Text | Google Scholar

Izhikevich, E. M. (2007). Dynamical systems in neuroscience. Cambridge, MA: MIT Press.

Google Scholar

Kopell, N., Ermentrout, G., Whittington, M. A., and Traub, R. D. (2000). Gamma rhythms and beta rhythms have different synchronization properties. Proc. Natl. Acad. Sci. 97, 1867–1872. doi:10.1073/pnas.97.4.1867

PubMed Abstract | CrossRef Full Text | Google Scholar

Laing, C. R. (2014). Derivation of a neural field model from a network of theta neurons. Phys. Rev. E 90, 010901. doi:10.1103/PhysRevE.90.010901

PubMed Abstract | CrossRef Full Text | Google Scholar

Laing, C. R. (2015). Exact neural fields incorporating gap junctions. SIAM J. Appl. Dyn. Syst. 14, 1899–1929. doi:10.1137/15m1011287

CrossRef Full Text | Google Scholar

Laing, C. R., and Chow, C. C. (2001). Stationary bumps in networks of spiking neurons. Neural Comput. 13, 1473–1494. doi:10.1162/089976601750264974

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, L., Barreto, E., and So, P. (2020). Synaptic diversity suppresses complex collective behavior in networks of theta neurons. Front. Comput. Neurosci. 14, 44. doi:10.3389/fncom.2020.00044

PubMed Abstract | CrossRef Full Text | Google Scholar

Luke, T. B., Barreto, E., and So, P. (2013). Complete classification of the macroscopic behavior of a heterogeneous network of theta neurons. Neural Comput. 25, 3207–3234. doi:10.1162/NECO_a_00525

PubMed Abstract | CrossRef Full Text | Google Scholar

Manoranjani, M., Gopal, R., Senthilkumar, D., and Chandrasekar, V. (2021). Role of phase-dependent influence function in the Winfree model of coupled oscillators. Phys. Rev. E 104, 064206. doi:10.1103/PhysRevE.104.064206

PubMed Abstract | CrossRef Full Text | Google Scholar

Mizuseki, K., and Buzsaki, G. (2014). Theta oscillations decrease spike synchrony in the hippocampus and entorhinal cortex. Philosophical Trans. R. Soc. B Biol. Sci. 369, 20120530. doi:10.1098/rstb.2012.0530

PubMed Abstract | CrossRef Full Text | Google Scholar

Mohanty, P., and Politi, A. (2006). A new approach to partial synchronization in globally coupled rotators. J. Phys. A Math. General 39, L415–L421. doi:10.1088/0305-4470/39/26/l01

CrossRef Full Text | Google Scholar

Montbrió, E., and Pazó, D. (2018). Kuramoto model for excitation-inhibition-based oscillations. Phys. Rev. Lett. 120, 244101. doi:10.1103/PhysRevLett.120.244101

PubMed Abstract | CrossRef Full Text | Google Scholar

Montbrió, E., and Pazó, D. (2020). Exact mean-field theory explains the dual role of electrical synapses in collective synchronization. Phys. Rev. Lett. 125, 248101. doi:10.1103/PhysRevLett.125.248101

PubMed Abstract | CrossRef Full Text | Google Scholar

Montbrió, E., Pazó, D., and Roxin, A. (2015). Macroscopic description for networks of spiking neurons. Phys. Rev. X 5, 021028. doi:10.1103/physrevx.5.021028

CrossRef Full Text | Google Scholar

Munyayev, V. O., Bolotov, M. I., Smirnov, L. A., Osipov, G. V., and Belykh, I. (2023). Cyclops states in repulsive Kuramoto networks: the role of higher-order coupling. Phys. Rev. Lett. 130, 107201. doi:10.1103/PhysRevLett.130.107201

PubMed Abstract | CrossRef Full Text | Google Scholar

O’Keeffe, K. P., and Strogatz, S. H. (2016). Dynamics of a population of oscillatory and excitable elements. Phys. Rev. E 93, 062203. doi:10.1103/PhysRevE.93.062203

PubMed Abstract | CrossRef Full Text | Google Scholar

Olmi, S., Politi, A., and Torcini, A. (2011). Collective chaos in pulse-coupled neural networks. Europhys. Lett. 92, 60007. doi:10.1209/0295-5075/92/60007

CrossRef Full Text | Google Scholar

Omelchenko, I., Omel’chenko, E., Hövel, P., and Schöll, E. (2013). When nonlocal coupling between oscillators becomes stronger: patched synchrony or multichimera states. Phys. Rev. Lett. 110, 224101. doi:10.1103/PhysRevLett.110.224101

PubMed Abstract | CrossRef Full Text | Google Scholar

Ott, E., and Antonsen, T. M. (2008). Low dimensional behavior of large systems of globally coupled oscillators. Chaos An Interdiscip. J. Nonlinear Sci. 18, 037113. doi:10.1063/1.2930766

PubMed Abstract | CrossRef Full Text | Google Scholar

Pazó, D., and Gallego, R. (2020). The Winfree model with non-infinitesimal phase-response curve: Ott–Antonsen theory. Chaos An Interdiscip. J. Nonlinear Sci. 30, 073139. doi:10.1063/5.0015131

CrossRef Full Text | Google Scholar

Pazó, D., and Montbrió, E. (2014). Low-dimensional dynamics of populations of pulse-coupled oscillators. Phys. Rev. X 4, 011009. doi:10.1103/physrevx.4.011009

CrossRef Full Text | Google Scholar

Pazó, D., and Montbrió, E. (2016). From quasiperiodic partial synchronization to collective chaos in populations of inhibitory neurons with delay. Phys. Rev. Lett. 116, 238101. doi:10.1103/PhysRevLett.116.238101

PubMed Abstract | CrossRef Full Text | Google Scholar

Pazó, D., Montbrió, E., and Gallego, R. (2019). The Winfree model with heterogeneous phase-response curves: analytical results. J. Phys. A Math. Theor. 52, 154001. doi:10.1088/1751-8121/ab0b4c

CrossRef Full Text | Google Scholar

Pietras, B., Cestnik, R., and Pikovsky, A. (2023). Exact finite-dimensional description for networks of globally coupled spiking neurons. Phys. Rev. E 107, 024315. doi:10.1103/PhysRevE.107.024315

PubMed Abstract | CrossRef Full Text | Google Scholar

Pietras, B., Devalle, F., Roxin, A., Daffertshofer, A., and Montbrió, E. (2019). Exact firing rate model reveals the differential effects of chemical versus electrical synapses in spiking networks. Phys. Rev. E 100, 042412. doi:10.1103/PhysRevE.100.042412

PubMed Abstract | CrossRef Full Text | Google Scholar

Pyragas, V., and Pyragas, K. (2022). Mean-field equations for neural populations with q-Gaussian heterogeneities. Phys. Rev. E 105, 044402. doi:10.1103/PhysRevE.105.044402

PubMed Abstract | CrossRef Full Text | Google Scholar

Pyragas, V., and Pyragas, K. (2023). Effect of cauchy noise on a network of quadratic integrate-and-fire neurons with non-cauchy heterogeneities. Phys. Lett. A 480, 128972. doi:10.1016/j.physleta.2023.128972

PubMed Abstract | CrossRef Full Text | Google Scholar

Ratas, I., and Pyragas, K. (2018). Macroscopic oscillations of a quadratic integrate-and-fire neuron network with global distributed-delay coupling. Phys. Rev. E 98, 052224. doi:10.1103/physreve.98.052224

CrossRef Full Text | Google Scholar

Schmidt, H., Avitabile, D., Montbrió, E., and Roxin, A. (2018). Network mechanisms underlying the role of oscillations in cognitive tasks. PLoS Comput. Biol. 14, e1006430. doi:10.1371/journal.pcbi.1006430

PubMed Abstract | CrossRef Full Text | Google Scholar

Schöll, E. (2016). Synchronization patterns and chimera states in complex networks: interplay of topology and dynamics. Eur. Phys. J. Special Top. 225, 891–919. doi:10.1140/epjst/e2016-02646-3

CrossRef Full Text | Google Scholar

Shilnikov, A., Gordon, R., and Belykh, I. (2008). Polyrhythmic synchronization in bursting networking motifs. Chaos An Interdiscip. J. Nonlinear Sci. 18, 037120. doi:10.1063/1.2959850

PubMed Abstract | CrossRef Full Text | Google Scholar

Skardal, P. S., Ott, E., and Restrepo, J. G. (2011). Cluster synchrony in systems of coupled phase oscillators with higher-order coupling. Phys. Rev. E 84, 036208. doi:10.1103/PhysRevE.84.036208

PubMed Abstract | CrossRef Full Text | Google Scholar

So, P., Luke, T. B., and Barreto, E. (2014). Networks of theta neurons with time-varying excitability: macroscopic chaos, multistability, and final-state uncertainty. Phys. D. Nonlinear Phenom. 267, 16–26. doi:10.1016/j.physd.2013.04.009

CrossRef Full Text | Google Scholar

Somers, D., and Kopell, N. (1993). Rapid synchronization through fast threshold modulation. Biol. Cybern. 68, 393–407. doi:10.1007/BF00198772

PubMed Abstract | CrossRef Full Text | Google Scholar

Taher, H., Torcini, A., and Olmi, S. (2020). Exact neural mass model for synaptic-based working memory. PLoS Comput. Biol. 16, e1008533. doi:10.1371/journal.pcbi.1008533

PubMed Abstract | CrossRef Full Text | Google Scholar

Terman, D., Kopell, N., and Bose, A. (1998). Dynamics of two mutually coupled slow inhibitory neurons. Phys. D. Nonlinear Phenom. 117, 241–275. doi:10.1016/s0167-2789(97)00312-6

CrossRef Full Text | Google Scholar

Van Vreeswijk, C., Abbott, L., and Bard Ermentrout, G. (1994). When inhibition not excitation synchronizes neural firing. J. Comput. Neurosci. 1, 313–321. doi:10.1007/BF00961879

PubMed Abstract | CrossRef Full Text | Google Scholar

Zillmer, R., Livi, R., Politi, A., and Torcini, A. (2007). Stability of the splay state in pulse-coupled networks. Phys. Rev. E 76, 046102. doi:10.1103/PhysRevE.76.046102

PubMed Abstract | CrossRef Full Text | Google Scholar

Zwillinger, D. (1992). The handbook of integration. New York: A K Peters/CRC Press.

Google Scholar