Coherent dynamics in soft-threshold integrate-and-fire networks

TL;DR

This paper investigates bifurcations and spatiotemporal patterns in stochastic integrate-and-fire neural networks using a mean-field approach, revealing mechanisms for oscillations, localized activity, and wave phenomena.

Contribution

It introduces a mean-field framework to analyze complex bifurcations and spatiotemporal dynamics in neural networks with delays and structured interactions.

Findings

Uniform oscillations arise via subcritical Hopf bifurcation.

Turing bifurcation leads to localized sustained activity.

Spatiotemporal patterns like standing and traveling waves are observed.

Abstract

We study bifurcations in networks of integrate-and-fire neurons with stochastic spike emission, focusing on the effects of the spatial and temporal structure of the synaptic interactions. Using a deterministic mean-field approximation of the population dynamics, we characterize spatial, temporal, and spatiotemporal patterns of macroscopic activity. In the mean-field theory, synaptic delays give rise to uniform oscillations across the population through a subcritical Hopf bifurcation of the stationary uniform equilibrium. With local excitation and long-range inhibition the network undergoes a Turing bifurcation, resulting in a localized area of sustained activity, or stationary bump. When the coupling has both delays, local inhibition, and long range excitation, the network undergoes a Turing-Hopf bifurcation leading to spatiotemporal dynamics, such as standing and traveling waves. When multiple instabilities are excited, we observe other complex spatiotemporal dynamics. We confirm all these predictions of the mean-field theory in simulations of the underlying stochastic model.

Equations

$$d{v}_i(t)=\frac{dt}{\tau }\left(-v_i(t)+E_i+\frac{1}{N}\sum _{j=1}^N\int J_{ij}(t-s-D)d{n}_j(s)\right)-d{n}_i(t^{+})\left(v_i(t)-r\right)$$

$${\partial }_t\bar{v}(x,t)=-\bar{v}(x,t)+E+{\int }_{-\infty }^{t-D}{\int }_{-\pi }^{\pi }J(x-y,t-s-D)\bar{f}(y,s)dyds-\bar{f}(x,t)\bar{v}(x,t)$$

