Continuation of localised coherent structures in nonlocal neural field equations

We study localised activity patterns in neural field equations posed on the Euclidean plane; such models are commonly used to describe the coarse-grained activity of large ensembles of cortical neurons in a spatially continuous way. We employ matrix-free Newton-Krylov solvers and perform numerical continuation of localised patterns directly on the integral form of the equation. This opens up the possibility to study systems whose synaptic kernel does not lead to an equivalent PDE formulation. We present a numerical bifurcation study of localised states and show that the proposed models support patterns of activity with varying spatial extent through the mechanism of homoclinic snaking. The regular organisation of these patterns is due to spatial interactions at a specific scale associated with the separation of excitation peaks in the chosen connectivity function. The results presented form a basis for the general study of localised cortical activity with inputs and, more specifically, for investigating the localised spread of orientation selective activity that has been observed in the primary visual cortex with local visual input.

1. Introduction. One of the most challenging research questions in neuroscience is understanding the relationship between spatially-structured cortical states and the underlying neural circuitry that supports them. A popular approach for analysing coarse-grained activity of large ensembles of neurons in the cortex is to model cortical space as a continuum. Since the pioneering work of Wilson and Cowan [58,59] and Amari [1,2], continuous neural field models have become a popular and effective tool in neuroscience. In such models, the large-scale activity of spatially-extended networks of neurons is described in terms of nonlinear integro-differential equations, whose associated integral kernels represent the spatial distribution of neuronal synaptic connections. The canonical Wilson-Cowan-Amari neural field equation [2,59] ∂ ∂t u(r, t) = −u(r, t) + Ω w(r, r )S u(r , t) dm(r ) (1.1) describes the evolution of the average membrane voltage potential of a neuronal population u(r, t), at position r ∈ Ω on the cortex and at time t. The nonlinear function S represents the neural firing rate, whereas the connectivity function w(r, r ) models how a population of neurons at position r on the cortex interacts with a population at position r . Frequently-used firing rate functions S include the Heaviside step function, piecewise-linear functions or smooth sigmoidal functions. Various choices are also possible for the connectivity function, which is often assumed to be translation invariant (that is, dependent on the Euclidean distance r − r ) and localised in space. The cortical domain Ω is usually a subset of R d , with d = 1 or, in more realistic models, d = 2. For a recent review on neural fields modeling, we refer to Bressloff [8].
Unlike spiking neural network models, continuous field models have the advantage that analytic techniques for partial differential equations (PDEs) can be adapted to study the formation of patterns and their dependence upon control parameters. Various types of coherent structures have been observed in neural field models, ranging from spatially and temporally periodic patterns to travelling waves and spiral waves [19,27,42]. Neural field equations have also successfully been used to model a wide range of neurobiological phenomena such as visual hallucinations [9,28], mechanisms for short term memory [44] and feature selectivity in the visual cortex [6,35].
A common strategy to derive analytical and numerical results for the nonlocal Eq. (1.1) is to assume translation invariance and exploit the freedom in the choice of the connectivity function w: if the Fourier Transform of w is a rational function, it is possible to derive a PDE formulation that is equivalent to the integral model [42]. Coherent structures supported by the original model can then be conveniently constructed and analyzed in the PDE framework. Indeed, previous studies have been carried out in cases where the synaptic kernel led to an equivalent PDE formulation [43,44]. To the best of the author's knowledge, there has been no attempt to propose efficient path-following methods for general connectivity functions w. This paper is motivated by the desire to develop numerical algorithms for Eq. (1.1) without relying on an equivalent PDE formulation. More precisely, we discuss how to solve Eq. (1.1) when Ω = R 2 , S is a smooth function of sigmoidal type and w has a generic Fourier Transform.
The main tools for our investigation are time simulation and numerical continuation. When a PDE model is available, we use standard techniques for both tasks, in line with what is typically done in several other works in this field [41][42][43]. However we show that, when the integral of Eq. (1.1) can be written as a convolution, it is convenient to employ a fast Fourier transform (FFT) for both time stepping and numerical continuation. Direct numerical simulations of Eq. (1.1) using FFTs have been performed before on full integral models [23,32]. In the present paper, we combine FFTs with Newton-Krylov solvers [38,39], thus opening up the possibility to perform numerical continuation directly on the integral model.
We concentrate our effort on the emergence and bifurcation structure of stationary localised patterns in planar neural field models of the form Eq. (1.1). Indeed, this type of solution is of great interest in models of prefrontal cortex, where localised states are believed to characterise working memory [18,33]. Recently, localised regions of activity have been observed in the cat primary visual cortex [17] when the animal is presented with localised-oriented input. Moreover, some reported drug-induced visual hallucinations have also been found to be spatially localised [54] indicating the existence of spatially localised regions of activity in the human primary visual cortex.
Localised states have been observed in a wide variety of nonlinear media [40]. The bifurcation structure of localised solutions has been studied extensively in the Swift-Hohenberg equation posed on the real line [11,12,16,25,60] and on the plane [3, 5, 45-47, 51, 52]. In this context, a well-known mechanism for the formation of localised states is homoclinic snaking: solutions with one or more bumps at the core emerge from the trivial homogeneous state and undergo a series of fold bifurcations, giving rise to a hierarchy of states with an increasing number of bumps. This scenario seems to be a common footprint of localised patterns, extending far beyond the prototypical Swift-Hohenberg equation [36,53] and to have also been found in nonlocal equations such as neural fields [21,29,43,44].
An example of homoclinic snaking is given in Fig. 1.1, where we show a bifurcation diagram for the integral model (1.1) posed on the real line with where b and µ control the decay of the synaptic kernel and the slope of the sigmoidal firing rate respectively, while θ is a threshold value. As µ varies, the trivial steady state u = 0 bifurcates at a subcritical Reversible Hopf bifurcation, from which a branch of periodic solutions and a branch of localised states originate. Beyond the fold L p on the periodic branch, a stable periodic state, shown in the inset (a), coexists with the trivial state. The branch of localised states features solutions with an odd . . number of bumps and snakes for µ ∈ [µ 1 , µ 2 ]; in this interval, localised solutions with different spatial extent coexist and are stable (see insets (b)-(d)). We note the existence of a counterpart even-numbered-bump solution branch along with so-called ladder branches connecting the odd and even branches (not shown). For the same connectivity function used here, snaking has been shown to occur in terms of the parameter b [44]. Elsewhere, Elvin et al. [26] used the Hamiltonian structure of the steady states of the model (1.1) with (1.2)-(1.3) and developed numerical techniques to find homoclinic orbits of the system. Snaking has also been shown to occur with a Mexican-hat connectivity function [21] and for a wizard-hat connectivity function [29]. In the latter study, normal form theory for a Reversible Hopf bifurcation was applied to prove the existence of localised solutions and a comprehensive parameter study was carried out with numerical continuation in terms of two parameters controlling the nonlinearity and a third controlling the shape of the connectivity function.
For the neural field equation posed on the Euclidean plane, various types of spatially localised two-dimensional states have been found including radially symmetric solutions [10,30,32,43,44,55,57], rings [23,50], hexagonal patches [30,43] along with more complex breathing and travelling states [22,31,50]. When working in the Euclidean plane, it is still possible to derive an equivalent PDE for suitable choices of the connectivity function w or of its Fourier transform w [43]. Following this approach, normal form theory has recently been applied to prove the existence of localised solution branches in both the Euclidean plane and on the Hyperbolic disk with a wizard-hat connectivity function [30]. When these solutions were pathfollowed using numerical continuation, it was found that the branches do not undergo snaking-type behaviour. However, for the connectivity (1.2), snaking was shown to occur for branches of radially-symmetric solutions [43]. Furthermore, the existence of D6-symmetric and D3-symmetric localised states were found at isolated parameter values.
In Fig. 1  with radially-symmetric connectivity function w(r, r ) := w( r − r ) = e −b r−r (b sin r − r + cos r − r ), (1.4) and sigmoidal firing rate function given by Eq. (1.3); various combination of initial conditions and control parameters lead to three different steady states. We note that this models supports localised states, such as the radially-symmetric spot of panel (c) or the hexagonal pattern of panel (f). Furthermore, changes in the slope of the sigmoidal firing rate affect the stability properties of the solutions, leading to other localised states or to domain-covering patterns such as the one in panel (i). In the present paper we will focus on localised planar patterns (with various symmetry properties) that coexist with the trivial state u = 0 and with fully periodic states, similar to the one shown in Fig. 1.1 for the 1D case.
A key result of the present paper is that homoclinic snaking occurs in planar neural field models for non-radial patterns and that the choice of the connectivity function has a considerable impact on the snaking structure, as well as on the stability properties and the selection of the localised states. In two spatial dimensions, periodic and localised solution branches bifurcate from the trivial state at a Turing instability (as opposed to a Reversible Hopf in 1D) and snake irregularly, in a similar fashion to what is found for the planar Swift-Hohenberg equation [46].
The outline of the paper is as follows. In Sec. 2, we present the different models that we study, and then, in Sec. 3, we describe the numerical methods that are used to analyze each model. Our homoclinic snaking results of localised states are presented in Sec. 4.

