Coarse-grained dynamics of an activity bump in a neural field model

We study a stochastic nonlocal PDE, arising in the context of modelling spatially distributed neural activity, which is capable of sustaining stationary and moving spatially-localized ``activity bumps''. This system is known to undergo a pitchfork bifurcation in bump speed as a parameter (the strength of adaptation) is changed; yet increasing the noise intensity effectively slowed the motion of the bump. Here we revisit the system from the point of view of describing the high-dimensional stochastic dynamics in terms of the effective dynamics of a single scalar"coarse"variable. We show that such a reduced description in the form of an effective Langevin equation characterized by a double-well potential is quantitatively successful. The effective potential can be extracted using short, appropriately-initialized bursts of direct simulation. We demonstrate this approach in terms of (a) an experience-based"intelligent"choice of the coarse observable and (b) an observable obtained through data-mining direct simulation results, using a diffusion map approach.


Introduction
Pattern formation at large spatial scales in the cortex is a field of much recent interest [5,10,15,31]. Neural dynamics are intrinsically noisy, so it is natural to study the effects of noise on models for such pattern formation. In this paper we revisit a noisy model for pattern formation in a neural system [27] but this time concentrating on describing its dynamics in a low-dimensional way. We will study the partial integrodifferential system on the domain −π ≤ x ≤ π with periodic boundary conditions, originally analyzed by Laing and Longtin [27] as a simple model for pattern formation in a neural field. The variable u(x, t) describes the neural activity at position x and time t. The function J describes the spatial coupling within the network, with J(x − y) being the strength of coupling between neurons at position x and neurons at position y. The function f (I) describes the firing rate of a single neuron with input I. The quantities A, I and τ are constant. The value of A determines the strength of adaptation, I is a background current applied to all neurons and τ is the adaptation time constant. Equation (2) represents the effects of spike frequency adaptation. Without (2) systems of this form have been used to model orientation tuning in the visual system, working memory and the head direction system [15,5]. Variants including (2) have been subsequently explored [10,11,31]. In our numerical implementation we spatially discretise the domain with M equally spaced points, resulting in the 2M ODEs for i = 1, . . . , M, where J ij = J(2π|i − j|/M). As in Laing and Longtin [27] we set I = −0.1, τ = 5 and use f (x) = (1 + tanh(10x))/2 and J(x) = 0.05 + 0.24 cos x, with M = 100. We add noise to the system by adding a Gaussian white noise term ξ i (t) to each of (3), with ξ i (t) = 0 and ξ i (t)ξ j (s) = 2ην ij δ(t − s), where ν ij = 0 if i = j and 1 if i = j. Initially we set the parameters A = 0.17 and η = 10 −4 ; we will study the effect of varying them below.
At the deterministic limit η = 0, for these parameter values the system dynamics ultimately exhibit a single stable activity bump that moves around the (periodic) domain in one direction (determined by the initial condition), at a constant velocity. If η = 0, the bump motion becomes stochastic, exhibiting occasional switches in direction. This effect, studied by Laing and Longtin [27], is illustrated in figure 1; they concentrated on the net distance travelled during a finite time interval. For small η the bump moved in only one direction, so that the net distance travelled was quite far (effectively proportional to the time); for η greater than a certain value, switches in direction lead to a much smaller net distance travelled. This phenomenon was referred to as noise-induced stabilisation.
Laing and Longtin [27] showed that for a coupling function of the form used here, the system (1)-(2) underwent a pitchfork bifurcation in speed as A was increased. They then added Gaussian white noise to the normal form of a pitchfork bifurcation, and qualitatively reproduced the noise-induced stabilisation. However, the actual derivation of a "noisy normal form" from the original system with noise was not attempted. In this paper we will demonstrate numerically that the noisy discretised system (3)- (4) undergoes what we will characterize as an effective pitchfork bifurcation in speed as A is increased, and investigate the effects of varying noise intensity. This computerassisted bifurcation analysis will be performed assuming that the dynamics of the highdimensional stochastic system (3)-(4) can be effectively described by the dynamics of a single, scalar coarse-grained observable (variable). The evolution of this variable can be determined by running many short, appropriately-initialized simulations of (3)-(4), in the spirit of the "equation-free" framework [23,26].
The value of our scalar variable V (defined below) can be interpreted as the position of a "particle" moving in a potential, subject to noise. Beyond the effective pitchfork bifurcation the potential is double-welled; before it, it is single-welled. Each of the two wells corresponds to persistent motion in one direction (left or right); in the case of a single well we have a stationary (on average) activity bump. Since the system is isotropic, we expect the double well to be symmetric. Although we cannot analytically derive the effective potential and noise in our hypothesized Langevin description for V , we will show how they can be estimated using appropriately initialized short bursts of simulation of (3)-(4). Once the effective potential is approximated, we can locate its maxima and minima (which take the place of unstable and stable fixed points for the particle, respectively) and estimate average transition times between minima (corresponding to bump direction switches). By determining how the potential changes as A and η are varied we can also numerically obtain a bifurcation diagram for its extrema, which takes the place of the traditional steady-state bifurcation diagram for the deterministic problem; identifying transition states (potential maxima, corresponding to unstable steady states) would be extremely difficult if not impossible using long simulations alone. An interesting point is that we can quantify the relationship between the noise intensity in the original system, η, and the effective noise intensity, D(V ), to which the particle in the potential is subject.

