Ghosts of bump attractors in stochastic neural fields: Bottlenecks and extinction

We study the effects of additive noise on stationary bump solutions to spatially extended neural fields near a saddle-node bifurcation. The integral terms of these evolution equations have a weight kernel describing synaptic interactions between neurons at different locations of the network. Excited regions of the neural field correspond to parts of the domain whose fraction of active neurons exceeds a sharp threshold of a firing rate nonlinearity. For sufficiently low firing threshold, a stable bump coexists with an unstable bump and a homogeneous quiescent state. As the threshold is increased, the stable and unstable branch of bump solutions annihilate in a saddle node bifurcation. Near this criticality, we derive a quadratic amplitude equation that describes the slow evolution of the even mode (bump contractions) as it depends on the distance from the bifurcation. Beyond the bifurcation, bumps eventually become extinct, and the time it takes for this to occur increases for systems nearer the bifurcation. When noise is incorporated, a stochastic amplitude equation for the even mode can be derived, which can be analyzed to reveal bump extinction time both below and above the saddle-node.

1. Introduction. Continuum neural fields are a well-accepted model of spatiotemporal neuronal activity evolving within in vitro and in vivo brain tissue [6,11]. Wilson and Cowan initially introduced these nonlocal integrodifferential equations to model activity of neuronal populations in terms of mean firing rates [54]. While they discount the intricate dynamics of neuronal spiking, these models can qualitatively capture a wide range of phenomena such as propagating activity waves observed in disinhibited slice preparations [27,44,45,47]. Neural field models exhibit a wide variety of spatiotemporal dynamics including traveling waves, Turing patterns, stationary pulses, breathers, and spiral waves [14,17,19,31,38]. A distinct advantage of utilizing these continuum equations to model large-scale neural activity is that many analytical methods for studying their behavior can be adapted from nonlinear partial differential equations (PDEs) [6]. Recently, several authors have explored the impact of stochasticity on spatiotemporal patterns in neural fields [8,30,35] by employing techniques originally used to study stochastic front propagation in reaction-diffusion systems [48]. Typically, the approach is to perturb about a linearly stable solution of the deterministic system, under the assumption of weak noise. However, some recent efforts have been aimed at understanding the impact of noise on patterns near bifurcations [30,36].
In this work, we are particularly interested in how noise interacts with stationary pulse (bump) solutions near a saddle-node bifurcation at which a branch of stable bumps and a branch of unstable bumps annihilate [1]. Bumps are commonly utilized as a model of persistent and tuned neural activity underlying spatial working memory [21,55]. This activity tends to last for a few seconds, after which it is extinguished, to allow for subsequent memories to be formed [23]. One possible way to terminate these sustained activity patterns is by transiently synchronizing the spiking patterns of excitatory neurons that participate in the signal [25]. Another proposed mechanism for terminating persistent activity is a strong and brief global inhibitory signal, which would drive the system from the stable bump state to a stable uniform quiescent state [10]. In terms of neural field and spiking models, this can be thought of as momentarily raising the firing threshold of the system, temporarily driving it beyond the saddle-node bifurcation from which the stable bump emerges.
We focus on a scalar neural field model that supports stationary bump solutions for appropriate choices of parameters and constituent functions [1,11]: where u(x, t) is the total synaptic input arriving to location x and time t, and w(x−y) describes the strength (amplitude) and polarity (sign) of synaptic connections from neurons at location y to neurons at location x. We assume w(x) is an even-symmetric function w(x) = w(−x) with a bounded integral Ω w(x)dx over the spatial domain x ∈ Ω = (−x ∞ , x ∞ ). The nonlinearity f (u) is a firing rate function, which we take to be the sigmoid [54] f (u) = 1 1 + e −η(u−θ) , (1.2) and we also find it useful to take the high gain limit η → ∞, in which case: allowing for analytical tractability in several of our calculations. It is important to note that (1.1) neglects several known features of neuronal networks including spike rate adaptation [26], propagation delays [29], synaptic depression [33], and refractoriness [16]. Thus, we assume we are focusing on a network where these effects are weak enough as to not impact our main results.
Amari was the first to analyze (1.1) in detail, showing that when f (u) is defined to be a Heaviside function (1.3), the network supports stable stationary bump solutions when the weight function w(x) is a lateral inhibitory (Mexican hat) distribution satisfying: x 0 ]; and (iv) w(x) has a unique minimum on [0, x ∞ ) at x = x 1 with x 1 > x 0 and w(x) strictly increasing on (x 1 , x ∞ ) [1]. Based on restrictions (i)-(iv), Amari made use of the integral of the weight function to prove some of the main results of his seminal work. For instance, it is clear that W (0) = 0 and W (x) = −W (−x) based on the above assumptions. Moreover, there will be a single maximum of the function W (x) on the interval (0, x ∞ ) given at , due to conditions (i) and (ii), and w(x 0 ) = 0. When θ < W (x 0 ) there are two bump solutions: one stable and one unstable (up to translation symmetry), and when θ > W (x 0 ) there are no bump solutions to (1.1). When θ = θ c ≡ W (x 0 ), there is a single marginally stable bump solution. It is at this point that the two branches (stable and unstable) of bump solutions meet and annihilate in a saddle-node bifurcation ( Fig. 2.1). Dynamics of (1.1) for values of θ beyond this saddle-node bifurcation evolve to quasi-stationary solutions resembling the ghost of the bump at θ c , lasting for a period of time inversely related to |θ − θ c | [51]. A principled exploration of these dynamics (section 2) is one of the primary goals of this paper. When θ is below the critical threshold θc at the saddle-node, there are two stationary bump solutions to (1.1): one stable as and one unstable au. When θ > θc, there are zero equilibria, but the dynamics of (1.1) are slow in the bottleneck near Uc(x).
As mentioned, the neural field equation (1.1) in the absence of noise has been analyzed extensively [1,11,17]. We expand upon these previous studies by also exploring the impact of noise on stationary bump solutions to (1.1) near a saddle-node bifurcations (section 3). Additive noise is incorporated, so that the evolution of the neural field is now described by the spatially extended Langevin equation [4,6,30,40]: where the term dW (x, t) is the increment of a spatially varying Wiener process with mean defined by dW (x, t) = 0 and correlations dW (x, t)dW (y, s) = C(x − y)δ(t − s)dtds and describes the amplitude of the noise, assumed to be weak ( 1). The function C(x − y) describes the spatial correlation in each noise increment between two points x, y ∈ Ω.
2. Slow bump extinction in the deterministic system. We begin by examining the dynamics of stationary bump solutions near a saddle-node bifurcation, where a stable and unstable branch of solutions annihilate. Our initial analysis focuses on the noise-free case W (x, t) ≡ 0, allowing us to derive an amplitude equation that approximates the evolution of the bump height. Linearization of bumps in (1.1) typically reveals that they are marginally stable to translating perturbations, so the overall stability is characterized by the stability to even perturbations that expand/contract the bump [17]. Our analysis will emphasize the region of parameter space near where bumps are marginally stable to even perturbations.
2.1. Existence and stability of bumps. We now briefly review existence and stability results for stationary bump solutions to the neural field equation (1.1). These results are analogous to those presented in [1,35,52]. For transparency, we focus on the case of a Heaviside firing rate function (1.3). This allows us to cast bump stability in terms of a finite dimensional set of equations, focusing on the evolution of the two edge interfaces of the bump [1,15]. Assuming a stationary solution u(x, t) = U (x), we find (1.1) requires Given a unimodal bump solution U (x), without loss of generality, we can fix the center and peak of the bump to be at the origin x = 0. In the case of even-symmetric bumps U (x) = U (−x) [1], we will have the conditions for the bump half-width a: U (x) > θ for x ∈ (−a, a), U (x) < θ for x ∈ Ω\[−a, a], and U (±a) = θ. In this case, (2.1) becomes By utilizing the integral function (1.4), we can write the even-symmetric solution To determine the half-width a, we require the threshold conditions U (±a) = θ of the solution (2.2) to yield Note that when θ < W max = max x W (x), there will be a stable and unstable bump solution to (1.1). When θ = θ c ≡ W max , there is a single marginally stable bump solution U c (x) to (1.1), as illustrated in Fig. 2.1B. Differentiating W (2a) by its argument yields W (2a c ) = w(2a c ) ≡ 0 as an implicit equation for the half-width a c at this criticality. Utilizing the notation of Amari condition (i), we have that a c = x 0 /2. Note, the relation w(2a c ) = 0 is explicitly solvable for a c for several typical lateral inhibitory type weight functions. For instance, in the case of the difference of Gaussians For the "wizard hat" w(x) = (1 − |x|)e −|x| on x ∈ (−∞, ∞) [12], we have a c = 1/2 and θ c = e −1 . For a cosine weight w(x) = cos(x) on the periodic domain x ∈ [−π, π] [35], we have a c = π/4 and θ c = 1.
To characterize the stability of bump solutions to (1.1), we will study the evolution of small smooth perturbations εψ(x, t) (ε 1) to stationary bumps U (x) by utilizing the Taylor expansion u(x, t) = U (x) + εψ(x, t) + O(ε 2 ). By plugging this expansion into (1.1) and truncating to O(ε), we can derive an equation whose solutions constitute the family of eigenfunctions associated with the linearization of (1.1) about the bump solution U (x). We begin by truncating (1.1) to O(ε) assuming u is given by the above expansion and that the nonlinearity f (u) is given by the Heaviside function (1.3), so and we can differentiate the Heaviside function, in the sense of distributions, by noting which we can rearrange to find One class of solutions, such thatψ(±a, t) = ψ(±a, 0) = 0, lies in the essential spectrum of the linear operator that defines (2.5).
In this case,ψ(x, t) =ψ(x, 0)e −t , so perturbations of this type do not contribute to any instabilities of the stationary bump U (x) [24]. Assuming separable solutions ψ(x, t) = b(t)ψ(x), we can characterize the remaining solutions to (2.5). In this case, Solutions to (2.6) that do not satisfy the condition ψ(±a) ≡ 0 can be separated into two classes: (i) odd ψ(a) = −ψ(−a) and (ii) even ψ(a) = ψ(−a). This is due to the fact that the equation (2.6) implies the function ψ(x) is fully specified by its values at x = ±a. Thus, we need only concern ourselves with these two points, yielding the two-dimensional linear system For odd solutions ψ(a) = −ψ(−a), the eigenvalue reflecting the fact that (1.1) is translationally symmetric, so bumps are marginally stable to perturbations that translate their position. Even solutions ψ(a) = ψ(−a) have associated eigenvalue Thus, when θ < θ c , the wide bump a s > a c will be linearly stable to expanding/contracting perturbations since w(2a s ) < 0 due to Amari's condition (ii) [1]. The narrow bump a u < a c is linearly unstable to such perturbations since w(2a u ) > 0 due to condition (i). When θ = θ c , we have w(2a c ) = 0 so that λ e = 0 and |U (±a c )| = w(0). In anticipation of our derivations of amplitude equations, we define the eigenfunctions at the criticality θ = θ c . Utilizing the fact that |U (±a c )| = w(0) and the linear system (2.7a), we have that the odd eigenfunction at the bifurcation is 8) and the even eigenfunction is Note, this specifies that ψ e (±a) = ψ o (a) = −ψ 0 (−a) = 1. Furthermore, we will find it useful to compute the derivatives which is odd (ψ e (−a c ) = −ψ e (a c )). Lastly, we note that we will utilize the fact that, for even symmetric functions,