$$v_0=\{\begin{array}{ll}{v}_Q=E& \text{if }E<1,\\ v_{+}=\left(J_0+\sqrt{J_0^2+4(E-J_0)}\right)/2& \text{if }{J}_0>2+2\sqrt{1-E}\text{ or }{J}_0<2,E>1,\\ v_{-}=\left(J_0-\sqrt{J_0^2+4(E-J_0)}\right)/2& \text{if }{J}_0>2+2\sqrt{1-E}\text{ and }{J}_0>2,E<1,\end{array}$$

$${\partial }_{t}w(x,t)=-2v_{0}w(x,t)+{\int }_{-\infty }^{t-D}{\int }_{-\pi }^{\pi }J(x-y,t-s-D)w(y,s)dyds$$

$${\lambda }_k=-2v_0+e^{-{\lambda }_{k}D}{\int }_0^{\infty }{\hat{J}}_k(s)e^{-{\lambda }_{k}s}ds,k\in \mathbb{Z}$$

$$J(x,t-D)=$$

$$J_{\mathrm{time}}(t-D)=$$

$$J_{\mathrm{space}}(x)=$$

$${\lambda }_k=-2v_0+{\hat{J}}_k{e}^{-{\lambda }_{k}D},$$

$$J_0=2+2\sqrt{1-E},$$

$$\frac{-1}{D}{\cos }^{-1}\left(\frac{2v_{+}}{J_0}\right)=J_0\sqrt{1-{\left(\frac{2v_{+}}{J_0}\right)}^2}$$

$$4v_0=J_1,J_1>4.$$

$$\frac{-1}{D}{\cos }^{-1}\left(\frac{4v_0}{J_1}\right)=\frac{J_1}{2}\sqrt{1-{\left(\frac{4v_0}{J_1}\right)}^2}.$$

$${\partial }_{t}v(x,t)=-v(x,t)+E-f(x,t)v(x,t)+{\int }_{-\infty }^{t-D}{\int }_{-\pi }^{\pi }J(x-y,t-s-D)f(y,s)dyds$$

$${\partial }_{t}w=\left(-1-f(v_0)-v_0{f}^{\prime }(v_0)\right)w+f^{\prime }(v_0){\int }_{-\infty }^{t-D}{\int }_{-\pi }^{\pi }J(x-y,t-s-D)w(y,s)dyds$$

$${\partial }_t{\hat{w}}_k=(-1-f(v_0)-v_0{f}^{\prime }(v_0)){\hat{w}}_k+f^{\prime }(v_0){\int }_{-\infty }^{t-D}{\hat{J}}_k(t-s-D){\hat{w}}_k(s)ds\text{ , }k\in \mathbb{Z}$$

$${\lambda }_k=-1-f(v_0)-v_0{f}^{\prime }(v_0)+f^{\prime }(v_0)e^{-{\lambda }_{k}D}{\int }_0^{\infty }{\hat{J}}_k(\tilde{s})e^{-{\lambda }_k\tilde{s}}d\tilde{s}$$

$${\lambda }_k=-2v_0+{\hat{J}}_k{e}^{-{\lambda }_{k}D}$$

$$0$$

$${\lambda }_{i,k}$$

$$J_0=2+2\sqrt{1-E}$$

$$\frac{-1}{D}{\cos }^{-1}\left(\frac{2v_0}{J_0}\right)=J_0\sqrt{1-{\left(\frac{2v_0}{J_0}\right)}^2}.$$

$$1\pm \frac{|J_0|}{J_0}\sqrt{1+\frac{4(E-J_0)}{J_0^2}}=\cos (-D{\lambda }_i)$$

$$\left|\frac{2v_0}{J_0}\right|=\left|\frac{J_0+\sqrt{J_0^2+4(E-J_0)}}{J_0}\right|<1$$

$$\pm \sqrt{J_0^2+4(E-J_0)}=J_1/2-J_0\text{ if }{J}_1>4$$

$$E=J_1^2/4+J_0(4-J_1)$$

$$\frac{-1}{D}{\cos }^{-1}\left(\frac{4v_0}{J_1}\right)=\frac{J_1}{2}\sqrt{1-{\left(\frac{4v_0}{J_1}\right)}^2}$$

$$J_{\mathrm{time}}(t-D)=\frac{J_0}{\tau }{e}^{-(t-D)/\tau }H(t)$$

$$J_{\mathrm{time}}(t-D)=\frac{J_0}{{\tau }^2}(t-D)e^{-(t-D)/\tau }H(t-D).$$

$${\lambda }_k=-2v_0+{\hat{J}}_k{e}^{-{\lambda }_{k}D}\frac{1}{(1+{\lambda }_k\tau )}$$

$$-\tau {\lambda }_{i,k}^2=-2v_0+{\hat{J}}_k\cos ({\lambda }_{i,k}D)$$

$${\lambda }_{i,k}={\hat{J}}_k\sin (-{\lambda }_{i,k}D)-2v_0{\lambda }_{i,k}\tau$$

Full text

{MSCcodes}

92B20, 92C20, 35Q70, 35Q92, 35B36, 60G55

1 Introduction

Spatially and temporally patterned behavior in large-scale neuronal activity are thought to play a role in a wide range of neurobiological phenomena. Some examples include geometric visual hallucinations [ermentroutMathematicalTheoryVisual1979], orientation tuning in the visual cortex [ben-yishaiTheoryOrientationTuning1995], and short term working memory [camperiModelVisuospatialWorking1998, LaingEtAllMultipleBumpsNeuronal2002]. These spatiotemporal dynamics are often understood using neural field equations, which are macroscopic descriptions of large networks’ activity [amariDynamicsPatternFormation1977, bressloffSpatiotemporalDynamicsContinuum2012, grossbergContourEnhancementShort1982, grossbergLearningEnergyentropyDependence1969, wilsonExcitatoryInhibitoryInteractions1972, wilson_mathematical_1973]. Although they have proven extremely useful in modeling neuronal dynamics, these macroscopic equations lack biophysical detail and can fail to capture spatiotemporal behaviors observed in microscopic spiking network models [byrneNextgenerationNeuralField2019]. The lack of a direct connection between scales obscures the relation between single-cell and population dynamics.

Integrate-and-fire (IF) models can capture a range of single-cell neural dynamics [jolivetQuantitativeSingleneuronModeling2008, naudFiringPatternsAdaptive2008, teeterGeneralizedLeakyIntegrateandfire2018]. In these models, the nonlinear dynamics of spiking are replaced with a simple fire-and-reset rule where the membrane voltage returns to a fixed value after the neuron emits a spike. IF models can be lower dimensional approximations of more bio-physical conductance based models such as the Hodgkin-Huxley model [abbottModelNeuronsHodgkinHuxley1990]. Networks of IF neurons can exhibit spatial, temporal, and spatiotemporal patterns. Oscillations can be induced by either synaptic delays or excitatory-inhibitory interactions [brunelFastGlobalOscillations1999, brunelDynamicsSparselyConnected2000, vanvreeswijkWhenInhibitionNot1994, whittingtonInhibitionbasedRhythmsExperimental2000]. The spatial profile of synaptic projections controls transitions between spatially incoherent and coherent activity [laingStationaryBumpsNetworks2001a, rosenbaumBalancedNetworksSpiking2014]. Despite the simplicity of the single-neuron dynamics, large networks of IF neurons are high dimensional and complex. Therefore, to analytically study temporal and spatial transitions in the network-level dynamics, there has been much interest in field theories that are both analytically tractable and directly connected to the microscopic single-neuron properties, e.g., [byrneNextgenerationNeuralField2019, montbrioMacroscopicDescriptionNetworks2015, schwalgerMindLastSpike2019].

Here, we study soft-threshold leaky IF (sLIF) networks using a mean-field theory [ockerDynamicsStochasticIntegratefire2023]. In particular, we explore oscillatory, spatial, and spatiotemporal instabilities in networks of sLIF neurons with temporal or spatial structure in their synaptic interactions. First, we introduce the sLIF network and mean-field approximation, and calculate a dispersion relation to identify various instabilities (Section 3.1). These instabilities, along with bifurcations of the resulting patterns, delineate regions in the parameter space phase diagram where there exist one or more qualitatively distinct steady state solutions of the mean-field theory. Each of these steady states correspond to a different type of behavior exhibited by the network.

Taken together, our findings reveal how different structures of connectivity and temporal dynamics give rise to a rich repertoire of network behaviors in a spatially extended sLIF network. In particular, we show that transmission delays and spatially structured excitation-inhibition in the synaptic connections support multiple kinds of activity states across different network configurations, as well as multi-stability within individual ones, and we further examine analyze the effect of more complex temporal kernels on the emergence of instabilities (Appendix 5.4). We identify and analyze four primary types of behavior and ten distinct regions of parameter space:

(Sections 3.1–3.2) Homogeneous steady states, constant in both time and space,

corresponding to globally active or quiescent network states,

• (i)

Quiescent region,

• (ii)

High activity region,

• (iii)

Quiescent and high activity region,

(Section 3.3) Homogeneous oscillatory states, corresponding to global oscillations in network activity,

• (iv)

Oscillatory region,

• (v)

Oscillatory and high activity region,

(Section 3.4) Stationary spatial patterns, in which only a portion of the network is above threshold, corresponding to spatially patterned activity,

• (vi)

Spatially patterned activity region,

• (vii)

Spatially patterned and high activity region,

• (viii)

Quiescent and spatially patterned activity region,

• (ix)

Spatially patterned, quiescent, and high activity region,

(Section 3.5) Spatiotemporal patterns, including standing waves, traveling waves, and oscillating Turing patterns,

• (x)

Standing and traveling wave region.

2 The model

We use a soft-threshold leaky integrate-and-fire neuron (sLIF) as our microscopic model to describe how individual neurons integrate incoming signals and generate action potentials (or spikes). We consider a network of $N$ sLIF neurons with sparse random connectivity on the ring (Fig. 1A). One-dimensional networks defined on periodic domains have been used to model feature selectivity in neuronal populations—for example, the sensitivity of neurons in the visual cortex to the orientation of visual stimuli [ben-yishaiTheoryOrientationTuning1995, hansel199813]. Neuron $i$ at location $x_{i}\in[-\pi,\pi)$ has membrane voltage $v_{i}(t)$ at time $t$. Let $n_{j}(t)$ be the counting process associated with neuron $j$, recording the total number of spikes that neuron $j$ has emitted up until time $t$. The increments of that process are $dn_{j}(t)$. After neuron $j$ emits a spike, its membrane potential is reset to the fixed value $r$. A postsynaptic neuron $i$ receives and integrates spikes from presynaptic neuron $j$ through the synaptic filter $J_{ij}(s-D)$, where $s$ is the time since the presynaptic spike and $D$ the synaptic delay (Fig. 1B).

We implement stochastic spike emission to model the variability in the membrane potentials at which neurons emit a spike [gerstnerPopulationDynamicsSpiking2000]. This soft threshold for spike emission distinguishes the sLIF model from classic LIF models, where neurons deterministically emit a spike at when the voltage hits a threshold value from below. The increments $dn_{j}(t)$ define a point process, $\dot{n}_{j}(t)$, also called the spike train. Its intensity is given by a non-negative function, $f(v_{j}(t))$. The intensity function is often chosen to be zero below a threshold, so there is no chance of emitting a spike at low membrane voltages. We take $f(v_{j})$ to increase unboundedly with $v_{j}$, so that if the voltage blows up the neuron will be guaranteed to emit a spike and reset.

Finally, each neuron has resting potential $E_{i}$, also incorporating any external input currents. Together, the membrane voltage of neuron $i$, $v_{i}(t)$, evolves according to the Itô stochastic differential equation

[TABLE]

The neuron integrates input spikes from the rest of the network through the synaptic filter, or coupling function, $J_{ij}(t)$. We show that with simple pulse coupling and no synaptic delay, $J_{ij}(t)=J_{0}\delta(t)$, the sLIF network can exhibit homogeneous spiking behavior (Fig. 1C). When $J$ has more temporal and spatial structure, the network exhibits spatiotemporal dynamics. For non-zero delay values, the sLIF network can exhibit homogeneous oscillatory activity (Fig. 1D). Additionally, with sufficient spatial modulation in $J$, the sLIF network exhibits spatially patterned activity (Fig. 1E). Finally, if our coupling function has both delays and spatial modulation, the network can exhibit spatiotemporal dynamics with patterns of traveling or standing waves (Fig. 1F, which also shows a spontaneous, noise-induced transition between those two patterns). To understand these dynamics, we develop a neural field approximation of (1).

2.1 Mean-field theory

We will study the following mean-field theory for (1),

[TABLE]

Here, $\bar{v}(x,t)$ is the mean membrane potential on the ring $x\in[-\pi,\pi)$, $\bar{f}:=f(\bar{v})$ is the mean-field approximation of the spike trains $\dot{n}$, $J\in L^{1}([\pi,\pi)\times\mathbb{R}_{+}$) is the synaptic filter, and $E$ is the mean resting potential. We have non-dimensionalized the model by measuring time relative to $\tau$ and the membrane potential relative to $r$, setting $\tau=1$ and $r=0$. The presence of the reset term $-\bar{f}(x,t)\bar{v}(x,t)$ differentiates this model from the classic Amari-Grossberg activity equations [amariDynamicsPatternFormation1977, grossbergLearningEnergyentropyDependence1969].

Rigorous mean-field limits for the population density of soft-threshold integrate-and-fire networks with delays, spatial structure, and disordered connectivity have each been studied separately, e.g., [chevallier_mean-field_2017, cormierHopfBifurcationMeanField2021, duarte_hydrodynamic_2015, jabinDenseNetworksIntegrateandfire2024, jabin_mean-field_2025, robertDynamicsRandomNeuronal2016]. To our knowledge, a mean-field limit for soft-threshold integrate-and-fire networks with all those factors combined has not yet been rigorously established. Here, we instead study a simple approximate mean-field theory for the expected voltage (2), which can be derived following the method of [ockerDynamicsStochasticIntegratefire2023]. Even in the continuum limit, $N\to\infty$, (2) will not be an exact description of the microscopic model; it results from truncating a cumulant hierarchy. We leverage its simplicity, however, for a qualitative understanding of the instabilities of stationary homogeneous states and the resulting activity patterns.

3 Results

We begin our study of the dynamics of (2) by characterizing its homogenous equilibria and their stability. Throughout, we will take a rectified linear intensity function, $f(v)=\lfloor v-1\rfloor_{+}$.

3.1 Homogeneous equilibria and dispersion relation

There are three possible homogeneous equilibria of (2): a quiescent state below threshold ($v_{Q}$) and two solutions above threshold ($v_{\pm}$). These are given by

[TABLE]

where $J_{0}=\int\int J(x-y,t-s)\,ds\,dy$.

The quiescent sub-threshold state $v_{Q}$ is always stable. Below threshold, the convolution and reset terms of (2) vanish, and the voltage converges to the resting potential $E$. The supra-threshold state $v_{+}(v_{-})$ is stable(unstable) in absence of delay and spatial modulation in the coupling. Therefore, without spatial or temporal structure in the coupling, the mean activity will converge to $v_{Q}$ or $v_{+}$ [ockerDynamicsStochasticIntegratefire2023]. The first three regions we identify in the phase diagram (Fig. 2A), where the stable homogeneous steady states exist, are:

1. i)