Reconstructing the potential
We choose as our coarse scalar variable, V (t), the instantaneous difference in position between the peak of u(x, t) and the peak of a(x, t). This choice comes from observations of the system dynamics: practically all system profiles during a long simulation are characterized by a single peak in each of the variable fields. Furthermore, the difference between these peaks is clearly related to the instantaneous speed of the bump: the peak of a normally lags the peak of u relative to the direction of bump motion. One can clearly see in figure 2 that periods of bump motion to the left (right) are characterized by effectively constant negative (positive) V (t). Factoring out the small amplitude noisy oscillations riding on the bump waveform, we identify the position of the peak of u(x, t) as the location c u that satisfies where x i = 2πi/M. (Effectively, we are finding the amount by which a simple sine wave must be shifted so that its inner product with u(x, t) is zero.) Similarly, c a satisfies Results for a typical simulation are shown in figure 2 and a typical probability distribution of V , for the same parameters, is shown in figure 3.
To establish the connection between V and the effective potential Φ(V ) we assume that V (t) satisfies an (unknown) Langevin equation. Then, the probability density P (V, t) satisfies a Fokker-Planck equation of the form The (unknown in closed form) quantities µ(V ) and D(V ) are estimated from ensembles of brief simulations initialized at V : where ∆V = V (t + ∆t) − V (t) and the average is taken over the realization ensemble. More accurate, maximum likelihood based estimation techniques [1] can also be used.    At stationarity we have P (V ) ∝ exp [−βΦ(V )] where Φ(V ) is the effective potential, related to µ(V ) and D(V ) by [13,35] Data for the estimation of µ(V 0 ) and D(V 0 ) for a given value of V 0 can be obtained in two ways. Firstly, we find occurrences of this value in a long simulation; we then track V (t) over a subsequent time interval (we use 14 time units) for each of these occurrences, and then average the appropriate quantities. This approach has been used previously to estimate drift and diffusion terms in Langevin equations [19]. Typical results of estimating µ(V ) and D(V ) using this approach are shown in figure 4. Performing this for a grid of V 0 values provides enough data for the numerical estimation of the integral in (9); the results for the estimates in figure 4 are shown in figure 5.
In the absence of a database from a long enough (equilibrium) simulation, a a better -more economical -way to find µ(V ) and D(V ), and thus Φ(V ), is to deliberately initialize the system with a given value of V 0 and then track V (t) over a subsequent time interval (we again use 14 time units), repeating and averaging as above. The advantage of this approach, as compared with the previous one, is that it need only be performed for as many V 0 values as necessary for a given accuracy in the integral evaluation. Furthermore, for each such V 0 value we can reinitialize as many independent runs as necessary for accurate estimation at will, without requiring a simulation so long that rare values of V 0 are revisited enough times. Results of the latter type of calculation are shown in figure 6, and figure 7 shows the reconstructed potential from the estimates in figure 6. The only significant difference in this case is that initializing the system with a particular value of V seems to result in a lower estimate of both D(V ) and Φ(V ) when |V | is large. Recall, however, that these regions are rarely visited in a simulation, and thus the discrepancy may be attributed to an insufficiently long database for the first approach.