Models.
In this section, we introduce several neural field models used in the paper: our starting points are the models introduced by Laing et al. [42,43], in which the cortical space Ω is assumed to be the Euclidean plane R 2 .
2.1. Integral model. The first model that we will consider is obtained from Eq. (1.1) assuming a translation-invariant, radially-symmetric kernel with sigmoidal firing rate (shown in Fig. 2.1(a)) radial connectivity function and external inhomogeneous input In the visual cortex regions of activation have been shown to have a Gaussian spread for radially-symmetric visual inputs [17], hence our choice for the function g(r). The connectivity function is plotted in Fig 3) was proposed in models of working memory as a description of synaptic connections in the prefrontal cortex [34,44]. The connectivity describes local excitation and longer-range connections that alternate between inhibition (w(r) < 0) and excitation (w(r) > 0). We argue that this type of connectivity pattern is also relevant to the study of patterns of activity in early visual areas like V1 where there is a characteristic length scale associated with the average orientation hypercolumn width. It has been shown in anatomical studies that the number of lateral connections decay with distance, that the number of excitatory connections peak each hypercolumn width and the number of inhibitory connections peak each half-hypercolumn width [13]. The net effect is alternating bands of inhibition and excitation that decay with distance. This is also consistent with auto-correlations computed for the orientation selectivity map [49] given that connections tend to be reinforced between neurons with similar orientation preference. Henceforth the model (2.1)-(2.3) will be referred to as the integral model (IM).

