The Dynamics of Neural Fields on Bounded Domains: An Interface Approach for Dirichlet Boundary Conditions

Continuum neural field equations model the large-scale spatio-temporal dynamics of interacting neurons on a cortical surface. They have been extensively studied, both analytically and numerically, on bounded as well as unbounded domains. Neural field models do not require the specification of boundary conditions. Relatively little attention has been paid to the imposition of neural activity on the boundary, or to its role in inducing patterned states. Here we redress this imbalance by studying neural field models of Amari type (posed on one- and two-dimensional bounded domains) with Dirichlet boundary conditions. The Amari model has a Heaviside nonlinearity that allows for a description of localised solutions of the neural field with an interface dynamics. We show how to generalise this reduced but exact description by deriving a normal velocity rule for an interface that encapsulates boundary effects. The linear stability analysis of localised states in the interface dynamics is used to understand how spatially extended patterns may develop in the absence and presence of boundary conditions. Theoretical results for pattern formation are shown to be in excellent agreement with simulations of the full neural field model. Furthermore, a numerical scheme for the interface dynamics is introduced and used to probe the way in which a Dirichlet boundary condition can limit the growth of labyrinthine structures.


Introduction
Neural field models are now widely recognised as a natural starting point for modelling the dynamics of cortical tissue. Since their initial inception in the 1970s by Wilson and Cowan [1,2], Amari [3,4], and Nunez [5], they have been extensively studied in idealised one-dimensional or planar settings, which are typically either infinite or isomorphic to a torus. This has facilitated both the mathematical and the numerical analyses of spatio-temporal patterns, and much has been learnt about localised states, global periodic patterns, and travelling waves. Indeed there are now a number of reviews summarising work to date, such as [6][7][8][9], and how neural field modelling has shed light on large-scale brain rhythms, geometric visual hallucinations, mechanisms for short term memory, motion perception, binocular rivalry, and anaesthesia, to list just a few of the more common application areas. For the most recent perspective on the development and use of neural field modelling we recommend the book by Coombes et al. [10], which also includes a tutorial review on the relevant mathematical methodologies (primarily drawn from functional analysis, Turing instability theory, applied nonlinear dynamics, perturbation theory, and scientific computation). This substantial body of knowledge is still expanding with further refinements of the original neural field models to include other important aspects of cortical neurobiology, including axonal delay [11], synaptic plasticity [12], and cortical folding [13], as well as rigorous mathematical results for existence and uniqueness of stationary solutions on bounded subsets of R n without regard to imposition of boundary conditions [14], and new numerical algorithms for their evolution and numerical bifurcation analysis [15,16]. Neural field models are typically expressed in the form of integro-differential equations, whose associated Cauchy problems do not require the specification of boundary conditions. The value attained by the activity variable at the boundary is determined by the initial condition and by the non-local synaptic input. However, very little work has been done on the enforcement of boundary conditions in neural fields, or on their effect on inducing patterned states. An exception to this statement is the work of Laing and Troy [17], who proposed an equivalent partial differential equation (PDE) formulation of the neural field equation. While boundary conditions must be specified in the PDE setting, they are often chosen to ensure the smooth decay of localised solutions rather than model any biophysical constraint. It is already appreciated that the continuum neural fields can be extended to include different properties that can strongly influence the spatio-temporal dynamics of waves and patterns. For example, heterogeneities may give rise to wave scattering [18] or even extinction [19]. The topic we address in this paper is to ponder the role that a boundary can have on spatio-temporal patterning. Given the historical success of analysing neural fields with a Heaviside firing rate, our first step in this direction will be taken within the so-called "Heaviside world" of Amari [20]. Amari's seminal work developed an approach for analysing localised solutions of neural field models posed on the real line, and has recently been extended to the planar case by Coombes et al. [21], albeit assuming that the synaptic connectivity can be expressed in terms of a linear combination of zeroth order modified Bessel functions of the second kind. This approach is not only able to describe localised stationary solutions, often called bumps in one dimension and spots in two dimensions, but also dynamically evolve states such as travelling pulses and their transients as well as spreading labyrinthine patterns. Since the Amari approach, in either one or two spatial dimensions, effectively tracks the boundary between a high and low state of neural activity, where the firing rate switches, we shall refer to it as an interface dynamics. Importantly it gives a reduced description of solutions to a neural field model without any approximation. On the down side the approach cannot be generalised to treat smooth firing rate functions, though simulations by many authors have shown that the behaviour of the Amari model is consistent with neural field models utilising a steep sigmoidal function. Here we show how the interface dynamics approach can be generalised to include the effects of Dirichlet boundary conditions for arbitrary choices of synaptic connectivity.
In Sect. 2 we introduce a simple scalar neural field model in the form of an integrodifferential equation defined on a finite domain, and discuss natural boundary conditions for neural tissue. Focussing on Dirichlet boundary conditions we develop the key mathematical idea in this paper. Namely that the re-formulation of the original scalar model in terms of the evolution of its gradient allows for an interface description that respects Dirichlet boundary conditions. To illustrate the effectiveness of this approach we first treat the example of localised states in a one-dimensional model in Sect. 3. This is a useful primer for the construction of an interface dynamics in a two-dimensional model, presented in Sect. 4. The first part of Sect. 4 also shows how to generalise the original treatment in [21], for infinite domains, to handle arbitrary choices of the synaptic connectivity function (removing the restriction to combinations of Bessel functions). Localised bump and spot solutions of the interface dynamics are explicitly constructed and their stability determined. In Sect. 5 we extend this approach to treat Dirichlet boundary conditions, and in Sect. 6 we show explicitly how this approach can be used to handle spots and their azimuthal instabilities. We work with standard Mexican hat synaptic connectivities, as well as piece-wise constant caricatures for which calculations simplify. All our theoretical results are found to be in excellent agreement with direct simulations of the original neural field model. We also develop a numerical scheme to evolve the interface dynamics and use this to highlight how a Dirichlet boundary condition can limit the growth of a spreading pattern arising from the azimuthal instability of a spot. Finally in Sect. 7 we discuss natural extensions of the work in this paper.