2.2.
Saddle-node bifurcation of bumps. Motivated by the above linear stability analysis, we now carry out a nonlinear analysis in the vicinity of the saddle-node bifurcation from which the stable and unstable branches of stationary bumps emanate. Specifically, we will perform a perturbation expansion about the bump solution U c (x) at the critical threshold value θ c . We therefore define θ = θ c + µε 2 , ε 1, so that µ is a bifurcation parameter determining the distance of θ from the saddle-node bifurcation point. As demonstrated above, the linear stability problem for U c (x) reveals two zero eigenvalues λ o = λ e = 0 associated with the odd ψ o and even ψ e eigenfunctions (2.8) and (2.9), respectively. Our analysis employs the ansatz: where τ = εt is a temporal rescaling that reflects the vicinity of the system to a saddle-node bifurcation associated with the even expanding/contracting eigenmode ψ e [51]. Similar expansions have been utilized in the analysis of bifurcations for spatial patterns in reaction-diffusion systems [3,49] and neural field models [5,28,53]. Upon plugging (2.10) into (1.1) and expanding in orders of ε, we find that at O(1), we simply have the stationary bump equation (2.1) at θ = θ c . Proceeding to O(ε), we find The right hand side of (2.11) vanishes due to the formula for the even (2.9) eigenfunction associated with the stability of the bump U c (x). At O(ε 2 ), we obtain an equation for higher order term u 2 : where L is the non-self-adjoint linear operator Both ψ o (x) and ψ e (x) lie in the nullspace N (L), as demonstrated in the previous section by identifying solutions to (2.3). Thus, the ψ o terms on the left hand side of (2.12) vanish. We can ensure a bounded solution to (2.12) exists by requiring that the right hand side be orthogonal to all elements of the nullspace of the adjoint operator L * . The adjoint is defined with respect to the L 2 inner product (2.14) Thus, we find defined in the sense of distributions under the L 2 inner product given in (2.14). It is straightforward to show that To show ϕ o , ϕ e ∈ N (L * ), we simply plug these formulas into (2.16) to find for j = o, e, which is true due to the fact that ψ o and ψ e lie in N (L). Thus, we will impose solvability of (2.12) by taking the inner product of both sides of the equation with respect to for j = o, e, where we have defined the convolution w * F = Ω w(x−y)F (y)dy. Due to odd-symmetry, terms of the form H (U c − θ c )ψ j , ψ k , j = k, vanish. In a similar way, Isolating the temporal derivatives A j in (2.17), we find that the amplitudes A j (j = o, e) satisfy the following fast-slow system of nonlinear differential equations (2.18b) With the system (2.18) in hand, we can determine the long term dynamics of the amplitudes as the bifurcation parameter µ is varied. We begin by computing the constituent components of the right hand sides, using properties of the eigenfunctions ψ o and ψ e . To start, we will compute the second derivative H (U c −θ c ), which appears in the coefficient of the quadratic term A 2 e . Differentiating the function H(U c (x) − θ c ) twice with respect to x, using the chain and product rule, we find the following formula where we have applied the identity (2.4) for the first derivative H (U −θ). Rearranging terms, we find that We can further specify the formula (2.19) by differentiating dH(Uc−θc) where δ (x − x 0 ) is defined, in the sense of distributions, for any smooth function F (x) by using integration-by-parts [32]: . Thus, we can at last write Computing the inner products in (2.18) then simply amounts to evaluating the integrals in the sense of distributions. First, we use (2.4) to note where we have utilized ψ e (±a c ) = 1 and w(2a c ) ≡ 0. Finally, we compute the quadratic terms using the identity (2.20), starting with and we note that individual terms under the integral from the sum defining (2.20) are for the terms involving the distributional derivative δ (x − x 0 ), whereas the terms The integrals in (2.25) are identical to those in (2.22), so it is straightforward to compute, using (2.23) and (2.24) that Thus, we can at last compute all the terms in (2.18), specifying that where we have noted the fact that w (2a c ) < 0 due to Amari's conditions (iii) and (iv) on the weight function w(x) [1]. Equation (2.26a) reflects the translational symmetry of the original neural field equation (1.1), so bumps are neutrally stable to translating perturbations ψ o regardless of the bifurcation parameter µ. On the other hand, as the bifurcation parameter µ is changed, the dynamics of the even eigenmode ψ e reflect the relative distance to the saddle-node bifurcation at which point bumps are marginally stable to expanding/contracting perturbations. When µ < 0, there are two fixed points of equation (2.26b) at A e = ±w(0) |µ/w (2a c )|, corresponding to the pair of emerging stationary bump solutions which are wider (+) and narrower (−) than the critical bump U c . As expected, the wide bump is linearly stable since a linearization of (2.26b) yields λ + = − |µ · w (2a c )|/w(0) < 0, and the narrow bump is linearly unstable since λ − = + |µ · w (2a c )|/w(0) > 0. Crossing through the subcritical saddle-node bifurcation, we find that for µ ≡ 0, there is a single fixed point A e ≡ 0, which is marginally stable, since λ 0 = 0.
Lastly, note when µ > 0, there are no fixed points of the differential equation (2.26b). However, starting at the initial condition A e (0) = 0 (correspondingly u(x, 0) = U c (x)), we find that the dynamics of the amplitude A e (τ ) are strongly determined by the ghost of the fixed point at A e = 0 [51]. Note in Fig. 2.2A that the transient bump retains a shape much like that of the critical bump for an appreciable period of time before extinguishing. Trajectories of the full system (1.1) evolve more slowly when the distance to the bifurcation |θ − θ c | = |µ|ε 2 is smaller. Solving for A e (τ ) in this specific case and reverting the the original time coordinate t = τ /ε, we find Thus, the residence time t b in the bottleneck, or neighborhood of the ghost of the fixed point A e = 0, is given by the amount of time it takes for A e (t) to traverse to some set value. Of course, this is dependent on the bifurcation parameter µ. For illustration, we examine how long it takes until A e (t b ) = −1. By explicitly focusing on the region where |A e (t b )| ≤ 1, we are roughly restricting to the time interval during which |u(x, t) − U c (x)| = O(ε), where we would expect the expansion (2.10) to be valid. Using the formula (2.27), it is straightforward to find that We compare this formula to the results of numerical simulations in Fig. 2.2B, utilizing the difference of Gaussians weight function w(x) = e −x 2 − Ae −x 2 /σ 2 on x ∈ (−∞, ∞).
Comparisons are made by noting that when A e (t b ) = −1, then u(x, t) ≈ U c (x) − εψ e (x), so that the peak of the activity profile will be   Fig. 2.2C,D that, as predicted, the time spent in the bottleneck increases as the amplitude of the small parameter ε is decreased. The attracting impact of the ghost is stronger when the parameters of the system lie closer to the bottleneck. For further comparison, we consider the case w(x) = cos(x) in Fig. 2.3. In this case the constituent functions a c = π/4, w(0) = 1, and w(2a c ) = −1. Furthermore, by setting µ = 1 the formulas for the amplitude (2.27) and residence time (2.28) simplify considerably to A e (t) = − tan(εt) and t b = π/[4ε].