Q: Only a quiescent state ($v_{Q}$, below threshold) exists and is stable. 2. ii)

H: Only a high activity state ($v_{+}$, above threshold) exists and is stable. 3. iii)

Q-H: Both the quiescent and high activity states are stable.

In the bistable region Q-H, delimited by saddle-node bifurcations (Section 3.2), the unstable equilibrium $v_{-}$ is the separatrix between $v_{Q}$ and $v_{+}$.

To identify the instabilities of a homogeneous state $v_{0}$, we perform a linear stability analysis of the mean-field approximation (2), which yields a dispersion relation characterizing the stability of $v_{0}$ with respect to perturbations in each spatial Fourier mode. Given the threshold-linear intensity function, $v_{0}$ above threshold ($v_{0}>1$), and recalling the delay $D$ in the synaptic coupling, the linearization of (2) about $v_{0}$ is

[TABLE]

By taking the spatial Fourier transform of (4) $\left(w(x,t)=\sum_{k\in\mathbb{Z}}\hat{w}_{k}(t)e^{ikx}\right)$ and making the ansatz $\hat{w}_{k}(t)=e^{\lambda_{k}t}$, we obtain the dispersion relation

[TABLE]

where $\lambda_{k}$ is the eigenvalue associated with the $k^{th}$ mode, $\hat{w}_{k}$. The stability depends on the structure of interactions between the neurons through the time-dependent Fourier coefficients of the coupling function $\hat{J}_{k}(t)=\int e^{-ikx}J(x,t)\,dx$. The dispersion relation for a general intensity function can be found in Appendix 5.2.

For simplicity, we assume delayed pulse and cosine coupling for the temporal and spatial structure of the coupling function:

[TABLE]

With delayed pulse coupling, the dispersion relation (5) simplifies to

[TABLE]

where $\hat{J}_{k}$ is the $k^{th}$ Fourier coefficient of the spatial coupling. Given cosine coupling (8), there are only three nonzero Fourier coefficients, $\hat{J}_{0}=J_{0}$ and $\hat{J}_{\pm 1}=J_{1}/2$. Therefore, all spatial modes other than $k=0$ and $k=1$ are always stable.

The analysis above extends to other coupling structures. Other $L^{2}$ spatial kernel functions can be studied in a similar way by substituting their Fourier coefficients into (5). For other choices of temporal kernels, the temporal integral in (5) will differ. We also analyze two other temporal kernels which model the rise and fall times of the action potential in addition to the delay (Appendix 5.4).

The dispersion relation (9) only pertains to equilibria above threshold. As discussed above, the only sub-threshold homogeneous equilibria $v_{Q}$ is always stable. Consequently, we restrict our analysis of instabilities to the two steady states above threshold ($v_{\pm}$). The location of primary instabilities are given by the dispersion relation (9) when Re$(\lambda_{k})=0$ for different modes $k$ and temporal frequencies Im$(\lambda_{k})$.

These instabilities correspond to four smooth bifurcations of $v_{0}$: a saddle-node bifurcation (Section 3.2), a Hopf bifurcation (Section 3.3), a Turing bifurcation (Section 3.4), and a Turing-Hopf bifurcation (Section 3.5). At threshold ($v_{0}=1$), non-smooth bifurcations can occur due to the non-smooth point of $f$, not described by the linearized analysis; we therefore study these numerically.

3.2 Homogeneous, static instabilities

The first instability of the homogeneous state occurs at $k=0$ with Im($\lambda_{0})=0$. Setting Re$(\lambda_{0})=0$ in the dispersion relation (9) gives the saddle-node bifurcation where $v_{\pm}$ collide and vanish,

[TABLE]

where $J_{0}>2$ and $E<1$. Beyond the bifurcation, at lower $J_{0}$ values, solutions tend toward the quiescent state $v_{Q}$. The saddle-node curve forms the boundary between the bistable regime (Q-H) and the quiescent regime (Q) (Fig. 2A). This bifurcation is independent of the delay and spatial coupling.

There is also a non-smooth saddle-node bifurcation at $E=1$ where $v_{-}$ and $v_{Q}$ collide and vanish. This bifurcation cannot be identified via the dispersion relation, which is derived from the linearized dynamics above threshold. The non-smooth saddle-node curve is the boundary between the bistable region (Q-H) and high activity region (H) (Fig. 2A). These bifurcations were identified in [ockerDynamicsStochasticIntegratefire2023]. If the intensity function $f(v)$ is not threshold-linear, another saddle-node bifurcation between high and low activity states can emerge [paliwalMetastabilityNetworksStochastic2025].

3.3 Delay-induced instability

Transmission delays can induce a Hopf bifurcation, at which homogeneous (bulk) oscillations of the network emerge. This oscillatory instability occurs at $k=0$ with Im($\lambda_{0})\neq 0$ when a pair of complex conjugate eigenvalues cross the imaginary axis ($\mathrm{Re}(\lambda_{0})=0$). From the dispersion relation (9), we find that $v_{+}$ undergoes a Hopf bifurcation when

[TABLE]

and $v_{+}>1$, $3J_{0}^{2}+4J_{0}>4E$, $D>0$, and $J_{0}<0$ (Appendix 5.3). Beyond the Hopf bifurcation, at lower $J_{0}$ values, there are homogeneous oscillatory solutions $v_{osc}$ (Fig. 2B, upper) if the delay is sufficiently large (Fig. 2C). Thus uniform oscillations require inhibitory coupling with sufficient delay and the reversal potential $E$ above threshold. However, an excessively large $E$ offsets the inhibition and oscillations do not occur.

We find, through numerical continuation of $v_{osc}$ (Appendix 5.5), that the Hopf bifurcation (11) is subcritical (Fig. 2D-F). Shortly after the minimum voltage on branch of unstable oscillations emerging from the bifurcation falls below threshold, there is a fold of large amplitude limit cycles, where the unstable branch meets a stable branch. The Hopf bifurcation and the fold of limit cycles bound a narrow bistable region of both $v_{+}$ and $v_{osc}$ (Fig. 2A, “O-H”: dotted line). As the magnitude of $J_{0}$ or $E$ decrease or $D$ increase from the their values at the fold, so does the amplitude and the period of the oscillation (Fig. 2D-F). As, $E$ further decreases, $v_{osc}$ vanishes at a supercritical non-smooth Hopf bifurcation at $E=1$, where $v_{+}$ drops below threshold (Fig. 2E).

To summarize, we have found with the presence of delays, there are two additional regions of the phase diagram (Fig. 2A) involving $v_{osc}$:

1. iv)

O: Only a homogeneous, large-amplitude oscillation ($v_{osc}$) around $v_{+}$ is stable. 2. v)

O-H: Both a homogeneous oscillation ($v_{osc}$) and high activity ($v_{+}$) state are stable.

The Hopf curve (11) is the boundary between the oscillatory (O) and bistable (O-H) regions, while the nearby fold of limit cycles is the boundary between the bistable (O-H) and high (H) activity regions. The non-smooth supercritical Hopf bifurcation at $E=1$ separates the oscillatory (O) and quiescent (Q) regions.

Oscillatory behavior in the sLIF network (Fig. 2B lower panel, orange) qualitatively matches the mean-field approximation (Fig. 2B lower panel, black), with some discrepancy expected due to the mean-field approximation. The maximum amplitude (Fig. 2D-F) of the sLIF network oscillations (grey circles) are generally lower than those of the mean-field approximation (solid black line). Moreover, the oscillatory region of the sLIF network is narrower: the onset of oscillation at the Hopf bifurcation occurs at a lower value of $E$ than predicted by the mean-field theory, while the non-smooth Hopf bifurcation remains at threshold (Fig. 2E). The sLIF simulation data are shown only for large amplitude oscillations, which we further investigate in Fig. 3.