A Neural Field Model with a Boundary Condition
Although single neuron models are able to predict dynamical activity of real neurons that have a wide variety of spiking behaviour [22,23], they are not well suited to describe the behaviour of tissue on the meso-or macro-cortical scale. To a first approximation the cortex is often viewed as being built from a dense reciprocally interconnected network of cortico-cortical axonal pathways [24]. These fibres make connections within the roughly 3 mm outer layer of the cerebrum. Given the large surface area of the (folded) cortex (∼800-1500 cm 2 ) and its small depth it is sensible to view it as a two-dimensional surface. Neural field modelling, on a line or a surface, is a very well-known framework for capturing the dynamics of cortex at this coarse level of description [10]. As well as being relevant to large-scale electroencephalogram (EEG) and magnetoencephalogram (MEG) neuroimaging studies [8], the understanding of epileptic seizures [25], visual hallucinations [26,27], and neural spiral waves [28,29], they have also been used to investigate localised states linked to short term working memory in the prefrontal cortex [30,31]. In the latter regard the idealised neural field model of Amari has proven especially advantageous [32]. This was originally posed on an infinite domain, without regard to the role of boundary conditions in shaping or creating patterns. However, the neural circuits of the neocortex are adapted to many different tasks, giving rise to functionally distinct areas such as the prefrontal cortex (for problem solving), motor association cortex (for coordination of complex movement), the primary sensory cortices (for vision, hearing, somatic sensation), Wernicke's area (language comprehension), Broca's Area (speech production and articulation), etc. Thus it would seem reasonable to parcellate their functional activity by the use of appropriate boundaries and boundary conditions. Previous work by Daunizeau et al. [33] on dynamic causal modelling for evoked responses using neural field equations has used Dirichlet boundary conditions. Here we extend the standard Amari model with the inclusion of a finite domain with an imposed Dirichlet boundary condition that clamps neural activity at the boundary to a specific value. Of course other choices are possible, though this one is a natural way to enforce a functional separation between cortical areas.
The scalar neural field model that we consider is given by where Ω is a planar domain Ω ⊆ R 2 , with x ∈ Ω and t ∈ R + . The variable u represents synaptic activity and the kernel w represents anatomical connectivity. For simplicity we shall only consider the case that this depends on Euclidean distance. The nonlinear function H represents the firing rate of the tissue and will be taken to be a Heaviside so that the parameter κ is interpreted as a firing threshold. We assume that a suitable initial condition is specified for (1), and we aim to impose on the corresponding solution u(x, t) the Dirichlet boundary condition where u BC is the prescribed boundary activity. For simplicity, we treat the case of homogeneous boundary conditions. It was the essential insight of Amari that the Heaviside choice allows the explicit construction of localised states (stationary bumps and travelling pulses) on infinite domains, as well as the construction of these on finite domains without a boundary condition. Our key observation that allows the extension of the Amari approach to handle the boundary condition (2) is to expose this constraint by writing the state of the system in terms of a line integral: Here Γ (x) denotes an arbitrary path that connects a point on the boundary to the point x within its interior, and z = ∇ x u ∈ R 2 . An evolution equation for z is easily constructed by differentiation of (1) to give with u given by (3).
We shall now consider Eqs. (3) and (4) as the neural field model of choice, and in the next sections develop the extension of the Amari interface dynamics. To set the scene we first consider a one-dimensional spatial model with a focus on stationary bump solutions.