Fourth-order PDE approximation.
In the cases when the Fourier transform of the synaptic kernel is a rational function, it is possible to derive an equivalent PDE formulation of Eq. (2.1) [42]. For simplicity, we will consider models without an external input g(r). If w(k) = P (k)/Q(k) with P and Q even functions in k where k = k for k ∈ R 2 , then the Fourier transform of Eq. (2.1) gives An inverse Fourier transform of the equation above leads to the desired PDE where L P and L Q are linear operators containing spatial derivatives of even order.
Since the Fourier transform of the connectivity function (2.3) does not have an analytic expression, the integral model does not admit an equivalent PDE. However, we can approximate w with a function whose Fourier transform is rational and then derive an approximate PDE for the integral model. We specify the approximate connectivity function w 4 (r) through its Hankel Transform where w 4 (k) ≈ w(k) is given by and the coefficients A, B, M are determined using a least-squares best fit algorithm (see Sec. 3.4.3 for further details). We compare the approximation with the original connectivity function in the physical and Fourier domains in Fig. 2.2(a) and (b). In physical space the two functions appear to be similar. The key qualitative difference is that in Fourier space w 4 (0) > 0, which is not consistent with the IM connectivity function, for which w(0) < 0. Biologically this means that w 4 represents a globally excitatory connectivity function, whereas w represents a globally inhibitory connectivity function. We will see in Sec. 4.2.2 that it is necessary to increase the order of the polynomials in the numerator and denominator of Eq. (2.5) in order to accurately capture the sign at k = 0.
Starting from the expression for w 4 , we derive the corresponding PDE, containing spatial derivatives up to the fourth order .
.  where the sigmoidal firing rate function is identical to the integral model case (2.2). In Eq. (2.6) we have denoted by ∆ the standard Euclidean Laplacian. Henceforth, this model will be referred to as PDE4.

Eight-order PDE approximation.
In order to obtain a more accurate representation of the connectivity function (2.3), we repeat the steps outlined in Sec. 2.2 with the following approximation where the values of A, B, C, D and M are determined using a least-squares best fit algorithm; see Sec.
where the sigmoidal firing rate function is identical to the integral model case. Henceforth this model will be denoted as PDE8.
3. Numerical methods. In this section, we review the numerical methods employed for the computation of localised states in IM, PDE4 and PDE8.
Similarly, we form vectors w, S(u), g ∈ R N 2 for the approximations to w, S(u) and g, respectively. Further, we introduce the discrete convolution, where we have denoted by F and F −1 the 2D Fourier Transform and its inverse, respectively. In summary, the discrete version of the evolution equation (2.1) is given byu This type of discretization has been applied before in direct numerical simulations of neural models (see, for instance [23,32]) even though it has not been used for numerical continuation. For smooth firing rate functions, the right-hand side can be evaluated accurately and efficiently using a Fast Fourier Transform (FFT) and its inverse (IFFT). In passing, we note that since the FFT of w can be performed and stored at the beginning of the computation, one function evaluation of the righthand side requires just one FFT and one IFFT. Furthermore, standard de-aliasing techniques can be applied to the convolution operator if required [14]. Once a stable steady-state of Eq. (3.2) is found via direct numerical simulation, it is possible to continue it in one of the control parameters using standard numerical continuation techniques. In previous studies of neural field equations, numerical continuation was performed on an equivalent PDE formulation of the integral system.
A key observation is that path following can be applied directly to IM (or to similar models), employing FFTs and Newton-Krylov solvers [38,39]. Such methods do not require the formation of a Jacobian matrix, but rely only on Jacobian-vector multiplications: for IM, this is conveniently done using just a single application of FFT and IFFT. We remark that Newton-Krylov methods are often used in conjunction with sparse systems, but the performance of FFTs and IFFTs makes them a favourable choice for IM, even though the system is full.
For numerical continuation of steady states of IM, we solve the system of algebraic equations whose associated Jacobian-vector product is given by where S (u) = diag(S (u 11 ), . . . , S (u N N )) ∈ R N 2 ×N 2 . We solve the system (3.3) using a Newton-GMRES solver implemented in Matlab and continue the solution with a secant method. Eigenvalue computations can also be performed using the Jacobian-vector products (3.4). Details of the numerical implementation and of the numerical parameters can be found in Sec. 3.4.1. Remark 3.1. The external input g guarantees that the system of algebraic equations (3.3) is not translation invariant, even when the problem is complemented with periodic boundary conditions. Unless otherwise stated, we will use a negligible external input for the IM, so that Newton iterations can be applied directly to Eq. (3.3). We point out that a similar result can be obtained without imposing any external input, by perturbing the right-hand side with a term containing one extra unknown and then closing the system with a suitable phase condition [3,15,46].