Amplitude equations for smooth nonlinearities.
Our nonlinear analysis in the case of Heaviside nonlinearities f (u) ≡ H(u − θ) made extensive use of the specific form of the distributional derivatives. Inner products with these functions lead to dynamical equations focused on a finite number of discrete points in space, rather than over the spatial continuum x ∈ Ω. Here, we show it is straightforward to extend this analysis to the case of arbitrary smooth nonlinearities f (u). There are several detailed analyses of stationary bumps in neural field with smooth firing rate, showing a similar bifurcation structure to that presented in Fig. 2.1: a stable and an unstable branch of bump solutions annihilate in a saddle-node bifurcation as the threshold of the firing rate function is increased. We refrain from such a detailed analysis here and refer the reader to these works [13,18,35,37,41,52]. Again, defining θ = θ c + µε 2 , ε 1, so µ determines the distance of θ from the bifurcation and on which side of θ c it lies. Following our previous analysis, we utilize the ansatz (2.10) and rescale time τ = ετ . In this case, ψ o (x) and ψ e (x) will still be odd and even eigenmodes associated with the linear stability of stationary bump solutions to (1.1). At the criticality θ ≡ θ c , their associated eigenvalues will be λ o = λ e ≡ 0, as in the case of Heaviside firing rates [52]. Expanding (1.1) in orders of ε using the ansatz (2.10) yields a similar amplitude equation to (2.18) at O(ε 2 ). Again, we apply solvability conditions to the equation for u 2 . After canceling odd terms and isolating the derivatives A j , we find the amplitudes A j satisfy the system: We can derive the coefficients in the system (2.29) by computing the inner products therein. To do so, we must choose a specific nonlinearity, such as the sigmoid (1.2), and a weight kernel. For illustration, we consider the cosine kernel w(x) = cos(x) on the ring x ∈ Ω = [−π, π] with periodic boundaries. As shown in previous studies, the bump solution U c (x) = A c cos x while the eigenmodes ψ o (x) = sin(x) and ψ e (x) = cos(x) [26,35,52]. Since Lψ j ≡ 0 for j = o, e, this means sin(x) = f (A c sin(y)) cos(y) 2 dy = 1. Furthermore, (2.32) Lastly, we can simplify the integrals in the quadratic term by again making use of the identity π −π w(x − y)ϕ e (y)dy = ψ e (x), so 3. Stochastic neural fields near the saddle-node. We now study the impact of stochastic forcing near the saddle-node bifurcation of bumps. Our analysis utilizes the spatially extended Langevin equation with additive noise (1.5). Guided by our analysis of the deterministic system (1.1), we will utilize an expansion in the small parameter ε, which determines the distance of the system from the saddle-node. To formally derive stochastic amplitude equations, we must specify the scaling of the noise amplitude as it relates to the small parameter ε, as this will determine the level of the perturbation hierarchy wherein the noise term dW will appear. We opt for the scaling = ε 5/2 , as this introduces a nontrivial interaction between the nonlinear amplitude equation for A e and the noise.
It is important to note that our derivations are only carried up to O(ε 2 ) in the hierarchy of the regular perturbation expansion in ε. Were we to continue this expansion further, we would likely find that the = ε 5/2 noise term does indeed shift the location of the bifurcation at higher order as in [2,30]. Thus, as the amplitude of noise is increased, the validity of the expansion we derive here will begin to break down, since the terms beyond O(ε 2 ) will have a more substantial effect on the dynamics. Hence, the results we derive in this section are valid for small noise levels only. An understanding of the effects of larger noise terms, employing scalings = ε p with p < 5/2, warrants further study which is beyond the scope of our current work.