One Spatial Dimension: A Primer
Prior to describing the analysis for a two-dimensional Amari neural field model with a Dirichlet boundary condition, we first consider the more tractable one-dimensional case. This illustrates the main components of our mathematical approach, as well as delivers new results about stable boundary induced bumps.
The one-dimensional version of (3) and (4) on the finite domain [−L, L] with an imposed boundary condition takes the explicit form with Here x ∈ [−L, L], t ∈ R + , and u BC denotes a constant boundary value imposed on the left end of the interval, namely u(−L) = u BC . In passing, we note that u(L) is determined once u(−L) is fixed, and some choices of u BC will result in an even bump u(x), for which u(L) = u(−L) = u BC . We now focus on a bump solution for which R(u) = {u(x) > κ} is a finite, connected open interval. The edges of the bump x i (t), i = 1, 2, are defined by a level-set condition that takes the form We shall refer to the two bump edges as the interface, as they naturally separate regions of high and low firing activity. The differentiation of the level-set condition (7) generates a rule for the evolution of the interface according tȯ Using the second fundamental theorem of calculus we obtain an expression for the interfacial velocitieṡ where A closed form expression for z(x, t) may also be found by integrating (5) using the method of variation of parameters to give where η(t) = e −t H (t), and H is the Heaviside step function. Equations (9)-(11) determine the interface dynamics for time-dependent spatially localised bump solutions that respect the Dirichlet boundary condition.
Since it is well known that the Amari model supports a stationary bump solution when the synaptic connectivity has a Mexican hat shape we now revisit this scenario and choose where b 1 , b 2 , c > 0. Moreover, we will focus on the case that the stationary bump is symmetric about the origin. In this case demanding that the interface velocity is equal to zero requires that the numerator in (9) vanish. The formula for ψ given by (10) will also become time independent, and if we denote this by P(x) then we have where we have set x 1 = − /2 and x 2 = /2 so that the bump width is given by = Hence, the bump width is determined implicitly by the single Eq. (13), and the bump shape, q(x), is calculated from (6) as To determine the stability of the bump solution we can follow the original approach of Amari and linearise the interface dynamics around the stationary values for x i . Alternatively we can follow the Evans function approach, reviewed in [34], which considers perturbations at all values of x (rather than just at the bump edges). Here we pursue the latter approach, though it is straightforward to check that the former approach gives the same answer.
To determine the linear stability of a bump we write u(x, t) = q(x) + e λtũ (x) whereũ 1. In this case the corresponding change to z is given by z(x, t) = dq(x)/ dx + e λtz (x), wherez(x) = ∂ũ(x)/∂x. Expanding (5) to first order gives For the Dirac-delta function occurring under the integral, we can use the formal identity and integrate (16) from −L to x and useũ(−L) = 0 to obtain Here From (18) we may generate two equations for the amplitudes (ũ(x 1 ),ũ(x 2 )) by setting x = x 1 and x = x 2 . This gives a linear system of equations that we can write in the form Requiring non-trivial solutions gives a formula for the spectrum as det[A − (λ + 1)I ] = 0, which yields Hence a bump solution will be stable provided Reλ ± < 0. In Fig. 1 we plot the bump width as a function of the threshold κ for a neural field posed on a finite domain Fig. 1(A) and, for the reformulated neural field with an imposed Dirichlet boundary condition u BC = 0 Fig. 1(B), using solid (dashed) lines for stable (unstable) solutions. For the former case we recover the expected Amari result, namely that there is coexistence between two bumps, the widest of which is stable. However, when we impose a Dirichlet boundary condition, four coexisting bumps are found for sufficiently large κ, and two of these bumps are stable. In other words, the Dirichlet boundary condition induces a new stable bump, whose active region occupies a large portion of the domain.

