We study localized patterns in an exact mean-field description of a spatially extended network of quadratic integrate-and-fire neurons. We investigate conditions for the existence and stability of localized solutions, so-called bumps, and give an analytic estimate for the parameter range, where these solutions exist in parameter space, when one or more microscopic network parameters are varied. We develop Galerkin methods for the model equations, which enable numerical bifurcation analysis of stationary and time-periodic spatially extended solutions. We study the emergence of patterns composed of multiple bumps, which are arranged in a snake-and-ladder bifurcation structure if a homogeneous or heterogeneous synaptic kernel is suitably chosen. Furthermore, we examine time-periodic, spatially localized solutions (oscillons) in the presence of external forcing, and in autonomous, recurrently coupled excitatory and inhibitory networks. In both cases, we observe period-doubling cascades leading to chaotic oscillations.

Spatially extended networks of spiking model neurons are capable of producing spatiotemporal patterns that are observed experimentally in neuronal tissues. An important tool in investigating such patterns is low-dimensional neural field models, which describe the macroscopic dynamics of such networks. Many neural field models are derived heuristically and do not fully describe the dynamics of the underlying network of spiking neurons. We utilize a recently derived mean-field description for networks of quadratic integrate-and-fire (QIF) neurons, which yield an accurate description of the mean firing rate and the mean membrane potential of the network. Contrary to other neural field models, this model only contains nonlinearities of the mean-field variables up to quadratic order, which are amenable to the development of Galerkin methods for the numerical approximation of the problem. This allows us to study the bifurcation structure of stationary localized solutions and of time-varying localized solutions up to the emergence of chaos.

Localized states in neuronal networks, so-called bumps, are related to working memory1,2 and feature selectivity,3 whereby neurons encoding similar stimuli or features show an increased firing rate for the duration of the related cognitive task. Neural fields are well-known coarse-grained models of spatiotemporal neuronal activity,4–6 capable of reproducing dynamic phenomena found experimentally, such as traveling waves, temporal oscillations, and spatially localized states.7,8 A challenge faced in the derivation of neural field models is to establish an accurate mean-field description of the spiking dynamics of the underlying microscopic neural network. Classical neural field models recover the microscopic dynamics only in the limit of slow synapses,9 and the derivation of neural mass or neural field models from networks of spiking neurons is still an active area of research.10–18 In addition, in neural fields, the network firing rate is not an emergent quantity but rather the result of a modeling choice.

Some limitations can be overcome if the microscopic model is a heterogeneous network of synaptically coupled θ or QIF neurons, subject to random, Cauchy-distributed background currents. Recently, it has been shown that heterogeneous networks of θ- and QIF neurons admit an exact mean-field description,19,20 which has later been expanded to spatially extended networks.21–24 In the thermodynamic limit, the network admits an exact mean-field description in terms of the network mean firing rate and voltage20 or in terms of a complex-valued order parameter.21,25

We study a network of n quadratic integrate-and-fire neurons,

V˙i=Vi2+ηi+Jsi,i=1,,n,
(1)

where Vi is the membrane potential of the ith neuron, ηi an intrinsic current, si the synaptic input, and J>0 a global coupling parameter. The ith neuron emits a spike when Vi reaches the firing threshold Vθ, and Vi is reset immediately to Vr. Following Ref. 20, we distribute {ηi:i=1,,n} according to a Lorentzian distribution using the formula ηi=η+Δtan[π/2(2jn1)/(n+1)], where η is the center and Δ is the half-width of the Lorentzian distribution, respectively. An important difference in the model considered in the present paper is that neurons are distributed in space, in a domain Ω=(L/2,L/2], with L1, at evenly spaced positions {xi=iδxL/2:i1,,n}, with δx=L/n. We associate with each lattice point xi a random component of the vector {ηj}, without repetitions. The synaptic current received by a neuron is determined by the synaptic footprint as follows:

si(t)=1nj=1nw(xi,xj)k:tjk<tta(tt)δ(ttjk)dt,
(2)

where w(x,y) models the synaptic coupling strength between neurons from position y to position x in the network and tjk is the emission time of the kth spike of the jth neuron. The kernel a(t) represents synaptic activation in response to incoming spikes, e.g., exponential synapses a(t)=et/τs/τs with synaptic time scale τs.26,27 In the mean-field description, we will let τs0, i.e., a(t)=δ(t). We note in passing that, as demonstrated for leaky integrate-and-fire neurons,28 the mean-field derivation and the limit τs0 do not commute.

A neural field model that describes without approximation of the average firing rate r(x,t) and the average membrane potential v(x,t) of the spatially extended networks presented above has been developed recently,23 

tr=Δπ+2rv,tv=v2+η+Jwrπ2r2,xΩ.
(3)

This neural field model inherits the coupling parameter J and the parameters η and Δ from the microscopic, spiking network. The mean-field description is exact in the limit n and Vθ=Vr. The spatial coupling, or synaptic footprint, is given by the integral operator,

[wr](x)=Ωw(x,y)r(y)dy,xΩ.
(4)

For the concrete calculations presented below, we will assume Ω=R or Ω=(L/2,L/2]S with L1 (a ring with large width). We use Ω=R for the theoretical framework in Sec. II and the Hermite–Galerkin method and use Ω=(L/2,L/2]S for the Fourier–Galerkin method and the numerical integration of both the macroscopic and microscopic model equations, as periodicity is enforced in this setting. We will study the model with a variety of kernels but, unless stated otherwise, we assume, with a small abuse of notation, w(x,y)=w(|xy|) and

w(x)=e|x|14e|x|/2,
(5)

when Ω=R, or take the 2L-periodic extension of w when Ω=(L/2,L/2]S. Hence, our default synaptic kernel will depend on the distance between two points in Ω and will have long-range inhibition and short-range excitation. With these choices, wr is a convolution and Ωw(y)dy=1.

This neural field model is related to mean-field descriptions of networks of theta neurons19,21,24,25 and was obtained using the Ott–Antonsen ansatz.29 It retains the transient dynamics of the microscopic network, including spike synchrony, and has, therefore, a richer dynamic repertoire than purely rate-based models.30 The derivation of such neural field models is analogous to mean-field approaches for spatially extended networks of phase-coupled oscillators.31 An example of localized solutions in this model is shown in Fig. 1(a), alongside a numerical simulation of the microscopic system of spiking neurons in Fig. 1(b).

FIG. 1.

The formation of a stationary localized solution (bump) in (a) the QIF neural field model and (b) the corresponding network of spiking neurons. The bump solution is induced in the model with synaptic kernel (5), by applying a localized transient current I(x,t)5 if (x,t)[2.5,2.5]×[0,5] and I(x,t)0, otherwise. In (b) we show a rastergram of 300 of 105 neurons used in the simulation of the underlying spiking network. Parameters: Δ=2, J=152, and η=10.

FIG. 1.

