{'fraction_non_alphanumeric': 0.045021723662897195, 'fraction_numerical': 0.01439665743221763, 'mean_word_length': 4.261070404587448, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 20, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Scalar Reduction of a Neural Field Model with Spike Frequency Adaptation

We study a deterministic version of a one- and two-dimensional attractor neural network model of hippocampal activity first studied by Itskov et al 2011. We analyze the dynamics of the system on the ring and torus domain with an even periodized weight matrix, assum- ing weak and slow spike frequency adaptation and a weak stationary input current. On these domains, we find transitions from spatially localized stationary solutions ("bumps") to (periodically modulated) solutions ("sloshers"), as well as constant and non-constant velocity traveling bumps depending on the relative strength of external input current and adaptation. The weak and slow adaptation allows for a reduction of the system from a distributed partial integro-differential equation to a system of scalar Volterra integro-differential equations describing the movement of the centroid of the bump solution. Using this reduction, we show that on both domains, sloshing solutions arise through an Andronov-Hopf bifurcation and derive a normal form for the Hopf bifurcation on the ring. We also show existence and stability of constant velocity solutions on both domains using Evans functions. In contrast to existing studies, we assume a general weight matrix of Mexican-hat type in addition to a smooth firing rate function.


Introduction
Spatially coherent activity states exist during normal brain function including mammalian path integration, head direction tracking, visual hallucination, working memory, spatial object location, and object orientation [4,5,12]. Neural field models (also called continuous attractor neural networks) are one way to understand the mechanism underlying such spatially coherent phenomena [17,15,16]. In neural recordings and field models, these spatio-temporal dynamics manifest as traveling waves, spirals, or single or multiple localized "bumps" or "pulses" of activity [4].
Extensive literature exists on the analysis of these behaviors. In particular, [32] show the existence and stability of traveling bumps using multiplelayer neural fields. Several other studies use one of or a combination of short term depression and spike frequency adaptation. In [24], the authors show that traveling pulses exist in a model with synaptic depression and adaptation when synaptic depression is sufficiently weak. For stronger synaptic depression, the traveling pulse ceases to exist via a saddle-node bifurcation. In [15] the authors show that spontaneous motion of a bump solution exists for a neural field with only spike frequency adaptation, and in a similar neural field model with only short term synaptic depression. The authors in [28] show the existence of a traveling pulse solution in a neural field model with spike frequency adaptation. The previous two studies also show the existence of traveling wavefronts in their respective neural field models.
In addition to the analysis of traveling bumps or wavefronts, rich oscillatory solutions of neural fields are also possible. For example, with spatially localized input current and spike frequency adaptation, a bump solution may oscillate in diameter (breathers) [3,11,12,13], which may play a role in generating epileptiform activity [13] and the processing of sensory stimuli [12]. There also exist studies of a combination of traveling and breathing pulses in an inhibitory-excitatory neural field [10]. In addition to breathers, there exist pulse-emitting neural fields [24,23], oscillatory wavefronts [2,3], and spiral waves [23].
Despite this large body of literature, the analyses often require particular assumptions. For example, the existence of "sloshing" solutions -bump solutions that oscillate periodically in the centroid -that arise through a Hopf bifurcation is known under certain assumptions. In early work, sloshers are shown to exist numerically using a rate model with a threshold nonlinearity [18]. In recent work, the authors of [6] show the existence of a Hopf bifurcation with a cosine kernel and a particular choice of smooth firing rate function. In [9], Folias computes a normal form for the Hopf bifurcation using a general kernel, but for a Heaviside firing rate function.
Proving existence of other phenomena also require special assumptions. In [8], the authors consider a neural field model on the real line with synaptic depression and prove the existence of a traveling pulse without a Heaviside assumption, but use the particular choice of a normalized exponential kernel. In [22], the authors use a center manifold reduction to analyze the existence of moving bump solutions. They allow the firing rate to be sigmoidal or a Heaviside, but require a cosine kernel. Similar assumptions are made in [], where they assume a hyperbolic tangent firing rate function and a cosine kernel.
The most general of such studies, [28], considers a neural field model on the real line with spike frequency adaptation and a singular perturbation approach to construct a constant velocity traveling pulse on the real line with a general firing rate function and a general kernel. However, the existence of other phenomena are not shown.
In this paper, we introduce a method to analyze the dynamics of a neural field model on a one-and two-dimensional domain with periodic boundary conditions and assume a smooth firing rate and an even, periodic kernel. Using our method, with standard numerical and analytical dynamical systems tools, we show existence and stability of traveling pulse solutions and oscillatory dynamics. In particular, we analyze sloshing solutions on the ring and torus.
The neural field we consider in this paper is defined as where the parameter ε is small, 0 < ε 1, and x, y ∈ R n . For n = 2, the kernel function K is an even function in the sense that, K(−x, y) = K(x, −y) = K(x, y), and doubly periodic in the sense that K(x + 2nπ, y + 2mπ) = K(x, y), for any integers n, m. The function w represents heterogeneity of neural connections, and q, g, β are constants. For convenience, we will denote the domain Ω = [−π, π) m , with m = 1, 2. Thus in onedimension the domain is a ring and in two-dimensions a torus. The variable z(x, t) represents linear adaptation [28] and I(x) an external input to the network. External inputs represent persistent stimuli that can be used to entrain the bump and move it to a specific location ( [1]. We have chosen to make both the timescale of adaptation and its magnitude to be small. While there is good biological justification for the former assumption as there are many forms of slow adaptation ( [20] section 7.4), the assumption that the adaptation is small is less biological. For the existence of traveling waves, adaptation need not be small ( [28]), but in order to study how the adaptation interacts with stimuli, we need both the adaptation and the stimuli to be the same order of magnitude. The effects of large stimuli to general neural field models are not easy to analyze, so that by treating them as perturbations, we are able to consider the effects in a great deal of detail. Thus, one can regard this assumption as a starting point for the continuation of these phenomena to large amplitude stimuli and adaptation.
Our goal in this paper is to analyze Equations (1), (2) when ε is small. When ε = 0, there is a stable "bump" attractor, u 0 (x), in the scalar neural field (1), i.e., a local stationary peak of u(x, t) centered at x = 0. The bump attractor satisfies where u 0 is nonconstant and even.
Although we allow for a general even, doubly periodic kernel in one-and two-dimensions and a general smooth threshold nonlinearity f , we make particular choices for numerical simulations. We choose f as where r = 15, u th = 0.25. In the one-dimensional case, we choose the kernel to be K(x) = A + B cos(x) with A = −0.5, B = 3 unless stated otherwise.
In the two-dimensional case, we form the Mexican-hat function, where r ≡ r(x, y, n, m) = (x + 2πn) 2 + (y + 2πm) 2 . We make the the functionK periodic in two dimensions using the definition (r(x, y, n, m)).
The parameters here are where σ e = 2, and σ i = 3. For numerical simulations, we find it sufficient to replace the infinite sum with a finite sum from n, m = −5 to n, m = 5. This is because the functionK(r) is a decaying exponential and therefore negligible for large r. For example, if a bump solution remains close to the origin, contributions from terms a distance of 10π (i.e., n or m=5) are negligible because exp −(10π) 2 ≈ 2 × 10 −429 .
To analyze particular dynamics in more detail, we numerically compute the period kernel above, then take the Fourier truncation of this doubly periodic kernel, K(x, y) = k 0 + k 1 (cos x + cos y) + k 2 cos x cos y.
We now outline the organization of the paper, as follows: We reduce Equations (1),(2) to a set of integro-differential equations for the centroid of the bump solution on the ring and torus. We study bifurcations of these equations using numerical and analytical techniques to show existence and stability of constant velocity traveling bumps and sloshing bumps. Depending on the parameter values g, q, these traveling bumps may traverse the domain periodically or exhibit chaos. Next we turn to the torus domain and perform similar analyses: we study bifurcations of these equations using numerical and analytical techniques to show existence and stability of constant velocity traveling bumps. In addition to the sloshing solutions found in the one-dimensional model, we also find several types of traveling bumps and modulated traveling bumps that densely fill the torus. We also find chaotic motion in some cases. We conclude with a discussion and some contrasts to previous analyses. We remark that all figure generation code and relevant data files with documentation is available on GitHub at https://github.com/youngmp/park_and_ermentrout_2017

Derivation of the Phase Equation
We start with Equations (1), (2). Let τ = εt be a slow timescale and assume that both z and u depend only on (x, τ ). In this case, we can integrate equation (2) to obtain: Since we are mainly interested in long term behavior, we can ignore the first exponentially decaying term. With these assumptions, we obtain the following scalar integro-differential equation: We will assume u(x, τ ) = U (x, τ, ε) and expand U as a power series in ε to get an approximate solution. Thus, Substituting this power series into (3), we get (with a bit of re-arrangement): where The equation for U 0 (x, τ ) is the equation for the bump solution and since it is translation invariant, we see that where θ(τ ) is a τ −dependent phase shift of the bump. Our goal, then is to determine the dynamics of θ(τ ). Figure 1 shows typical examples of the stationary bump U 0 (x) for one-and two-dimensions. On the ring, we denote the centroid by θ, while we denote the centroid of the bump on the torus by (θ 1 , θ 2 ). Our phase model (Equation (9)) describes shifts in the centroid.
Before continuing with the perturbation calculation, we establish a few preliminaries. We define the compact linear operator and establish several properties of it. Recall that the bump, u 0 (x) satisfies By making a change of variables and noting that all functions are periodic in x (that is, periodic in each of the components of x), then u 0 (x) satisfies Recalling that the domain is Ω = [−π, π) m , with m = 1, 2, we center u 0 at the origin. Thus, u 0 (x) is an even periodic function of x, component-wise.
Let ∂ i u(x) denote the partial derivative of u(x) along the x i direction where x = (x 1 , x 2 ). If we differentiate (6) along one of the axes, we see that and changing variables again, we rewrite this as so we see that L 0 ∂ i u 0 (x) = 0. In other words, the linear operator, L 0 has an m−dimensional nullspace spanned by the principle directional derivatives of u 0 (x). With the natural inner product By multiplying equation (7) by f (u 0 (x)), we see that the nullspace of L * is spanned by v * i (x) = f (u 0 (x))∂ i u 0 (x). Since u 0 (x) is an even periodic function, componentwise, we note that ∂ 1 u 0 (x) is even in x 2 and odd in x 1 where x = (x 1 , x 2 ); v * 1 (x) has the same property, while ∂ 2 u 0 (x), v * 2 (x) are even in x 1 and odd in x 2 . These properties imply the ∂ i u 0 (x), v * k (x) = 0 when i = k. We also have Finally, the Fredholm alternative holds for L 0 . That is, for any continuous periodic function b(x), has a bounded solution if and only if v * i (x), b(x) = 0 for i = 1, . . . , m [21]. With these technical issues aside, we turn to equation (5), which we can rewrite as Writing θ(τ ) = (θ 1 (τ ), θ 2 (τ )) and applying the m conditions for the Fredholm alternative, we arrive at We note that because of the symmetry of u 0 (x), the functions, H i (θ) have a similar symmetry which we will exploit in the analysis of Equation (8).
The derivation here has been fairly general and holds in any dimension although we will focus only on one-and two-dimensional bumps in this model. Since both W i and J i have no explicit time dependence and act mainly as heterogeneities, we will ignore W i and focus on J i which is conveniently parameterized by q. Figure 2 shows the functions H(θ), J(θ) in the onedimensional case, while Figure 3 shows the functions H i (θ), J i (θ) in the two-dimensional case. We have reduced the problem of the bump dynamics to slow timescale phase shifts of the bump solution, represented as an integro-differential equation. For simplicity and convenience, we ignore transients by changing the Figure 3: Numerically computed functions for the two-dimensional phase model A: H 1 , B: H 2 , C: J 1 , D: J 2 . Note that H 2 (x, y) = H 1 (y, x). We always choose I(x) = u 0 (x). That is, we always use the steady-state bump as the pinning function. With this choice, J i = −H i in 1-and 2-dimensions. (8) from [0, τ ] to (−∞, τ ], and apply the change of variables ξ = τ − s, yielding

limits of integration in Equation
With a trivial change of notation, we arrive at the equations We study stability properties and bifurcations in this form. Note that H i is implicitly a function of the kernel K.
To facilitate calculations, we first prove the following statements: . In particular, H 1 is odd in the first coordinate and even in the second coordinate.
3. If the input current I(x) is defined as the steady-state bump solution, For the first statement, fix θ 1 , θ 2 and consider the sum H 1 (−θ 1 , θ 2 ) + H(θ 1 , θ 2 ). By definition this sum is the sum of integrals Given x 2 , and for the sake of clarity, consider the temporary function Φ( These lines follow by the even assumption on each coordinate oa:simplifyf the bump solution u 0 . The remaining terms in the integrand, ∂ 1 u 0 (x), and f (u 0 (x)), are odd and even in x 1 , respectively. Thus, the integrand is odd in x 1 and the integral evaluates to zero for each x 2 (and indeed, for each θ 1 , θ 2 ). It follows that H 1 (−θ 1 , θ 2 ) + H(θ 1 , θ 2 ) = 0, i.e., that the first coordinate is odd.
To show that the second coordinate is even, we use a similar argument. Again, fix θ 1 , θ 2 and consider the sum H 1 (θ 1 , −θ 2 )−H(θ 1 , θ 2 ). By definition, this sum is the sum of integrals Given x 1 , we redefine our temporary function Φ as Φ( Again, these lines follow by the even assmption on each coordinate of the bump solution u 0 . The integrand term ∂ 1 u 0 (x) is even in the second coordinate as is the term f (u 0 (x)). Thus, the integrand is odd in x 2 and the integral evaluates to zero for each x 1 (and indeed, for each θ 1 , θ 2 ). It follows that H 1 (θ 1 , −θ 2 ) − H(θ 1 , θ 2 ) = 0, i.e., that the second coordinate is even.
We have shown that H 1 is an odd function that is odd in the first coordinate and even in the second coordinate. The proof of H 2 being an odd function that is even in the first coordinate and odd in the second follows using the same arguments, or by using the second statement, which we prove next.
To prove the second statement, we proceed by definition.
The steady-state bump solution is invariant under reflections about the unit line, and due to the radial symmetry of the bump solution, its partial derivatives are related by Next we relabel the coordinates and flip the order of integration Then we flip the coordinates of u 0 , and the resulting integral is by definition H 2 (θ 2 , θ 1 ): To prove the third statement, suppose that a functionĥ on a periodic two-dimensional domain [0, 2π] × [0, 2π] is odd in the first coordinate and even in the second so thatĥ(x 1 , x 2 ) =ĥ(x 1 , −x 2 ) = −ĥ(−x 1 , −x 2 ). In particular, it follows that for a given value x 2 , This integral property holds whenĥ is even in the first coordinate and odd in the second with a similar argument.
If we choose I(x) to be the steady-state bump, then Then taking the sum H i (θ) + J i (θ) yields For a given θ, the term [u 0 (x+θ)+u 0 (x−θ)] in the integrand is even in both coordinates. The remaining term, f (u 0 (x))∂ i u 0 (x), when i = 1 (i = 2), is odd (even) in the first coordinate and even (odd) in the second. Therefore, when i = 1 (i = 2), the integrand is an odd function in the first (second) coordinate and the integral evaluates to zero. It follows trivially that This property remains true on the ring using the same argument. These statements will come in useful in the sections to follow. We now proceed with an analysis of the reduced equations on the ring domain.
Remark We have assumed linear adaptation in our derivation of the reduced model, but, this is not necessary. We could replace Equation (2) by where M (u) is an arbitrary monotonically increasing continuously differentiable function. In this case, we find The new version of H i has exactly the same properties as the linear case since M (u 0 (x)) is an even function and its derivative with respect to x is an odd function.

The Ring Domain
In this section, we choose the domain Ω to be the ring. First, we thoroughly analyze the full neural field model through a bifurcation analysis. We then turn to Equation (9) on the ring and perform the same bifurcation analysis and through analytical study.

Equivalent Neural Field Model on the Ring
To classify the bifurcations of the full neural field model on the ring, we transform the equations to an equivalent 6-dimensional system of ODEs, allowing us to use dynamical systems software and techniques to analyze the model. Recall that for numerical simulations on the ring, we choose a cosine kernel K(x) = A + B cos(x). This technique and choice of kernel is the same as that used in [25], where as part of the study they analyze a rate model similar to the model in the current study, but in contrast, the adaptation and input current terms are input directly to the firing rate function. They provide sufficient detail with regards to transforming their rate model to a system of ODEs, but as the details differ from our model, we include the derivation of our model here (in particular they include a phase lag between the peak of the bump activity u and the peak of the adaptation activity z which results in slightly different equations). Note that with this choice of kernel, the bump solution is also sinusoidal and without loss of generality takes the form u 0 (x) = C + D cos(x). For simplicity we choose J(x) = u 0 (x). We are now ready to transform the equations.
Since the functions u(x, t), z(x, t) are periodic in x, we expand them in a Fourier series, u(x, t) =â 0 (t) + ∞ n=1â n (t) cos nx +b n (t) sin nx, z(x, t) =ĉ 0 (t) + ∞ n=1ĉ n (t) cos nx +d n (t) sin nx. and plug into equations (1), (2). First, a direct substitution into the dynamics of u yieldŝ a 0 + ∞ n=1â n cos(nx) +b n sin(nx) We have used the elementary trigonometric identity cos(x−y) = cos(x) cos(y)+ sin(x) sin(y) to separate the kernel into multiple integrals. A direct substitution into the dynamics of z yieldŝ n (t) cos nx +b n (t) sin nx Next, we group like terms in the Fourier basis, starting with the Fourier coefficients of u: We repeat this grouping for the Fourier coefficients of z: The pattern is clear at this point: The coefficients of all Fourier modes greater than 1 satisfy where a i are placeholders for the Fourier coefficients of u and b i are placeholders for the Fourier coefficients of z. Through an elementary stability analysis, all solutions to these equations decay to zero so they are unnecessary to consider. We proceed with the remaining nontrivial terms, Note that these are still the first two Fourier modes, but we have dropped the hat notation and relabeled the coefficients. Using this notation, we have the system Note that we do not need to explicitly write the full Fourier series of f or extract any of its coefficients. In fact, the Fourier modes of f (u) greater than 1 vanish as we will now show. Consider the Fourier series of f (u(x, t)): This expansion exists because f is bounded and integrable on [0, 2π]. We now evaluate each integral, Ω f (u(y, t))dy, Ω cos(y)f (u(y, t))dy, and Ω sin(y)f (u(y, t))dy in turn. First, α n cos(ny) + β n sin(ny)dy α n cos(ny) + β n sin(ny) dy α n Ω cos(y) cos(ny)dy + β n Ω cos(y) sin(ny)dy = πα 1 and finally, Thus, the nonlinearity f in the integrand only appears in the dynamics of the first few Fourier coefficients. At each time step in the numerics, we compute the integrals Ω f (u(y, t))dy, Ω cos(y)f (u(y, t))dy, and Ω sin(y)f (u(y, t))dy using Riemann integration at each time step as it is more straightforward than extracting the necessary Fourier coefficients.
We focus our numerical studies on the coefficients a 1 and a 2 because they produce the most salient features of the bump solution (the a 0 coefficient changes as a function of time, but only up to order O(ε), while the b i terms represent aggregate behavior of the adaptation variable z). By following the fixed points and oscillatory behavior in a 1 and a 2 , we produce a bifurcation diagram of this system in Figure 4. Figure 4 shows that there are three main solution types: the pinned or stationary bump, the sloshing bump, and the traveling bump, which traverses the ring at some finite speed. In addition there are small regions of bistability between the sloshing bump and the traveling bump.

Phase Model on the Ring
We now turn to the analysis of the phase dynamics on the ring. The analysis to follow depends on proving the following statements: depends on the parameters A, B. Recall the scalar version of the functions H i (θ) and v * i (x): Because u 0 is even, it follows that f (u 0 ) is even, u 0 is odd, and therefore u * is odd. Noting that where the function u * (y)u 0 (y) is odd, the first statement follows.
For the second statement, we follow the definitions to arrive at The function f is an increasing sigmoidal, thus f > 0. In addition, Next, we prove the third statement. With the kernel choice K(x) = A + B cos(x), the steady-state bump solution is some shifted multiple of cosine, u 0 (x) = C + D cos(x), where C, D implicitly depend on the kernel parameters A, B. Plugging this u 0 into H(θ) yields Let h(y) = f (C + D cos(y))[D sin(y)], which is an odd function. Recalling that cos(y + x) = cos(y) cos(x) − sin(y) sin(x), H simplifies to Because h(y) is odd, some integrals cancel, and we are left with where A = D π −π h(y) sin(y)dy. From these statements, it follows that dJ dθ θ=θ < 0 and J(θ) = 0 where θ represents a steady-state bump peak. WLOG, we let θ = 0 because we generally choose the center of the steady-state bump to be the origin.
Remark: For the more general adaptation (c.f. above), as long as M (u) is differentiable and monotonically increasing, we still have that H (0) > 0.

Equivalent Phase Model on the Ring
We next show that there are really only two relevant parameters. We can rescale time to obtain where we have re-used τ as the now scaled time βτ . Next, divide by µβ to obtain withq = q µβ andĝ = g µβ . This rearrangement shows that making adaptation slower by decreasing β is equivalent to increasing the rescaled parameters ĝ andq. For analytic calculations, we will often reference this equation without the hats on the parameters.
For numerical studies of bifurcations in this system, we let H(θ) = A sin θ, from which J follows immediately (see Equation (10) and statement 2 above). We once more abuse notation and absorb A intoĝ and intô q, then drop the hats. So we will now study To numerically integrate this phase equation, we rewrite this differential equation as a system of three equations by exploiting basic differentiation properties of integrals. To begin, we use a trigonometric identity to rewrite the integral in the right-hand side With the change of variables s = τ − s, S, C become By differentiating, we rewrite S and C as ODEs: We have transformed a single integro-differential equation into a system of three ODEs, simplifying the numerics considerably: The bifurcation diagram in Figure 5 summarizes the dynamics of the phase model on the ring. On the left panel, we fix a parameter value q = 0.5 and as we vary the parameter g, the system transitions from steadystate to sloshing solutions, then to a co-existence of large-amplitude and relatively small amplitude sloshing solutions, and eventually to a steady traveling pulse. On the right panel, we find that the parameter space is separated into several regions. In particular, for q ≥ 0 arbitrarily small, there exists a traveling bump for some nonzero g.
In the following sections, we analyze the existence of these bifurcations including the Hopf bifurcation leading to sloshing solutions, and the saddlenode bifurcation leads to the constant-velocity traveling bump.

Constant Velocity Bump Solution on the Ring
To show the existence of a constant velocity traveling bump solution, we require that q = 0 and g > 0. For the first part of this analysis, we do not require the kernel to take a particular form. We only require the kernel to be even and admit a steady-state bump solution to Equation (1). We make a traveling bump ansatz, θ(τ ) = ντ , where ν corresponds to the traveling bump velocity. We first determine the existence and stability of the zero velocity bump solution. Plugging the ansatz into Equation (11) yields where the last line follows by the oddness of H. Because H(0) = 0, ν = 0 is a solution. To determine the stability of the zero velocity solution, we consider a small perturbation, θ(τ ) = ντ + εψ. By plugging this perturbation into Equation (11), we extract the dynamics of the perturbed variable ψ, Assuming ψ(τ ) = e λτ and ν = 0, we obtain the stability equation, We integrate the right-hand side and rearrange to yield Thus, either λ = 0, or λ = −1+gH (0). Moreover, the zero velocity solution becomes unstable when g > 1/H (0). In general, we may view the relationship between g and ν by rearranging Equation (12) into the function We show examples of Γ in Figure 6. In the left panel, the relationship between the adaptation strength g and bump velocity ν is straightforward: as g increases, there is some critical value ν where a nonzero velocity traveling bump exists. However, the choice of kernel may change the shape of Γ, and therefore change the relationship between g and ν, as well as the stability of traveling bump solutions. For example, a kernel of the form K(x) = a + b cos(x) + c cos(2x) results in an H function of the form Using this H function to plot Γ results in the right panel of Figure 6. The branch with negative slope represents another traveling bump solution. We now show that if Γ (ν) < 0, then the traveling bump with velocity ν is unstable. Recall again Equation (13). Assuming ψ(τ ) = e λτ and ν = 0, we obtain the stability equation, To prove the statement, we seek to show that lim λ→∞ f (λ) > 0 and lim λ→0 f (λ) < 0. Then by continuity of f , there exists a positive root to Equation (15). We take each limit in turn, starting with the limit as λ → ∞.
Thus, the positive λ limit is positive. For the other limit, we rearrange Equation (14) into where D(ν) = ∞ 0 e −βs H(νs)ds, and differentiate with respect to ν to obtain Solving for Γ (ν) yields Note that D(ν) > 0 at least within a neighborhood of ν = 0 since H(0) = 0 and H (0) > 0. In addition, D (ν) = ∞ 0 e −s H (νs)sds. Using the hypothesis that Γ (ν) < 0, we have the inequality We use this fact in the next limit Thus, the zero λ limit is negative. Because f (0) is negative, and f (λ) is positive for asymptotically large values of λ, there exists a positive root λ of f (λ) by continuity. It follows that branches of Γ(ν) with negative slope indicate an unstable traveling bump at least within a neighborhood of ν = 0. For the next part of this analysis, we show how to compute a formula for the velocity of the traveling bump when the kernel is K(x) = A + B cos(x). With this kernel, the H function is proportional to sin(x), and Equation (12) becomes explicitly computable. Computing the integral results in a formula for the nontrivial bump velocity ν, Equation (16) corresponds to the branches of a pitchfork bifurcation in the velocity of the traveling bump. We show a particular example of a constantvelocity traveling bump in Figure 7A. We note that any odd H will lead to a pitchfork bifurcation to a traveling bump. In particular, it is trivial to derive the following bifurcation equation: This equation tells us that the pitchfork bifurcation is super-critical if H (0) < 0 and sub-critical otherwise.

Andronov-Hopf Bifurcation on the Ring
We now prove the existence of a Hopf bifurcation. For this analysis, we do not require H or J to take a particular form. However, we do require that H and J be sufficiently differentiable, along with the properties H(0) = J(0) = 0, H (0) > 0,J (0) < 0, H odd, and g, q > 0. Consider again the simplified phase model, Equation (11). Let us fix q and absorb the parameter into J. We write J and H as Taylor expansions, Then to first order, Letting θ = e λt and rearranging the resulting equation yields or equivalently, Since j 1 < 0 and h 1 > 0, there exists a Hopf bifurcation when This bifurcation leads to oscillations in the peak of the bump solution. We show a particular example of this oscillatory behavior in Figure 7B.

Normal Form for the Hopf Bifurcation on the Ring
We wish to analyze the bifurcation to a sloshing pulse for the general integral equation: as g increases. For simplicity, we will assume J(θ) is an odd periodic function (as is the case for H(θ) and through suitable rescaling of g, q, we will assume: We also assume q > 0 so that θ = 0 is stable without adaptation. If, we use H(θ) = J(θ) = sin(θ), then j 3 = h 3 = −(1/6). The linearization about θ = 0 has the form: which has the general solution, e λτ . After some simplification, we find that so there is an imaginary eigenvalue, i √ q := iω when g = 1 + q ≡ g 0 , so we expect a Hopf bifurcation will occur. All nonlinearities are odd, so we can assume the multiple timescale expansion where δ is the amplitude of the bifurcating solution, ζ = τ is a "fast" time, and ξ = δ 2 τ is a "slow" time. We detail the remaining steps of the normal form analysis in Appendix A and jump to the conclusion, where To get the actual normal form, we divide (18) by α, to obtain: If we assume that j 3 = h 3 as would be the case if the input was the bump, itself, then We compare the normal form calculation to the numerics in Figure 8. We use XPPAUTO to compute the numerical bifurcation diagram. As expected, the normal form approximation is quite accurate near the bifurcation.

Non-Constant Velocity Bump Solution on the Ring
When adaptation is made even stronger, the solution breaks free from the oscillating state and travels across the periodic domain ( Figure 7C. The onset is shown numerically in the bifurcation diagrams of Figures 4 and 5, purple LP2). Due to the pinning term, the velocity of the bump is nonconstant.
Remark. As q → 0, we see in figure 5B that all the twoparameter curves converge to the point g * which is the point of onset of the traveling bumps with no input stimulus.

Chaos on the Ring
With g, q > 0, there exists a small parameter range in which the neural field exhibits chaotic movement about the ring. Examples of this behavior are shown in Figure 9. In both panels, the initial conditions differ by 1e-7. The solutions in each panel remain nearly identical for a long time (we have truncated a significant portion of the simulation). This section completes our analysis of the one-dimensional case. We have found a good match between the phase-reduced equations and the full neural model. For a fixed amplitude of the external input, we find a transition from a stationary bump to "sloshers", and, finally to bumps that move completely around the ring, modulated traveling bumps. In the sections to follow, we repeat the analytical and numerical analysis for the two-dimensional domain.

Torus Domain
In this section, we define the domain Ω as the torus, or the square [−π, π] × [−π, π] with periodic boundary conditions. We seek to analyze the full neural field model on this domain using the same bifurcation analysis performed in the one-dimensional case. We begin by considering simplifications that allow us to use standard bifurcation analysis tools like XPPAUTO.

Approximation of the Neural Field Model on the Torus
In order to numerically investigate the full neural field equation (1-2), we need to either discretize space in two-dimensions or use an approximation of the kernel that is degenerate. (Since the integral operator is compact, it can be approximated to arbitrary precision by a degenerate integral operator ;see section 2.8 [29]). Thus to study the dynamics of the full neural field model on a two-dimensional domain, we take a Fourier truncation of the kernel to make the integral in Equation (3) separable. This truncation allows us to rewrite the infinite dimensional system as a finite system of ODEs and use traditional dynamical systems tools like XPPAUTO to analyze the system. To begin, take the truncated Fourier approximation to the kernel, K(x, y) = k 00 + k 10 cos(x) + k 01 cos(y) + k 11 cos(x) cos(y), (19) and plug it into Equation (3): After expanding the kernel using standard trigonometric identities, we derive the time-varying solutions u(x, t) = a 00 (t) + a 10 (t) cos(x 1 ) + a 01 (t) cos(x 2 ) where the coefficients satisfy where i = 0, 1 and j = 0, 1. The dummy variables ξ, ζ represent each of the pairs (a ij , E ij ), (b ij , F ij ), and (c i , G i ). The time-varying functions p ij , r ij , s i are defined as p 00 (t) = Ω f (u(y, t))dy p 01 = Ω cos(y 2 )f (u(y, t))dy p 10 (t) = Ω cos(y 1 )f (u(y, t))dy p 11 = Ω cos(y 1 ) cos(y 2 )f (u(y, t))dy r 01 (t) = Ω sin(y 2 )f (u(y, t))dy r 10 = Ω sin(y 1 )f (u(y, t))dy r 11 (t) = Ω sin(y 1 ) sin(y 2 )f (u(y, t))dy s 1 (t) = Ω cos(y 1 ) sin(y 2 )f (u(y, t))dy s 2 = Ω sin(y 1 ) cos(y 2 )f (u(y, t))dy, and the coefficients u ij are taken from the truncated Fourier series of the steady-state solution, u 0 (x, y) = u 00 + u 10 cos(x) + u 01 cos(y) + u 11 cos(x) cos(y).
The coefficient values of the kernel and steady-state bump solutions are shown in Tables 2 and 3, respectively. Details on how we approximate the spatial integrals of p ij , r ij , s i are in Appendix C.
We show the bifurcation diagram and salient solutions of this system in Figure 10. The only bifurcations classified by AUTO are a subcritical Hopf bifurcation (HB, orange) and limit points (limit points occur at each change in stability of periodic solutions). The Hopf bifurcation leads to small amplitude unstable solutions. We attribute the existence of a subcritical Hopf bifurcation to the coarse discretization of the spatial domain in the parameters of Equation (21). Although not shown here, a finer discretization of the spatial domain using 200 intervals results in a qualitatively supercritical Hopf bifurcation.
To summarize, we find the usual types of oscillatory solutions in this truncated neural field model as we found in the neural field model. A stable limit cycle of this system is shown in the first bottom left panel (A, which corresponds to the initial conditions taken from the point A in the bifurcation diagram). A large-sloshing solution exists for slightly larger g (B,C), and eventually, for sufficiently large g, there exist only traveling bump solutions D. Additional non-periodic attractors are shown in panels E-G. The attractors shown in this figure are simply those with the largest basins of attraction. Generally, starting random initial conditions with g anywhere in the range 1.55 < g < 1.95 (the gray shaded area labeled F in the main plot) results in solutions that qualitatively match panel F. The same holds for the shaded areas E and F, with their corresponding panels. Indeed, there exist several attractors not shown in this figure that are more difficult to find numerically. However, the focus of this study is not the thorough classification of attractors in the truncated neural field model, so we move on to the analysis of the phase model on the torus.

Approximations of the Phase Model on the Torus
We now turn to the analysis of the phase dynamics in two-dimensions and begin by reducing the number of parameters with the same rescaling used to obtain Equation (11) in the one dimensional case, and recall that Details on the numerical integration of the integro-differential equation, Equation (22), are shown in Appendix C. In panels E-G, we show solutions in parameter regimes without stable periodic attractors. These solutions are displayed in a relatively short time window after integrating for long times and travel from light to dark. In panel E (g = 1.4), we integrate for t = 5000 time units and show the last 30% of the data. In panel F (g = 1.7), we integrate for t = 8000 time units and show the last 10% of the data. In panel G (g = 2.4), we integrate for t = 8000 time units and show the last 6% of the data. We initialize the solutions of panels E-F using standard normally distributed random variables. Parameter ε = 0.01.
To facilitate the study of existence and stability of solutions, we consider two approximations to H i to be used in Equation (22): the first is a highaccuracy Fourier series of H i , and the second is a low-accuracy Fourier series of H i . We detail these approximations in turn.
For the accurate Fourier series approximation of H i , we use one of two Figure 11: Cartoon of the parameter space of the approximate neural field model on the torus (Equation (21)). The most salient solutions are shown. Solutions advance in time from light to dark, thin to thin. For sufficiently small g or sufficiently large q, the bump solution tends to a stationary solution. By increasing g or decreasing q to g = 1.2, q = 0.1, the centroid of the bump solution oscillates about the origin. For larger g, say g = 4, q = 0.5, the solution begins to traverse chaotically about the domain. When q = 0, there exists a constant velocity traveling bump solution for g sufficiently large, e.g., g = 3. Parameter ε = 0.01. equivalent forms where, due to the odd (even) property of the first (second) coordinate, the coefficients have the property that a n,±m = −a −n,±m . We can then rewrite this Fourier series into the equivalent form, whereâ nm = 4a nm . This equivalent form makes integrals much easier to compute. We use both forms interchangeably as we see fit, and abuse notation in Equation (24) by removing the hats from the coefficients. We find that 30 Fourier coefficients provides a sufficiently good approximation for simulations on a 64 × 64 domain (the error is on the order of 1e-7). For the low-accuracy Fourier series, we consider a more substantial truncation of the interaction function using only 3 Fourier coefficients. While this truncation is drastic, it allows us to analyze Equation (22) more rigorously. We derive the 3 term H i -function starting with the same Fourier truncation of the kernel as above, which leads to the same steady-state bump solution, u 0 (x, y) = u 00 + u 10 cos(x) + u 01 cos(y) + u 11 cos(x) cos(y), which in turn leads to a truncated H i function, H F 1 (θ 1 , θ 2 ) = sin(θ 1 )(h 10 + h 11 cos(θ 2 )), where h 10 = 4u 10 Ω sin 2 (x)(u 10 + 2 cos(y)u 11 )dxdy, h 11 = 8u 11 Ω cos(y) sin 2 (x)(u 10 + 2 cos(y)u 11 )dxdy.

For simplicity, analysis of this H i function uses the simpler form
where we have absorbed h 10 into the parameter g of Equation (22), and relabeled h 11 /h 10 as b. Naturally, it follows that H F 2 (x, y) = H F 1 (y, x) and H F i = −J F i . Using the truncated interaction function H F 1 enables us to use traditional dynamical systems tools and techniques to identify qualitative dynamics of Equation (22) through a bifurcation analysis.
We show the many bifurcations and salient solutions of this system in Figure 12. We find that there exists a stable sloshing bump solution that arises from a Hopf bifurcation (solution A, bifurcation HB). Due to the symmetry of the system, there is also an unstable sloshing solution in an axial direction (G) that arises from the same Hopf bifurcation. For slightly larger parameter values, there is bistability of large-sloshing solutions (H and B), and an even larger-sloshing solution (F) that loses stability through a torus bifurcation (TR). For this choice of q = 0.1, the solutions are chaotic for parameter values between the first torus bifurcation (TR1) and the first period doubling bifurcation (PD1). The multitude of period doubling bifurcations beyond this point represents the onset of chaotic behavior of the system due to the Fourier truncation, the error of which is proportional to g. This error is more apparent beyond limit point LP2 where there are no more traveling bump solutions, which qualitatively disagrees with the original system where traveling bump solutions exist for even relatively large g.
The most salient bifurcations are captured in the two-parameter bifurcation in Figure 13. There are several qualitative similarities to the two parameter bifurcation diagram of the phase model on the ring. In particular the transition from the stationary bump to the sloshing bump, and from sloshing to large-sloshing. However, for the two-dimensional domain, there are much larger regions of chaotic solutions.  Figure 13. Panels A-E show stable attractors. In panels F-H, we show solutions in parameter regimes without stable periodic attractors. These solutions are displayed in a relatively short time window after integrating for long times and travel from light to dark. In panel F (g = 1.05), we integrate for t = 500 time units and show the last 9% of the data. In panel F (g = 1.5), we integrate for t = 500 time units and show the last 7% of the data. In panel G (g = 2.15), we integrate for t = 500 time units and show the last 5% of the data. We initialize the solutions of panels F-H using standard normally distributed random variables. Parameter ε = 0.01. In the following sections, we study the dynamics of the original phase model and the truncated phase model and repeat most of the analysis as completed in the ring domain. In particular, using a combination of numerical and analytical methods, we analyze the existence and stability of traveling bump solutions, and the existence of a Hopf bifurcation.

Constant Velocity Bump Solution on the Torus
In this section, we analyze the existence and stability of constant velocity bump solutions on the torus for q = 0, the only case in which there can be constant velocity traveling bumps. Figure 14

Existence
To show existence of solutions in the axial directions, we only need to show existence of the solution θ 1 (τ ) = ντ and θ 2 (τ ) = 0. We plug this ansatz into (22) and rearrange to yield The analysis of this equation is identical to the one-dimensional case, Equation (14). By varying ν from zero, we can find the values of g where there are solutions. Values of g for which ν is nonzero imply there exists a traveling bump solution. In the case of the truncated H function, we compute this integral explicitly to derive the velocity ν as a function of adaptation strength g: To determine the critical value for the existence of axial constant velocity bump solutions, we take the limit lim ν→0 Γ(ν): To show existence of non-axial solutions, we use the ansatz θ 1 (τ ) = ν 1 τ and θ 2 (τ ) = ν 2 τ where ν 1 , ν 2 = 0. There exist non-axial traveling bump solutions if ν 1 , ν 2 simultaneously satisfy 0 = −ν 1 + gG(ν 1 , ν 2 ), where We can not compute the velocities ν 1 , ν 2 explicitly as a function of g, but we can exploit the Fourier series of H i to compute G explicitly, allowing us to use XPPAUTO to follow the velocities as a function of the adaptation parameter g. The existence of traveling solutions using the accurate Fourier series is shown in Figure 15, and the existence of traveling solutions using the truncated Fourier series is shown in Figure 16.
In these figures, we find that the truncated model ( Figure 16) exhibits a similar set of traveling bump solutions as the full phase model ( Figure  Figure 16: Existence of traveling bump solutions using the truncated interaction function H F i (q = 0, b = 0.8). Left: After a first critical value (g * ) of the bifurcation parameter g (red plane), there exist traveling bumps in the axial directions. After a second critical value (g * * ) (marked by a green plane), off-diagonal solutions form and continue to persist for large g. The dark to light color gradient and thin to thick thickness gradient corresponds to increasing values of g. Right: The projection of the branches on the left onto the g = 4 plane. A given point on one of these branches marks the magnitude and direction of a traveling bump. If necessary, one can approximate the parameter value g of this traveling bump by looking at the thickness and color of the chosen point and looking back at the branches in the left panel. 15). In particular, each system at a critical value g * , bifurcates into two axial solutions and one diagonal solution. For larger g, the system bifurcates again at another critical value g * * , giving rise to two non-axial, non-diagonal constant velocity directions. Indeed, negative velocity solutions exist, but as these solutions are symmetric up to multiples of a 90-degree rotation about the g-axis, we only show the positive directions. The "mixed" solutions that branch off for g > g * * are, in general, not rationally related so that the resulting traveling bumps will densely cover the torus. As such quasiperiodic solutions are often not structurally stable, we expect to see complex and possibly chaotic behavior when q > 0. Indeed, looking at Figure 13, we see that most of the complex behavior occurs for a small value of q and g sufficiently large.
Now that we have shown existence of traveling bump solutions, we proceed with a stability analysis.

Stability
We begin this section with stability of traveling bump solutions in the axial directions. We perturb off the axial solution, θ 1 (τ ) = ντ +εe λ 1 τ and θ 2 (τ ) = 0 + εe λ 2 τ , with Re(λ i ) > −1. The first order terms yield two independent eigenvalue problems which we combine with Equation (27) to yield two independent eigenvalue equations, Using these equations, we may determine stability of a traveling bump solution as a function of its velocity. We rephrase this problem and consider the independent scalar valued functions For a given Λ i , the zero level curves in (ν, λ) space determine stability properties of traveling bump solutions. We begin the analysis of these equations using the accurate Fourier series of H i and compute the integrals explicitly. The zero level set of the resulting function is shown in Figure 17. On the left panel, find that for any velocity, the x-direction is always stable. On the right, we find that for sufficiently small velocities, θ 2 (τ ) = 0 is a stable solution. Thus, constant velocity traveling solutions in this parameter regime will converge to the x-axis. Finally, for greater traveling bump velocities, the vertical direction loses stability, giving rise to non-axial solutions.
With the truncated H i , which we recall to be H F i (θ 1 , θ 2 ) = sin(θ 1 )(1 + b cos(θ 2 )), we may compute the equations Λ i = 0 explicitly as polynomials, where The coefficients of Equation (31) determine the stability of the horizontal traveling solution. Note that for ν sufficiently small, all coefficients are positive and the product c 1 c 2 dominates the coefficient c 0 . Thus, for small velocities, the coefficients have the properties c 1 c 2 > c 0 and c 2 , c 0 > 0, which implies stability by the Routh Hurwitz criterion. When ν * = ± √ 2b − 1, the coefficient c 0 is no longer positive and the stability condition fails.
We have found that horizontal traveling bump solutions lose stability at some critical velocity ν * , giving rise to non-axial traveling bump solutions. By symmetry, this argument holds for vertical traveling bump solutions: after the same critical ν * , constant velocity traveling bumps in the vertical direction lose stability and become non-axial solutions. Now that we understand the existence of and stability of axial constant velocity traveling solutions, we turn our attention to the stability of nonaxial traveling bump solutions.
To determine the stability of non-axial directions, we consider the solution, θ 1 (τ ) = ντ + φ 1 e λτ and θ 2 (τ ) = ντ + φ 2 e λτ . This ansatz results in the equations, By rewriting the integrals in the more compact form, − 1)ds, the problem reduces to finding an eigenvector (φ 1 , φ 2 ) T with corresponding eigenvalue −λ: This condition holds if and only if the determinant is zero. This determinant is the Evans function, and we use its roots to determine stability properties of the constant velocity solutions. Using the accurate Fourier series of H i , the integrals of the eigenvalue problem (32) are explicitly computable. Given a value g, it is straightforward to compute the contours of E = 0 using a standard contour plot routine. In panels A and B of Figure 18, we show the Evans function when g = 1.5, and g = 3, respectively.
The domain values where the real part of the Evans function is zero is shown as a black contour, while the domain values where the imaginary part is zero is shown in gray. Intersections of these contours show roots of the Evans function. Generally, there exists a root of the Evans function at the origin due to translation invariance of the underlying bump solution. Thus we ignore this root and consider only those nontrivial roots with real part sufficiently greater than −1. These nontrivial roots are marked with red dots.
We follow these roots using XPPAUTO and generate the bifurcation diagram shown in the right panel of Figure 18. The real part of the root remains negative for the range of g that we consider, thus the constant traveling bump solution remains stable for a large range of adaptation strengths.
We repeat the analysis of the Evans function using the truncated Fourier interaction function, H F . Once again, the integrals of the eigenvalue problem (32) are explicitly computable and we follow the roots of the Evans function using XPPAUTO in two parameters, b and g, the Fourier coefficient, and adaptation strength, respectively. The right panel of Figure 19 shows the result of this continuation: within the unstable region (marked in light blue), constant velocity solutions are unstable, as demonstrated by the lower inset showing θ 1 as a function of time. Because the instability arises through a Hopf bifurcation, the traveling bumps begin to "wobble". In the stable region, bump solutions travel with constant velocity, as demonstrated by the upper inset showing θ 1 as a function of time.
The left panels (A and B) of Figure 19 demonstrates the existence of a Hopf bifurcation on the boundary between stable and unstable regions. These panels correspond to points labeled A and B on the right panel. In each case, we find a complex conjugate pair of eigenvalues that cross the imaginary axis.
In this section, we analyzed the reduced neural field model with nonzero adaptation strength (g > 0) and no input current (q = 0). We now explore the dynamics arising from activating the time-invariant input current.

Hopf Bifurcation on the Torus
We have seen in Figures 13, 10 that for nonzero g and q, the system may produce a traveling bump solution that oscillates about the origin-centered input. For a fixed parameter value q, the origin is stable for g = 0, and with increasing g eventually becomes unstable through a Hopf bifurcation.
We study this phenomenon with Equation (22) Evaluating the integral and solving for λ yields Thus, as in the case of the ring, for a fixed q and given g sufficiently large, there exists a Hopf bifurcation at the critical value For the truncated interaction function H F i , the critical value is

Non-Constant Velocity Bump Solution on the Torus
As we have seen in earlier sections, stable oscillating solutions exist for particular choices of input current strength and adaptation on both the ring and torus. The similarities of solutions on the ring and torus continue as adaptation strength increases. On the ring, the oscillating solution gives way to a bump solution that travels around the ring with non-constant velocity. Similarly, with sufficiently large adaptation g, the bump solution on the torus also breaks free from the oscillating solution and traverses the domain with non-constant velocity. Figure 21 shows  There are plenty of other examples of these types of solutions (Figures 12,11,10) that are in fact chaotic. To demonstrate the existence of chaos numerically, we use the truncated phase model and a Poincaré section through cy = 0, as we find that generically the variable cy consistently crosses zero throughout simulations.

Chaos on the Torus
For a given g, we simulate the truncated phase model (Equation (26)) for t = 150000 time units and ignore the first 7000 time units to remove transients. By plotting the appropriate state variables, we are able to determine whether a system is aperiodic (and possibly chaotic) or periodic. The top panel of Figure 22 shows one example of one such plot, where for each g value we plot all cy values for the duration of the simulation. The black regions of Figure  22 correspond to the gray regions of Figure 12: the approximate range 0.85 < g < 1.1 corresponds to region F, the approximate range 1.18 < g < 1.61 corresponds to region G, and the approximate range g > 2.05 corresponds to region H, respectively.

Discussion
Our motivation for this work was to understand the behavior of the model presented in [19] where the authors showed that heterogeneities in a recurrent network with adaptation produced a seemingly randomly moving bump of activity. Similar moving bump dynamics was also found in a homogeneous bump model with adaptation in [14]; here the authors report only axially moving bumps with no external inputs. The neural field model considered in this paper (Equations (1),(2)) is also capable of producing a rich variety of solutions on the ring and torus. On the ring, the centroid of the bump solution exhibits sloshing and large-sloshing behaviors (for moderate strengths of input current and adaptation) that, for stronger adaptation, lead to non-constant velocity traveling bump solutions. With no input current and sufficient adaptation, the system generates a constant velocity traveling bump solution. We also observe chaotic solutions for a narrow range of adaptation strengths.
On the torus, the qualitative solutions are similar to those on the ring: Stationary bump solutions give rise to sloshing solutions (for moderate strengths of input current and adaptation), as well as non-constant velocity traveling bump solutions (with sufficient adaptation) and constant velocity traveling solutions (with no input current and sufficient adaptation). In this system, we do not see pulses that change in diameter with a fixed centroid (breathers) on the ring or torus.
Neural fields with nonsmooth firing rate functions (i.e., the Heaviside or rectifying nonlinearity) reproduce many of these qualitative behaviors. The existence and bifurcation of sloshing solutions on the ring are analyzed in [9,6], and constant velocity solutions are shown to exist on the ring ( [31,6]), the real line ( [24]), and the plane ( [15]). Nonconstant velocity bump solutions are shown to exist in [6]. However, there are no studies showing the existence of aperiodic attractors on the torus (assuming a deterministic system with an even kernel), or the existence of chaos on the ring.
In this study, we contributed to the analysis of the known behaviors by using a smooth firing rate function and a caveat of weak and slow adaptation. This assumption on the adaptation variable allowed us to reduce the neural field model, which is a distributed partial integro-differential equation, to a system of scalar delay integro-differential equations describing the centroid of the bump solution. Moreover, our only restriction on the kernel is of Mexican-hat type. Put together, these assumptions and our results are more general than what currently exists in literature.
In one spatial dimension, for example, we derived the normal form for the Hopf bifurcation in the one-dimensional neural field model and determined the conditions for super-and sub-criticality. Although normal form calculations exist for neural field models on the ring or the real line, our calculation allows for a general choice of kernel and a smooth firing rate function (as opposed to a particular choice of kernel or a non-smooth Heaviside firing rate function [6,9]).
As mentioned previously, existing studies require particular choices of kernels or Heaviside firing rate functions where a smooth firing rate function would be desirable ( [2] note that a smooth firing rate function allows for a straightforward normal form analysis). Although these assumptions are restrictive, these studies have advantages that the current study does not address. In particular, our analysis requires that adaptation is weak and slow and that the input current is weak. As a result, we can only study phenomena that evolve on a slow timescale. These weak and slow assumptions are well-suited for studying long-lasting sequences of spatially coherent activity in the absence of changing external stimuli [27,19], but may not be as well suited to study phenomena on a faster timescale, like the effects of weak modulatory interactions mediated by the reciprocal, long-range patchy connections in primary visual cortex [9].
Generally speaking, one might ask why we need the adaptation to be both slow and weak. For example, in [28], the adaptation was slow but not weak. One could imagine doing a perturbation analysis for a weak stimulus such as in [7] where a weak slowly moving stimulus is applied to a system that has a stable traveling bump. However, strong adaptation, will alway induce movement in a bump so that we can never pin the bump with weak inputs. Furthermore, by keeping the adaptation O(1), one needs to compute the adjoint solution to a two-variable traveling bump, a difficult task in one spatial dimension, and impossible (as far as we can tell) in two spatial dimensions. Thus, by working with weak inputs and weak/slow adaptation, we have hit a sweet spot from which many of the interesting dynamics emerges.
One type of behavior that has been observed in this class of models that does not occur in our analysis is the so-called breathing solutions [11]. Breathers are periodic solutions to the neural field equations that occur when the bump solution loses stability via a symmetric mode. In contrast, sloshers appear when there is a Hopf bifurcation to an anti-symmetric mode. In the breather case, the centroid of the bump does not change, so our reduced equations cannot detect such a bifurcation. In contrast, sloshers lead to modulation of the centroid and thus our analysis can capture that. To further explore this, we were able to induce bifurcation to a breathing solution in equations (1,2) but only when ε is sufficiently large. We find that it is possible to continue this bifurcation in ε and make ε quite small, but only if we increase both the strength of adaptation g and the heterogeneity, q such that εq, εg remain O(1). That is, breathers can only occur when the adaptation and input magnitudes are large compared to the rate of adaptation. Our analysis, therefore, cannot include the appearance of breathers.
The effects of noise on the phase equations is one possible direction for future study. Several cited papers analyze the movement of bump solutions in the presence of noise. Near the drift bifurcation for traveling bump solutions with sufficiently strong linear adaptation, it is possible to derive a stochastic amplitude equation when the adaptation strength and stochastic forcing are similar in magnitude [22]. Sufficiently far from the bifurcation, the stochastic forcing leads to diffusive wandering of the bump solution. Other studies analyze the diffusive behavior of solutions to neural field models and how pinning eliminates diffusive behavior [30]. In [25], the authors compute the effects of adding noise to the normal form of the pitchfork bifurcation. As mentioned previously, these studies assume either a particular firing rate function or kernel. The general case remains unexplored.

A Normal Form for the Hopf Bifurcation on the Ring
Recall that we analyze the bifurcation to sloshing pulses for the general integral equation, as g increases. For simplicity, we assume the expansions and q > 0. Based on the eigenvalue equation, we expect a Hopf bifurcation to occur. To analyze the Hopf bifurcation, we use a multiple time scale expansion. We assume that θ(τ ) is a function of a "fast" time ζ = τ and a "slow" time ξ = δ 2 τ where δ measures the amplitude of the bifurcating solution. As the nonlinearities are all odd, we can assume that g = g 0 + δ 2 g 2 , θ = δθ 1 (ζ, ξ) + δ 3 θ 3 (ζ, ξ) to order δ 3 . We develop a perturbation expansion to obtain the normal form. We expand this expression to order δ 2 to get:  We plug in all the expansions and find to first order that where z(ξ) is a complex function of ξ and c.c means complex conjugates. Our goal is to derive equations for z. To cubic order, we obtain: (Lθ 3 )(ζ) = z ξ e iωζ −1 + g 0 1 + 2iω − ω 2 + c.c where α = 1 − g 0 1 + 2iω − ω 2 = 2 1 + q (q + √ qi) γ 0 = g 2 iω 1 + iω = g 2 1 + q (q + √ qi)

B Computation of Functions H i and J i
To numerically integrate the phase models on the ring or torus, we require an approximation to the functions H i , and J i . These functions depend on and use lookup tables for the steady state bump u 0 (u0ss), the derivative of the firing rate evaluated at the steady state bump f (u 0 ) (df u0b), and the partial derivatives of the steady state bump, ∂u 0 /∂x, ∂u 0 /∂y (ux,uy). On the toroidal domain, each lookup table has N × N entries, where for the coefficients below, we choose N = 64.
To compute H i in Equation (9),we use the following procedure error, revealing that J i = −H i . Thus, for all phase computations involving J i , we use the same Fourier approximations for H i and J i . On the ring, the computations are virtually identical with the obvious exception of array shapes.

B.1 Fourier Approximations
After creating the lookup tables H1,H2, we perform a Fourier approximation to make numerical integration easier. The following function and corresponding coefficients and frequencies (Table 1) provide an excellent approximation to the lookup tables H1,H2. A basic error analysis shows that the supremum norm difference between the lookup tables H1,H2 and their Fourier approximations,H 1 , H 2 , is H1 − H 1 ∞ ≈ 3.354e−7.
H 1 (x, y) = − 26 k=1 a k N 2 sin(xn k + ym k ). (37) The coefficients in Table 1 are computed using Python with Numpy by taking the Fourier transform of the lookup tables H1,H2.

C Numerical Integration
In this section, we detail the various numerical methods used to evaluate the many integro-delay-differential equations of this paper.

C.1 Truncated Neural Field Model on the Torus
The integration of Equation (21) requires the approximation of several double integrals. In the interest of reducing computation time, we use Riemann integrals and a relatively coarse discretization of the spatial domain. For example, for a given time t, the coefficient p 10 (t) is approximated as cos(y m )f (u(x n , y m , t)) (2π) 2 N 2 .
Because a linear increase in N leads to a quadratic increase in the total number of operations, we keep N = 100, which is an acceptable compromise between speed and accuracy for this problem. All other double sums that appear in p ij , r ij , and s i are computed this way. When computing the bifurcation diagram using this system, we use XPPAUTO and the numerical options shown in Table 4. The most important options are Ntst and Dsmin. If Ntst is less than 1000, XPPAUTO is unreliable in determining the stability of periodic solutions. If Dsmin is too large, XPPAUTO will skip bifurcation points. 1.30119169839e-07 (-3,1) 10 0.00255677958311 Table 2: Fourier Coefficients of the steady-state coefficients. Plotting u 00 + 2u 10 cos(x) + 2u 01 cos(y) + 4u 11 cos(x) cos(y) gives a reasonable approximation to the numerically computed steady-state bump solution.  Table 3: Fourier Coefficients of the kernel. Plotting k 00 + 2k 10 cos(x) + 2k 01 cos(y) + 4k 11 cos(x) cos(y) gives a reasonable approximation to the original periodix kernel.

C.2 Delay Integro-Differential Equations
We implement the right hand side of the integro-differential in Equation The initial condition for this algorithm requires an array of M time steps. If the parameters are chosen such that a limit cycle exists, then we initialize in an arc: When plotting solutions, we disregard at least the first M entries of the solution vector.