3.2.
Fourth-order PDE approximation. In order to continue stationary localised solutions to PDE4 with D6 symmetry, we use polar coordinates and pose a boundary-value problem on the sector Ω 1/6 = { (r, θ) ∈ R 2 | 0 < r < R, 0 < θ < π/3 } with Neumann boundary conditions where the Laplacian operator ∆ is expressed in polar coordinates. We discretise the system above using finite differences in r and a Fourier collocation method in θ, with N r and N θ evenly spaced points, respectively, leading to a system of nonlinear algebraic equations of the form where I is the identity matrix. The Laplacian matrix L is formed explicitly, starting from differentiation matrices D r , D rr , D θθ for spatial derivatives with Neumann boundary conditions and then combining them with Kronecker products [3,46,56] where R = diag(1, r 2 , . . . , r Nr ) and I θ is the N θ -by-N θ identity matrix. For purely radial patterns, we adapt the boundary-value problem so as to contain only the radial direction r and impose Neumann boundary conditions. Numerical continuation of the system (3.6) is performed with a secant method. Further details on the numerical implementation can be found in Sec. 3.4.2.
3.3. Eight-order PDE approximation. For localised solutions of PDE8, we follow a similar approach to the on used for PDE4. In order to avoid the discretisation of 8th-order differential operators, we recast PDE8 as a system of two 4th-order PDEs and seek solutions to the following boundary-value problem Again, the discretisation of this system uses finite-differences in r and Fourier spectral collocation in θ. The example above shows also that, as we approximate w more accurately, the order of the underlying PDE increases and its numerical continuation becomes more demanding. A more convenient approach would be to use directly the integral form for the model, with connectivity function w 8 and proceed with a Newton-GMRES, solver. In this way, the computational cost would be the same as for IM.

Implementation and numerical parameters.
3.4.1. Integral model. Time simulations are carried out using a standard 4th order Runge-Kutta method with fixed step size. At each time step, the right-hand side of Eq. (3.2) is evaluated four times using an Nvidia Graphic Processing Unit (Tesla C2070). To compute the Discrete Fourier Transform we use the CUFFT library provided by Nvidia as part of its CUDA framework [48]. This software library allows us to easily exploit the parallelism of a GPU to obtain a fast implementation. The vector u is kept in the GPU global memory throughout a time step in order to avoid memory transfers and it is updated in parallel once the four stages of the Runge-Kutta scheme have been computed. Transfers between CPU and GPU memory only occur when the result of a time step needs to be saved to a file. Time simulations use a grid of 10 3 × 10 3 points and a stepsize of 0.5. Computation of each time step requires approximately 0.1 seconds.
Numerical continuation for IM is done in Matlab using a in-house secant continuation code which employs a Newton-GMRES method for the nonlinear solves. Unless otherwise stated, we used a grid of 2 10 × 2 10 points and fixed an absolute tolerance of 10 −3 for the nonlinear iterations. The Newton-GMRES solver uses the MATLAB in-built function gmres without preconditioners and with parameters restart = 20, tol = 10 −3 and maxit = 10. Initial guesses are obtained directly from the Runge-Kutta time stepper, interpolating with MATLAB's interp2 function where necessary. Stability computations are performed with Arnoldi iterations via MATLAB's eigs function, passing the Jacobian-vector product (3.4) and computing (with the default tolerance) the first 20 eigenvalues with the largest real part. Computations are performed on a MacPro with a 3 GHz Quad-Core Intel Xeon processor employing exclusively the CPU.
3.4.2. PDE4 and PDE8 models. Numerical continuation for the PDE models have been carried out with a secant code similar to the one used for IM, but using MATLAB's in-built function fsolve for the nonlinear iterations. Unless otherwise stated, we used 300 grid points in the radial directions and 20 in the angular direction. We use the Levenberg-Marquardt algorithm implemented in fsolve and set TolFun = 10 −6 . The sparse Jacobians of these problems are formed and passed directly to the solver. Initial guesses for the continuation have been obtained using the expressions given in Eqs. (A.1) and (A.2). Computations are performed on a standard laptop on a single core.