We observe the small bistable region O-H in the sLIF network with a slow ramp of the bifurcation parameter through the region, starting on either side. We show this process first with the mean-field approximation, where we know the location of the Hopf bifurcation analytically and the location of the fold of limit cycles from numerical continuation. We then repeat the process using the sLIF network and observe similar phenomena.

We initialize a simulation of the mean-field (Fig. 3A, grey) in the oscillatory region (O) at a stable limit cycle and slowly increase the parameter $J_{0}$, through the bistable region (Fig. 3A, green). As the ramp continues, the limit cycle (Fig. 3B, Fig. 3A blue) persists until it loses stability at the fold marking the end of bistable region. The network activity then transitions to the stable homogeneous state and the oscillation decays.

For the reverse ramp (Fig. 3C), we choose initial conditions above the bistable region near the stable homogeneous state. As $J_{0}$ is slowly decreased through the bistable region, the network stays close to the stable homogeneous solution (Fig. 3D), until it loses stability at the Hopf bifurcation. A temporal oscillation then begins to grow and converges to the stable limit cycle. The delayed loss of stability, where the system stays close to the unstable homogeneous state after the bifurcation point, is expected when simulating smooth dynamical systems near a Hopf bifurcation [izhikevichDynamicalSystemsNeuroscience2014]. It is due to the time it takes for the system to diverge from the unstable equilibrium.

We repeat both parameter ramps with the sLIF network. When increasing $J_{0}$ from the oscillatory region (Fig. 3E), we observe a stable oscillation persisting in the bistable region (Fig. 3F) and a similar loss of stability of the limit cycle, but at a slightly lower $J_{0}$ value than that of the mean-field approximation. Decreasing $J_{0}$ from the homogeneous regime, we again observe similar behavior to the mean-field (Fig. 3G) where the behavior stays close to the homogeneous solution in the bistable region (Fig. 3H), and then loses stability, but slightly lower $J_{0}$ values than the mean-field.

In the sLIF network, the loss of stability of the homogeneous solution occurs at a lower parameter value than where the limit cycle loses stability. There is thus a parameter range where the stationary and oscillatory solutions coexist. We also demonstrate this bistability at set parameter values where a stimulus switches the network between the homogeneous and oscillatory states in both the mean-field approximation (Fig. 3I) and the sLIF network (Fig. 3J).

3.4 Spatial instabilities

Next, we investigate the instabilities induced by spatial modulation of the coupling function. A spatial instability occurs at $k\neq 0$ with Im$(\lambda_{k})=0$ when Re$(\lambda_{k})=0$. Recall that the amount of the spatial modulation of the synaptic coupling is given by the parameter $J_{1}$. Setting $k=1$ and $\lambda_{1}=0$ in the dispersion relation (9) with cosine coupling (8) reveals a spatial instability at

[TABLE]

The spatial instability (12) takes one of two forms, depending on $J_{0}$. For $J_{0}<J_{1}/2$, it is a Turing bifurcation of $v_{+}$ since the $k=0$ mode remains stable, while for $J_{0}>J_{1}/2$, it is a secondary instability of $v_{-}$, whose $k=0$ mode is already unstable (Fig. 4A, solid and dashed red curves).

The spatial instability switches from $v_{+}$ to $v_{-}$ at $J_{0}=J_{1}/2$ and $E=J_{0}-J_{0}^{2}/4$, where the Turing and saddle-node curves meet tangentially at a codimension-2 bifurcation (Fig.  4A, red and black curves). At this point, the $k=0$ and $k=1$ Fourier coefficients of the coupling function are equal and the saddle-node (instability of the $k=0$ mode) and Turing (instability of the $k=1$ mode) bifurcations coincide.

The existence and location of the spatial instability curve depends on the amount of spatial modulation in the coupling function (Fig. 4B). If $J_{1}\leq 4$, no spatial instability occurs. This critical value emerges because at $J_{1}=4$, the spatial instability curve lies at the boundary of the regions of where $v_{\pm}$ are above threshold ($E=1$). If $J_{1}<4$, then (12) is only defined in regions where $v_{\pm}$ are below threshold and no longer correspond to equilibria. Therefore $J_{1}$ must exceed this critical value ($J_{1}>4$) for pattern formation to occur. Additionally, the Turing bifurcation of $v_{+}$ only occurs when the spatial coupling is either purely inhibitory ($J_{0}<-|J_{1}|$) or has local excitation and long range inhibition ($|J_{0}|<J_{1}$), as illustrated later in Figs. 5A and B (Section 3.5).

3.4.1 A Turing pattern

Beyond the spatial instability (12), the system supports a stationary periodic Turing pattern (Fig. 4C). The profile of the Turing pattern has portions above and below threshold, corresponding to the active and inactive portions of the network. We refer to this as the spatially patterned state (S). The sLIF network exhibits qualitatively similar spatial patterns (Fig. 4D,E).

The Turing bifurcation is subcritical with a nearby fold of Turing patterns. We used numerical continuation of the Turing pattern and determined its stability by numerically computing eigenvalues of the mean-field linearization about the Turing pattern (Appendix 5.5). An unstable branch of small-amplitude Turing patterns (Fig. 4F-H, dashed black) emerges backward when the active state $v_{+}$ (Fig. 4F-H, purple) destabilizes through the Turing bifurcation as $E$ or $J_{0}$ decrease or $J_{1}$ increases. The numerical continuation revealed that all completely supra-threshold Turing patterns are unstable. The minimum value of the unstable Turing pattern drops below threshold at $v=1$, it soon after meets a stable branch of Turing patterns at a fold. The stable branch extends off in both $J_{0}$ (Fig. 4G) and $J_{1}$ (Fig. 4H), but ends at a non-smooth Turing at $E=1$ when $v_{+}$ drops below threshold (Fig. 4F).

3.4.2 Multi-stability with the Turing pattern

Due to the subcritical nature of the Turing bifurcation, there exists a bistable region of both the homogeneous state and Turing pattern between the bifurcation point and the fold. Within this region, networks can support either spatially patterned or high activity (Fig. 4A, “S-H”). Additionally, we find that in a small neighborhood of the codimension-2 bifurcation of the Turing and saddle-node, the non-smooth Turing is also subcritical with a nearby fold extending into the quiescent region. This gives rise to two additional regions of multi-stability with the Turing pattern involving the quiescent state: a quiescent and spatially patterned region (“Q-S”) and a region of triple stability between the quiescent , spatially patterned, and high activity states (“S-Q-H”) when the subcritical regions of the two bifurcations overlap. A detailed characterization of the regions and their boundaries, through numerical continuation and tracking the folds, is provided in Appendix 5.6.

We have identified four additional regions of the phase diagram (Fig. 4A) involving the Turing pattern:

1. vi)

S: A stable spatially patterned solution exists. 2. vii)

S-H: Spatially patterned and High ($v_{+}$) activity solutions exist and are stable. 3. viii)

Q-S: Stable Quiescent ($v_{Q}$) and spatially patterned solutions exist. 4. ix)

S-Q-H: Spatially patterned, Quiescent ($v_{Q}$), and High ($v_{+}$) activity solutions exist and are stable.

The spiking network also exhibits multi-stability in these regions. We observe bistability between the spatially patterned and high activity states, as predicted by the subcritical Turing bifurcation in the mean-field approximation (Fig. 4I). Additionally, we observe the bistability between the quiescent and spatially patterned activity states (Fig. 4J), as predicted by the subcritical non-smooth Turing, and multi-stability between all three states (Fig. 4K).

3.5 Spatiotemporal instabilities

In the presence of both a delay and spatial modulation in the coupling, a spatiotemporal instability occurs and gives rise to spatiotemporal patterns. From the dispersion relation (9) we find the spatiotemporal instability in the first spatial mode ($k=1$) at nonzero frequency (Im($\lambda_{1})\neq 0$). It occurs when $v_{0}$ is above threshold and

[TABLE]

This instability is a Turing-Hopf bifurcation of $v_{+}$ if $J_{0}<J_{*}$ and a secondary instability of $v_{-}$ if $J_{0}>J_{*}$, where $J_{*}$ is the $J_{0}$ value of the codimension 2 point at which the Turing-Hopf curve and saddle-node curve intersect.

Four primary bifurcations occur throughout the phase space when $D\neq 0$ and $J_{1}\neq 0$: the saddle-node, Hopf, Turing, and Turing-Hopf bifurcations (Fig. 5A,B: blue, green, red, purple respectively). The locations of these instabilities for any resting potential above threshold ($E>1$), are topologically similar to Fig. 5A with different scaling dependent on $E$ and the delay. Similarly, the locations of the instabilities when $E<1$ are similar to Fig. 5B.