Two Spatial Dimensions: Infinite Domain
Before discussing the extension of Sect. 3 to a finite two-dimensional domain with an imposed boundary condition it is first instructive to consider the problem posed on R 2 . An interface description for this case was originally developed in [21], albeit for a special choice of synaptic connectivity kernel. By exploiting certain properties of the modified Bessel function of the second kind it was possible to reformulate integrals over two-dimensional domains in terms of one-dimensional line integrals. This allowed the interface dynamics to be expressed solely in terms of the shape of the active region in the tissue, namely a one-dimensional closed curve. Here we extend this approach to a far more general class of synaptic connectivity kernels, which include combinations of radially symmetric Gaussian functions (12).
We consider the integro-differential equation given by (1) with Ω = R 2 . We decompose the domain Ω by writing Ω = Ω + ∪ ∂Ω + ∪ Ω − where ∂Ω + represents the level-set which separates Ω + (excited) and Ω − (quiescent) regions. These domains are given explicitly by Ω + = {x|u(x) > κ}, Ω − = {x|u(x) < κ}, and ∂Ω + = {x|u(x) = κ}. We shall assume that ∂Ω + is a closed contour (or a finite set of disconnected closed contours). In Fig. 2 we show a direct numerical simulation of the full space-time model to illustrate that a synaptic connectivity function that is a radially symmetric difference of Gaussians can support a spreading labyrinthine pattern. Similar patterns have previously been reported and discussed in [21,35] for both Heaviside and steep sigmoidal firing rate functions. A description of the numerical scheme used to evolve the full space-time model is given in Appendix 1. Differentiation of the level-set condition u(∂Ω + (t), t) = κ gives the normal velocity rule: where we have introduced the normal vector n = −∇ x u/|∇ x u| along ∂Ω + (t). We will now show that c n can be expressed solely in terms of integrals along ∂Ω + (t).
Let us first consider the denominator in (21). The temporal integration of (4), using variation of parameters, gives where η(t) = e −t H (t), z 0 (x) = ∇ x u(x, 0) denotes gradient information at t = 0, and The term ∇ x ψ in (22) can be constructed as a line integral using the integral vector identity: Thus the denominator in (21) can be expressed solely in terms of a line integral around the contour ∂Ω + (t). The representation of the numerator in (21) in terms of a line integral rather than a double integral is more challenging. In previous work we have shown that this can be achieved for the special case that the weight kernel is constructed from a linear combination of zeroth order modified Bessel functions of the second kind [21]. In Appendix 2 we show that a line-integral representation can be constructed for a far more general class of anatomical connectivity patterns, making use of the divergence theorem. Using this result the numerator of (21) can be written where  Fig. 2 and s is a parametrisation for points on the contour γ ∈ ∂Ω + . Here, Hence the normal velocity rule (21) can be expressed solely in terms of onedimensional line integrals involving the shape of the active region Ω + (which is prescribed by ∂Ω + ). This is a substantial reduction in description as compared to the full space-time model, yet is exact.
As an example of the approach above let us consider a difference of Gaussians with w(r) given by (12). A simple calculation for this choice shows that In Fig. 3 we show a numerical simulation prescribed by the interface method, with initial data equivalent to that from the full space-time simulation shown in Fig. 2. The excellent agreement between the two figures is easily observed. The full details of our numerical scheme for implementing the interface dynamics are given in Appendix 3.