Bifurcations
A traditional bifurcation diagram for the deterministic problem with respect to a parameter such as A would involve tracing branches of steady states and constant shape travelling waves (which can also be reduced to solutions of a fixed point problem) using pseudo-arclength and branch switching techniques (e.g. [2]). In the coarse-grained stochastic case, it is natural to trace instead the zeros of the drift µ(V ) as the same parameter A is varied, using the same standard bifurcation codes. It is worth noting that for the case of state-dependent noise, the local maxima of the probability density do not exactly correspond to zeros of µ(V ) and one needs instead to find fixed points of the effective potential given in (9) (zeros of the right hand side of (9) differentiated with respect to V ). For the scalar case, one can implement secant-type iterative methods to converge to the zeros of the appropriate function using function estimates only; for the multivariable case matrix-free iterative techniques like Newton-Krylov GMRES can be used (see [22,33]). In the cases we study µ(V ) is well approximated by a cubic function; we take advantage of this simplification by estimating µ(V ) for just four values of V , uniquely defining this cubic whose zeros we can then easily find. Such a "noisy" bifurcation diagram is shown in figure 8, where we vary A. Stability can be determined from the local slope, µ ′ (V ), at the fixed point, and is indicated in the figure. Parameters of a different nature, such as the noise intensity or colour, can also be varied in this context. Figure 9 shows the noise-induced "effective pitchfork" bifurcation that occurs as η is varied. The location of these effective bifurcations was verified by running long simulations in their vicinity (results not shown). Note that if we assume in advance that the potential Φ(V ) is symmetric (as the underlying system is), the cubic function representing µ(V ) would have no quadratic or constant terms, so that only two evaluations of µ(V ) would be needed to approximate the required cubic.

Switching times
One of the most important statistics of our problem is the average time between changes of direction, or, more generally, the statistics of these switching times for the bump. Based on our effective Fokker-Planck model and Kramers' theory, once we have βΦ(V ) we can also estimate the average time between switches in direction [13], as where V min is the value of V at which Φ has a minimum, D = [D(0) + D(V min )]/2 and ∆Φ = Φ(0) − Φ(V min ). Measurements of D and Φ for η = 10 −4 and A = 0.17 give τ ≈ 3 × 10 3 . A typical distribution of waiting times from a long simulation is shown in figure 10. The mean for the data shown is 3.5 × 10 3 , in excellent agreement with our Kramers' approximation. Along the same lines, we can quantify how the average switching time, τ , depends on parameters. Typical results are shown in figure 11. Here, the effective potential is not estimated from a long simulation database, but rather "on demand" by initializing the system at particular values of V , performing short time simulation bursts, and processing their results as before. A clear advantage of this approach, as opposed to running the system for long enough to measure sufficient occurrences of switches, is that large values of τ can be inferred from a reasonable number of short simulations. See figure 11, for example, where we can "measure" switching times greater than 10 6 .
For a simulation of a given, fixed duration T , these curves can be used to infer the parameter values for which we expect to see direction changes. If τ > T , we will typically see the bump move in one direction only, whereas if τ < T we will see switching, leading to a drop in the absolute value of the distance travelled during the time T . For example, from figure 11 (b) we see that for simulations of length T = 10 3 , η must be greater than approximately 10 −3 to observe any switching of direction.