The formation of a stationary localized solution (bump) in (a) the QIF neural field model and (b) the corresponding network of spiking neurons. The bump solution is induced in the model with synaptic kernel (5), by applying a localized transient current I(x,t)5 if (x,t)[2.5,2.5]×[0,5] and I(x,t)0, otherwise. In (b) we show a rastergram of 300 of 105 neurons used in the simulation of the underlying spiking network. Parameters: Δ=2, J=152, and η=10.

Close modal

The main aim of the present paper is to study spatiotemporal localized patterns supported by this model, such as the one presented in Fig. 1. Our investigation will be primarily numerical, and, therefore, we will also introduce several numerical schemes for the approximation of the QIF neural field model. The paper is structured as follows: in Sec. II, we discuss analytical methods to study stationary solutions and their bifurcations; in Sec. III, we introduce the numerical methods used to perform numerical bifurcation analysis of stationary and time-periodic localized structures, which are presented in Secs. IV and V, respectively; and we make a few concluding remarks in Sec. VI.

Stationary states of Eq. (3) are determined by the conditions tr=tv=0. Bounded solutions with r(x)>0 satisfy

0=Δ24π2r2+η+Jwrπ2r2,v=Δ2πr.
(6)

The model supports both uniform and non-uniform steady states, which we discuss below in further detail.

Solutions to (6) depend in general on x. Spatially uniform solutions, for which Jwr=Jr, satisfy the quartic equation,

r4Jπ2r3ηπ2r2Δ24π4=0,
(7)

which has the following four solutions:

r1,2=J4π2+12S±2pS2q/S,r3,4=J4π212S±2pS+2q/S,
(8)

where p, q, and S are given by

p=ηπ238J2π4,q=J38π612Jηπ4,
(9)
S=23p+13(Q+R0/Q),
(10)

with

Q=(12(R1+R124R03))1/3,R0=η2π43Δ2π4,R1=2η3π6274J2Δ2π818ηΔ2π6,
(11)

respectively. Physically relevant solutions are positive and real, and an inspection of the equations above reveals that r4 must be discarded, and the system admits either one or three homogeneous steady states. At sufficiently small (large) η, only one stable fixed point exists, represented by r3 (r1); also, there exists an interval in parameter space where the stable solutions r1, r3 coexist with r2, which is unstable. The conclusions presented above justify the bifurcation diagram found in Refs. 20 and 23 and reported in Fig. 2(a).

Loci of saddle-node bifurcations in the (η,J)-plane can be found by setting dη/dr=0 in the first equation in (6) which combined with (7) yield a parameterization in r,

ηsn=π2r23Δ24π2r2,Jsn=2π2r+Δ22π2r3,
(12)

or, more explicitly,

Jsn=2π2ηsn±ηsn23Δ2+2π2Δ2(ηsn±ηsn23Δ2)3/2,
(13)

where ± denote two bifurcation branches of saddle-node bifurcation which collide at a cusp

(ηc,Jc)=(3Δ,4π323Δ).
(14)

A first step toward the construction of heterogeneous steady states is the determination of Turing bifurcations, which mark points in parameter space where a spatially uniform solution becomes unstable to spatially periodic patterns. We remark that it is known that spatially extended networks of QIF or θ neurons display this instability,23,24 and here, we present an analytic determination of the loci of such bifurcations in parameter space. Turing bifurcations of a homogeneous steady state (r,v) can be identified by linear stability analysis of the model equations in Fourier space, which results in the following eigenvalue problem:

λ(k)(r~v~)=(2v2rJw^(k)2π2r2v)(r~v~):=A^(k)(r~v~),

where w^(k) is the Fourier transform of the connectivity kernel. A sufficient condition for a Turing bifurcation is the existence of a critical wavenumber kc>0 for which detA(kc)=0, which yields