3.4.3.
Least-squares data fitting. Before comparing the connectivity functions w 4 (r) and w 8 (r) with w(r), it was necessary to tune the parameters in the definitions of w 4 and w 8 . In each case we performed a nonlinear least-squares optimization of the parameters using the lsqcurvefit function in Matlab. For w 4 the objective was to minimize the L 2 -norm of the difference between w 4 and w whilst varying the parameters A, B and M in Eq. (2.5), where w is computed numerically using the Hankel transform at 300 points. Similarly, for w 8 , the L 2 -norm of the difference between w 8 and w was minimised whilst varying the paramaters A, B, M , C, D in Eq. (2.7). For reference, the L 2 -norms of w and w are given above the panels in Fig. 2.2; the norm of the difference between the two functions plotted in each panel is also given. We note the largest Fourier mode of the connectivity dictates the location of bifurcations in terms of the sigmoid parameters θ and µ. By minimising the difference between the connectivity functions in Fourier space, we expect to find similar behaviour for each connectivity over the same parameter ranges. On the other hand, if one minimises the difference between the connectivity functions in physical space and the amplitudes of the largest Fourier modes are not matched, bifurcations occur in different parameter ranges in each model and a direct comparison cannot be made.

Numerical results.
4.1. Convergence of the Newton-GMRES solver. Since Newton-GMRES methods with pseudospectral evaluation of the right-hand side have not been used before for integral neural field models, we report briefly on our solver. To test convergence, we perturbed a localised steady state of IM to obtain an initial guess (panel (b) of Fig. 4.1) and converge back to the original solution (panel (a) of Fig. 4.1) using our Newton-GMRES solver. In panel (c) we plot the relative residuals of each iteration, showing that we achieve convergence within a few nonlinear iterations. Similar convergence plots (not shown) are obtained for the numerical continuation, albeit solutions in that case are achieved with fewer iterations, owing to the more accurate initial guess provided by the secant predictor scheme. The experiment is repeated for various values of N : the convergence diagrams are indistinguishable from the one reported in panel (c). The wall time for the numerical experiment scales linearly with the number of unknowns N 2 , as reported in panel (d). We remark that, even without using any GPU acceleration and without enforcing explicit parallelisation in the CPU, the Newton-GMRES solver finds a solution to a full problem with 1,048,576 unknowns in less than 40 seconds. We remark that translation invariance is removed in the IM by the negligible external input G 0 while in PDE4 and PDE8 this is achieved by the boundary conditions of the problems (3.5) and (3.7), so we choose G 0 = 0. In the present section we focus on the no-(or equivalent negligible-) input case which should be well understood before the addition of an input. In Sec. 5.2 we will provide examples of the model behaviour with input and discuss the implications.
The bifurcation points in the diagrams presented in this section will be labelled as follows. F represents a fold bifurcation and P a spatial-symmetry-breaking bifurcation from a radial state. Superscripts indicate the symmetry properties of the bifurcation, where R represents a bifurcation on a branch with radially symmetric solutions and D6 represents a bifurcation on a branch of D6-symmetric solutions The labels l and r in the subscripts for fold bifurcations indicate whether the fold occurs on the left or right of the snaking structure. The indices n in the subscripts indicate the ordering moving up the snaking structure. On radial branches n corresponds to the number of rings around a central spot solution, for example, on the branch between F R l1 and F R r1 there is one ring around a central spot. On D6-branches there are n(n + 1)/2 × 6 additional spots glued around a central spot, for example, on the branch between F D6 l1 and F D6 r1 the solution has a total of 7 spots.