White versus coloured noise
In [27] the authors investigated the effects of adding coloured noise to (3)-(4). They found that increasing the correlation time of the coloured noise made the bump less likely to move. This was rationalized in terms of "frozen" noise (i.e. spatial inhomogeneity) which was likely to "pin" the bump, preventing it from moving. Increasing the correlation time of the coloured noise was thought of as interpolating between Gaussian white noise (with delta function autocorrelation in time) to frozen, spatially structured, noise. Here we investigate the effect again, quantifying the influence of noise colour on the location of the underlying bifurcation.
As in [27] we add noise to the system by adding a term η i (t) to each of (3), with η i (t) = 0 and { η i (t)η j (s) } = 2ǫν ij e −|t−s|/λ , where ν ij = 0 if i = j and 1 if i = j. The notation {· · ·} indicates averaging over the initial distribution of η(0) values, taken from a Gaussian with zero mean and variance 2ǫ. We keep ǫ = 10 −4 and vary the correlation time λ. The results are shown in figure 12 where we have plotted zeros of µ(V ) as a function of A for λ = 1 and λ = 100. The figure was computed using the short runs initialized at prescribed V values. The results are consistent with those found in [27], where it was seen that coloured noise was more effective at slowing the bump than Gaussian white noise of the same power, and that the effect was stronger for longer correlation times. The results presented here show that this effect can be rationalized in terms of shifting a bifurcation point. Results for several other values of λ indicate that as λ is increased, the bifurcation moves to higher values of A (results not shown).

Diffusion maps and the data-based detection of coarse observables
The previous results relied on our ability to choose a scalar variable, V , whose value correlates with the state of the high-dimensional system (after brief initial transients have equilibrated). However, we are often faced with dynamical systems for which choosing such a variable is far from obvious. We then need a systematic procedure for determining, from the results of a simulation, a good low-dimensional representation of the system state. Such data-mining, so-called "manifold learning" techniques have Figure 13. Schematic of datapoints (filled circles) in R 2 that lie along a curve (solid line). Euclidean distance is a good measure of the separation between points A, B, and C along the curve but is a poor measure of the separation of points D, E, and F along the curve.
been a focus of intense research in recent years [8,9,14]; we will use here the recently developed diffusion map approach [28,29]. A large data ensemble from direct simulations of our system can be represented as a cloud of points in a high-dimensional space (here, R 200 ); see the schematic in figure 13. When two such points (activity profiles) are very close to each other (i.e. when their Euclidean distance is small, as in the case of points A, B and C in the schematic), we can consider this distance as representative of the "intrinsic similarity" of the two configurations -in some sense, of how easy it is for the system dynamics to cause a transition from one configuration to the other. When, however, this Euclidean distance is larger than some threshold (as in the case of points D, E and F in the schematic) the Euclidean distance stops being a good measure of the "intrinsic similarity" between configurations -the "effort to transition" from E to F, measured by the arclength between them, is clearly less than the effort to transition from E to D, even though the Euclidean distances DE and EF are similar. The main idea underpinning diffusion maps is to perform a random walk on a graph in which data points are vertices, and connection strengths between data points are given by a Gaussian kernel of the form K (x, y) = exp (− x − y 2 /σ 2 ).
When the points are farther away than a cutoff (controlled by the parameter σ) the vertices are effectively disconnected; when they are very close, the strength of the connection is controlled by their (small) Euclidean distance. A random walk on such a graph gives rise to a diffusion distance -starting from the same source point and diffusing on the graph for some time, we then look for equal density contours; points on such a contour are equally easy to access from the source point, and are therefore at equal diffusion distance from it, even though their Euclidean distances from the source may vary substantially. The procedure also provides a set of transformed coordinates; the Euclidean distance in these new coordinates is a true measure of intrinsic datapoint similarity. The procedure can be thought of as a nonlinear generalization of Principal Component Analysis [21].
In our case, we start from a long simulation (of 30,000 time units) during which we sample the u(x), a(x) profile every 8 time units, giving N = 3,750 data points.
(From now on, we use the same parameter values as those in figure 1.) Using the procedure briefly outlined in the Appendix we process these data forming a Markov matrix whose leading eigenvectors provide useful "reduction coordinates" for our datain particular, the second eigenvector Φ 2 of this matrix is our computer-assisted candidate coarse observable (the first eigenvector, corresponding to the eigenvalue 1, is trivial). Figure 14 shows the simulation data projected on the first two "reduction coordinates" -the first two nontrivial eigenvectors Φ 2 and Φ 3 of the appropriate Markov matrix. The right panel of figure 14 shows a three-dimensional diffusion map, in terms of the first three nontrivial reduction coordinates. The coordinates of the i th datapoint in this map are (Φ 4 ) . Simulation points in the diffusion map are coloured according to their associated value of the empirical coordinate V -it is clear that points are ordered along the curve by their values of this known reaction coordinate. This sorting of highdimensional data vectors is obtained in an "automated" fashion by the diffusion map calculation.
All data instances clearly collapse on a one-dimensional curve. Figure 15 indicates five representative data points (at values of Φ 2 = −1.8, −1, 0, 1, and 1.8) on this curve. The insets, showing brief space-time plots initialized at the corresponding data points, clearly indicate that right-moving states reside on one end of this curve, left moving states on the opposite end, and transition states between the two reside in the interior of the curve. It is clear that the curve is one-to-one with the variable Φ 2 (for every value of Φ 2 there exists one point on the curve); this suggests that Φ 2 would be a good scalar observable for our system. Figure 2 (third row) shows the evolution of this data-based observable during our system dynamics; it is clearly capable of capturing the direction switches in the system in a manner comparable with the "experience-based" coordinate V . Another possible scalar observable we considered is the arc-length distance along the curve connecting components of the top two eigenvectors. The variation in this observable is comparable to that of Φ 2 and also accurately describes the direction of travel of the bump during the simulation (results not shown). Figure 16 plots the simulation data in terms of the two coarse coordinates: the experience-based V and the data-based Φ 2 . At first sight, a one-to-one correspondence between the two observables is apparent, suggesting that our experience-based coordinate was indeed a good choice.
The processing of the database from a long, equilibrium run simulation that gave us an effective potential in terms of the observable V can now also be performed in terms of the computer-assisted observable Φ 2 ; the resulting two-well effective potential and the corresponding effective Langevin description are practically isomorphic to the V -based description. The effective potential computed in terms of the computer-assisted observable Φ 2 is shown in figure 17.

