Snakes and ladders in an inhomogeneous neural field model

Continuous neural field models with inhomogeneous synaptic connectivities are known to support traveling fronts as well as stable bumps of localized activity. We analyze stationary localized structures in a neural field model with periodic modulation of the synaptic connectivity kernel and find that they are arranged in a snakes-and-ladders bifurcation structure. In the case of Heaviside firing rates, we construct analytically symmetric and asymmetric states and hence derive closed-form expressions for the corresponding bifurcation diagrams. We show that the ideas proposed by Beck and co-workers to analyze snaking solutions to the Swift-Hohenberg equation remain valid for the neural field model, even though the corresponding spatial-dynamical formulation is non-autonomous. We investigate how the modulation amplitude affects the bifurcation structure and compare numerical calculations for steep sigmoidal firing rates with analytic predictions valid in the Heaviside limit.

1. Introduction. Continuous neural field models are a common tool to investigate large-scale activity of neuronal ensembles. Since the seminal work of Wilson and Cowan [57,58] and Amari [1,2], these nonlocal models have helped understanding the emergence of spatial and spatio-temporal coherent structures in various experimental observations. Stationary spatially-extended patterns have been found in visual hallucinations [23,9], while stationary localized structures, commonly referred to as bumps [14], are related to short term (working) memory [28] and feature selectivity in the visual cortex [6,30]. Traveling waves of neural activity are relevant for information processing [24] and can be evoked in vitro in slice preparations of cortical [60], thalamic [36] or hippocampal [47] tissue by electric stimulation (for recent reviews see [49,54]). Furthermore, traveling waves have also been observed in vivo in the form of spreading depression in neurological disorders such as migraine [41].
The simplest neural field models are (systems of) integro-differential equations posed on the real line or on the plane. The corresponding nonlocal terms feature a synaptic kernel, a function that models the neural connectivity at a macroscopic scale. For mathematical convenience, neural field models are often chosen to be translationally invariant, that is, the synaptic kernel depends on the Euclidean distance between points on the domain. This assumption reflects well experiments in which cortical slices are pharmacologically prepared. However, in vivo experiments by Hubel and Wiesel [31,32,33,34] revealed that a complex microstructure is present in several areas of mammalian cortex. In order to model this microstructure, Bressloff [7] incorporated a spatially-periodic modulation of the synaptic kernel into a one-dimensional neural field model. The translational invariance is thus broken, leading to slower traveling waves and, for sufficiently large modulation amplitudes, to propagation failure (a similar effect is also caused by inhomogeneities in the input [10]). In the present article we show how inhomogeneities in the synaptic connectivity can give rise to a multitude of stable stationary bumps which are organized in parameter space via a characteristic snaking bifurcation structure.
The formation and bifurcation structure of stationary localized patterns has been studied extensively in partial differential equations (PDEs) posed on domains with one [59,21,11,12,13,22,5], two [52,53,43,42,3,46] and three spatial dimensions [56,4,44]. Most analytical studies focus on the Swift-Hohenberg equation (or one of its variants) posed on the real line: stationary localized solutions to the PDE connect a homogeneous (background) state to a patterned state at the core of the domain, hence they can be interpreted as homoclinic connections in the corresponding spatialdynamical system. In a suitable region of parameter space, close to the so-called Maxwell point, there exist infinitely many homoclinic connections, corresponding to PDE solutions with varying spatial extent. Localized states with different symmetries belong to intertwined solution branches that snake between two (or more) limits and are connected by branches of asymmetric states. This bifurcation structure was called snakes and ladders by Burke and Knobloch [12]. It is known that snaking localized structures arise also in systems with nonlocal terms. For instance, snaking bumps are supported by neural field models posed on the real line [38,40,17,26,25] and on the plane [51], as well as the Swift-Hohenberg equation with nonlocal terms [48]. In neural field models, the choice of the kernel has an impact on the bifurcation structure [51], hence it is interesting to study how inhomogeneities affect the existence and stability of localized states, an investigation that has been carried out very recently by Kao et al. in the context of the Swift-Hohenberg equation [35].
In the present paper, we study a neural field model with a synaptic kernel featuring a tunable harmonic inhomogeneity [55,16]. As pointed out by Schmidt et al. [55], the inhomogeneity gives rise to stable bumps which do not exist in the homogeneous case. We will show here that the synaptic modulation is also responsible for the snaking behavior of such solutions.
A characteristic of neural field models is that they can be conveniently analyzed in the limit of Heaviside firing rates: for the model under consideration, bumps can be constructed analytically, hence, following Beck et al. [5], we can derive closed form expressions for the snaking bifurcation curves. In addition, we show that the Heaviside limit provides a good approximation to the case of steep sigmoidal firing rates, for which the theory by Beck et al. can not be directly applied.
This article is structured as follows: in Section 2 we present the neural field model and discuss stability of stationary solutions. In Section 3 we show numerical simulations of the model in the case of steep sigmoidal firing rates, for which an equivalent PDE formulation is available. In Section 4 we move to the Heaviside firing rate limit, for which we discuss the construction and stability of generic localized steady states. In Sections 5 and 6 we calculate explicitly periodic and localized steady states and infer the relative bifurcation diagrams. In Section 7 we discuss how the bifurcation structure is affected by changes in the modulation amplitude. We conclude the paper in Section 8.
2. The integral model. We consider a neural field model of the Amari type, posed on the real line, where u is the synaptic potential, W the synaptic connectivity and f a nonlinear function for the conversion of the synaptic potential into a firing rate. In general, both W and f depend upon a set of control parameters, which have been omitted here for simplicity.
Several studies of neural field models assume translation invariance in the model (see [8,15] and references therein), therefore the synaptic strength W depends solely on the Euclidean distance between x and y, that is W (x, y) = w(|x − y|). A neural field of this type is said to be homogeneous.
A simple way to incorporate an inhomogeneous microstructure is to multiply the homogeneous kernel w by a periodic function A(y) that modulates the synaptic connectivity and thus breaks translational invariance. Following Bressloff [7] we choose A(y) to be a simple harmonic function W (x, y; a, ε) = w(|x − y|)A(y; a, ε) Here, a is the amplitude of the modulation and 2πε its wavelength. With this choice, the neural field model is invariant with respect to transformations In this paper we study stationary localized states of system (2.1) with inhomogeneous kernels (2.2). The firing rate f will be either a Heaviside function f (u; h) = H(u − h), where h is the firing threshold, or a sigmoidal firing rate with ν 1. In the limit ν → ∞, the sigmoidal firing rate (2.4) recovers the Heaviside case. As we shall see, a Heaviside firing rate will be more convenient for analytical calculations, whereas a steep smooth firing rate will be employed for numerical computations.
Stationary solutions to the system (2.1)-(2.2) satisfy Linear stability is studied posing u(x, t) = q(x) + e λt v(x), with v 1, and linearizing the right-hand side of (2.1). This leads to the following nonlocal eigenvalue problem where we have formally denoted by f the derivative of f with respect to u. This linear stability analysis is standard in the study of localized solutions in neural field models [14].
3. PDE formulation for smooth firing rates. When f is a smooth sigmoid, it is advantageous to reformulate the nonlocal problem (2.1) as a local PDE, more suitable for direct numerical simulation and numerical continuation. Following [40,20,16,18,39], we take the Fourier transform of (2.1), with kernel expressed by (2.2) whereŵ(ξ) = (2π) −1 /(ξ 2 +1). Multiplying the previous equation by ξ 2 +1 and taking the inverse Fourier transform we obtain where the dependence of u on x and t has been omitted for simplicity. Once complemented with suitable initial and boundary conditions, the equation above constitutes an equivalent PDE formulation of the model problem. Steady states are solutions to and linear stability is inferred via the generalized eigenvalue problem In passing we note that time simulations of (3.1) and stability calculations (3.3) can be carried out numerically without forming a discretization for (1 − ∂ 2 x ) −1 (see [18]). In   [16], whereas in the present paper we focus on the existence and bifurcation structure of localized states.
The simulations in Figure 3.1 are compatible with a snakes-and-ladders bifurcation structure and, owing to the spatial modulation, we expect to find stable localized states that are spatially in-phase with A and centered around its local minima and maxima. In Figure 3.2 we fix h and perturb a localized steady state with abrupt phase slips in the kernel modulation. More precisely we set where t i = 10 + iκ and χ [ti,ti+1] is the indicator function with support [t i , t i+1 ]. After four phase slips, we return to the original spatial inhomogeneity A(x) = 1+a cos(x/ε), which is kept constant thereafter. Perturbations with κ = 10 elicit a localized steady state that is symmetric with respect to the axis x = −2π, whereas shorter phase slips, with κ = 5, give rise to states that are symmetric with respect to the axis x = −π.
In local models supporting localized states, symmetries of the PDE are reflected in the bifurcation structure: each snaking branch includes solutions with the same symmetry and intertwined branches are connected by ladders of asymmetric solutions. In one-dimensional snaking systems with spatial reversibility, localized states can be interpreted from a spatial-dynamical systems viewpoint and symmetries of the PDE correspond to reversers of the spatial-dynamical system [5,45]. Following this approach, we recast (3.2) as a first-order non-autonomous Hamiltonian system in x d dx where we posed (U 1 , U 2 ) = (q, q x ). Localized steady states of the nonlocal model correspond to bounded solutions to (3.5) that decay exponentially as |x| → ∞. System (3.5) is reversible: for each n ∈ Z, we consider the following autonomous extension We say that a stationary state q is even-(odd-)symmetric if Rq = q and n is even (odd), that is, q(x) is symmetric with respect to the axis x = nπε and n is even (odd). Conversely, we say that a solution is asymmetric if Rq = q. The stationary profiles plotted in Figures 3.2(a) and 3.2(b) correspond to an even-and odd-symmetric solution, respectively. The spatial-dynamic formulation developed in [5,45] for the Swift-Hohenberg equation allows to predict snaking branches of localized patterns from the bifurcation structure of fronts connecting the trivial (background) state to the core state. We can not directly apply this theory to our case, in that system (3.5) is non-autonomous and (0, 0) is not an equilibrium. However, we shall see that the ideas presented in those papers remain valid: in the limit of Heaviside firing rate, which gives rise to a nonsmooth spatial-dynamic formulation, we are able to compute explicit expressions for connecting orbits and, hence, for the snaking bifurcation diagram, which we partially present in Figure 3.3; for sigmoidal firing rates we will adopt numerical continuation and compute snaking bifurcation branches solving the boundary-value problem (3.2) and the associated stability problem (3.3).
4. Steady states for Heaviside firing rate. In the case of Heaviside firing rate, localized steady states with two threshold crossings can be constructed explicitly for the inhomogeneous model and their stability can be inferred solving a simple 2by-2 eigenvalue problem. To each steady state q with firing threshold h, we associate an active region B = { x ∈ R | q(x) > h }, that is, a subset of the real line in which q is above threshold. In the case of Heaviside firing rate, this implies that H(q(x)) ≡ 1 if x ∈ B and 0 otherwise, hence Equation (2.5) can be rewritten as where x 1 and x 2 are the lower and upper boundary of B, which is assumed to be compact. If the threshold crossings x 1,2 are known, then (4.1) yields the profile of the stationary solution. The boundaries x 1 and x 2 can be determined as functions of system parameters by enforcing the threshold crossing conditions q( If f is the Heaviside function, the nonlocal eigenvalue problem (2.6) is written as where Lx = πε for periodic solutions and Lx = 30 for localized solutions. The trivial steady state u(x) ≡ 0 coexists with the fully periodic state for h ∈ 0, 1 − aε 2 /(1 + ε 2 ) . A snaking branch of even-symmetric localized solutions emanates from B 0 . As we ascend the snaking diagram, more bumps are formed. Example solutions are plotted in the panels. For reference, we also plot the activity threshold u(x) ≡ h (dashed magenta). There exist (not shown) a snaking branch of localized odd-symmetric solutions, as well as ladder branches connecting the two snaking branches.
where δ denotes the Dirac delta function and we have formally made use of the identity Evaluating equation (4.2) at x = x 1,2 we find that the stability of a steady state of the inhomogeneous integral problem with Heaviside firing rate is governed by a 2-by-2 eigenvalue problem whose eigenpairs {(λ k , ξ k )} k=1,2 can be found explicitly. Here, the eigenvector ξ k has In particular, we find For each eigenpair (λ k , ξ k ), the corresponding eigenfunction v k of (4.2) is then given by In the following sections we will apply this framework to both periodic and localized solutions in the Heaviside limit.
Remark 4.1 (Number of threshold crossings). The framework presented here can be extended to patterns with an arbitrary number of threshold crossings; however, throughout this paper we will restrict analytic calculations to solutions that have only two threshold crossings, or to spatially-periodic patterns with two threshold crossings per period. The linear stability analysis outlined here is valid for small perturbations v that have the same number of threshold crossings of q.
Remark 4.2 (Stability of solutions with no threshold crossing). Solutions that do not cross threshold are linearly stable, in that the eigenvalue problem (4.2) gives a single eigenvalue λ = −1.
5. Homogeneous and spatially periodic solutions for Heaviside firing rates. We now begin exploring steady state solutions to the integral model (2.1) with inhomogeneous kernel (2.2) and Heaviside firing rate f (u; h) = H(u − h). If the kernel is homogeneous, a straightforward computation shows that localized solutions exist and are linearly unstable. These patterns are organized in parameter space with a non-snaking bifurcation diagram: we integrate (4.1) with a = 0, x 1,2 = ±L/2 and obtain where h = (1 − e −L )/2. Using (4.5) we find λ 1,2 ≥ 0. We plot these solutions and their bifurcation diagram in Figure 5.1. From now on, we will concentrate on the more interesting case a > 0.
Owing to the inhomogeneity, the only spatially-homogeneous solution is the trivial from which we deduce 0 = κ < h. The trivial solution is linearly stable for strictly positive h (see Remark 4.2).
Spatially-periodic states are also supported by the integral model. In A we show that 2πε-periodic solutions satisfy  In other words, if we seek for a stationary 2πε-periodic solution, then we may pass from an integral equation posed on R to a reduced integral formulation posed on the interval [−πε, πε], provided that we use the amended kernel w instead of w. In passing, we note that similar conditions for periodic solutions can be derived for generic exponential kernels. We now specialize the problem (5.1)-(5.3) to the case of Heaviside firing rate f (u; h) = H(u − h), construct 2πε-periodic stationary solutions and explore their bifurcation structure. The simplest type of stationary periodic state of the model is the above-threshold solution q at , that is, a solution that lies above threshold h for all x ∈ R. We then formulate the following problem: Problem 5.1 (Above-threshold periodic solutions). For fixed h, a, ε ∈ R + , find a smooth 2πε-periodic function q at such that An explicit solution q at can be computed in closed form for the specific kernel (2.2), yielding q at (x) = 1 + aε 2 1 + ε 2 cos for h ∈ 0, 1 − aε 2 /(1 + ε 2 ) . Since there are no threshold crossings, q at is stable in this interval of h for all values of a and ε. In Figure 3.3, we show an example of q at for a = 0.3, ε = 1 (solution label 5).
We now turn to the more interesting case of periodic solutions that cross threshold. The simplest of such cross-threshold states, q ct , are solutions that attain the value h exactly twice in [−πε, πε), as shown in Figure 5.2(a). More precisely, we derive cross-threshold solutions as follows: Problem 5.2 (Cross-threshold periodic solutions). For fixed h, a, ε ∈ R + , find an even 2πε-periodic smooth function q ct and a number L ∈ (0, 2πε) such that q ct (L/2) = h, In analogy with [5], we call the equation above a bifurcation equation for periodic solutions q ct . Explicit formulae for the solution profile q ct and the corresponding bifurcation equation are given in B.
The stability of a stationary profile (q ct , L) is found in a similar fashion to what was done for stationary states in Section 4, with the original kernel w replaced by the amended kernelw. We find where x 1,2 = ∓L/2. Evaluating the equation above at x = x 1,2 yields the pair of eigenvalues where we have made use of the fact that |q ct (x)| and A(x) are even. We are now ready to study the bifurcation structure of periodic solutions in greater detail.
In the Heaviside limit we use Equations (5.5)-(5.8) which allow us to compute the solution profile, its activity region L and its stability as a function of h. The resulting bifurcation diagrams are shown in Figure 5.2(b). The main continuation parameter is h and we set ε = 1, a ∈ {0.3, 0.7, 1}: for small values of a the trivial state q 0 coexists with the above threshold solution q at for 0 < h < 0, 1 − aε 2 /(1 + ε 2 ) . At the grazing point h = 1 − aε 2 /(1 + ε 2 ), the above threshold solution becomes tangent to u(x) ≡ h.
The branches of q 0 and q at are connected by a branch of cross-threshold solutions which are initially unstable. As we increase a, two saddle node bifurcations emerge on the cross-threshold branch, at a cusp, and there exists an interval of h in which q 0 , q ct and q at coexist and are stable. As a is further increased, only one saddle node persists and we have an extended bistability region. We refer the reader to Section 7 for a more detailed study of the two-parameter bifurcation diagram.
We can also study the case of continuous sigmoidal firing rates (2.4) using standard numerical bifurcation analysis techniques: we find steady states q solving (3.2) with Neumann boundary conditions and we continue the solution in parameter space with pseudo-arclength continuation [29] using the secant code developed in [51]. A comparison between bifurcation diagrams for Heaviside and sigmoidal firing rates is presented in Figure 5.3. The solution branches are in good agreement, with the exception of the fold points, as it can be seen in the insets.

Construction and bifurcation structure of localized solutions for
Heaviside firing rates. Localized steady states are solutions to (2.1) which decay to zero as |x| → ∞ and for which the activity region B is a finite disjoint union of bounded connected intervals [2,25]. In of the PDE model (3.1) posed on a large finite domain with Neumann boundary conditions and steep sigmoidal firing rate with ν = 50. The parameters are chosen such that the trivial solution q 0 and the above-threshold periodic solution q at are supported in the Heaviside firing rate case. As expected, stable localized patterns are found in this region.
In this section, we construct such patterns analytically and study their stability. As it was done in Section 5, we will perform analytical or semi-analytical calculations in the Heaviside limit, whereas we will employ numerical continuation for sigmoidal firing rates.
As seen in Section 4, a generic bump q b with active region B = (x 1 , x 2 ) ⊂ R satisfies, in the Heaviside limit, w(|x − y|)A(y; a, ε) dy. (6.1) Without loss of generality, we pose x 1,2 = x 0 ∓ L/2. We note that if L = 0 then q b coincides with the trivial solution. In analogy with the periodic case, we find a localized solution as follows: q b and scalars x 0 ∈ R, L ∈ R + , such that Remark 6.2. In the problem above we do not enforce explicitly asymptotic conditions for q b , since they are implied by (6.4) for our particular choice of w and A. Indeed, let A * (a, ε) = max x∈R |A(x; a, ε)|, then In Section 3 we discussed symmetric and asymmetric solutions in the context of spatial-dynamical systems of the PDE associated with the integral model. Equivalently, a solution is symmetric if q b (x − x 0 ) = q b (x 0 − x) and asymmetric otherwise. Problem 6.1 does not provide a direct way to distinguish between symmetric and asymmetric states, but it can be reformulated so as to avoid this limitation. Each solution (q b , x 0 , L) to Problem 6.1 is such that q b (x 0 − L/2) = q b (x 0 + L/2), which can be written as Crucially, (6.5) holds if either Ψ sym = 0 or Ψ asym = 0, so we are now ready to construct symmetric and asymmetric localized solutions as follows: 13 Problem 6.3 (Symmetric and asymmetric localized solutions). For fixed h, a, ∈ R + , find a smooth nonnegative function q b and scalars x 0 ∈ R, L ∈ R + , such that Ψ sym (x 0 ; ε) = 0, (or Ψ asym (L; ε) = 0) (6.8) In symmetric states, the symmetry condition (6.6) fixes the value of x 0 ; more precisely we have x 0 = nπε for n ∈ Z, therefore we distinguish between even-and oddsymmetric solutions, depending on the value of n. On the other hand, in asymmetric states the width L is fixed by the asymmetry condition (6.7) and x 0 is not restricted to assume discrete values.
For our choice of the connectivity function w and modulation A we derive closedform expressions for symmetric and asymmetric localized states.
For the profile of symmetric solutions we find where the auxiliary functions Θ 1 and Θ 2 are given by In the above expressions we posed Φ = arctan −1 and we exploited the fact that sin(x 0 / ) = 0. Similarly, for asymmetric solutions we obtain with auxiliary functions Λ 1 and Λ 2 given by Examples trivial solution q 0 and the periodic above-threshold q at solution coexist. As expected, localized solutions are in-phase with the inhomogeneity A. where Θ 2 and Λ 2 are auxiliary functions defined above. In the bifurcation function I sym the value of x 0 is fixed by the condition Ψ sym (x 0 ) = 0, hence cos(x 0 / ) = ±1. Similarly, L is fixed in the expression of I asym and its value is determined by Ψ asym (L) = 0. Following [5], we notice that the bifurcation equation (6.12) is a parametrization of snaking branches of even-and odd-symmetric solutions, whereas equation (6.13) is a parametrization of ladder branches of asymmetric solutions: indeed both x 0 and L depend on h, as they solve Problem 6.3. In this case, however, the bifurcation equations are available in closed form so we can proceed directly to plot snakes and ladders. In Figure 6.2, we fix a and ε, construct localized solutions and plot their bifurcation diagrams as loci of points on the (L, h)-plane that satisfy the bifurcation equations. In particular we use I sym (L; 0, a, ε) and I sym (L; πε, a, ε) to plot representative branches of even-and odd-symmetric solutions, respectively. As expected, in the limit for large L, I sym is well approximated by a cosinusoidal function. On the other hand, ladders are found using I asym (x 0 ; L, a, ε), where L satisfies the asymmetry condition Ψ asym (L; ε) = 0.
The stability problem of a localized state (q b , x 0 , L) is determined following the scheme outlined in Section 2: we use Equations  Figure 6.2 (see Equation (6.14)). Right: eigenfunctions of selected states, in the proximity of a saddle node and a pitchfork on the snaking branch.
x 1,2 = x 0 ∓ L/2. For symmetric solutions we find (6.14) In Figure 6.3 we plot eigenvalues λ 1,2 , along the even-symmetric snaking branch for n = 0. The results show that solutions on this branch undergo a sequence of saddlenodes and pitchfork bifurcations, as indicated by the corresponding eigenfunctions. Similar results (not shown) are found for odd-symmetric states. For asymmetric solutions we obtain Here we have made use of the fact that, with our choice of the synaptic kernel, we have which is found by differentiating (6). Further, we note that By using (6.13), (6.17) and (6.18) we see that the eigenvalues λ 1,2 are such that λ 1 > 0 and λ 2 ≤ 0. As a consequence, all asymmetric solutions are linearly unstable. at which The snake-and-ladder bifurcation structure derived here for Heaviside firing rates is also found in the case of steep sigmoidal firing rates: in particular, we have performed numerical continuation for the firing rate function (2.4) with ν = 50 and found an analogous bifurcation diagram (not shown).
7. Changes in the modulation amplitude. The framework developed in the previous Sections can be employed to study two-parameter bifurcation diagrams. So far, we have fixed the parameters a, ε and used h as our main continuation parameter. It is interesting to explore how variations in secondary parameters affect the snaking branches. In [35], the authors explore variations in the spatial scale of the heterogeneity for the Swift-Hohenberg equation. Here, we concentrate on the amplitude a of the heterogeneity A(x) for the integral neural field model. Following the previous sections, we study the Heaviside case analytically and then present numerical simulations for the steep sigmoid case. We begin by considering Heaviside firing rate and outlining the region of parameter space where the trivial steady state q 0 and the above-threshold periodic solution q at coexist and are stable, that is, we follow the grazing point In Figure 7.1 we present a two-parameter bifurcation diagram and indicate with a green line the locus of grazing points (7.1): q at and q 0 coexist and are stable if (a, h) is below the green line. Next, we compute the snaking limits, for large L, as functions of h and a. We use the bifurcation equation for symmetric localized states, Equation (6.12), and find in the limit for large L the following snaking limits These curves are plotted in Figure 7.1 (solid blue lines). Further, we compute the loci of saddle-node bifurcations of the cross-threshold solutions q ct (which are labeled SN 1 and SN 2 in Figure 5.2) by solving for (a, h) the following system The loci of saddle-node bifurcations are plotted with dark magenta lines in Figure 7.1. The area between these two curves identifies a region in which q ct , q at and q 0 coexist and are stable. In passing we note that the curve for SN 2 intersects the curve for the grazing point B 1 at a = 1.
We found a snake-and-ladder bifurcation structure, as discussed in Section 6, in a wedge delimited by the lines h 1 and h 2 for a 0.57 (dark blue area in Figure 7.1).
Bifurcation diagram for a = 0.8 ν = 50, ε = 1. The snaking branch is composed by solutions with two (green), six (blue) or more (magenta) threshold crossings. The snaking structure reflects these three types of solutions and their occurrence is predicted adequately by the two-parameter bifurcation analysis for Heaviside case (reference intervals are reported on top of the bifurcation diagram). Stability is not indicated and a second intertwined branch of odd localized states is also found (not shown). Snaking branches in this region are formed of solutions with exactly two threshold crossings at x 0 ∓ L/2. However, there exist snaking branches of solutions with more threshold crossings. An example is given for the steep sigmoidal case for a = 0.6 in Figure 7.3: the snaking branch collides with neighbouring branches of solutions with multiple crossings and give rise to an intricate bifurcation structure.
In order to understand the occurrence of such curves we return to the Heaviside case and concentrate on the even-and odd-symmetric solutions featuring a threshold crossing followed by a threshold tangency at a local minimum (for an example with large L, see pattern 1 in Figure 7.1). More precisely, we denote by x * the point with largest absolute value at which q b attains a local minimum and solve for (a, h, x * ) the system where q b is given by Equation (6.10). We follow solutions to the system above as L varies in a given range and show the corresponding loci of solutions in the (a, h)plane in Figure 7.1: the dashed curve T w contains solutions to (7.2)-(7.4) with a wide active domain (L varies approximately between 45 and 57), whereas T n corresponds to solutions with a narrow active domain (L varies approximately between 0.7 and 12). Even though T w and T n are not loci of bifurcations, they are indicative of regions of parameter space where solutions with multiple threshold crossings may occur. In Figure 7.3 we show a bifurcation diagram for a = 0.8 for the steep sigmoid: the snaking branch is composed by solutions with two (green), six (blue) or more (magenta) threshold crossings. The snaking structure reflects these three types of solutions and their occurrence is predicted adequately by the two-parameter bifurcation analysis for Heaviside case (reference intervals are reported on top of the bifurcation diagram of Figure 7.3). Stable and unstable branches alternate in the usual manner and an intertwined branch of localized odd solutions exists as well (not shown). A similar scenario, with an even wider snaking diagram, is found for a = 1.2 (see Figure 7.4): for large modulation amplitudes, the bifurcation diagram also contains cross-threshold solutions, but this time their occurrence is marked by the grazing point B 1 and the saddle nodes SN 2 (see also the bifurcation diagram for a = 1 in Figure 5.2(b)).

Conclusions.
In the present paper we have studied the existence and bifurcation structure of stationary localized solutions to a neural field model with inhomogeneous synaptic kernel. For Heaviside firing rates, we computed localized as well as spatially-periodic solutions and we followed them in parameter space. We recovered the classical snakes and ladders structure that is found in the one-dimensional Swift-Hohenberg equation as well as previous studies in neural field models: for our model, however, both solutions and bifurcation equations are found analytically. Since linear stability can also be inferred with a simple calculation, it is possible to draw the snaking bifurcation diagrams analytically or semi-analytically (using elementary quadrature rules for the integrals).
Interestingly, we found that the interpretation of the snake and ladder structure proposed by Beck and co-workers [5] and extended by Makrides and Sandstede [45] is valid for the specific inhomogeneous case presented here, for both Heaviside and sigmoidal firing rates: it seems plausible that their framework could be extended to tackle the corresponding non-autonomous spatial-dynamical formulation (3.5).
With reference to the particular system presented here, we found that a harmonic modulation with an O(1) spatial wavelength promotes the formation of snaking localized bumps and we note that these structures are driven entirely by the inhomogeneity: in the translation-invariant case, a = 0, the system supports localized fronts belonging to a non-snaking branch (a scenario that is also found in the homogeneous Swift-Hohenberg equation [37,3]) We also remark that, in a wide region of parameter space, a ≤ 1, the kernel is purely excitatory, yet snaking stable bumps are supported. When a is further increased and the kernel becomes excitatory-inhibitory (a > 1), the snaking limits become wider and involve solutions with multiple threshold crossings. We note that with a modulated but translation-invariant kernel, A(x − y; a, ), the integral over the resulting kernel, W (x) would be monotonically increasing and would then prevent the formation of stable bumps for a < 1 [40]. The inhomogeneity is thus a key ingredient to produce stable solutions in the absence of inhibition when a < 1.
The analytical methods presented in this paper could be useful in the future to study time-periodic spatially-localized structures (often termed oscillons). A simple mechanism to obtain oscillatory instabilities in neural field models is by introducing linear adaptation [50]. This modification seems amenable to study oscillons, since localized bumps of the extended system can be constructed in the same way presented in this paper, yet the corresponding stability problem changes slightly and may lead to a Hopf bifurcation of the localized steady states. This approach has recently been used by Folias and Ermentrout [27] and Coombes and co-workers [18] in two component models supporting breathers and other spatio-temporal patterns. Another possible extension is to study the effect of spatial modulation in planar neural field models, in which case one could build upon the interface method developed in [19] for homogeneous planar neural fields. Here, we have made use of the fact that q(x) and A(x) are 2πε-periodic functions. Since we obtain the reduced formulation (5.1) with amended kernel (5.3).
Appendix B. Explicit solutions for cross-threshold solutions q(x). An explicit solution for equation (5.6) with Heaviside nonlinearity and kernel (2.2) is found by carrying out a direct integration, which gives