Two Spatial Dimensions: Dirichlet Boundary Condition
Using the notation of Sect. 4 we now show how to extend the one-dimensional approach of Sect. 3 to develop an interface dynamics for planar Amari models on a bounded domain with Dirichlet boundary conditions. For a single active region the dynamics for z(x, t) is given by (4), which can be written succinctly as with ψ given by (23), or in terms of ∂Ω + (t), by (26). Using (3) the level-set condition is Using the identity and differentiating (30) with respect to t, we obtain the normal velocity rule Here the normal vector is given by n = −z/|z| along the contour ∂Ω + . Using (29) we may write the numerator in the normal velocity rule (32) as where ζ : ∂Ω + (t) → ∂Ω is a mapping from points on the contour ∂Ω + (t) to points on the boundary ∂Ω. Hence, using the formulae for z and ψ from Sect. 4, namely Eqs. (22), (24), and (26), then all of the terms in the normal velocity rule (32) may be expressed as onedimensional line integrals. This yields the interface dynamics for Dirichlet boundary conditions, and once again we see that it is a reduced yet exact alternative formulation to the full space-time model. In contrast to the interface dynamics on an infinite domain one needs only develop further numerical algorithms for computing the line integral in (3). The numerical method for implementing the interface dynamics can be based upon that, for an infinite domain, with a specific choice for the paths Γ defining this integral. Each of the paths Γ connects a point x in the interior of the domain to a point on the boundary, and we set ζ (∂Ω + (t)) to be the endpoint of Γ (∂Ω + (t)) (see Appendix 3 for details on the numerical scheme).
Note that we do not have to numerically integrate along this path (to determine the normal velocity), and that we need only to determine the values of ψ(x, t) at the two endpoints.  Figure 4 shows a direct numerical simulation computed using the evolution of the gradient z = ∇ x u as well as the corresponding interface dynamics. We see excellent agreement between the two approaches. The obvious advantage of the interface dynamics is that one need only evolve the shape of the active region to fully reconstruct the full space-time dynamics using (3) and (24). We see from Fig. 4 that the main effect of the Dirichlet boundary condition is to limit the spread of a labyrinthine structure and ultimately induce a highly structured stationary pattern, as expected.

Spots in a Circular Domain: Dirichlet Boundary Condition
Given the large amount of historical interest in spot solutions of neural field models on infinite domains, and those on finite domains without incorporating the role of boundary conditions [35][36][37][38][39], it is worthwhile to revisit this specific class of solutions on a finite disc with an imposed Dirichlet boundary condition. We shall consider radially symmetric synaptic connectivity kernels and a disc of radius D with a spot (circularly symmetric) solution of radius R. In this case u(r, t) = q(r) with r = |r| for all t, and q(D) = u BC , with q(R) = κ and q(r) > κ for r < R and q(r) < κ for R < r < D. We shall denote the corresponding stationary field for ψ by ψ(r), and this is conveniently constructed from (26).