Stochastic amplitude equation for bumps.
Motivated by our quantitative analysis in the noise-free case, we rescale time in the stochastic term of (1.5) using τ = εt, so u(y, t))dy dt + ε 2 dŴ (x, τ ), (3.1) where dŴ (x, τ ) := √ εdW (x, ε −1 τ ) is a rescaled version of the Wiener process dW that is independent of ε [22]. We then apply the ansatz (2.10) once again and take Heaviside firing rate functions (1.3), thus finding (2.11) at O(ε). The O(ε) equation is satisfied due to the fact that ψ e ∈ N (L), where L is the linear operator given by (2.13). Finally, proceeding to O(ε 2 ), we find As before, the ψ o terms on the left vanish since Lψ o ≡ 0, and we ensure a bounded solution to (3.2) exists by requiring the inhomogeneous part is orthogonal to ϕ o , ϕ e ∈ N (L * ), where L * is the adjoint linear operator given by (2.15). Taking inner products yields for j = o, e. Isolating temporal derivatives, we find the amplitudes A o (t) and A e (τ ) obey the following pair of nonlinear stochastic differential equations Utilizing the formulas for H (U c − θ c ) (2.4) and H (U c − θ c ) (2.19) we derived in the previous section, we can simplify the expressions in (3.4). Additionally, we make use of the fact that  where we have converted the noise term in (3.5a) back to the original time coordinate: dW o (t) = dŴ o (εt)/ √ ε [22]. Note that in equation (3.5a), we essentially recover the diffusion approximation of the translating mode of the bump A o (t) 2 = εD o t, which is analyzed in [35]. Equation (3.5b) is a stochastic amplitude equation, so that the noise term dW e is projected onto the direction of the neutrally stable even perturbation ψ e .