Formulaically, the Turing-Hopf curve (13) is similar to the Hopf curve (11). This is because they are both instabilities with non-zero frequencies, but at different spatial modes. Despite this, the Turing-Hopf curve increasingly resembles a reflection of the Turing curve across $J_{1}=0$ as the delay increases. This becomes exact in the limit $D\to\infty$.

Beyond the Turing-Hopf bifurcation, we observe the emergence of standing and traveling waves (Fig. 5E, F) induced by the $O(2)$ symmetry of $J$ to translations and reflections in $x$. The sLIF exhibits both these patterns (Fig. 1G). In that simulation, a standing wave spontaneously transitions into a traveling wave, suggesting metastability in the sLIF network. Bistability of the standing and traveling waves with other states of the network remains to be investigated, and could be done using numerical continuation of solutions implementing both delays and spatial modulation. This is the final region of the phase diagram we describe:

1. x)

SW-TW: Standing and traveling waves exist and are stable.

In regions beyond multiple primary instabilities, both the mean-field approximation and sLIF network exhibit a variety of additional spatiotemporal patterns, some likely arising from secondary or higher-order bifurcations. Beyond both the Turing and Hopf bifurcations, we observe a Turing-Hopf pattern (Fig. 5H, I). Beyond the Turing-Hopf bifurcation, we observe additional dynamics, likely due to an additional instability (Fig. 5P-Q). Finally, beyond both the Turing-Hopf and Hopf bifurcations, we observe a variety of mixed mode oscillations (Fig. 5J-S). Investigating these secondary instabilities and the rich variety of dynamical patterns to which they lead remains a direction for future research.

Remark 3.1**.**

*The classic neural field equation $\partial_{t}x=-x+f(E+J\ast x)$, on the ring, has similar transitions from a homogeneous equilibrium to oscillatory and Turing patterns in the presence of delays and spatially modulated coupling [roxinRoleDelaysShaping2005]; see also [amariDynamicsPatternFormation1977, ben-yishaiTheoryOrientationTuning1995, wilsonExcitatoryInhibitoryInteractions1972, wilson_mathematical_1973]. The locations of those bifurcations differ in the classic rate and sLIF networks, however. For example, in the classic neural field equation the Turing instability of the homogenous active state occurs at $J_{1}=2$; in the sLIF mean-field dynamics we studied, the Turing instability occurs at $2J_{1}=J_{0}+\sqrt{J_{0}^{2}+4(E-J_{0})}$ (12). Secondary bifurcations and codimension-2 points also differ between the two models. *

4 Discussion

Patterns in neural activity are classically modeled using neural field equations such as the Wilson-Cowan and Amari-Grossburg models [amariDynamicsPatternFormation1977, grossbergLearningEnergyentropyDependence1969, wilson_mathematical_1973]. These can be derived as explicit mean-field theories for highly simplified microscopic models or by making strong assumptions like a separation of timescales between neural and synaptic dynamics [buiceFieldtheoreticApproachFluctuation2007, ginzburgTheoryCorrelationsStochastic1994, ohiraMasterequationApproachStochastic1993, pintoQuantitativePopulationModel1996]. This complicates their relation to biophysical microscopic models. This discrepancy is one motivation for the development of next-generation neural field theories from specific microscopic models [byrneNextgenerationNeuralField2019, montbrioMacroscopicDescriptionNetworks2015, schwalgerMindLastSpike2019].

Soft-threshold integrate-and-fire networks replace the nonlinear dynamics of spike emission with a probabilistic spike-and-reset rule, but other biophysical detail can be directly incorporated. This family of models is often studied with a population density approach, which exposes rigorous mean-field limits [demasiHydrodynamicLimitInteracting2015, gerstnerPopulationDynamicsSpiking2000, jabin_mean-field_2025]. The mean-field limit can also be exposed through the joint density functional of the network’s state trajectories, rather than the population density [ockerDynamicsStochasticIntegratefire2023, robertDynamicsRandomNeuronal2016]. The population density approach has recently been extended to networks with either spatial connectivity or delays, exposing Hopf or Turing bifurcations separately [cormierHopfBifurcationMeanField2021, dumontOscillationsFullyConnected2023, dumontPatternFormationSpiking2024, jabinDenseNetworksIntegrateandfire2024].

Here, we instead used a deterministic approximation to study the dynamics of soft-threshold leaky integrate-and-fire (sLIF) networks with both synaptic delays and spatial connectivity. This approximation has a similar form to the classic Amari-Grossberg equations with an additional term due to the voltage resets after action potentials. We found oscillatory, spatial, and spatiotemporal instabilities that generate coherent activity patterns such as bulk oscillations, stationary Turing patterns, standing and traveling waves. We identified various multi-stable regions of different network states due to the sub-criticality of the Hopf and Turing bifurcations. In these regions, the network can support both patterned behavior and a homogeneous state (quiescent or active), and can be switched from one state to the other by global or local perturbations. We confirmed all these predictions of the mean-field dynamics in simulations of the underlying microscopic stochastic system.

From a mathematical point of view, all of the primary bifurcations we described could be studied using a parameter dependent infinite-dimensional center manifold reduction [avitabile2020local, haragus2010local, vanderbauwhede1992center, veltz13]. Due to the lack of other positive spectrum of the homogeneous equilibrium at the bifurcation, invariant foliation results would give that stability on the manifold implies local stability in the full system.

The sLIF networks we studied are members of a broader family of soft-threshold integrate-and-fire networks with richer voltage dynamics. Soft-threshold integrate-and-fire networks with nonlinear voltage dynamics can have oscillatory solutions [cormierHopfBifurcationMeanField2021]. Adaptive exponential integrate-and-fire neurons can capture a broad range of single-neuron spike patterns [naudFiringPatternsAdaptive2008, touboulDynamicsBifurcationsAdaptive2008]. Adaptive nonlinear soft-threshold integrate-and-fire networks are also amenable to the same type of mean-field approximation we used here. Including excitatory and inhibitory cell types and their relative spatial coupling profiles and synaptic timescales may also give rise to new spatiotemporal dynamics [huangCircuitModelsLowDimensional2019].

Furthermore, it would be of interest to consider these patterns on unbounded domains as well as infinite dimensions. Near the Turing instabilities found here, we expect a family of (subcritical) periodic waves to bifurcate. The subcriticality and bistability in the bounded-domain problem indicate that it should be possible, in the unbounded domain, to construct fronts connecting the periodic state to the quiescent state, spatial patterns, and other complex modulated waves [faye2015modulated, faye2018center].

Next-generation neural field theories also aim to explicitly describe fluctuations, correlations, and synchrony [byrneNextgenerationNeuralField2019]. The neural field equations we study can be straightforwardly extended to include these through a fluctuation expansion of the density functional of the underlying microscopic model [bressloffStochasticNeuralField2010, buiceFieldtheoreticApproachFluctuation2007, ockerDynamicsStochasticIntegratefire2023]. For the sLIF networks studied here, there are three possible types of fluctuation correction: 1) corrections to the mean-field rate $f(\bar{v})$ due to fluctuations in the voltage, 2) a correction to the mean-field voltage dynamics due to the spike-voltage covariance, and 3) variability in the synaptic field due to correlated or finite-size fluctuations in the activity. These have each been studied in non-spatial networks [ockerDynamicsStochasticIntegratefire2023, paliwalMetastabilityNetworksStochastic2025, schmutzFiniteSizeNeuronalPopulation2023]. How these different fluctuations interact with neural nonlinearities to shape spatiotemporal patterns in sLIF networks remains to be investigated.

Declarations

Authorship:

All authors have made substantial intellectual contributions to the study conception, execution, and design of the work. All authors have read and approved the final manuscript. In addition, the following contributions occurred: Conceptualization: Ryan Goh, Gabriel Koch Ocker; Formal analysis and investigation: Lauren Forbes, Jared Grossman, Ryan Goh; Writing - original draft preparation: Lauren Forbes; Writing - review and editing: Lauren Forbes, Ryan Goh, Gabriel Koch Ocker; Supervision: Montie Avery, Ryan Goh, Gabriel Koch Ocker.

Conflicts of interest:

The authors declare no conflicts of interest.

Data and code availability:

Code to reproduce the results of this manuscript can be found at https://github.com/lcforbes4/sLIFspatioTemporal

Funding:

The research of the authors was partially supported by NSF-DMS 2307650 (RG), and by a grant from the Allen Institute Mindscope Phase 4 program (GKO).

5 Appendix

5.1 Further reading

We used largely standard methods for linear stability analysis of partial differential equations and numerical continuation, as described in [coombesNeuralFieldsTheory2014, bressloffSpatiotemporalDynamicsContinuum2012, uecker2021numerical]. In the following sections, we describe these computations in detail.

5.2 General dispersion relation

To derive the dispersion relationship, we consider the mean-field approximation with a general coupling function $J(x,t-D)$ where $J(x,t)\in L^{1}((-\pi,\pi)\times(0,\infty))$ along with the intensity function $f(v(x,t))$,

[TABLE]

with periodic spatial domain, $x\in(-\pi,\pi]$.

Suppose the stationary uniform solution $v_{0}$, satisfying $v_{0}=E-v_{0}f(v_{0})+(J\ast f(v_{0}))(x,t-D)$, lies away from any non-differentiable points of $f$. Writing $v(x,t)=v_{0}+w(x,t)$ with $w(x,t)$ sufficiently small, linearizing (14) about $v_{0}$ yields

[TABLE]

Taking the Fourier transform of (15) in space results in the following system of differential equations indexed by the wave number $k$

[TABLE]

where $\hat{w}_{k}(t)=\int e^{-ikx}w(x,t)\,dx$ and $\hat{J}_{k}(t)=\int e^{-ikx}J(x,t)\,dx$ are the time-dependent $k^{th}$ Fourier coefficients of $w$ and $J$. Finally, making the ansatz $\hat{w}_{k}=e^{\lambda_{k}t}$ and the substitution $\tilde{s}=t-s-D$, we get the dispersion relation

[TABLE]

With the additional assumption of a threshold linear intensity function of the form $f(v)=\lfloor v-1\rfloor_{+}$, (17) simplifies to (5).

5.3 Instabilities with delayed pulse coupling

To derive the dispersion relation for delayed pulse coupling, we assume $\hat{J}_{k}(t)=\hat{J}_{k}\delta(t)$, $f(v)=[v-1]_{+}$, and $v_{0}>1$. The dispersion relation (17) then simplifies to

[TABLE]

Separating the real and imaginary components of (18) and assuming that $Re(\lambda_{k})=0$ yields

[TABLE]

where $\lambda_{i,k}:=Im(\lambda_{k})$. The eigenvalue equation can also be formulated with the Lambert W function to achieve the same results [corlessLambertWFunction1996]. From now on, we assume the Fourier coefficients of cosine spatial coupling (8) which are $\hat{J}_{0}=J_{0}$, $\hat{J}_{\pm 1}=J_{1}/2$, and $\hat{J}_{k}=0\,\forall k\notin\{0,\pm 1\}$. We next identify the locations of four instabilities of the stationary uniform equilibria $v_{\pm}$ given in (3).

I. Saddle-node Bifurcation ($k=0,\lambda_{i,0}=0$): This is an instability of the uniform spatial mode with zero temporal frequency. The imaginary part of the eigenvalue is zero, so (20) is trivial. By substituting the form of $v_{\pm}$ given in (3), (19) simplifies to $J_{0}=2\pm 2\sqrt{1-E}$. If $J_{0}=2-2\sqrt{1-E}$, the values of $v_{\pm}$ are both below the threshold in $f$, an so neither equilibrium exists. Therefore the instability only occurs along the curve

[TABLE]

II. Hopf Bifurcation ($k=0,\lambda_{i,0}\neq 0$): This is an instability of the uniform spatial mode with non-zero temporal frequency. We solve (19), yielding $\lambda_{i}=\frac{-1}{D}\cos^{-1}(\frac{2v_{0}}{J_{0}})$, using the principle branch of $\cos^{-1}$. Substituting this into (20) gives the Hopf curve as an implicit function of the parameters:

[TABLE]

We next show that only $v_{+}$ can undergo a Hopf Bifurcation and it requries $3J_{0}^{2}+4J_{0}>4E$ and $J_{0}<0$. Using (3), the real part of the dispersion relation (19) simplifies to

[TABLE]

Due to the range of cosine, (23) is only well defined for the minus case (instability of $v_{-}$) if $J_{0}>0$. (23) is well defined for the plus case (instability of $v_{+}$) if $J_{0}<0$. But since the principle branch of arccosine is a non-negative function, the Hopf curve given in (22) is only well defined for $J_{0}<0$. Therefore, we conclude that only $v_{+}$ can undergo a Hopf bifurcation.

Finally, to satisfy the domain of arccosine, we have the condition

[TABLE]

which can be rewritten as the set of inequalities $3J_{0}^{2}+4J_{0}>4E$, $v_{+}>1$, and $J_{0}<0$.

Using implicit differentiation, one can also show the transverse crossing condition of a Hopf bifurcation is satisfied along the curve (22). As the delay is increased, the Hopf curve increases in $J_{0}$ and approaches the curve $3J_{0}^{2}+4J_{0}=4E$, which is the boundary of the region on which the Hopf curve (22) is well-defined.

III. Turing Bifurcation ($k=1,\lambda_{i,1}=0$): This is an instability of a non-uniform spatial mode with zero temporal frequency. Similar to the saddle-node bifurcation case, since the imaginary part of the eigenvalue is zero, (20) is trivial and (19) simplifies to $2v_{0}=J_{1}/2$. Note that $v_{0}>1$ only if $J_{1}>4$. Therefore, there is a spatial instability of $v_{\pm}$ at

[TABLE]

This can be rewritten as

[TABLE]

which corresponds to the Turing Bifurcation of $v_{+}$ if $J_{0}<J_{1}/2$ and a secondary instability of $v_{-}$ if $J_{0}>J_{1}/2$.

IV. Turing-Hopf ($k=1,\lambda_{i,1}\neq 0$): This is an instability of a non-uniform spatial mode with non-zero temporal frequency. Similar to the Hopf bifurcation case, we substitute $\lambda_{i,1}=\frac{-1}{D}\cos^{-1}(\frac{2v_{0}}{J_{0}})$ from (19) into (20). This gives the curve of a spatiotemporal instability of $v_{0}$

[TABLE]

if $v_{0}>1$ as well, which occurs if $\frac{J_{1}}{4}\cos(D\lambda_{i,1})>1$ from (19).

Note that this curve intersects the saddle-node curve in a codimension-2 bifurcation if $J_{0}=\frac{J_{1}}{2}\cos(-D\lambda_{i,1})$. The instability (27) is a Turing-Hopf bifurcation of $v_{+}$ if $J_{0}<\frac{J_{1}}{2}\cos(-D\lambda_{i,1})$ and a secondary instability of $v_{-}$ if $J_{0}>\frac{J_{1}}{2}\cos(-D\lambda_{i,1})$.

5.4 More realistic synaptic responses

Synaptic potentials are characterized by three time scales: the latency or transmission delay (the time between the presynaptic spike emission and the postsynaptic response), as well as the rise time and the decay time [destexheKineticModelsSynaptic1998]. We consider two additional forms for the temporal coupling which model the rise and decay times of the synaptic potential in addition to the delay. The first is the delayed exponential

[TABLE]

where $H$ is the Heaviside step function. This kernel implements a decay time of the action potential in addition to the delay (Fig. 6A). The additional decay time smooths the postsynaptic response to incoming spikes compared to pulse coupling (Fig. 6B). The locations of the four primary instabilities are very similar to the pulsed coupling case (Fig. 6C and D). See Appendix 5.4.1 for the derivation of the dispersion relations and instability curves.

The delayed alpha function (Fig. 6E) is another option for the temporal profile which also includes the rise time of the synaptic potential. It is defined by

[TABLE]

This synaptic kernel further smooths the change in voltage of the postsynaptic neuron upon the arrival of a spike (Fig. 6F). Again, the instability diagrams are qualitatively similar to the two previous cases of temporal kernels (Fig. 6G and H). See Appendix 5.4.2 for the derivation of the dispersion relations and instability curve. Notably, a nonzero delay is not required for the existence of a Hopf bifurcation and homogeneous oscillatory behavior with this temporal profile. The rise time of the alpha function acts as an effective delay in the transmission of action potentials between neurons.

5.4.1 Instabilities with delayed exponential coupling