Computing with diffusion map coordinates
When we attempt to directly construct the effective potential in terms of the diffusion map coordinate Φ 2 (i.e., not from processing a long run database), it becomes necessary to initialize short simulation bursts consistent with particular values of Φ 2 as well as observe the Φ 2 values corresponding to results of detailed simulation. The diffusion map calculation that leads to the identification of the coordinate Φ 2 automatically provides its value on each of the data points used in the calculation (see the Appendix). It therefore becomes important to establish computational protocols that routinely allow the translation between physical coordinates and diffusion map coordinates. This translation must be possible in both directions: to initialize a short system run we need to translate a value of Φ 2 to system initial conditions (living in R 200 ); processing of the short simulation bursts involves extracting the Φ 2 values of the resulting states. These translation operations are termed "lifting" and "restriction" in the equation-free framework [18,23]. We now briefly discuss possible implementations of such protocols, starting with the restriction: obtaining the Φ 2 values of new system configurations, not in the original database, that result from our short simulation bursts.
For points in the original dataset (used to assemble the neighbourhood matrix K) For a particular point x i in the dataset we have For a new datapoint x (new) the same formula should be valid. However, while the right hand side is unknown, the left hand side may be computed for the new datapoint, which gives the following formula for the extension of the diffusion map where M s is a generalized kernel which extends to new sample points and is given by Equation (13)   by eigenspace interpolation instead of repeated eigendecomposition of the Markov matrix augmented by the new data vectors. The expectations appearing in the denominator of (14) are computed from the empirical data using Use of the Nyström formula (13) is central to the efficient computation of diffusion map coordinates associated with a high-dimensional vector (restriction) and also the preparation of a data vector x (lifting) with desired diffusion map coordinate values. Equation (13) is used to first compute the eigenvector components Ψ The Nyström formula allows calculation of the diffusion map eigenvector components Φ new j associated with a new data vector x (new) . A full eigendecomposition is typically performed first for a representative subset of simulation datapoints, identifying the relevant eigenvectors. The Nyström formula is then used to perform the restriction operation in (13) which amounts to interpolation in diffusion map space.
Our protocol for lifting from prescribed diffusion map coordinates to consistent system states using stochastic optimization is presented in figure 18. The main step in this lifting procedure, which produces u(x) and a(x) solution profiles with specified diffusion map coordinates, is the minimization of a quadratic objective function given by where λ OBJ,2 is a weighting parameter that controls the shape of the objective away from its minimum at Φ The implicit dependence of the eigenvector components Φ (i) j (associated with a data vector x (i) ) on x (i) and all other data vectors used in the construction of the neighbourhood matrix K makes this optimization problem challenging. We used here, for simplicity, the method of Simulated Annealing (SA) [24,32] to minimize the objective function defined in (17), and identify a data vector x (targ) with the target diffusion map coordinate Φ targ 2 . Table 1 compares drift and diffusion coefficient values (at Φ 2 = −0.5) estimated by multiple short simulation bursts, initialized using this lifting protocol, with values computed from a long time simulation. The agreement is very good.
Given these protocols, the entire effective potential and, more generally, all the bifurcation/switching time computations performed above with our empirical V coarse variable can be repeated with the Φ 2 variable. "Coarse-grained smoothness" in the diffusion map coordinates can be exploited to guide the efficient exploration of the effective potential surface [17].  Method Lifting −3.39 × 10 −3 3.20 × 10 −3