PDE4 results and comparison with IM.
We discuss both radiallyand D6-symmetric localised solutions of PDE4 with connectivity function defined via Eq. (2.5). Other solution branches with different symmetry properties do exist but these are only discussed for IM in Sec. 4.3. An unstable radial spot solution bifurcates from the trivial state u = 0 at a Turing instability with µ-value to the right of the range shown in the subsequent diagrams. It is this unstable radial spot branch that . .  appears in the bottom right of Fig. 4.2(a) and (b). Figure 4.2(b) shows the snaking structure for radial and D6 branches. We first focus on the radial branch, detail of which is shown in panel (a). A radial spot branch enters the diagram in the bottom-right-hand corner and undergoes a fold at F R l0 . The radial spot solution existing on the branch segment between F R l0 and F R r0 is plotted in panel (a1). This solution is stable between F R l0 and P D6 . At P D6 on the radial branch a D6 instability produces the bifurcating branch of D6-symmetric solutions that leaves panel (a) in the bottom-left-hand corner. Beyond P D6 the radial branch is unstable and undergoes a further fold F R r0 . After another fold F R l1 a ring has formed around the radial spot. The spot with ring solution existing on the branch segment between F R l1 and F R r1 is shown in panel (a2). The branch remains unstable and undergoes a series of further folds (F R l2 , F R r2 , F R l3 , etc) adding additional rings as is shown in panel (b).
The D6-symmetric branch that bifurcates from the radial branch at P D6 also undergoes a series of fold bifurcations as shown in Fig. 4.2(b). In this case a series of additional spots are added to the central spot in a configuration that preserves the D6-symmetry. Planar plots in panels (b1), (b2) and (b3) show the stable 7spot, 19-  that can be found between F D6 r2 and F D6 l3 , there exists a 25-spot solution for which two spots have been glued on the long edge of each of the six sides of the solution shown in panel (b2).
The same radially-and D6-symmetric branches shown in the previous section have been computed for IM. Here we test the accuracy of PDE4 both qualitatively in terms of the types of solutions produced and their bifurcations, and quantitatively in terms of the parameter ranges for which the different solution types persist. We are also interested to see whether the relative ranges of existence for different types of solution is consistent between PDE4 and IM. IM . The first major point to make is that in terms of the types of solution encountered, the series of bifurcations encountered and the stability of each branch segment, there is an exact agreement between PDE4 and IM. Furthermore, the quantitative agreement on the radial branch is good up until the fold point F R r1 . Above F R r1 , the radial branch for PDE4 makes a large excursion away from the IM branch and the branches remain well separated as the snaking continues; see panel (a). Similarly, for the D6 branch the level of agreement is good up to F D6 r1 above which PDE4 branch deviates and remains well separated from the IM branch; see panel (b). The range of existence in µ for each stable branch segment on the D6 branch is greatly under estimated by PDE4.
We now highlight a key qualitative difference between the bifurcation diagrams for PDE4 and IM. In IM, there is a range of µ ∈ [2.5, 3.0] for which the stable branch segments corresponding to 7-spot, 19-spot and 37-spot all overlap. This is not the case for PDE4, in particular, the branch segments corresponding to stable 7-spot and 37-spot solutions between the fold-pairs (F D6 l1 , F D6 r1 ) and (F D6 l3 , F D6 r3 ) do not overlap.
IM r r r r j This can be seen by the fact that F D6 r3 occurs at a smaller µ-value (indicated by the first thin vertical line in Fig. 4.2(b)) than F D6 l1 (indicated by the second thin vertical line). This organisation of the solutions in parameter space is qualitatively inconsistent with IM. Figure 4.4(a) shows branches of both radially-and D6-symmetric solutions of PDE8. Globally the bifurcation diagram is the same as that of PDE4 in terms of the types of solution observed and the sequence of bifurcation encountered. There is an important difference between PDE4 and PDE8 in terms of the organisation in parameter space of the solution branches. For PDE4, there is no overlap in parameter ranges for which the 7-and 37-spot branches are stable; see the two vertical lines in Fig. 4.2 and note that F D6 r3 occurs before F D6 l1 in this case. In the PDE8 case, as shown in Fig. 4.4(a), there is an overlap in the parameter ranges as indicated by the grey shaded region. This organisation of the solution branches in parameter space is now consistent with the full model as shown in panel (c). Indeed PDE8 provides better agreement with IM; the branches remain close for both the radial and the D6 branches as we move up the snaking structure as can be seen in panels (b) and (c). We note that when compared with IM, the upper section of the radial branch occurs at smaller values of µ for PDE4 and at larger values of µ for PDE8; compare Fig. 4.2(a) with Fig. 4.4(b).

4.3.
Snaking of D2, D3 and D4 patterns in IM. In this section we discuss several patterns that possess neither radial nor D6 symmetry. For non-radial patterns discussed so far in this article, see Fig. 4.2(b1)-(b3), the individual spots lie on a regular hexagonal lattice. However, there exist an infinite number of stable configurations that conform to the same lattice spacing but without the full D6 symmetry; here we present two such examples. We also present one further example of a stable configuration that does not conform to the regular hexagonal lattice.  stability (not shown). The panels (a1)-(a5) show stable solutions on the first five full excursions in µ of the snaking structure. The existence of three-spot (see panel (a1)) and twelve-spot (see panel (a3)) solutions for PDE4 with connectivity given by Eq. (2.5) was shown in [43]. Here we have shown that these solutions exist in IM and that they form part of a larger snaking structure.

Two-spot pattern (D2).
Similarly, there is an unstable two-spot solution that connects to the trivial state u = 0 at the Turing instability (not shown). Figure 4.6(a) shows that this solution also undergoes a sequence of fold bifurcations giving rise to larger D2-symmetric patterns. We note that the spacing between the spots in these patterns still conforms to the regular hexagonal lattice. The panels (a1)-(a5) show solutions on the first five stable branch segments moving up the snaking structure; we note that the pattern (a3) is on an intermediate branch that does not make a full excursion in µ. Fig. 4.7(a2) a stable configuration that lies on a square lattice but with a spacing between the spot peaks that is double that of the hexagonal-lattice solutions encountered thus far. The solution is formed by five spots that interact at the first excitatory peak away from 0 in the connectivity function (2.3) as shown in Fig. 2.1(d). The configuration of four spots forming a square with an additional spot in the center is typically referred to as a quincunx pattern that is found, for example, on dice and dominoes. As shown in Fig. 4.6(a) these solutions exist on an isola in parameter space where other branches, see panels (a1) and (a3), are unstable. Although the pattern does not undergo snaking-type behaviour it may be possible to construct larger patterns on the double-spaced square lattice.