To identify the locations of instabilities with the delayed exponential temporal profile of the coupling, ie $\hat{J}_{k}(t)=\hat{J}_{k}\tau^{-1}e^{-(t)/\tau}H(t)$, we assume $f(v)=[v-1]_{+}$ and $v_{0}>1$. The dispersion relation (17) then simplifies to

[TABLE]

Separating the real and imaginary components of (30) with the assumption $Re(\lambda_{k})=0$ results in

[TABLE]

Note that when $\lambda_{i,k}=0$, (30) is equivalent to (18). Therefore the instabilities with zero temporal frequency occur along the same curves as the delayed pulse case. The saddle-node bifurcation in the delayed exponential case is then also given by (21) and the Turing bifurcation is given by (26).

To identify the locations of instabilities with non-zero temporal frequencies, we first rewrite (31b) as $\lambda_{i}=\hat{J}_{k}\sin(-\lambda_{i,k}D)/(1+2v_{0}\tau)$ and substitute it into (31a) to get

[TABLE]

Then, by using the identity $\sin^{2}(x)=1-\cos^{2}(x)$ and defining $y:=\cos(\lambda_{i,k}D)$, we can transform the equation above into the quadratic

[TABLE]

We solve this quadratic for $y$ and thus for $\lambda_{i,k}=D^{-1}\arccos(y),$

[TABLE]

Plugging this $\lambda_{i}$, which is in terms of the parameters $\hat{J}_{k},E,\tau,D$, into the imaginary part of our dispersion relation (31b) gives

[TABLE]

where $A=1+4v_{0}\tau+4v_{0}^{2}\tau^{2}$. This is the Hopf curve when $k=0$ and $\hat{J}_{k}=J_{0}$ or the Turing-Hopf curve when $k=1$ and $\hat{J}_{k}=J_{1}/2$.

Next, we will look at which parameter values and which equilibria $v_{0}$ the curve (33) is well defined for. First, we note that the range of the principle branch of arccosine is non-negative, so (33) is not defined for $\hat{J}_{k}>0$.

To satisfy the domain of arccosine in (32) to be defined, we need

[TABLE]

When $k=0$, the minus case does not satisfy this inequality and the plus case can be rewritten as the set of conditions. $J_{0}<-2$, $v_{+}>1$, and $4E<4J_{0}+3J_{0}^{2}$. Therefore the Hopf and Turing-Hopf bifurcations are both instabilities of $v_{+}$.

5.4.2 Instabilities with alpha function coupling

To identify the dispersion relation and primary instabilities of the homogeneous solution with the delayed alpha function chosen for the temporal coupling profile, we assume $v_{0}>1$, $f(v)=[v-1]_{+}$ and $\hat{J}_{k}(t)=\frac{\hat{J}_{k}}{\tau^{2}}te^{-t/\tau}H(t)$. The dispersion relation (17) becomes

[TABLE]

Separating the real and imaginary components of (34) with the assumption $Re(\lambda_{k})=0$,

[TABLE]

When $\lambda_{i,k}=0$, as in the precious cases, the saddle-node and Turing bifurcations are given by (21) and (26) respectively.

Next we identify the other two types of instabilities which have non-zero temporal frequency. From the real part of the dispersion relation (35),

[TABLE]

and substitute it into (36),

[TABLE]

Then, using the identity $\sin^{2}x=1-\cos^{2}x$ and defining $y:=\cos(\lambda_{i,k}D)$, we have the equation

[TABLE]

which can be rewritten as cubic in $y$

[TABLE]

We can explicitly solve for the roots of the cubic, one of which is always real.

Finally, substituting $\lambda_{i,k}=\frac{1}{D}\arccos(y)$ into the real equation (35) gives an implicit equation of the instability in terms of the system parameters,

[TABLE]

This is the Hopf curve when $k=0$ and $\hat{J}_{k}=J_{0}$ and a Turing-Hopf curve when $k=1$ and $\hat{J}_{k}=J_{1}/2$.

5.5 Numerical methods

5.5.1 Simulations

Simulations of the mean-field approximation (2) were done using the forward Euler method with a step size of $dt=0.001$. The ring $(0,2\pi]$ was descritized as $x_{i}=-\pi+\frac{2\pi}{N}i$, $i=1,...,N$, and $N=100$ ($N=1000$ for MF simulations of spatiotemporal dynamics in Fig. 5).

Simulations of the sLIF network (1) were done using the stochastic forward Euler method with step size $dt=0.001$ with $N=1000$. The connectivity matrix $\mathbf{J}\in\mathbb{R}^{N\times N}$ was randomly sparse, with connection probability $p=0.5$. The connections between neuron pairs $(i,j)$ were independently sampled as Bernoulli random variables with probability $p=0.5$. Each non-zero connection then has synaptic weight given by $J_{ij}=\frac{2\pi}{pN}(J_{0}+J_{1}\cos(x_{i}-x_{j}))$. Self-connections are excluded by setting $J_{ii}=0$. We generated spikes by sampling $dn_{i}$ as Bernoulli random variables at each time step with success probabilities $f(v_{i})dt$. Finally, synaptic delays were incorporated with a time-shift of $N_{\text{delay}}=\lfloor D/dt\rfloor$ to the synaptic input.

The initial conditions must be defined for $t\in[-D,0]$ due to the presence of the delay. Most patterns were initialized with constant initial conditions or the stationary form $A+B\cos(x)$. For the traveling wave pattern seen in Fig. 5F, we used an initial condition of the form $A+B\cos(x+tv)$ with $v\neq 0$ to induce wave propagation. Details for specific initial conditions are found in the code repository https://github.com/lcforbes4/sLIFspatioTemporal.

5.5.2 Inducing network transitions in bistable regions

To switch from the high activity state to the oscillatory state (Fig. 3 I-J), a positive pulse with short duration was applied to the entire network. To turn off the oscillatory state and return to the high activity state, a negative pulse was applied during the peak of the oscillations to suppress oscillatory activity. The timing, amplitude, and duration of the negative pulse were fine tuned to the particular simulation and set of parameters to achieve the return to the high activity state. See figure caption for specific pulse parameters.

In all three simulations (Fig. 8A–C), a spatially dependent drive of the form $-2\cos(x)$, applied for a duration of 1, was used to induce a transition to spatially patterned network activity. To transition from the spatially patterned activity state to the high activity state (Fig. 8, A and C), a pulse with amplitude -2 was applied to the quiescent middle third of the network (neurons $N/3:2N/3$). To transition from the spatially patterned activity state to the quiescent state (Fig 8 B), a pulse with amplitude -2 was applied to the active outer two thirds of the network (neurons $0:N/3$ and $2N/3:N$).

5.5.3 Continuation of oscillatory solutions

Continuation of oscillatory solutions was performed using a one-dimensional, spatially homogeneous mean-field reduction (2). This approximation assumes zero spatial modulation ($J_{1}=0$) in the coupling function, such that each neuron in the network evolves according to the spatially independent delay-ODE. We used the Matlab package DDE-Biftool, which was developed for bifurcation analysis of delay differential equations ([engelborghsNumericalBifurcationAnalysis2002, sieberDDEBIFTOOLManualBifurcation2016]). First, we computed and continued a branch of equilibria along with the spectrum of its linearization to identify the location of the Hopf bifurcation. Then we initialized a small-amplitude orbit near the Hopf bifurcation to construct and continue a branch of oscillatory solutions. Stability of periodic orbits was determined by simultaneously computing Floquet multipliers of the linearized periodic map. We found the Hopf bifurcation to be subcritical, with a branch of unstable oscillations emerging backwards out of the bifurcation point as $J_{0}$ and $E$ decrease and $D$ increases.

5.5.4 Continuation of Turing patterns

The periodic Turing patterns computed in Section 3.4 are steady-state solutions of the mean-field approximation (2) with spatial connectivity ($J_{1}\neq 0$) but without delay ($D=0$). To compute these solutions numerically in Matlab, we implemented a Newton-GMRES FFT method [rankinContinuationLocalizedCoherent2014]. We discretized (2) using a Fourier spectral method with $N=2^{12}$ modes on the periodic domain $[0,2\pi)$. We denote the discretized time-independent version of (2), as the non-linear system $F_{1}(v)=0$ where $v\in\mathbb{R}^{N}$ is the spatially discretized steady-state solution.

To obtain a good initial guess for the Turing pattern, we ran a simulation initialized with a small perturbation (of order $10^{-1}$) of the spatially homogeneous steady state. To fix the continuous translational symmetry on the ring, we appended a phase condition on our system,