Construction
An implicit equation for the radius of the bump is obtained after setting the normal velocity to zero. Using (32) and (33) this yields For the specific choice of a difference of Gaussians given by (12) we may write this in the form and Q(θ ) = √ R 2 + r 2 − 2Rr cos θ . Although (36) is in closed form it is a challenge to perform the integral analytically. Thus it is also of interest to consider synaptic connectivity kernels for which more explicit progress can be made. A case in point is that of piece-wise constant functions.
Let us first consider a top-hat connectivity defined by In this case it is easier to construct ψ(r) directly from (23) as For the top-hat shape (37) we may split the above integral as Introducing the area A + (r, σ ) as where A + (R, σ ) = κ, the self-consistent equation for a spot (34) takes the form Following the work of Herrmann et al. [40], we now show how to evaluate the integral (40) using simple geometric ideas. For example, the area A + (R, σ ) can be calculated in terms of the area of overlap of two circles, one of centre 0 and radius |r| = R, and the other of centre r and radius σ subject to the constraint r = R. Using the results from Appendix 4 we have A + (r, σ ) = A(R, φ 0 (r, σ )) + A (σ, φ 1 (r, σ ) with R > D − σ . Another natural piece-wise constant choice is the piece-wise constant Mexican hat shape given by Using a similar argument to that for the top-hat connectivity we find that with R > D − σ 1 .

Stability
The stability of spots without boundary conditions has been treated by several authors, and see [38] for a recent overview. Here we extend this approach to treat a finite domain with an imposed Dirichlet boundary condition following very similar arguments to those presented in Sect. 3.
To determine the linear stability of a spot we write u(r, t) = q(r) + e λt cos(mθ )ũ(r) whereũ 1 and m ∈ N. In this case the corresponding change to z is given by z(r, t) = ∇ r q(r) + e λt cos(mθ )z(r), wherez(r) = ∇ rũ (r). Expanding (4) to first order gives (45) where |r − r | = √ r 2 + r 2 − 2rr cos θ . Using properties of the Dirac-delta distribution we find Since the term in square brackets in (46) is radially symmetric we may integrate in the radial direction usingũ(D) = 0 to obtain Setting r = R in (47) and demanding non-trivial solutions gives an equation for the eigenvalues λ in the form Thus a spot solution will be stable provided λ m < 0 for all m ∈ N where λ m is a zero of E m (λ). Once again the choice of a piece-wise constant connectivity function considerably simplifies further calculations. For example for the top-hat function given by (37) it is simple to show that and 2π 0 dθ cos(mθ )w r − r | r =r=R = 2 where θ * is the smaller of the two roots of the equation R √ 2(1 − cos θ) = σ for θ ∈ [0, 2π). Equation (50) allows for the explicit evaluation of (48) for a piece-wise constant synaptic connectivity.
Using the above analysis we find that, for the smooth Mexican hat function, given by (12), that for large domains a wide and narrow spot can coexist for a sufficiently low value of the threshold κ. Moreover, the narrow spots are always unstable (to modes with m = 0, reflecting uniform changes of size), whilst the wider spots can develop instabilities to modes with m ≥ 2. We note that the mode with m = 1 is always expected to exist due to rotational invariance (and would give rise to a zero eigenvalue for all parameter values). This is entirely consistent with previous results for Mexican hat connectivities on domains where no boundary condition is used, as reviewed in [38]. However, on a finite size disc and with an imposed Dirichlet boundary condition further spots can be induced, with sizes commensurate to that of the radius of the disc. Both of these scenarios are summarised with the use of Fig. 5. Qualitatively similar behaviour is found for the piece-wise constant Mexican hat function given by (43) (not shown). Interestingly for the simple top-hat connectivity, given by (37), we find similar results for existence, though without azimuthal instabilities to modes m ≥ 2.

Discussion
In this paper we have revisited the seminal work of Amari on neural fields and shown how to incorporate Dirichlet boundary conditions. We have built on the previous work of Coombes et al. [21] to develop an interface dynamics approach for the evolution of closed curves defining pattern boundaries. Compared to the full space-time model with imposed Dirichlet boundary conditions the interface dynamics is reduced, yet requires no approximations. The interface framework has been illustrated in a number of settings in both one-and two-dimensions, with a focus on localised states and their instabilities. In all cases we have highlighted the excellent correspondence between results obtained from numerical simulations of the full space-time model and the interface approach. Moreover, we have also emphasised that for piece-wise constant synaptic connectivities the interface approach becomes quasi-analytical, in that many of the terms required for the computation of the normal velocity of the interface can be calculated by hand rather than have to be found numerically. For spreading patterns that may arise from the azimuthal instability of a localised spot, the main effect of a Dirichlet boundary condition has been to limit the growth of the pattern. This was entirely expected, although the precise shape of the resulting stationary pattern is of course hard to predict without simulation. However, the induction of other branches of localised states in a neural field model on a disc was more surprising, even though all near the boundary proved to be stable. It should be noted that the imposition of different boundary conditions may effect the spatio-temporal evolution of a pattern and the conditions for its dynamics instability. For the sake of computational simplicity, the value attained by the activity variable at the boundary was chosen to be a constant (u BC = 0) throughout this paper. However, a full analysis which also treats space-and time-dependent boundary conditions can be readily developed for the direct numerical simulations, as well as for the equivalent interface description. There are a number of natural extensions of the approach that we have presented here to treat other, more biophysically rich, neural field models which we outline below.
Although we have focussed exclusively on the Amari model with a Heaviside firing rate, direct numerical simulations (not shown) readily confirm that boundary induced patterns can also be seen in models with a smooth sigmoidal firing rate. Thus it would also be of interest to extend the elegant functional analytic treatment of localised states in bounded domains by Faugeras et al. [14] to incorporate imposed boundary conditions. Some form of spike frequency adaptation (SFA) is often included in neural field models to mimic a negative feedback process to diminish sustained firing. This can cause a travelling front to transition to a travelling pulse [41], or subserve the generation of planar spiral waves [28]. If this current is linear, as is often the case [42], or itself described by dynamics involving a Heaviside switch, as in [43], then the interface approach presented here can be generalised. Given previous work on equivalent PDE models on bounded domains with SFA that analyses spiral wave behaviour, the treatment of spiral waves from an interface perspective would be an advance as it is not limited to synaptic connectivities with a rational Fourier structure [44]. Another natural extension of the work in this paper is to neural fields on feature spaces. For example, in the primary visual cortex (V1), cells respond preferentially to lines and edges of a particular orientation. A standard neural field model that links points at r and r (in the plane) with a weight w(r|r ), should be replaced by a more general form such as w(r|r ) = w(r, θ|r , θ ), where θ (θ ) would represent an orientation preference at r (r ). This model has recently been studied using a neural field dynamics with a Heaviside firing rate [45], and is thus ripe for a further analysis using an interface approach. Finally it is worth pointing out the rather pertinent difference between the flat models we have discussed here and the well-known folded characteristic of real cortex, with its sulci and gyri. Fortunately there is no substantial difficulty in formulating neural field models on curved surfaces, though to date there has been surprisingly little analysis of spatiotemporal pattern formation in this context. The exception to this rule is the simulation studies of Bojak et al. [46], and the recent work of Sato et al. for growing brains [13].
One obvious caveat to all of the above is that the interface approach is restricted to Amari style models with a Heaviside firing rate. Nonetheless the qualitative similarities between Amari models and those with a steep sigmoidal firing rate are well known. In summary the treatment of neural fields with boundary conditions is a relatively unexplored area of mathematical neuroscience whose further study should pay dividends for the understanding of neuroimaging data, and in particular waves of activity in functionally identified and folded cortices.

Appendix 1: Numerical Scheme for the Full Space-Time Model
The numerical simulation of the full space-time model (1) without boundary conditions was performed by discretising the domain on an N -by-N tensor grid, and using a Nyström scheme for the spatial discretisation [47]. It is known that the major costs of this scheme is in the evaluation of the Hammerstein operator occurring on the righthand side of (1), which on the grid outlined above requires N 4 operations. Owing to the convolutional structure of the operator, it is possible to decrease considerably the computational cost of each function evaluation by performing a pseudo-spectral evaluation of the convolution, using a Fast Fourier Transform (FFT), followed by an inverse Fast Fourier Transform (IFFT). This reduces the number of operations to O(N 2 log N 2 ) and allows one to simulate neural fields and compute equilibria efficiently. We refer the reader to [16,21] for further details. In these calculations we set N = 2 9 and used Matlab's in-built ode45 routine, with standard tolerance settings.
In simulations where the state variable is not periodic, such as the ones where we enforced boundary conditions (Sect. 5) we used a standard matrix-vector multiplication to evaluate the integral operator. A full matrix was precomputed and stored during the initialisation phase of the time stepping and the grid-size was limited to N = 2 7 points, owing to memory constraints. In this setting, the interface dynamic approach becomes a viable alternative to the full spatial simulation.

Appendix 2: Expressing ψ in Terms of Contour Integrals
In this appendix, we derive the identities (26) and (27). This allows us to represent the double integral for the non-local input ψ(x, t) given by (23) as an equivalent line integral. We recall divergence theorem for a generic vector field F on a domain B with boundary ∂B, where n is the unit normal vector on ∂B. We consider a rotationally symmetric twodimensional synaptic weight kernel w(x) = w(r) which satisfies R 2 dxw(x) = K, for some finite constant K, and we introduce a function g(x) : R 2 → R such that Now considering a function ϕ(r) : R + → R which satisfies the condition lim r→∞ rϕ(r) = 0, the vector field can be written using polar coordinates, that is, F = ϕ(r)(cos θ, sin θ) = ϕ(r)x/|x| with x = r(cos θ, sin θ). Transforming the expressions K and g into polar coordinates, integrating Eq. (52), and using the divergence theorem, yields where the line integral is described over a circle of radius R → ∞. Therefore, the weight kernel can be written in the form The integration of (56) yields Using the above results means that (23) can be evaluated as Here γ ∈ ∂Ω + , and the integration over the Dirac-delta function gives C = 1 if x is within Ω + , C = 0 if x is outside Ω + , and C = 1/2 if x is on the boundary of Ω + .

Appendix 3: Numerical Scheme for the Interface Dynamics
Time stepping for interface dynamics requires a novel integration scheme. We present here the implementation used in our numerical experiments, which were found to be in agreement with the full spatio-temporal simulation, and we defer a numerical analytical study of its properties to a later date. The method is formed of four constitutive parts: a scheme for approximating a closed curve (the interface ∂Ω + (t)), a scheme to approximate the instantaneous normal velocity of the interface, a scheme to propagate the contour according to the normal velocity, and a strategy to remesh or postprocess the contour, if needed.
Closed contours: We chose a periodic parametrisation where ξ 1 , ξ 2 are smooth and 2π -periodic in s for all t, and we approximated ξ 1 (s, ·), ξ 2 (s, ·) spectrally, using evenly spaced points in s. Using FFTs, we could approximate quickly and accurately the normal and tangent vectors to ∂Ω + at each point s, and each time t. We also parametrise the domain boundary ∂Ω as follows: where η i are continuous and 2π -periodic. We choose the paths Γ to be straight lines connecting x to its closest point on the boundary, and we define ζ using the endpoints of Γ , Normal velocity: To simplify the discussion, let us consider the case where boundary conditions are not imposed. We compute the normal velocity by approximating the numerator and denominator of (21). The denominator is updated at each time step using (22) and (24): this requires the normal to ∂Ω + at time t, as well as the full history of ∂Ω + (t) in the interval [0, t], in view of the integral (22). However, since η(t) = e −t H (t), we found that retaining in memory and using only the last 20-50 time steps was not detrimental for the accuracy of the solution; we use the trapezium rule for both the integration along the contour in (24) (which is therefore spectrally accurate) and for the integral over t in (22). The numerator of (21) is computed using (26), for which a further integration along ∂Ω + is performed. The integrand is singular, and can be treated as in [48]. A similar strategy is used for (32) and (33), in the case of a bounded domain with Dirichlet boundary conditions. In all cases this step is by far the most time consuming of the algorithm, due to the large number of integrals which need to be evaluated at each step.
Position update: The contour is propagated in the normal direction, using the velocity computed at each point of the contour, c n (s, t). To this end we use a simple Euler update ∂Ω + (t + t) = ∂Ω + (t) + c n t to find the new contour, given that c n is also computed with O( t) accuracy in time. Other choices are obviously possible, but require more function evaluations and more expensive quadrature rules. A stepsize of 0.05 or less has been used in our simulations.
Remeshing and postprocessing: The updated contour ∂Ω + (t + t) leads to a new parametrisation, that is, to an update of the functions ξ 1 (s, ·), ξ 2 (s, ·). Since we need a uniform distribution of the nodes with respect to the variable s, we redistribute points using standard interpolation [49]. As the pattern grows or shrinks, points are added or removed so as to keep the arclength between consecutive points approximately constant.

Appendix 4: Geometric Formulae for a Piece-Wise Constant Kernel
Consider a portion of a disk whose upper boundary is a (circular) arc and whose lower boundary is a chord making a central angle φ 0 < π, illustrated as the shaded region in Fig. 6(A).

(61)
Acknowledgements AG acknowledges the support from The University of Nottingham and Ministry of National Education in Turkey. We thank the anonymous referees for their helpful comments on our manuscript.
Funding SC was supported by the European Commission through the FP7 Marie Curie Initial Training Network 289146, NETT: Neural Engineering Transformative Technologies.
Availability of data and materials Please contact author for data requests.

Competing interests
The authors declare that they have no competing interests.
Authors' contributions AG, DA and SC contributed equally. All authors read and approved the final manuscript.

Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.