Discussion.
5.1. Summary. This paper explores patterns of localised activity in the neural field equation posed on the Euclidean plane with a smooth firing rate function. The choice of connectivity function is an important factor in determining whether, the localised behaviour found is restricted to individual spots, or whether multiple interacting spots can form coherent localised patterns. In [30] localised states were studied in a model with a radially-symmetric wizard-hat connectivity function describing local excitation and lateral inhibitions in the Euclidean plane. When spot solutions were tracked using numerical continuation no snaking behaviour was observed, i.e. the only steady states found consisted of a single spot. The radially symmetric connectivity function studied here and shown in Fig. 2.1(d)) features local excitation, lateral inhibition and long-range bands of excitation that decay with distance. Indeed, the distance between excitation peaks fixes a spatial scale that allows for regular spatial interactions and the formation of larger patterns of activity. In [43] it was shown that multiple-spot patterns could be obtained all with the common property of the peaks lying on a regular hexagonal lattice with spacing determined by the connectivity function. One of the main aims in the present article was to show how these solutions are connected in parameter space and how patterns with varying spatial extent grow via the mechanism of homoclinic snaking.
The results presented in [43] relied on working with an approximated connectivity function (see Fig. 2.2(a) and (b)) that allowed for solutions the full integral neural field equation to be studied in an equivalent fourth-order PDE . The initial parts of the results section are concerned with the relative agreement between solutions to the integral model and equivalent PDE formulations with approximated connectivity functions. We pursue the problem numerically by investigating the level of agreement in terms of an entire bifurcation diagram rather than comparing individual solutions at fixed parameter values. We compared both radially symmetric and D6-symmetric solution branches and found that the qualitative difference of the zero-mode in the Fourier domain for the approximated connectivity used in the fourth-order model (see Fig. 2.2(b)) led to significant discrepancies in the existence ranges of solution branches, in particular for solutions with a larger spatial extent. We demonstrated that the improved approximation of the connectivity function shown in Fig. 2.2(c) and (d), leading to an eighth-order PDE, provides a better agreement across the full bifurcation diagram. In particular, the eigth-order model captures the key feature of there being a specific parameter range in which multiple solutions coexist, each solution with a different spatial extent. It was not possible to capture this feature with the fourth-order approximation. We conclude that, although equivalent PDE formulations have proved to be a useful tool for the study of neural fields, it is important to ensure close agreement between the connectivity functions in the Fourier domain. Increasing the order of the PDEs used allows for improvements in this agreement. We believe that, while converting the integral formulation to higher-order PDEs could be useful in analytical studies, numerical calculations of these systems should be approached without resorting to PDE formulation where possible. In passing, we point out that the methodology proposed here for the integral model is applicable to inhomogeneous synaptic kernels, provided that the convolution structure of the integral is preserved. Furthermore, we point out that here we have used the standard Newton-GMRES method mainly for its simplicity, but more sophisticated choices are also possible [39].
Having investigated radial and D6 solutions with PDE formulations and compared the results with the full integral model, we proceeded to give an account of other types of solutions that, when path-followed with numerical continuation lead to patterns with different underlying symmetry properties. We worked with the full integral model and showed that patterns with D2 and D3 symmetry also give rise to snaking behaviour, generating spatial patterns with variable spatial extent. All the solutions of this type, including those shown earlier with D6-symmetry, have the common feature of the individual spots lying on a regular hexagonal lattice. Furthermore, these solutions of variable spatial extent exist within roughly the same ranges of the parameter µ. As the snaking diagram is ascended, new spots are glued to long edges of the pattern in a regular fashion. We note the existence of an arbitrary number of intermediate branches not shown in our bifurcation diagrams where symmetries can be broken via the (simultaneous) addition or subtraction of one (or more spots). We expect these solutions to lie on intermediate branches that are stable over ranges smaller than the full excursions in the main snaking structure.