Metastability and bump extinction.
To analyze the one-dimensional nonlinear SDE (3.5b), we further rescale the equation by setting A := |w (2ac)| w(0) 2 A e : where m := |w (2ac)| w(0) 2 µ. Thus, the effective diffusion coefficient of the rescaled noise term is Ŵ (τ ) 2 = Dτ = w (2a c ) 2 (C(0) + C(2a c ))τ / 2w(0) 4 . Note the rescaled equation (3.6) has an effective potential [42,51]: stochastic system (3.6) will eventually escape to the limit A → −∞ as τ → ∞. Such trajectories were observed in the noise-free system in the case m > 0, as demonstrated in Fig. 2.2 of the previous section. However, we show here that noise qualitatively alters the dynamics of the system, so its state will not remain in the vicinity of the stable attractor (at A = √ m) when m < 0. As before, we study the problem of bump extinction using the stochastic amplitude equation (3.6) in the case m > 0. We show that the noise decreases the average amount of time until an extinction event will occur. For clarity, we assume the initial condition A(0) = 0 (correspondingly u(x, 0) = U c (x)). We take the bottleneck to be the region A e ∈ [−1, 1], which in the rescaled variable is A ∈ [−|w (2a c )|/w(0) 2 , |w (2a c )|/w(0) 2 ]. The residence time τ b in the bottleneck is given by the amount of time it takes for A to escape this region. We can determine the statistics of τ b by considering it as a first passage time problem.
Let p(A, τ ) be the probability density for the stochastic process A(τ ) given the initial condition A(0) = A 0 . Then the corresponding Fokker-Planck equation is given For the non-generic case m = 0, the timescale of departure scales algebraically [50]. When m > 0, noise simply modulates the flows of the deterministic equationȦ = −m − A 2 , leading to an average speed-up in the departure from the bottleneck. In general, we consider solving the first passage time problem as an escape from the domain (−α, ∞) where α := |w (2ac)| w(0) 2 (equivalently where A e = −1) [22]. To do so, we impose an absorbing boundary condition at −α: p(−α, τ ) = 0. Now let T (A) denote the stochastic first passage time for which (3.6) first reaches the point −α, given it started at A ∈ (−α, ∞). The first passage time distribution is related to the survival probability that the system has not yet reached −α: which is S(τ ) := Pr(t > T (A)), so the first passage time density is [22] Substituting for the expression for ∂p/∂τ using the Fokker-Planck equation ( and V (x) is the potential function (3.7). Explicit expressions for the integral (3.11) can be found in some special cases [42,50]. For our purposes, we simply integrate (3.11) numerically to generate theoretical relationships between the mean first passage time and model parameters. For comparison, we focus on the case the weight function w(x) = cos(x) and the correlations C(x) = cos(x), so that U c (x) = √ 2 cos(x), a c = π 4 , w(0) = 1, w (2a c ) = −1, C(0) = 1, and C(2a c ) = 0. Therefore, α = 1, m = µ, D = 1/2. This allows us to write the formula (3.11) at A = 0 as Lastly, note that by rescaling time t = ετ , we have that the mean first passage time in units of t will bet b = T (0)/ε. We compare our theory (3.12) with the results of numerical simulations of the full stochastic neural field (1.5) in Fig. 3.2B. Note there is some discrepancy between our numerical simulations and theory as m is decreased. One of the primary reasons for this deviation is likely because of the moderate level of noise (ε = 0.6) used in comparison to the small parameter assumption (ε 1) using in the theory we have developed. Any minor mismatch will be exacerbated by the fact that mean first passage times for escape problems depend exponentially on parameters like noise amplitude and well depth, as in (3.12). Nonetheless, the theory does provide a rough estimate of the mean first passage times for smaller values of the parameter m.
4. Discussion. We have developed a weakly nonlinear analysis for saddle-node bifurcations of bumps in deterministic and stochastic neural field equations. While most of our analysis has focused upon Heaviside firing rate functions, we have also demonstrated the techniques can easily be extended to arbitrary smooth nonlinearities. In the vicinity of the saddle-node, the dynamics of bump expansion/contraction can be described by a quadratic amplitude equation. For deterministic neural fields, this low dimensional approximation can be used to approximate the trajectory and lifetime of bumps as they slowly extinguish. To do so, we focused on the initial time epoch in the bottleneck surrounding the ghost of the critical bump U c (x). In stochastic neural fields with appropriate noise scaling, a stochastic amplitude equation for the even mode of the bump can be derived. Importantly, we must choose the noise amplitude to scale as = ε 5/2 , in order for the noise term to appear in the stochastic version of the quadratic amplitude equation. We then cast the lifetime of the bump in terms of a mean first passage time problem of the reduced system, which is valid for the noise scaling we have chosen.
Our work extends a variety of recent studies that have derived low-dimensional nonlinear approximations of neural field pattern dynamics in the vicinity of bifurcations [5,7,20,30,35,36]. As in our work, most of these previous studies derived approximations where the location of the bifurcation was unaffected by noise terms. On the other hand, Hutt et al. showed that noise can in fact shift the position of Turing bifurcations in neural fields, and the amplitude of the bifurcation threshold shift was proportional to the noise variance [30]. Were we to have carried the hierarchy out to higher order, we would have found such a shift in the case we studied. Note, it was necessary in our work to apply a specific noise scaling (ε 5/2 ), as compared to the distance from criticality (ε 2 ), in order for the noise to simply appear as a modification of the even mode amplitude equation. Were we to have selected noise of larger amplitude, this could have induced bifurcation shifts at lower order, analogous to that found in [30]. Another potential future direction would be to consider the impact of axonal propagation delays [29] on the dynamics close to the saddle-node. As demonstrated in this work, the neural field (1.1) is quite sensitive to small perturbations near criticality, so delays may alter the duration of the bottleneck or even shift the saddle-node bifurcation point. In previous work [36], we derived amplitude equations describing propagation-generating drift bifurcations that arise when linear adaptation is incorporated into (1.1). We anticipate that similar analyses might be performed on networks with synaptic depression [33], although piecewise smooth methods may be necessary in the case of Heaviside firing rates [34]. Lastly, we note there have been recent efforts to systematically derive macroscopic descriptions of neural activity from recurrently coupled spiking networks [9,39,43]. Such alternative descriptions can also capture the form of steady state solutions like stationary bumps [39,46]. Bumps in these models do exhibit saddle-node bifurcations similar to those observed in the