Discussion and Conclusion
In this paper we have revisited a stochastic spatio-temporal pattern forming system originally used to model activity in the cortex. Previous work [27] investigated the effects of changing two parameters (spike frequency adaptation strength and noise intensity) on the dynamics of a spatially-localized "bump" of neural activity. Here, we have concentrated on deriving and using an effective low-dimensional description of the system in terms of a single scalar coarse variable. This variable was assumed to satisfy an (unknown) effective Langevin equation, so that its probability density satisfies an effective Fokker-Planck equation. By appropriate processing of the results of a long simulation, or -more efficiently -by appropriately initializing short bursts of simulation of the full system, we estimated the drift and diffusion functions that appear in the Fokker-Planck equation, and thus determined an "effective potential" for the system on demand. In the second part of the paper we showed how a similar analysis can be performed using a variable that is extracted by data-mining a sufficiently long system simulation in an automated fashion, using the diffusion map approach [28,29]. Such an approach could prove particularly useful for other systems, since in general a good coordinate (or coordinates) for a low-dimensional description of a complex system may not be known, nor easy to guess. The results obtained using this new variable showed good correspondence with those obtained in the first part of the paper. Considering more general systems, it is known that in certain limiting cases (e.g. weak coupling) it is possible to explicitly obtain accurate reduced descriptions (e.g. phase oscillator equations) for large coupled neural systems [16,20]. Clearly, if such equations can be derived, they are much easier to use than the approach presented here; our approach is intended for cases where we believe the reduction is possible, but cannot be explicitly performed. One of the most challenging tests for a reduction -as we move away from the conditions where we can guarantee its validity analytically -is to test, on line, whether it is accurate, and -importantly -whether a different level reduction, with more (or even possibly with fewer) variables is in order. In our case, this would correspond to devising tests to suggest, on line, that more than one coarse-grained variable must take part in our effective Langevin equation. The slight "thickness" of the line in figure 16 could be an indication that, for the parameter values used, more than one variable may be necessary in our model. Devising such tests -in effect, testing the hypothesis that the data are locally consistent with a particular model, e.g. a scalar Langevin equation -is the subject of ongoing research across several disciplines (statistics, financial mathematics). Integrating such tools in a multiscale simulation framework is only just starting. Even though the model discussed here is rather abstract (from a computational neuroscience point of view) it adds to the body of work demonstrating that it is possible to obtain useful low-dimensional descriptions of complex systems, both in neural modelling [26] and elsewhere [13,23,35]. These descriptions not only provide insight into the fundamental dynamics, but enable one to simulate and analyse such systems in an efficient manner.
Coarse-graining large scale, faithful neural network computations constitutes an important subject of intense current research (see, for example, the probability density approach [7,25,30]). Multiscale, coarse-graining numerical algorithms such as the one demonstrated here have an important role to play in elucidating the types of behavior possible for such models and their parametric dependence.