5.2.
Localised patterns with input. The numerical bifurcation analysis presented in this paper opens up the possibility to investigate the spread of cortical activity with inputs. An important principle in studying the neural field equation is that weak inputs to the equations should drive the system to states that are already solutions to the underlying equations without input. The bifurcation study presented here allows us to identify the types of solution that we may expect to encounter for localised inputs and the relevant parameter regime in which they occur. Two criteria need to be satisfied when identifying a suitable operating regime, 1) the model should only produce the trivial homogeneous state u = 0 before an input is introduced and 2) when an input is introduced, the model should be driven to one of the underlying non-trivial solutions. We have shown that, for the full integral model, there is an accumulation of fold bifurcations at around µ = 2.4 representing the first point for which localised patterns can be observed. Both criteria are satisfied when the model is operated just before these fold points. Introducing a weak input the sys- tem can be driven to states that have a spatial extent corresponding to that of the input. Figure 5.1(a) shows the profile of the input used in a series of simulations that were initiated with the u = 0 state plus a small random perturbation across the entire spatial domain. Figures 5.1(b), (c) and (d) show the result of three 150 time unit simulations, each with different spatial extent for the input. In each case, the trivial state u = 0 is no longer stable and the system naturally selects one of the solutions described in the early bifurcation analysis. As the spatial extent of the input increases, the size of the pattern selected by the model increases and this is an important consequence of the corresponding solutions existing as part of a snaking structure. The computational framework presented in this article will allow for the relationship between model inputs and spatially localised patterns to be investigated in future work.

Patterns on other lattices and parametric forcing.
The multi-spot solutions described in this article all have the common feature of the activated peaks falling onto a regular hexagonal lattice. Indeed, this has been found to be the default way for the radial symmetry to be broken in pattern forming systems, notably the archetypal Swift-Hohenberg equation [46]. We also investigated whether it is possible to find other states that do not conform to the regular hexagonal lattice. In Fig. 4.7 an example of a solution with D4 symmetry was shown that consists of five spots interacting at double the standard separation between excitation peaks. We found this solution to exist on an isola and not undergo snaking so as to find larger patterns tiling the plane with the same spacing. We also attempted to converge solutions with D4 symmetry that have the regular spacing between peaks. A suitable initial condition to find such solutions is given by Eq. (A.3) chosen such that there is a depression at (x, y) = (0, 0) and the surrounding square-lattice pattern decays away from the origin, see Fig. 5.2(b). In an appropriate parameter range these states appear to converge to stable patterns. However, we found the patterns to be weakly unstable, finally converging to a pattern on a hexagonal lattice after a long transient. Figure 5.2(a) shows a time-course of the L2 norm from a simulation with the initial condition shown in panel (b). The model reaches the apparently stable D4 configuration after approximately 30 time units (panel (c)) before finally converging after 230 time units to the D2-symmetric pattern (panel (d)) that was previously identified in Fig. 4.6(a2). It would be possible to stabilise such weakly unstable solutions with the use of parametric forcing, by introducing small modulations that encourage interactions on a fixed lattice in cortical space. In the neural field equations this is typically referred to as inhomogeneous neural media; travelling waves, travelling fronts, periodic patterns and pulsating fronts have been studied with such modulations [7,20,24]. The study of localised states in this context would be an interesting future direction. Furthermore, stable quasi-periodic patterns have been obtained in the Swift-Hohenberg equation through parametric forcing [37]. The question of introducing orientation-preference tilings on square and hexagonal lattices was addressed in [4]. However, one is restricted in the number of orientations that can be equally represented with such tilings, (two and three, respectively). Parametric forcing on a quasi-periodic lattice is of particular interest in the neural field equations as this could allow for near-continuous representations of features in a model without an abstracted feature space.
6. Conclusions. The organisation in parameter space of localised structures consisting of multiple spots has been revealed for the first time in planar neural field equations. In order to find such behaviour, one must choose a connectivity function with an excitatory peak away from the origin that fixes a regular spatial scale of interactions between spots. As localised solutions are path-followed using numerical continuation we find that these structures grow in a series of fold bifurcations through the mechanism of homoclinic snaking that has been well-studied in the Swift-Hohenberg equation. A numerical strategy has been proposed to perform a numerical bifurcation analysis without resorting to a PDE formulation, but taking advantage of matrix-free Newton-Krylov nonlinear solvers combined with a pseudospectral evaluation of the right-hand side. The novel application of these methods to the neural field equations allowed for numerical continuation to be applied to the full integral form of the model. Previous studies in 2D have relied exclusively on PDE approximations of the connectivity functions; here we demonstrated that these approximations can give a very close agreement with the full integral model if a sufficiently high-order approximation is taken. The numerical schemes presented here will allow for future studies of the neural field equations to use connectivity functions defined either directly in the real domain or the Fourier domain without recourse to PDE methods, provided that the sigmoidal firing rate be smooth and that the integral formulation can be expressed as a convolution (this extends also to inhomogeneous firing rates).
The neural field studied in the present paper can be considered as a model of the visual cortex and the localised patterns studied without inputs can be related to visual hallucinations that can be localised in the visual field [54]. Furthermore, we have shown that the localised states computed in our bifurcation analysis are exactly the types of solutions selected by the model in the presence of weak inputs. The persistence of these localised structures in the presence of a model input is new. This future direction will be of particular interest for the study of localised patterns of activity that have been observed in the primary visual cortex [17] with localised visual input.
with A = 2 and L = 65; subsequent panels (c) and (d) show transient states after 200 and 300 time units, respectively. The L 2 -norm is plotted over this time course in panel (a).