rT=1πηTw^(kc)2(2w^(kc)±12ηT2w^(kc)2(2w^(kc))2(2+w^(kc))Δ22w^(kc)
(15)

and

JT=1w^(kc)(Δ22π2rT3+2π2rT).
(16)

Combining (15) and (16) results in an equation for the loci of the Turing bifurcation in the (η,J)-plane. As kc0, the resulting equation recovers (13), since w^(kc)1. This analytic result agrees well with the numerical calculations of these loci, which will be presented further below.

After studying uniform and spatially periodic steady states, we move to the construction of localized steady states supported by the models. One strategy to study localized stationary states in nonlinear models posed on R is to construct solutions to boundary-value problems derived from the model’s steady state equations.32–36 With this approach, localized steady states correspond to homoclinic orbits of a dynamical system in which x plays the role of time (hence the term spatial dynamics).

In this section, we make some preliminary considerations on the spatial dynamics of steady state solutions to (3), although we do not explicitly study the associated spatial-dynamical system, as we will construct our solutions numerically in the following sections. Using the positions,

u(r)=Δ24π2r2η+π2r2=v(r)2η+π2r2,
(17)
f(u(r))=Jr,
(18)

the steady state equation (6) is recast as

0=u+wf(u).
(19)

We note that Eq. (19) is formally equivalent to the Amari steady state equation;6 in Amari’s theory, u represents the voltage, whereas in this case u combines the steady state’s voltage and rate and scales as uv2η for small r and uπ2r2η for large r, respectively.

Importantly, the identification with the Amari equation allows us to use spatial dynamics to characterize localized steady state solutions.37–40 The Fourier transform of w is of the form w^(k)=Π1(k2)/Π2(k2), with Πi being a polynomial of order i, if w is the bi-exponential kernel (5). Hence, the integral kernel can be regarded as Green’s function of a differential operator. In particular, the bi-exponential kernel (5) leads to the differential equation,

u54u+14u14f(u)+74[f(u)]=0,
(20)

where a prime denotes differentiation with respect to x. The equation above can be cast as a 4D, first-order spatial-dynamical system in the vector (u,u,u,u), which we omit here for brevity. To construct localized solutions to (3) we proceed in the same spirit as32–37,39,40 to each homogeneous steady state of (3) corresponds to one value rj in (8), and hence one value uj in (17), and one constant solution (uj,0,0,0) to (20); in addition, there exists a region in parameter space, where r1 and r3 coexist and are stable [see also Fig. 2(a)]. A localized steady state of (3) is identified with a bounded, sufficiently regular function u:RR which satisfies (20) with boundary conditions

limx(u(x),u(x),u(x),u(x))=(u1,0,0,0),limx+(u(x),u(x),u(x),u(x))=(u3,0,0,0).

Furthermore, we note that the quantity,

H(u,u,u,u,x)=uu12(u)258(u)2+18u214uf(z)dz+74x[f(u)](z)u(z)dz,

is conserved in the sense that, if (20) holds, then

ddxH(u(x),u(x),u(x),u(x),x)=0.

Therefore, we expect to construct a localized stationary state in a region of parameter space where H(u1,0,0,0,0)=H(u3,0,0,0,0). With a slight abuse of notation, we write this condition in terms of the variable r, as H(r1)=H(r3), where H is given by

H(r)=2ηr32Jr2+43r3.
(21)

In analogy with the literature mentioned above, we called Maxwell points the values on the (η,J)-plane where the condition H(r1)=H(r3) is met. We display the Maxwell point for our standard parameter set in Fig. 2(b), and we plot the locus of Maxwell points and the bistability region in Fig. 2(c).

FIG. 2.

(a) Solution branches of uniform solutions around the bistable regime with saddle-node bifurcations (S). (b) The associated conserved quantity H, as defined in (21), showing a Maxwell point (M) at ηM9.69. (c) A plot of the bistable region of homogeneous states (shaded) on the (J,η)-plane delimited by saddle-node bifurcations (S), emanating from a cusp (C). We also show the locus of Maxwell points (M) on the plane. Parameters: J=152 and Δ=2.

FIG. 2.

(a) Solution branches of uniform solutions around the bistable regime with saddle-node bifurcations (S). (b) The associated conserved quantity H, as defined in (21), showing a Maxwell point (M) at ηM9.69. (c) A plot of the bistable region of homogeneous states (shaded) on the (J,η)-plane delimited by saddle-node bifurcations (S), emanating from a cusp (C). We also show the locus of Maxwell points (M) on the plane. Parameters: J=152 and Δ=2.

Close modal

As anticipated in Secs. I and II, stationary states beyond onset are computed numerically; hence, we present, in this section, several numerical schemes used in the upcoming computations. In preparation for presenting the schemes, we rewrite the model as an ordinary differential equation on a function space. To simplify the notation, we apply, in this section, the scaling rr/π, JπJ to (3) and obtain

r˙=Δ+2rv,v˙=η+v2r2+Wr,
(22)

where W is the integral operator defined as (Wr)(x)=J(wr)(x)=JΩw(x,y)r(y)dy. In the system above, we assume r,v:RLρ2(Ω,C) (also denoted by Lρ2(Ω)), that is, at each time t, r(t) and v(t) belong to a weighted Lebesgue space of complex-valued functions defined on Ω, with inner product,

f,gρ=Ωf(x)g(x)ρ(x)dx,

and norm fρ=f,f1/2. Note that the subscript ρ will be omitted when ρ(x)1. We assume that, once complemented with initial conditions, system (22) defines a well-posed Cauchy problem on Lρ2(Ω)×Lρ2(Ω).

Galerkin schemes are derived by introducing a complete orthogonal basis {φi:iN} for the weighted space Lρ2(Ω) and seeking an approximation in the n-dimensional subspace spanned by {φi:iΛn}, where Λn is an index set with n elements, as follows:

rn(x,t)=iΛnRi(t)φi(x),vn(x,t)=iΛnVi(t)φi(x).

A Galerkin scheme for (22) is then given by

φi,r˙n+Δ+2rnvnρ=0,φi,v˙n+η+vn2rn2+Wrnρ=0,iΛn,

that is,

R˙i=αiΔ+2j,kΛnγijkRjVk,V˙i=αiη+jΛnβijRj+2j,kΛnγijk(VjVkRjRk),

for iΛn, with coefficients given by

αi=φi,1ρ,βij=φi,Wφjρ,γijk=φi,φjφkρ.

1. Fourier–Galerkin scheme

When Ω=(L/2,L/2]S, the functions r(t) and v(t) are L-periodic. Therefore, we choose the Fourier basis φj(x)=exp(ij2πx/L), jZ, which is a complete orthogonal basis for L2(Ω). The index set for this case is Λn={n/2,,n/21} with n even. Exploiting the trigonometric properties of the Fourier basis, we obtain

αi={Lif i = 0, 0otherwise, γijk={Lif i+j+k=0,0otherwise.

In passing, we note that βij can also be expressed compactly, in terms of the Fourier coefficients of the kernel w, if the operator W is convolutional. In addition, requiring r and v to be real-valued implies (Ri,Vi)=(Ri,Vi). We call this method the Fourier–Galerkin scheme.

2. Hermite–Galerkin scheme

When Ω=R, a natural basis for the Galerkin scheme is given by the Hermite polynomials,

φj(x)=Hj(x)=(1)jexp(x2)djdxjexp(x2),

which are a complete orthogonal set for Lρ2(Ω,R) with weight ρ(x)=exp(x2). For this scheme, Λn={0,,n1}. To avoid problems with the numerical evaluations of φj for large |x|, we derive an alternative scheme, which uses inner products with weight ρ(x)1, as the Fourier–Galerkin scheme. We seek a solution to (22) in the form

r=R0+r~,v=V0+v~,

with R0 and V0 constant in x and r~,v~L2(Ω,R). This leads to the system

R˙0=Δ+2R0V0,V˙0=η+V02R02+JR0,r~˙=2R0v~+2V0r~+2r~v~,v~˙=2V0v~2R0r~+v~2r~2+Wr~,

in which the homogeneous background dynamics for (R0,V0) is decoupled from (r~,v~) and follows the spatially clamped QIF mean field.20 Since the Hermite functions,

φj(x)=exp(x2/2)Hj1(x),jN>0,

are an orthogonal set for L2(Ω,R), an approximation to r~,v~ is sought in the space spanned by φj, with jΛn={1,,n}, giving the scheme

R˙0=Δ+2R0V0,V˙0=η+V02R02+JR0,R˙i=2jΛn(R0Vj+V0Rj)+2j,kΛnγijkRjVk,V˙i=jΛn[2V0Vj+(βij2R0)Rj]+j,kΛnγijk(VjVkRjRk)

for iΛn. We call this method the Hermite–Galerkin scheme.

A Fourier collocation scheme can be derived in the case Ω=(L/2,L/2]S. This method, which has been used in the past for Amari neural field models40,41 and the QIF neural field model,23 represents (rn,vn) by its values at the grid points xj=L+2Lj/n, jΛn={1,,n},

R˙i=Δ+2RiVi,V˙i=η+Vi2+Ri2+(Wrn)i,

and evaluates (Wrn)i either with a quadrature rule or, more efficiently, with a pseudospectral evaluation if W is convolutional.

To the best of our knowledge, the methods presented above are novel, and we leave the analysis of the numerical properties of these schemes to a separate publication. The calculations presented here have been tested against event-driven simulations of a large network of spiking neurons. We employ our schemes as follows: the Fourier collocation scheme with n=5000 is generally used for time simulations to obtain accurate initial guesses for the continuation. However, we observed that time-periodic orbits are reproduced with a similar accuracy by the Hermite–Galerkin scheme with just n=50 modes; hence, we select this scheme to continue periodic orbits. Finally, we use the Fourier–Galerkin scheme with n=200 for bifurcation analysis of steady states on large domains, when solutions are non-localized.

We compare the results of the QIF neural field model with the dynamics of the underlying network of spiking neurons. We integrate Eq. (1) using the Euler method, with time step Δt=104. The domain is chosen to be Ω=(L/2,L/2] with periodic boundary conditions. We choose L=50 and 5×105 model neurons, which ensures a good correspondence to the neural field model (δx=104); with L=50, the relative error of the normalization of the synaptic kernel, Ωw(y)dy=1, is less than 105. We note here that the microscopic description only matches the mean-field description if δx is much smaller than the characteristic length scale of the synaptic kernel w. The synaptic input to each neuron is computed with Eq. (2). We follow Ref. 20 in computing the synaptic integration across a time window with aτ(t)=Θ(τt)/τ, τ=103 and in setting Vθ=Vr=100 with a refractory period of 2/Vi once the ith neuron has exceeded Vθ. The latter approximates the limit Vθ=Vr. The refractory period is rounded to the nearest multiple of Δt, and after the refractory period, Vi is set to Vi. Rasterplots are generated using a subset of 103 randomly chosen model neurons.

We use the numerical schemes presented in Sec. III to study the bifurcation structure of stationary localized solution to the QIF mean field model. We initially study the model with our default excitatory-inhibitory kernel (5) and then show that a snaking bifurcation scenario is supported when the kernel is switched to a homogeneous oscillatory kernel, or to a kernel with harmonic heterogeneities, similarly to what is found for Amari neural field models.

We set w as in (5), generate a stationary localized solution by numerically integrating the model equations in time, and then implement the Fourier–Galerkin scheme to continue the localized solutions in η, using AUTO.

In Fig. 3, we show the bifurcation diagram of localized solutions. Across a range of parameters, these occur as a pair of one wide, stable solution and one narrow, unstable solution. The solution branch connects to the branch of uniform solutions at points where Turing bifurcations occur, which also give rise to a branch of periodic solutions. Using the Fourier basis, it can be shown that the stable solution branch approaches the Maxwell point asymptotically, and solutions grow wider, which resemble two (stationary) interacting wave fronts. Because of the periodic boundary conditions, the solution branch grows larger again and forms another stable/unstable solution pair of locally low activity (not shown). The latter could be regarded as stationary versions of traveling anti-pulses reported in Refs. 42 and 43.

FIG. 3.

(a) Bifurcation diagram in η of localized solutions (blue), periodic solutions (green), and uniform solutions (black). (b) Exemplary profiles of unstable narrow (1), stable narrow (2), unstable wide (3), and stable periodic (4) localized solutions, close to the Maxwell point [vertical line in (a)]. Parameters: Δ=2 and J=15Δ.

FIG. 3.

(a) Bifurcation diagram in η of localized solutions (blue), periodic solutions (green), and uniform solutions (black). (b) Exemplary profiles of unstable narrow (1), stable narrow (2), unstable wide (3), and stable periodic (4) localized solutions, close to the Maxwell point [vertical line in (a)]. Parameters: Δ=2 and J=15Δ.

Close modal

Because stable solutions are of particular interest, we present a two-parameter bifurcation diagram [Fig. 4(a)] of the saddle-node bifurcations that delimit the branch of stable solutions. As expected, the locus of saddle-node bifurcations of localized states encloses the Maxwell point. In addition, the loci of saddle-node bifurcations of localized and uniform steady states meet at two separate cusps, as shown in Fig. 4(b).

FIG. 4.

(a) The parameter space in which stable localized solutions exist is delimited by loci of saddle-node bifurcations. They can be approximated by the saddle-node bifurcations of spatially uniform solutions and the Maxwell point. (b) Inset of (a), showing additionally the loci of Turing bifurcations and saddle-node bifurcations of periodic solutions. The saddle-node bifurcations of the bump solutions form a cusp, where the Turing bifurcation changes from supercritical to subcritical. Parameter: Δ=2.

FIG. 4.

(a) The parameter space in which stable localized solutions exist is delimited by loci of saddle-node bifurcations. They can be approximated by the saddle-node bifurcations of spatially uniform solutions and the Maxwell point. (b) Inset of (a), showing additionally the loci of Turing bifurcations and saddle-node bifurcations of periodic solutions. The saddle-node bifurcations of the bump solutions form a cusp, where the Turing bifurcation changes from supercritical to subcritical. Parameter: Δ=2.

Close modal

The bifurcation behavior of localized solutions described above is robust to changes in coupling parameters but, as we shall see below, it is strongly affected by changes in the kernel.

Homoclinic snaking is a phenomenon that describes the formation of multiple, coexisting localized solutions in spatially extended models. Steady states are arranged in branches of intertwined snaking bifurcation diagrams, connected via ladders.33,35,44–46 Adopting the spatial-dynamics approach outlined above, localized solutions are interpreted as homoclinic orbits to a fixed point. Snaking solution branches correspond to symmetries of the problems, which are broken along the ladder branches.34,36 This scenario is not limited to partial differential equations, but have also been studied in the nonlocal Swift–Hohenberg equation,47 as well as in neural field models.37,38,40,41,48–50

In the simplest setting, localized snaking solutions are found in regions of parameter space where there is bistability between a stationary homogeneous state and a periodic state. In nonlocal neural fields, homoclinic snaking has been observed with the following homogeneous damped-oscillatory kernel:39 

w(x)=1+b24beb|x|(bsin|x|+cosx),
(23)

which we now adopt also for the QIF neural field model. This kernel leads to a subcritical Turing bifurcation of the lower stable branch of uniform solutions, from which an unstable branch of spatially periodic solutions emerges. This branch undergoes a saddle-node bifurcation, where spatially periodic solutions become stable. Eventually, the branch connects to the upper stable branch of uniform solutions (see Fig. 5).

FIG. 5.

(a) The damped-oscillatory kernel (23). (b) Fourier transform of the kernel, which shows a maximum at a non-zero wavenumber kc. (c) As η is varied, the branch of homogeneous steady states (black) undergoes a Turing bifurcation, from which a branch of periodic solutions with wavenumber kc emerges (green). Parameters: J=39, Δ=1, and b=0.4.

FIG. 5.

(a) The damped-oscillatory kernel (23). (b) Fourier transform of the kernel, which shows a maximum at a non-zero wavenumber kc. (c) As η is varied, the branch of homogeneous steady states (black) undergoes a Turing bifurcation, from which a branch of periodic solutions with wavenumber kc emerges (green). Parameters: J=39, Δ=1, and b=0.4.

Close modal

As anticipated, spatially localized snaking solutions are found in this region of parameter space, and they are arranged in a typical snakes-and-ladders bifurcation structure, which is displayed in Fig. 6.

FIG. 6.

(a) Snakes-and-ladders bifurcation scenario. (b) Representative solutions for branches of symmetric (1,3) and asymmetric solutions (2). Parameters: J=39, Δ=1, and b=0.4.

FIG. 6.

(a) Snakes-and-ladders bifurcation scenario. (b) Representative solutions for branches of symmetric (1,3) and asymmetric solutions (2). Parameters: J=39, Δ=1, and b=0.4.

Close modal

It is known that snaking bifurcation scenarios can be triggered by heterogeneities in the underlying evolution equations. Examples discussed in the literature include the Swift–Hohenberg,51 Amari,40 and Ginzburg–Landau52 equations. In neural field models, heterogeneities are naturally introduced via harmonic perturbations of a homogeneous (distance-dependent) kernel, which break the translational invariance of the problem.53–56 In Ref. 40, we have shown that the following kernel leads to snaking in the Amari model,

w(x,y)=12e|xy|(1+acos(ky)),
(24)

and we, therefore, investigate the effect of this kernel on the QIF neural field model.

In the absence of spatial forcing (a=0), a system with exponential connectivity does not yield stable localized solutions (see Fig. 7). In the presence of modulation, we find snaking branches that oscillate around the branch obtained for a=0 (see Fig. 8). Furthermore, for small values of a, the snaking width increases proportionally to the value of a (not shown). These findings indicate that the snaking phenomenon in the QIF neural field model is entirely determined by the kernel choice, as in the Amari case.

FIG. 7.

(a) Bifurcation diagram of spatially uniform solutions and of localized solutions generated with the exponential kernel (24) with a=0. The lack of lateral inhibition results in the entire branch of solutions being unstable. (b) Representative solutions. Parameters: J=15Δ and Δ=2.

FIG. 7.

(a) Bifurcation diagram of spatially uniform solutions and of localized solutions generated with the exponential kernel (24) with a=0. The lack of lateral inhibition results in the entire branch of solutions being unstable. (b) Representative solutions. Parameters: J=15Δ and Δ=2.

Close modal
FIG. 8.

Snaking induced by spatial periodic modulation of the exponential connectivity kernel (24). (a) Bifurcation diagram obtained for a=0.1 (orange and purple branches, ladders not shown) compared with the one obtained for a=0 (blue). (b) Representative solutions. Parameters: J=15Δ and Δ=2.

FIG. 8.

Snaking induced by spatial periodic modulation of the exponential connectivity kernel (24). (a) Bifurcation diagram obtained for a=0.1 (orange and purple branches, ladders not shown) compared with the one obtained for a=0 (blue). (b) Representative solutions. Parameters: J=15Δ and Δ=2.

Close modal

Various nonlinear models, including chemical, fluid-dynamical, and particle systems, support time-periodic, spatially localized states termed oscillons (see Ref. 57 and references therein). A comprehensive theory for the existence and bifurcation structure of such solutions is the subject of experimental, numerical, and analytical investigations. We study oscillons in the QIF neural field model in the two main settings where they are observed in other media: (i) a non-autonomous setting, whereby oscillons emerge as the medium is subject to a homogeneous, exogenous, time-periodic forcing and (ii) an autonomous setting, whereby oscillons emerge spontaneously as one of the model parameters is varied.

We setup the QIF neural field model subject to a time-dependent, homogeneous, sinusoidal forcing with frequency ω,

tr=Δπ+2rv,tv=v2+Jwrπ2r2+η+Asin(ωt),

and cast it in the following, equivalent autonomous model formulation to perform numerical bifurcation analysis:

tr=Δπ+2rv,tv=v2+Jwrπ2r2+η+Aξ,ξ˙=ξ+ωζ(ξ2+ζ2)ξ,ζ˙=ζωξ(ξ2+ζ2)ζ.
(25)

Note that the numerical framework proposed here is applicable also if the forcing is heterogeneous.

In this setting, we expect oscillons to emerge without bifurcation from a localized steady state of the QIF neural field model with A=0, upon imposing a small-amplitude forcing, A1. We, therefore, select the default kernel (5), set η=10, for which the model with A=0 supports one stable (wide) and one unstable (narrow) bump (see Fig. 3), and continue time-periodic solutions to (25) in A>0 for ω=4, close to the network’s resonant frequency.30 

One stable and one unstable branch of oscillons emerge from A=0, as shown in Fig. 9, and connect at a saddle-node bifurcation. The stable branch undergoes a sequence of period-doubling bifurcations leading to chaos, and examples of a period-doubled solution and a chaotic solution are shown in Fig. 9, demonstrating the correspondence between the QIF neural field model and the spiking network model.

FIG. 9.

Maximum values of R1(t) plotted against the amplitude of sinusoidal forcing. At A3.8, a period-doubling bifurcation occurs, which is the starting point of a period-doubling cascade leading to chaos at A4.6. A period-doubled solution (A=4) and a chaotic solution (A=4.6) are shown. Parameters: Δ=2, J=15Δ, η=10, and ω=4.

FIG. 9.

Maximum values of R1(t) plotted against the amplitude of sinusoidal forcing. At A3.8, a period-doubling bifurcation occurs, which is the starting point of a period-doubling cascade leading to chaos at A4.6. A period-doubled solution (A=4) and a chaotic solution (A=4.6) are shown. Parameters: Δ=2, J=15Δ, η=10, and ω=4.

Close modal

In a recent study, we have investigated the effect of periodic forcing on a population of excitatory spiking neurons,30 whose solutions correspond to the spatially uniform states of the present model. In that context, it was shown that a sufficiently large forcing amplitude is able to suppress homogeneous oscillations. Here, we report that the same statement holds true for forced oscillons: no localized time-periodic solution is found to the right of the saddle-node bifurcation in Fig. 9, where the attractor is a spatially homogeneous, time-periodic state, which can be found by continuing in A the low-activity uniform steady state (not shown).

In the second scenario, oscillons occur in autonomous systems. Direct numerical simulations of reaction diffusion systems display oscillons in the proximity of codimension-two Turing–Hopf bifurcations of the homogeneous steady state.58,59 Oscillons in these systems have typically been observed as large-amplitude structures; hence,they are conjectured to form via a subcritical Hopf bifurcation of a heterogeneous, spatially localized steady state. This conjecture, however, has not yet been confirmed by numerical bifurcation analysis which, in contrast to direct numerical simulations, allows to track both stable and unstable states.

Here,we employ the Hermite–Galerkin scheme to study the formation of oscillons in the QIF neural field model. As mentioned above, a necessary ingredient for oscillons is the presence of oscillatory bifurcations. These bifurcations are precluded in one-population networks of QIF neurons but, as we shall see, are possible in two-population models; therefore,we turn our attention to the following network of coupled excitatory and inhibitory populations:

r˙e=Δπ+2reve,v˙e=ve2+ηe+JewereJiτiwiriπ2re2,τi2r˙i=Δπ+2τirivi,τiv˙i=vi2+ηi+JewereJiτiwiriπ2τi2ri2.

The subscripts e and indicate whether a variable or parameter refers to the excitatory or inhibitory population, respectively: the two populations have, for simplicity, the same heterogeneity parameter Δ, but they have possibly different membrane time constants and average background currents. In single-population mean fields, excitation and inhibition are artificially lumped into a single excitatory-inhibitory kernel [see for instance (5), (23), and (24)], whereas in the new, more realistic model, the kernels are separate,

we(x)=e|x|,wi(x)=14e|x|/2.
(26)

The connectivity parameters are chosen to be Je=Ji=J to recover a similar setting used in the lumped model. In Fig. 10(a), we show the bifurcation diagram of localized solutions using ηe as the bifurcation parameter. The bifurcation structure is similar to the lumped model, with the exception that the range of parameters for which stable solutions exist is narrower. This computation confirms that stationary bumps are supported by the two-population network. In order to hunt for oscillons, we continue the solution for ηe=10 in the parameter τi: the bump becomes unstable at a subcritical Hopf bifurcation at τi1.14, restabilizes at a saddle-node bifurcation, and undergoes a sequence of saddle-node bifurcations leading to a torus bifurcation (i.e., generalized Hopf bifurcation). The branch eventually restabilizes at a further saddle-node, leading to a period-doubling cascade which initiates around τi1.26 and to chaos at τi>1.29 [Fig. 10(b)].

FIG. 10.

(a) Bifurcation diagram of bump solutions in the E-I network (blue) and spatially uniform solutions (black) for τi=1. (b) Bifurcation diagram of emerging limit cycles [showing maxima of R1(t)] using τi as the bifurcation parameter. H, Hopf bifurcation; S, saddle-node bifurcation; T, torus bifurcation; and P, period-doubling bifurcation. Parameters: ηe=ηi=10, Je=Ji=15Δ, and Δ=2.

FIG. 10.

(a) Bifurcation diagram of bump solutions in the E-I network (blue) and spatially uniform solutions (black) for τi=1. (b) Bifurcation diagram of emerging limit cycles [showing maxima of R1(t)] using τi as the bifurcation parameter. H, Hopf bifurcation; S, saddle-node bifurcation; T, torus bifurcation; and P, period-doubling bifurcation. Parameters: ηe=ηi=10, Je=Ji=15Δ, and Δ=2.

Close modal

In Fig. 10, we also show numerical examples of a stable period-doubled solution at τi=1.28 and a chaotic solution at τi=1.295. We do not observe oscillons beyond τi=1.3. Chaotic solutions can also be reproduced in the spiking network model (see Fig. 11).

FIG. 11.

Chaotic solution in a spiking network. (a) Rastergram of excitatory population. (b) Rastergram of inhibitory population. (c) Excitatory and inhibitory spike rates averaged on interval 1<x<1 and sliding window in t (width 103). (d) Phase portrait of spike rates in (c). Parameters: ηe=ηi=10, Je=Ji=15Δ, Δ=2, and τi=1.295. These parameters correspond to point (4) in Fig. 10.

FIG. 11.

Chaotic solution in a spiking network. (a) Rastergram of excitatory population. (b) Rastergram of inhibitory population. (c) Excitatory and inhibitory spike rates averaged on interval 1<x<1 and sliding window in t (width 103). (d) Phase portrait of spike rates in (c). Parameters: ηe=ηi=10, Je=Ji=15Δ, Δ=2, and τi=1.295. These parameters correspond to point (4) in Fig. 10.

Close modal

We introduced a framework to study localized solutions in a neural field model that was recently derived as an exact representation of the mean field dynamics of networks of spiking neurons. Although this model does not permit closed-form solutions such as the Amari model with Heaviside firing rates, we show that it is possible to give an analytical estimate for the range of model parameters for which stable localized solutions exist. The structure of the QIF neural field model permits the straightforward use of Galerkin methods, which unlike the Amari model has a linear nonlocal term.

We have demonstrated that stationary equations can be transformed into a formulation that is equivalent to the stationary Amari model, provided an effective firing rate function is defined. The significance of such a firing rate is chiefly mathematical: the neural field possesses a rate variable, which is combined with the voltage variable in the effective firing rate; however, this transformation allows to map out patterned steady states of the QIF neural field model using the same toolkit available for the Amari formulation. In both models, localized solutions emerge subcritically from a branch of homogeneous steady states, which then restabilize at a saddle-node bifurcation. In the Amari model, this behavior is parametrized by a firing threshold, whereas here we use the average excitability of the network to map out solutions. However, there is a correspondence between the excitability of the model used here and the firing threshold in the Amari model, in the sense that an increase in the firing threshold in the latter corresponds to a decrease in the excitability in the former. In addition, techniques developed for piecewise-linear firing rate functions in the Amari model60 could be adapted to work for steady states in the QIF neural field model, using the correspondence described above. Furthermore, all branches of stationary solutions computed in this paper, including the snaking branches, also occur in standard rate-based models. The crucial difference lies in the transient dynamics of the two models, which makes the model considered here dynamically richer and more realistic.

The development of a Galerkin method opened up the possibility to study oscillons using numerical bifurcation analysis. We focused here on sinusoidal forcing of bump solutions, which is a proxy of oscillations ubiquitous in neuronal systems. In previous work,30 the neural mass version of this model was studied in terms of its response to oscillatory forcing in various frequency bands, and the present paper makes this exploration feasible also in the spatially extended model. We leave this exploration to a future publication.

In coupled networks of excitatory and inhibitory populations, a small change in the inhibitory membrane time scale can have a significant effect on the existence and dynamics of bump solutions and can elicit oscillons. This was demonstrated for instantaneous synapses, and it remains to be seen how the dynamics change when synaptic delays are introduced to the model. Interestingly, oscillatory solutions which undergo torus bifurcations have been observed in spatially extended networks of excitatory and inhibitory neurons with conductance-based dynamics.61 Another natural extension would be to examine coupled multi-layer neural field models,62 which are known to give rise to localized bump solutions when neither layer does in isolation.63 

The Galerkin numerical methods derived in this paper can be applied directly to more general spatially extended models of QIF networks, such as the ones mentioned above. For instance, adding a synaptic variable can be accounted for with an additional Galerkin expansion and n scalar variables per additional evolution equation.

Single population, QIF neural mass models with chemical as well as electrical synapses have recently been developed,64 and it was found that oscillations originate at Hopf bifurcations. Spatially extended versions of this model would then have the possibility of forming oscillons with a single population, although it is not clear whether Hopf bifurcations of bumps will occur near Hopf bifurcations of homogeneous states, which are the ones mapped in Ref. 64.

Understanding how slow–fast temporal scales are generated by the discrete network is an open question, which has recently been addressed in networks of sparsely coupled networks of QIF neurons.65 Employing our numerical methodology to these macroscopic mean fields is also possible, and one could study how such slow–fast phenomena occur in more realistic, spatially extended networks.

H.S. acknowledges financial support from the Spanish Ministry of Economics and Competitiveness through the María de Maeztu Program for Units of Excellence in R&D (No. MDM-2014-0445) and Grant (No. MTM2015-71509-C2-1-R) and from the German Research Council [DFG (No. KN 588/7-1) within priority program “Computational Connectomics” (SPP 2041)].

1.
A.
Compte
,
N.
Brunel
,
P. S.
Goldman-Rakic
, and
X.-J.
Wang
, “
Synaptic mechanisms and network dynamics underlying spatial working memory in a cortical network model
,”
Cereb. Cortex
10
,
910
923
(
2000
).
2.
K.
Wimmer
,
D. Q.
Nykamp
,
C.
Constantinidis
, and
A.
Compte
, “
Bump attractor dynamics in prefrontal cortex explains behavioral precision in spatial working memory
,”
Nat. Neurosci.
17
,
431
439
(
2014
).
3.
S.
Kim
,
H.
Rouault
,
S.
Druckmann
, and
V.
Jayaraman
, “
Ring attractor dynamics in the drosophila brain
,”
Science
356
,
849
853
(
2017
).
4.
H. R.
Wilson
and
J. D.
Cowan
, “
A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue
,”
Kybernetik
13
,
55
80
(
1973
).
5.
P. L.
Nunez
, “
The brain wave equation: A model for the EEG
,”
Math. Biosci.
21
,
279
297
(
1974
).
6.
S.
Amari
, “
Dynamics of pattern formation in lateral-inhibition type neural fields
,”
Biol. Cybern.
27
,
77
87
(
1977
).
7.
S.
Coombes
, “
Waves, bumps, and patterns in neural field theories
,”
Biol. Cybern.
93
,
91
108
(
2005
).
8.
P. C.
Bressloff
, “
Spatiotemporal dynamics of continuum neural fields
,”
J. Phys. A
45
,
033001
(
2012
).
9.
G. B.
Ermentrout
, “
Reduction of conductance-based models with slow synapses to neural nets
,”
Neural Comput.
6
,
679
695
(
1994
).
10.
S.
Ostojic
and
N.
Brunel
, “
From spiking neuron models to linear-nonlinear models
,”
PLoS Comput. Biol.
7
,
e1001056
(
2011
).
11.
E. S.
Schaffer
,
S.
Ostojic
, and
L. F.
Abbott
, “
A complex-valued firing-rate model that approximates the dynamics of spiking networks
,”
PLoS Comput. Biol.
9
,
e1003301
(
2013
).
12.
M. A.
Buice
and
C. C.
Chow
, “
Dynamic finite size effects in spiking neural networks
,”
PLoS Comput. Biol.
9
,
e1002872
(
2013
).
13.
S.
Visser
and
S. A. V.
Gils
, “
Lumping Izhikevich neurons
,”
EPJ Nonlinear Biomed.
2
,
6
(
2014
).
14.
M.
Mattia
, “Low-dimensional firing rate dynamics of spiking neuron networks,” arXiv:160908855 (2016).
15.
T.
Schwalger
,
M.
Deger
, and
W.
Gerstner
, “
Towards a theory of cortical columns: From spiking neurons to interacting neural populations of finite size
,”
PLoS Comput. Biol.
13
,
e1005507
(
2017
).
16.
M.
Augustin
,
J.
Ladenbauer
,
F.
Baumann
, and
K.
Obermayer
, “
Low-dimensional spike rate models derived from networks of adaptive integrate-and-fire neurons: Comparison and implementation
,”
PLoS Comput. Biol.
13
,
e1005545
(
2017
).
17.
Y.
Park
and
G. B.
Ermentrout
, “
A multiple timescales approach to bridging spiking- and population-level dynamics
,”
Chaos
28
,
083123
(
2018
).
18.
S.-W.
Qiu
and
C. C.
Chow
, “
Finite-size effects for spiking neural networks with spatially dependent coupling
,”
Phys. Rev. E
98
,
062414
(
2018
).
19.
T. B.
Luke
,
E.
Barreto
, and
P.
So
, “
Complete classification of the macroscopic behavior of a heterogeneous network of theta neurons
,”
Neural Comput.
25
,
3207
3234
(
2013
).
20.
E.
Montbrió
,
D.
Pazó
, and
A.
Roxin
, “
Macroscopic description for networks of spiking neurons
,”
Phys. Rev. X
5
,
021028
(
2015
).
21.
C. R.
Laing
, “
Derivation of a neural field model from a network of theta neurons
,”
Phys. Rev. E
90
,
010901
(
2014
).
22.
C. R.
Laing
, “
Exact neural fields incorporating gap junctions
,”
SIAM J. Appl. Dyn. Syst.
14
,
1899
1929
(
2015
).
23.
J. M.
Esnaola-Acebes
,
A.
Roxin
,
D.
Avitabile
, and
E.
Montbrió
, “
Synchrony-induced modes of oscillation of a neural field model
,”
Phys. Rev. E
96
,
052407
(
2017
).
24.
A.
Byrne
,
D.
Avitabile
, and
S.
Coombes
, “
Next-generation neural field model: The evolution of synchrony within patterns and waves
,”
Phys. Rev. E
99
,
012313
(
2019
).
25.
S.
Coombes
and
À.
Byrne
, “Next generation neural mass models,” in Nonlinear Dynamics in Computational Neuroscience (
Springer
,
2019
).
26.
F.
Devalle
,
A.
Roxin
, and
E.
Montbrió
, “
Firing rate equations require a spike synchrony mechanism to correctly describe fast oscillations in inhibitory networks
,”
PLoS Comput. Biol.
13
,
e1005881
(
2017
).
27.
I.
Ratas
and
K.
Pyragas
, “
Macroscopic oscillations of a quadratic integrate-and-fire neuron network with global distributed-delay coupling
,”
Phys. Rev. E
98
,
052224
(
2018
).
28.
R.
Zillmer
,
R.
Livi
,
A.
Politi
, and
A.
Torcini
, “
Stability of the splay state in pulse-coupled networks
,”
Phys. Rev. E
76
,
046102
(
2007
).
29.
E.
Ott
and
T. M.
Antonsen
, “
Low dimensional behavior of large systems of globally coupled oscillators
,”
Chaos
18
,
037113
(
2008
).
30.
H.
Schmidt
,
D.
Avitabile
,
E.
Montbrió
, and
A.
Roxin
, “
Network mechanisms underlying the role of oscillations in cognitive tasks
,”
PLoS Comput. Biol.
14
,
e1006430
(
2018
).
31.
Y.
Kawamura
, “
From the Kuramoto-Sakaguchi model to the Kuramoto-Sivashinsky equation
,”
Phys. Rev. E
89
,
010901
(
2014
).
32.
A. R.
Champneys
, “
Homoclinic orbits in reversible systems and their applications in mechanics, fluids and optics
,”
Physica D
112
,
158
186
(
1998
).
33.
J.
Burke
and
E.
Knobloch
, “
Homoclinic snaking: Structure and stability
,”
Chaos
17
,
037102
(
2007
).
34.
J.
Burke
and
E.
Knobloch
, “
Snakes and ladders: Localized states in the Swift–Hohenberg equation
,”
Phys. Lett. A
360
,
681
688
(
2007
).
35.
J.
Knobloch
and
T.
Wagenknecht
, “
Homoclinic snaking near a heteroclinic cycle in reversible systems
,”
Physica D
206
,
82
93
(
2005
).
36.
M.
Beck
,
J.
Knobloch
,
D. J. B.
Lloyd
,
B.
Sandstede
, and
T.
Wagenknecht
, “
Snakes, ladders, and isolas of localised patterns
,”
SIAM J. Math. Anal.
41
,
936
972
(
2009
).
37.
C. R.
Laing
,
W. C.
Troy
,
B.
Gutkin
, and
G. B.
Ermentrout
, “
Multiple bumps in a neuronal model of working memory
,”
SIAM J. Appl. Math.
63
,
62
97
(
2002
).
38.
S.
Coombes
,
G. J.
Lord
, and
M. R.
Owen
, “
Waves and bumps in neuronal networks with axo-dendritic synaptic interactions
,”
Physica D
178
,
219
241
(
2003
).
39.
A. J.
Elvin
,
C. R.
Laing
,
R. I.
McLachlan
, and
M. G.
Roberts
, “
Exploiting the Hamiltonian structure of a neural field model
,”
Physica D
239
,
537
546
(
2010
).
40.
D.
Avitabile
and
H.
Schmidt
, “
Snakes and ladders in an inhomogeneous neural field model
,”
Physica D
294
,
24
36
(
2015
).
41.
J.
Rankin
,
D.
Avitabile
,
J.
Baladron
,
G.
Faye
, and
D. J. B.
Lloyd
, “
Continuation of localised coherent structures in nonlocal neural field equations
,”
SIAM J. Sci. Comput.
36
,
B70
B93
(
2013
).
42.
C. R.
Laing
and
S.
Coombes
, “
The importance of different timings of excitatory and inhibitory pathways in neural field models
,”
Network
17
,
151
172
(
2006
).
43.
H. G.
Meijer
and
S.
Coombes
, “
Travelling waves in models of neural tissue: From localised structures to periodic waves
,”
EPJ Nonlinear Biomed. Phys.
2
,
3
(
2014
).
44.
A.
Champneys
, “
Homoclinic orbits in reversible systems and their applications in mechanics, fluids and optics
,”
Physica D
112
,
158
186
(
1998
), .
45.
D. J. B.
Lloyd
,
B.
Sandstede
,
D.
Avitabile
, and
A. R.
Champneys
, “
Localized hexagon patterns of the planar Swift–Hohenberg equation
,”
SIAM J. Appl. Dyn. Syst.
7
,
1049
1100
(
2008
).
46.
D.
Avitabile
,
D. J. B.
Lloyd
,
J.
Burke
,
E.
Knobloch
, and
B.
Sandstede
, “
To snake or not to snake in the planar Swift-Hohenberg equation
,”
SIAM J. Appl. Dyn. Syst.
9
,
704
733
(
2010
).
47.
D.
Morgan
and
J. H. P.
Dawes
, “
The Swift–Hohenberg equation with a nonlocal nonlinearity
,”
Physica D
270
,
60
80
(
2014
).
48.
C. R.
Laing
and
W. C.
Troy
, “
PDE methods for nonlocal models
,”
SIAM J. Appl. Dyn. Syst.
2
,
487
516
(
2003
).
49.
G.
Faye
,
J.
Rankin
, and
P.
Chossat
, “
Localized states in an unbounded neural field equation with smooth firing rate function: A multi-parameter analysis
,”
J. Math. Biol.
66
,
1303
1338
(
2013
).
50.
G.
Faye
,
J.
Rankin
, and
D. J. B.
Lloyd
, “
Localized radial bumps of a neural field equation on the Euclidean plane and the Poincaré disk
,”
Nonlinearity
26
,
437
(
2013
).
51.
H.-C.
Kao
,
C.
Beaume
, and
E.
Knobloch
, “
Spatial localization in heterogeneous systems
,”
Phys. Rev. E
89
,
012903
(
2014
).
52.
B. C.
Ponedel
and
E.
Knobloch
, “
Forced snaking: Localized structures in the real Ginzburg-Landau equation with spatially periodic parametric forcing
,”
Eur. Phys. J. Spec. Top.
225
,
2549
2561
(
2016
).
53.
P. C.
Bressloff
, “
Traveling fronts and wave propagation failure in an inhomogeneous neural network
,”
Physica D
155
,
83
100
(
2001
).
54.
Z. P.
Kilpatrick
,
S. E.
Folias
, and
P. C.
Bressloff
, “
Traveling pulses and wave propagation failure in inhomogeneous neural media
,”
SIAM J. Appl. Dyn. Syst.
7
,
161
185
(
2008
).
55.
H.
Schmidt
,
A.
Hutt
, and
L.
Schimansky-Geier
, “
Wave fronts in inhomogeneous neural field models
,”
Physica D
238
,
1101
1112
(
2009
).
56.
C.
Coombes
and
C. R.
Laing
, “
Pulsating fronts in periodically modulated neural field models
,”
Phys. Rev. E
83
,
011912
(
2011
).
57.
E.
Knobloch
, “
Spatially localized structures in dissipative systems: Open problems
,”
Nonlinearity
21
,
T45
T60
(
2008
).
58.
V.
Vanag
and
I.
Epstein
, “
Stationary and oscillatory localized patterns, and subcritical bifurcations
,”
Phys. Rev. Lett.
92
,
128301
(
2004
).
59.
V. K.
Vanag
and
I. R.
Epstein
, “
Localized patterns in reaction-diffusion systems
,”
Chaos
17
,
037110
(
2007
).
60.
S.
Coombes
and
H.
Schmidt
, “
Neural fields with sigmoidal firing rates: Approximate solutions
,”
Discrete Contin. Dyn. Syst. A
28
,
1369
1379
(
2010
).
61.
S. E.
Folias
and
G. B.
Ermentrout
, “
Spatially localized synchronous oscillations in synaptically coupled neuronal networks: Conductance-based models and discrete maps
,”
SIAM J. Appl. Dyn. Syst.
9
,
1019
1060
(
2010
).
62.
P. C.
Bressloff
and
S. R.
Carroll
, “
Laminar neural field model of laterally propagating waves of orientation selectivity
,”
PLoS Comput. Biol.
11
,
e1004545
(
2015
).
63.
S. E.
Folias
and
G. B.
Ermentrout
, “
New patterns of activity in a pair of interacting excitatory-inhibitory neural fields
,”
Phys. Rev. Lett.
107
,
228103
(
2011
).
64.
B.
Pietras
,
F.
Devalle
,
A.
Roxin
,
A.
Daffertshofer
, and
E.
Montbrió
, “
Exact firing rate model reveals the differential effects of chemical versus electrical synapses in spiking networks
,”
Phys. Rev. E
100
,
042412
(
2019
).
65.
H.
Bi
,
M.
Segneri
,
M. d.
Volo
, and
A.
Torcini
, “Coexistence of fast and slow gamma oscillations in one population of inhibitory spiking neurons,” arXiv:1907.00230 (2019).