[TABLE]

which pins the solution, and modified $F_{1}$ with an additional drift term $c\partial_{x}v$. The imposition of this condition requires an additional unknown or dummy variable $c$. We then solved the resulting augmented system $F(v,c)=\begin{bmatrix}F_{1}(v)+c\partial_{x}v\\ F_{2}(v)\end{bmatrix}=0$ for $v$ and $c$ using a Newton–GMRES method, where vector products in the Jacobian were computed spectrally and GMRES without preconditioning was used to compute the Newton step.

Once the initial Turing pattern was obtained, we performed numerical continuation in a single bifurcation parameter $\beta$ using secant continuation. To do this, we introduced the augmented variable $z=[v,c,\beta]$ and appended the orthogonality condition

[TABLE]

where $z_{0}=[v_{0},c_{0},\beta_{0}]$ is the previously computed solution in the continuation and $sec\in\mathbb{R}^{N+2}$ is the normalized secant vector between the two most recent solutions. The variable $c$ stays uniformly close to zero throughout the continuation. The complete augmented system used for secant continuation is then

[TABLE]

At each step of the continuation, we computed the spectral stability of the linearized operator near the origin with leading eigenvalues of Jacobian of the original system $F_{1}$ (without the phase and orthogonality conditions) and Matlab’s eigs function. In Fig. 7, we plot the real part of the eigenvalues with $Re(\lambda)>-5$ along the entire continuation curve in the parameter $E$. The continuation starts at the Turing bifurcation and initially has an eigenvalue with $Re(\lambda)>0$, indicating instability and likely a subcritical bifurcation. When the largest eigenvalue drops below zero (which occurs near $E=2.90$), this marks the location of a fold of Turing patterns and the continuation’s switch onto the stable branch at the upper fold. The continuation then follows the stable branch until the lower fold where a real part of an eigenvalue becomes positive ($E=0.22$) and solutions are unstable until terminating at the non-smooth Turing bifurcation. The eigenvalue with $Re(\lambda)=0$, which remains present throughout the continuation, is a consequence of the translational invariance of solutions on the ring. The locations of the two fold curves in the phase diagram in Fig. 8D were identified by determining the location of the fold in multiple single-parameter continuations of the Turing pattern at varying parameters, not through continuation of the fold itself.

5.6 Tracking the fold points to identify the boundaries of multi-stable regions

The regions Q-S and S-Q-H only exist in a small neighborhood near the co-dimension 2 bifurcation of the saddle-node and Turing bifurcations (Fig. 8D). The boundaries of these and the surrounding regions are defined by the various instabilities which occur in this neighborhood: the subcritical Turing bifurcation (Fig. 8D, solid red) and nearby fold of Turing patterns (Fig. 8D, dashed green), the saddle-node bifurcation (Fig. 8D, solid black), and the subcritical non-smooth Turing (Fig. 8D, solid orange) with a nearby fold of the Turing patterns (Fig. 8D, dashed blue). We refer to this second fold as the ‘lower’ fold since it occurs at lower $E$ and $J_{0}$ values than the ‘upper’ fold. Since both of the spatial instabilities are subcritical in this neighborhood, the boundaries of the stable Turing pattern region are determined by the folds.

To track the fold points and understand the interaction of various instabilities in this region, we computed bifurcation diagrams across a range of parameter values (Fig. 8E-L). These bifurcation diagrams are slices of the phase diagram (Fig. 8D) with an additional dimension showing $v$. The upper row (Fig. 8E-H), which shows continuation of the Turing pattern in the parameter $E$, are vertical slices of Fig. 8A at different $J_{0}$ values. The lower row (Fig. 8I-L), showing continuation in $J_{0}$, are horizontal slices of Fig. 8A at various $E$ values.

In each of these bifurcation diagrams, the instabilities’ locations are marked with a circle whose color matches the corresponding curve in the phase diagram (Fig. 8D). For example, consider Fig. 8E, where five different instabilities (uniquely colored circles) are visible. The saddle-node bifurcation, where $v_{+}$ and $v_{-}$ collide and vanish, is marked with a black circle. The Turing bifurcation, where a unstable Turing pattern emerges from $v_{+}$, is shown with a red circle. Both a non-smooth Turing and saddle-node bifurcation occur at the same point, shown with an orange circle, where $v_{-}$, $v_{Q}$, as well as an unstable Turing pattern emerge as $E$ is decreased. Lastly, there are two folds of the Turing pattern: the upper fold (green circle), connected to the Turing bifurcation via an unstable branch of Turing patterns, and the lower fold (blue circle), connected to the non-smooth Turing via an unstable branch. The stable branch connecting the two folds determines where there is stable spatially patterned activity.

All three regions of multi-stability with the Turing pattern are present when $J_{0}=4$ (Fig. 8E). First, the S-H region, where only the Turing pattern and $v_{+}$ are stable, lies between the non-smooth Turing bifurcation and the upper fold. The Q-S region, where there is bistability of $v_{Q}$ and the Turing pattern, lies between the lower fold and the Turing bifurcation. The last region of multi-stability between three solutions (the Turing pattern, $v_{Q}$, and $v_{+}$) lies between the Turing and non-smooth Turing bifurcations.

As $J_{0}$ increases (Fig. 8F, $J_{0}=5$), the S-H and Q-S regions vanish and S-Q-H is the only remaining multi-stable region with the Turing pattern. The two fold points (green and blue) have moved closer together and passed by the locations of the spatial instability (red) and the non-smooth Turing (orange). Therefore the stable branch of Turing patterns between them is smaller and lies entirely in the region where both $v_{+}$ and $v_{Q}$ are also stable.

In addition to the loss of two of the bistable regions at $J_{0}=4$, the location of the spatial instability is now on the $v_{-}$ branch after having just passed through the location of the saddle-node bifurcation (black circle) in a co-dimension 2 bifurcation. Despite being a secondary instability, the emerging branch of unstable Turing patterns leads to a stable branch through a fold (Fig. 8F, green). We mark this secondary instability of $v_{-}$ with an open red circle (largely obscured by the black circle) to distinguish it from the Turing bifurcation of $v_{+}$ (marked with a solid red circle in the other panels).

If $J_{0}$ is further increased (Fig. 8G, $J_{0}=5.75$), the distance between the two folds continues to decrease and thus the size of the S-Q-H region decreases. They continue to approach each other until they collide and vanish in a codimension-2 cusp point (Fig. 8D, empty square) along with the stable branch between them, ending the S-Q-H region. Beyond the cusp point, an unstable branch of Turing patterns remains, connecting the spatial instability to the non-smooth bifurcation (Fig. 8H, $J_{0}=6.5$).

We can also see the emergence and disappearance of these bistable regions as the resting potential passes through threshold by looking at the continuation of the Turing pattern in the parameter $J_{0}$, at differing $E$ values. When above threshold ($E>1$, Fig. 8I), there is no lower fold and thus the only bistable region with spatially patterned activity is S-H. The stable branch of Turing patterns exists for all $J_{0}$ values below the upper fold extending off to $-\infty$ and slowly decreasing in amplitude as $J_{0}$ is decreased. The stable branch also decreases in amplitude as $E$ is decreased until at the threshold ($E=1$, Fig. 8J), it no longer extends to $-\infty$ and instead meets $v_{Q}$ at the non-smooth supercritical Turing bifurcation (Fig. 8J, $J_{0}\approx-1.3$). This is where the non-smooth Turing bifurcation becomes subcritical in the bifurcation parameter $E$ and the lower fold emerges. Recall there are multiple transitions between homogeneous solutions at threshold. For $J_{0}<2$, the above threshold state $v_{+}$ meets the below threshold state $v_{Q}$ (i.e. $v_{+}=v_{Q}=1$). For $J_{0}>2$, there is the non-smooth saddle bifurcation where both $v_{-}$ and $v_{Q}$ emerge. Thus the homogeneous solution at $E=1$ is stable for $J_{0}<2$ and half-stable (dash-dotted) for $J_{0}>2$. Then below threshold ($E<1$, Fig. 8K), there is a lower fold and an unstable branch of Turing patterns that extends off to $+\infty$ in $J_{0}$. Consequently, there are now Q-S and S-Q-H regions. As $E$ decreases, the two folds move closer (Fig. 8L) leaving only the S-Q-H region. The folds eventually meet in the codimension-2 cusp point (Fig. 8A, blue dashed, green dashed, empty square) ending the branch of stable Turing patterns and spatially patterned activity region.