Neural ‘Bubble’ Dynamics Revisited

In this paper, we revisit the work of John G Taylor on neural ‘bubble’ dynamics in two-dimensional neural field models. This builds on original work of Amari in a one-dimensional setting and makes use of the fact that mathematical treatments are much simpler when the firing rate function is chosen to be a Heaviside. In this case, the dynamics of an excited or active region, defining a ‘bubble’, reduce to the dynamics of the boundary. The focus of John’s work was on the properties of radially symmetric ‘bubbles’, including existence and radial stability, with applications to the theory of topographic map formation in self-organising neural networks. As well as reviewing John’s work in this area, we also include some recent results that treat more general classes of perturbations.

PCB I first met John in 1985 when he interviewed me for a PhD. position in string theory at King's College London. I had drifted away from academia at the time and was thus grateful that he decided to offer me the position. I was immediately struck by his sharpness, his larger-thanlife personality, and his enthusiasm for mathematics and science. Although that enthusiasm would sometimes lead John to follow rather controversial paths, it was inspiring for a young graduate student. When I started my PhD, John had a large group of students working on horrendous perturbation calculations in supergravity. I still remember the long line of students disappearing down the corridor waiting to discuss their latest results. However, all of these students graduated during my first year so I was essentially John's only student, which meant that we worked closely together. My thesis work formed part of a book on 'Finite Superstrings' [49] co-authored with John and one of his long-standing collaborators Alvarez Restuccia from Venezuela. After graduating, I joined the Long-range Research Lab at the Hirst Research Center in London, which carried out fundamental research for GEC-Marconi. I was asked to develop a research programme in neural networks, which coincidentally was a research area that John had embraced in the 1970s and returned to in the late 1980s. John became a research consultant for GEC-Marconi and this enabled us to continue collaborating. Unfortunately, John and I lost contact when I joined Loughborough University in 1993. In fact, I only saw John one more time, which was as a guest speaker at his 'retirement' colloquium. Of course, John remained as active as ever following his retirement, in particular, pursuing his lifetime fascination with the problem of consciousness.

Introduction
In this paper, we revisit a problem that John considered in 1999 [48], namely the existence and stability of radially symmetric spots in two-dimensional (2D) neural fields, also known as 'neural bubble dynamics'. Neural fields represent the large-scale dynamics of spatially structured networks of neurons in terms of nonlinear integrodifferential equations, whose associated integral kernels represent the spatial distribution of neuronal synaptic connections. One type of solution that emerges in the presence of nonlocal synaptic inhibition (lateral inhibition) is a stationary spot solution, also known as an activity bump. Such bumps are typically coexistent with a stable low-activity state (bistability) so that an initial stimulus is needed in order to transition from the low-activity state to the bump. However, the bump persists after removal of the stimulus, so that the bump represents a persistent spatially localised activity state [12,24,33]. One of the reasons why persistent activity bumps are of interest is that they are thought to arise in cortical circuits performing certain spatial working memory tasks. Working memory involves cortical ''memory neurons'' that establish a representation of a stimulus that persists after the stimulus is removed. A typical experiment is a delayed response task, in which a primate is required to retain information regarding the location of a sensory cue across a delay period between the stimulus and behavioural response. Physiological recordings in prefrontal cortex have shown that spatially localised groups of neurons fire during the recall task and then stop firing once the task has finished [22]. The stimulus response of a cell is characterised by a smooth tuning curve that is peaked at a preferred spatial cue and varies from cell to cell. At the network level, the memory of cue location is stored as an activity bump. Persistent activity bumps occur in a number of other systems that encode directional or spatial information, including head direction cells in thalamus and basal ganglia [46] and place cells in the hippocampus [38].
Prior to John's study [48], almost all analyses of bumps had been restricted to 1D networks. Wilson and Cowan established the existence of 1D bumps numerically [53,54], and Amari subsequently constructed an exact bump solution for a 1D scalar neural field equation with a Heaviside firing rate function [1]. He also showed how the stability of the bump depends on whether or not perturbations of the bump boundary (threshold crossing points) grow or decay. Interestingly, Amari also developed an analysis of 2D bumps. However, this work only appeared in a 1978 book that was never translated from the Japanese [3] and has only just been summarised in English [5]. John generalised Amari's 1D analysis by deriving conditions for the existence of radially symmetric 2D bumps and by determining the stability of the bumps with respect to radially symmetric perturbations of the circular bump boundary [48], see also [50]. This analysis was later extended to the case of nonradially symmetric perturbations using Fourier methods and properties of Bessel functions [16,20,21,25,39]. In related work, Laing and Troy [32] introduced PDE methods to study symmetry breaking of rotationally symmetric bumps and the formation of multiple bump solutions. However, such PDE methods can only be applied to specific forms of weight distribution. In particular, they break down if the weight distribution has compact support.
In section 'Dynamics of a 'Bubble' Boundary', we review the work of John (in an appropriately revised notation) on 2D bumps, as well as describe how to treat azimuthal instabilities, and the generation of breathing bumps in neural field models with adaptation and external drive. Next in section 'Self-Organising Neural Field Theory' we cover John's work on how 2D bump activity can drive topographic map formation in self-organising neural networks, also presented in [48], and describe a further set of techniques for treating nonradially symmetric perturbations. We end with a brief discussion in section 'Discussion'.

Dynamics of a 'Bubble' Boundary
Consider the following neural field equation for a 2D sheet of neural tissue s ouðr; tÞ ot ¼ Àuðr; tÞ þ Z R 2 wðjr À r 0 jÞf ðuðr 0 ; tÞÞdr 0 ð1Þ where r = (r, h) and r 0 = (r 0 , h 0 ). The neural field u(r, t) represents the local activity of a population of neurons at position r, s is a synaptic or membrane time constant (which we set to unity), w is the distribution of synaptic weights, and f denotes an output firing rate function. The weight distribution is assumed to depend on the Euclidean distance |r -r 0 | between interacting neurons at r and r 0 . A common choice for the firing rate function is a bounded, monotonic function such as a sigmoid. Here we follow Amari [1] and Taylor [48] and take this to be Heaviside function with threshold j such that f(u) = H(u -j). Below we review John's extension of Amari's 1D analysis to radially symmetric localised solutions in 2D. For a more detailed discussion of methods for analysing the existence and stability of bumps in neural fields, see recent reviews by the authors [9,14]. First let us rewrite the original model (1) in the form u t = -u ? w, where wðr; tÞ ¼ Z Bðr 0 ;tÞ dr 0 wðjr À r 0 jÞ; ð2Þ and Bðr; tÞ ¼ frjuðr; tÞ ! jg (which defines the 'bubble'). For radially symmetric spot solutions of radius DðtÞ then we have that w(r, t) = w(r, t) with Here DðtÞ is defined by the level set condition uðDðtÞ; tÞ ¼ j. Time-independent spot solutions such that DðtÞ ¼ D (namely with a constant radius) such that u(r, t) = U(r) with lim r!1 UðrÞ ¼ 0 and UðrÞ?j for r7D are given by Differentiating the level set condition with respect to time gives an equation for the velocity of the spot boundary in the form Using (1) we may derive an ODE for v ¼ ouðr; tÞ=orj r¼D as Hence, we may generate a system of two exact nonlinear ODEs for ðD; vÞ to describe the evolution of the (radially symmetric) spot: where The steady state of (7) is determined by the solution of WðDÞ ¼ j and v ¼ W 0 ðDÞ. A linear stability analysis around the steady-state solution shows that there are two eigenvalues, one with value -1 (that can not lead to any instability) and the other with a value k, where Hence, only the solution with W 0 ðDÞ\0 is stable (after realising that v is the gradient of the spot at its boundary and is negative).
In the 'Appendix', we show that U(r) can be written in the computationally useful form where b w is the 2D Fourier transform of w. For the sake of illustration, consider a wizard hat weight distribution given by a combination of modified Bessel functions of the second kind As shown in Fig. 1, this can generate a weight distribution with short-range excitation and long-range inhibition.
Using the fact that the corresponding 2D Fourier transform of K 0 (s r) is Hðq; sÞ ¼ ðq 2 þ s 2 Þ À1 , we have b wðqÞ ¼ 2 3p ðHðq; 1Þ À Hðq; 2Þ À AðHðq; 1=rÞ À Hðq; 2=rÞÞÞ: Thus, the integral (10) can be evaluated explicitly using the identity where I m is the modified Bessel function of the first kind of order m. Thus, the stationary bump U(r) has the form IðD; r; 1Þ À IðD; r; 2Þ À AðI ðD; r; 1=rÞ ½ À IðD; r; 2=rÞÞ: The bump radius may then be computed by finding the roots D of the equation j ¼ WðDÞ with ÀAðrI 1 ðD=rÞK 0 ðD=rÞ À r 2 I 1 ð2D=rÞK 0 ð2DrÞÞ : ð15Þ Note that the threshold condition is a necessary but not sufficient condition for proving existence of a 2D bump. One also has to check that there are no other threshold crossing points. In the case of a Mexican or wizard hat weight distribution, there is typically a maximum of two spot solutions as illustrated in Fig. 2 for w given by Eq. (11). Using (9), we find that the narrower spot is always unstable as found in 1D. However, the stable upper branch can develop instabilities as the threshold is decreased leading to the formation of solutions that break the rotational symmetry [32,39], see below.

Azimuthal Instabilities
It can be misleading to extrapolate results about bumps in 1D to spots in 2D. For example, it was known to Amari [1] in his seminal work on 1D models with Mexican hat weight distribution that bump solutions come in pairs, and that it is only the wider of the two that is stable. Extrapolating this to 2D, and focusing on radially symmetric spots, one would conclude a similar state of affairs for radial perturbations to the stationary spot, as done by John [48]. However, in the 2D setting one must also pay careful attention to azimuthal instabilities. This has been appreciated by a number of authors in recent years [20,32,39], though anticipated much earlier by Amari [3]. Indeed azimuthal instabilities can destabilise spots on the wide branch (that are stable to radial perturbations) and lead to the generation of intricate labyrinthine structures [16,39].
In order to determine linear stability of a bump solution U(r), substitute u(r, t) = U(r) ? p (r)e kt into Eq. (1) and expand to first order in p using equation (4). This gives the eigenvalue equation where a 0 = (a, / 0 ). We can now formulate the stability problem in terms of finding the spectrum of a compact linear operator acting on continuous, bounded functions w(r, /) defined on the disc of radius r D. One class of solution to Eq. (16) consists of functions p(r) that vanish on the boundary, w(a,/) = 0 for all /, such that k = -1. This belongs to the essential spectrum, which does not contribute to any instabilities. The discrete spectrum is determined by setting r ¼ a ðD; /Þ in Eq. (16): where we have simplified the argument of w(r) using Equation (17) can be solved in terms of Fourier eigenmodes, that is, pðD; /Þ ¼ P n ð/Þ ¼ c n e in/ þ c n e Àin/ with corresponding eigenvalue k n satisfying Note that k n is real since (after rescaling /) i.e., the integrand is odd-symmetric about p/2. Hence,  Fig. 2. It can be seen that the nth-order boundary perturbation has D n symmetry, meaning the resulting solution has the n reflectional and rotational symmetries of the dihedral group D n . The perturbations also have a simple geometric interpretation. For example, n = 0 corresponds to a uniform expansion or contraction of the spot, whereas n = 1 corresponds to a uniform shift of the spot. Since the n = 1 mode represents pure shifts of the spot solution, we expect k 1 = 0 from translation symmetry. In order to verify this, we evaluate the integral appearing in equation (21) using Bessel functions, along similar lines to the evaluation of U(r), Eq. (10). That is, Moreover, differentiating Eq. (10) with respect to r gives U 0 ðDÞ ¼ À2pD Hence, the eigenvalue (21) can be rewritten as It immediately follows that k 1 = 0. Hence, the 2D spot is linearly stable if k n \ 0 for all n = 1. For the weight distribution (11), the points of azimuthal stability as determined by the conditions k n = 0 can be calculated as [16]: where (A 1 , A 2 , A 3 , A 4 , a 1 , a 2 , a 3 , a 4 ) = (1, -1, -A, A, 1, 2, 1/r, 2/r). In Fig. 2 we plot the set of points for n = 0, 2 , …, 9 as determined by equation (27) to illustrate how the upper branch of rotationally symmetric spots can become unstable as the threshold is decreased, leading to the formation of solutions that breaks continuous rotational symmetry. The discrete rotational symmetry D n of a bifurcating solution reflects the order n of the dominant eigenvalue k n at bifurcation. Interestingly, if linear adaptation is included in the neural field model, then these nonrotationally symmetric solutions can undergo a secondary instability leading to the formation of a rotating wave [32,39]. Sufficiently strong adaptation can also destabilise a bump leading to a travelling spot [15,16]. An example of the shape of an expanding labyrinthine structure that can emerge when a spot destabilises to an n = 3 mode is shown in Fig. 3.

Stimulus-Induced Breathers
So far we have focused on activity bumps that persist in the absence of external stimuli due to the combined action of local recurrent excitation and lateral inhibition. We now describe some interesting instabilities that arise in the case of nonpersistent bumps [20,21]. For the sake of illustration, consider a 2D excitatory neural field with linear adaptation and an external input I: wðjr À r 0 jÞHðuðr 0 ; tÞ À jÞdr 0 À bvðr; tÞ þ IðrÞ 1 ovðr; tÞ ot ¼ Àvðr; tÞ þ uðr; tÞ; ð28Þ where v(r, t) represents a local negative feedback mechanism, such as spike-rate adaptation, with e, b determining the relative rate and strength of feedback [41]. Suppose that the inhomogeneous input is a radially symmetric Gaussian centred about the origin, IðrÞ ¼ I 0 e Àr 2 =r 2 s . In the absence of an input, the resulting excitatory network supports travelling waves rather than stationary bumps. On the other hand, for sufficiently strong input amplitude I 0 , the network supports a radially symmetric bump centred about the input. Such a bump is not persistent, since if the input is removed then the bump disappears as well. The basic problem we want to address is what happens to the stability of the bump as the input amplitude is slowly decreased.
The analysis of the existence and stability of radially symmetric 2D bumps proceeds as above with minor changes. First, the threshold condition for the existence of a bump becomes Second, the linear stability of the bump is determined by the pair of eigenvalues k = k n ± associated with the Fourier modes [c n e in/ ? c n e -in/ ]e k t , where [20] k AE n ¼ and l n ðDÞ ¼ D It follows that stability of a radially symmetric bump require K n [ 0 and C n \1 for all n C 0. Given the form of K n , this reduces to the following stability conditions: If the first condition is violated as some parameter is varied, then there is a saddle-node bifurcation, whereas a breakdown of the second condition signals a Hopf bifurcation. In the latter case, the bump instability leads to the formation of a breather. In the particular case of an excitatory network (such as obtained by setting A = 0 in Eq. (11)), with w(r) C 0 for all r C 0, we have so that l n B l 0 for all n. Since l 1 = 0, it follows that l 0 C 0 and, hence, a purely excitatory neural field cannot support stable radially symmetric bumps. In this case, we expect any instability to involve the growth of radially symmetric perturbations, and hence, the resulting breather will be radially symmetric. On the other hand, if there is a Mexican hat weight distribution, then nonradially symmetric breather and rotating waves can occur [21,39]. One way to induce a Hopf instability of a bump is to reduce the amplitude I 0 of the Gaussian input; this modifies both the pulse-width D and the slope of the bump at threshold, jU 0 ðDÞj. Interestingly, as the input amplitude is further reduced, the breather can undergo a Here we show the shape of such a structure at a fixed moment in time. As the pattern evolves, it fills more of the spatial domain. The black lines denote the borders between high and low firing regions. This example illustrates the type of pattern that one can expect to see beyond an n = 3 azimuthal instability. For further examples and discussion, see [16] secondary instability such that it now acts as an oscillating core that emits circular target waves. An example of such a periodic wave emitter is shown in Fig. 4. Thus, a spatially localised stationary input provides a mechanism for the formation of a network pacemaker oscillator. For a recent discussion of breathers in the absence of localised input, see [15].

Self-Organising Neural Field Theory
One of the striking features of the visual system is that the visual world is mapped on to the cortical surface in a topographic manner. This means that neighbouring points in a visual image evoke activity in neighbouring regions of visual cortex. Superimposed upon this topographic map are additional maps reflecting the fact that neurons respond preferentially to stimuli with particular features such as ocular dominance and orientation [27]. In recent years, much information has accumulated regarding the twodimensional distribution of both orientation preference and ocular dominance columns using optical imaging techniques [6,7]. These experimental studies indicate that there is an underlying periodicity in the microstructure of V1 with a period of approximately 1 mm (in cats and primates). The fundamental domain of this tiling of the cortical plane is the hypercolumn, which contains two sets of orientation preferences h [ [0,p) per eye, organised around a pair of orientation singularities or pinwheels [37]. It is generally accepted that the preference of cortical neurons for particular stimulus features such as orientation and ocular dominance arises primarily from the spatial arrangement of convergent feedforward afferents from the LGN (or from other layers of cortex). Although molecular cues and axon guidance are involved in the initial stages of cortical map formation, the resulting maps are rather crude and some form of spontaneous or stimulus-driven activity appears to be necessary for the subsequent refinement of these maps through the pruning of initially exuberant axonal arborisations.
A large number of models have been proposed that describe activity-dependent development as a self-organising Hebbian process (see the review of Swindale [44]). In the case of correlation-based developmental models [35], the statistical structure of input correlations provides a mechanism for spontaneously breaking some underlying symmetry of the neuronal receptive fields leading to the emergence of feature selectivity. When such correlations are combined with intracortical interactions, there is a simultaneous breaking of translation symmetry across cortex leading to the formation of a spatially periodic cortical feature map. Correlation-based models are essentially linear, so that considerable insight into the developmental process can be obtained by solving an associated eigenvalue problem [36]. One of the possible limitations of this class of model is that a regular topographic map is assumed already to exist before feature-based columns begin to develop. In order to model the joint development of topography and cortical feature maps, it appears necessary to introduce some form of nonlinear competition for activation [30,52], neurotrophic factors [17], or a combination of the two [51].
An alternative mathematical formulation of topographic map formation has been developed by Amari using the theory of self-organising neural fields [2,4,45]. At the simplest level, this model can be formulated in terms of a two-layer network corresponding respectively to the lateral geniculate nucleus (LGN) of the thalamus and the primary visual cortex (see Fig. 5). The cortical or output layer behaves pretty much as in the single-layer model given by equation (1). However, the inputs I are now determined by activity in the LGN or input layer together with the strength of the feedforward connections from the LGN to cortex. These connections are modifiable by experience via Hebbian learning. Initially, an activity pattern in the form of a bump occurs at some random location within the LGN and elicits a corresponding activity bump in cortex whose location depends on the given pattern of feedforward connections. Over timescales much slower than the relaxation time s for cortical dynamics, inputs arising from many different locations within LGN (and hence the visual image) occur. The resulting set of cortical responses to this set of inputs then Fig. 4 Sequence of snapshots of a 2D breather acting as a periodic pulse emitter in the case of a 2D excitatory neural field with linear adaptation and exponential weight function w(r) = e -r /(2 p). Parameters are b ¼ 4; j ¼ 0:2; ¼ 0:1 and I 0 = 0.2. Lighter colours indicate higher activity [20] determines how the feedforward weights are modified. In the one-dimensional case, Amari [45] derived conditions for the existence and stability of a continuous topographic map between the two layers. Moreover, he showed that under certain circumstances, the continuous topographic map can undergo a pattern-forming instability that spontaneously breaks continuous translation symmetry, and the map becomes partitioned into discretised blocks; it has been suggested that these blocks could be a precursor for the columnar microstructure of cortex [45]. Given that cortical columns tend to be associated with stimulus features such as ocular dominance and orientation, this raises the interesting question whether or not such features could also emerge through the spontaneous symmetry breaking of self-organising neural fields. In the same paper on bubble dynamics [48], John extended Amari's theory to 2D networks, again focusing on existence and stability with respect to radially symmetric perturbations. Stability with respect to a more general class of perturbations was subsequently analysed by Bressloff [8], who showed that stimulus preference maps could indeed emerge through the spontaneous symmetry breaking of self-organising neural fields. This analytical result was consistent with a number of numerical studies [19,56]. In this section, we outline the basic steps in the analysis of 2D self-organising neural fields.

Existence and Stability of Bumps
Let s 2 X 2 denote a point in the LGN layer and r 2 X 1 a point in the cortical layer, as shown in Fig. 5. Each point s labels a distinct input pattern which we denote by IðrjsÞ; r 2 X 1 . Let u(r|s) denote the corresponding cortical activity induced by the input I(r|s), which is assumed to converge to a stable equilibrium given by the solution of the integral equation wðjr À r 0 jÞHðuðr 0 jsÞ À jÞdr 0 þ IðrjsÞ: Now suppose that once the network has reached equilibrium, a new input pattern is presented to the network from a different point s 0 , and that this process is iterated multiple times so that the network samples over the space of inputs. How the network responds to the given set of input patterns then determines how these input patterns are themselves modified by experience. Adaptation is introduced into the model by assuming that the input patterns evolve according to the dynamical equation gðjs À s 0 jÞHðuðrjs 0 ; tÞ À jÞds 0 : Here g ) s where s is the relaxation time for cortical dynamics. Hence, uðrjs; tÞ; r 2 X 1 denotes the equilibrium cortical state given by Eq. (35) in response to the input pattern I(r|s, t) for fixed s, t. The integral kernel g incorporates the effect of statistical correlations between activity patterns at different points in LGN. Thus, increasing the strength of an input from s also leads to an increase in the strength of an input at s 0 provided that s and s 0 are sufficiently close. Inputs from well-separated points in LGN are anticorrelated. These features can be included by taking g to be of the form for constants c, ĉ. Equations (35) and (36) are the basic neural field equations for topographic map formation [2,45]. Rescaling the LGN and cortical coordinates appropriately, that is, ignoring the effects of cortical magnification, one can then look for homogeneous steady-state solutions of the form u(r|s) = U(|r -s |) and IðrjsÞ ¼ Iðjr À sjÞ, where U is a radially symmetric stationary pulse solution of width a such that wðjr À r 0 jÞHðUðjr 0 À sjÞ À jÞdr 0 þ Z X 2 gðjs À s 0 jÞHðUðjr À s 0 jÞ À jÞds 0 ¼ Z 2p 0 Z a 0 ½wðjr À r 0 jÞ þ gðjr À r 0 jÞr 0 dr 0 dh; ð38Þ cortex r'

I(r|s)
LGN s r I(r'|s) gðjs À s 0 jÞHðUðjr À s 0 jÞ À jÞds 0 : Such a solution represents a continuous topographic map in which an input from s 2 X 2 is mapped to the centre of cortical output activity at the corresponding point s 2 X 1 . Stability of the bump (on the fast time-scale of cortical dynamics) can be analysed along identical lines to section 'Azimuthal Instabilities', after setting b = 0 and taking w ? w ? g. The corresponding eigenvalues are k n ¼ À1þ c À1 m n with |U 0 (a)| = c and Differentiating Eq. (4) with respect to r and setting r = a shows that m 1 = c and, hence, k 1 = 0. The existence of a zero eigenvalue reflects the underlying translation symmetry of the two-layer network, which implies that the bump is marginally stable with respect to uniform shifts in space. Thus, the bump will be stable provided that m n \ c for all n = 1. In the following, we assume that there exists a unique stable bump. We now investigate the stability of the two-dimensional topographic map (on the slow timescale of the input or weights dynamics) by linearising Eqs. (35) and (36) about the homogeneous radially symmetric solution given by Eqs. (38) and (39). First, introducing the perturbations uðrjs; tÞ ¼ Uðjr À sjÞ þ pðrjs; tÞ; Iðrjs; tÞ ¼ Iðjr À sjÞ þ qðrjs; tÞ; ð41Þ and expanding to first order in p, q leads to the linear equations (on setting g = 1) gðjs À s 0 jÞH 0 ðUðjr À s 0 jÞpðrjs 0 ; tÞds 0 ; ð42Þ and pðrjs; tÞ ¼ qðrjs; tÞ þ Z R 2 wðjr À r 0 jÞH 0 ðUðjr 0 À sjÞ Â pðr 0 js; tÞdr 0 : The neural field perturbations p h (r, s) can be related to perturbations of the circular boundary of the activity bump. Let us write the perturbed threshold condition in the form u(r|s, t) = j at r = s ? (a ? q(h, n, t))e h . This yields j ¼ Uða þ qðh; n; tÞÞ þ pðs þ ða þ qðh; s; tÞÞe h js; tÞ to the boundary perturbation considered in section 'Azimuthal Instabilities'. Equations (48) and (50) Equations (54) and (55) can be analysed further by introducing the Fourier series This leads to the discrete set of equations

Calculation of Eigenmodes
Determining the stability of the two-dimensional topographic map has reduced to the problem of finding the eigenvalues of the infinite-dimensional matrix G nn 0 (k) for n; n 0 2 Z. It is not possible to do this analytically for general input kernel g. However, an explicit solution can be obtained in the limiting case of wide Gaussian inputs such that r ) a in Eq. (37). We can then carry out a perturbation expansion in a 2 /r 2 by writing gð2a sinðh=2ÞÞ ¼ ce Àa 2 sin 2 ðh=2Þ=r 2 Àĉ % c Àĉ À ca 2 2r 2 ð1 À cosðhÞÞ þ Oða 4 =r 4 Þ: ð61Þ Keeping only lowest order terms, we find that G nn 0 ðkÞ % aðc ÀĉÞ Similarly, substituting Eq. (61) into (40) and gives to lowest order l 0 % m 0 À 2paðc ÀĉÞ; l 1 % m 1 À pac a 2 2r 2 ; l n % m n ; n [ 1: There are two classes of solution to Eq. (66). If b Á b P ¼ 0 then k = -1 and the topographic map is stable with respect to excitation of the corresponding eigenmodes. On the other hand, if b Á b P 6 ¼ 0 then b P ¼ b Ã (up to a constant multiplicative factor). Substituting into the Fourier series (56), the resulting eigenmode is of the form P h ðkÞ ¼ Pðk; h À uÞ with k ¼ ðk; uÞ, ðÀ1Þ n J 2n ðkaÞ c À l 2n cosð2nhÞ " À 2 X n ! 1 ðÀ1Þ n J 2nÀ1 ðkaÞ c À l 2nÀ1 sinðð2n À 1ÞhÞ where C is an arbitrary amplitude. The corresponding eigenvalue is k = k(k) with 2paðc ÀĉÞ C À l n J n ðkaÞ 2 : For the sake of illustration, suppose that l n \ c for all n 2 Z. This is plausible given Eq. (65) and the conditions on l n . Equation (69) implies that if c\ĉ then the topographic map is stable since k(k) \ 0 for all k. On the other hand, if c [ĉ such that k(k c ) = max kk (k) [ 0 then the topographic map is unstable and the fastest growing eigenmodes have the critical wavenumber k c . Recall that c = m 1 . It then follows from Eq. (65) that l 1 & c and the dominant contribution to the sum in Eq. (69) will arise from the n = 1 term, at least away from the zeros of J 1 (ka). Hence, k c is approximately given by the point at which the first-order Bessel function attains its global maximum, that is, |J 1 (ak c )| = max k |J 1 (ak)|. Numerically one finds that k c & 3/a. Note that the eigenvalues k(k), k = 0, have an infinite degeneracy that reflects the additional rotation symmetry of the system. That is, all eigenmodes P h (k) with |k| = k have the same eigenvalue. It follows that the pattern-forming instability will be dominated by some linear combination of eigenmodes lying on the critical circle |k| = k c : where k i ¼ ðk c ; u i Þ and z i is a complex amplitude. Suppose that each eigenmode can be approximated by the first three terms of Eq. (68) so that Pðk c ; hÞ % C J 0 ðk c aÞ c À l 0 þ 2J 1 ðk c aÞ c À l 1 sinðhÞ À 2J 2 ðk c aÞ c À l 2 cosð2hÞ The first term generates an expansion of the bump, the second term generates a uniform shift of the bump, and the third term generates an elongation of the bump (see Fig. 2). In general, we expect the eigenmode (71) to be dominated by the first harmonic term sinðhÞ, since l 1 & c. However, if l 2 & c as well, then there could also be a significant contribution from the term cosð2hÞ. Thus, the spontaneous symmetry breaking mechanism has the potential for generating elongated or oriented bumps. Moreover, since each eigenmode in the sum (70) then represents an elongation in the direction u i or p=2 þ u i (depending on the sign of its associated coefficient zðrÞ ¼ z i e ik i Ár þ z Ã i e ik i Ár ), it follows that there is some complicated variation in the preferred bump orientation as r varies across the cortex.
In other words, the two-dimensional network has both translation and rotation/reflection symmetries. An isotropic and homogeneous equilibrium solution of the form uðr 2 jrÞ ¼ Uðjr 2 À rjÞ; Iðr 2 jrÞ ¼ Iðjr 2 À rjÞ then explicitly breaks the symmetry group from Eð2Þ Â Eð2Þ ! Eð2Þ with E(2) having the group elements T d = T d , d, R n = R n,n and R j = R j,j . As we have shown above, the homogeneous equilibrium solution can undergo a pattern-forming instability that spontaneously breaks the remaining Euclidean symmetry. The Euclidean symmetry of the full Eqs. (35) and (36) also manifests itself in a reduced form after linearising about the homogeneous solution. That is, the linear equations (48) and (50) are equivariant with respect to the socalled shift-twist action of the Euclidean group E(2) on the space R 2 9 S 1 [10,11]: It can be seen that the rotation operation comprises a translation or shift of the angle h to h ? s, together with a rotation or twist of the position vector r by the angle s. This is illustrated in Fig. 6. One of the consequences of the underlying Euclidean symmetry is that the associated eigenfunctions form irreducible representations of the shift-twist group action [10,11]. This explains why the eigenmodes P h (k) have the basic structure given by equation (68), with the angular variable h coupled to the direction of the wavevector k.
The emergence of elongated cortical activity bumps reflects a corresponding elongation in the receptive field properties of cortical cells. Hence, the spontaneous symmetry breaking of the topographic map between LGN and visual cortex provides a possible mechanism for the development of orientation preference columns in primary visual cortex. The variation in orientation with cortical position then determines the corresponding orientation preference map. Interestingly, a recent statistical analysis of orientation preference maps in primates indicates that there are correlations between the direction of the topographic axis joining pairs of columns with similar orientation preferences and their common orientation [34]. Thus, the orientation preference map exhibits a form of rotational shift-twist symmetry as predicted from our analysis of two-dimensional topographic maps. Numerical simulations of a feature-based dynamical spin model has led to the suggestion that such a symmetry could help to stabilise the emerging orientation preference map with its associated set of pinwheels [34]. As previously shown by Wolf and Geisel [55], in the absence of such a coupling, the pinwheels typically annihilate in pairs. Hence, in order to maintain pinwheels, either development has to be stopped or one has to introduce inhomogeneities that trap the pinwheels.

Discussion
In this review, we have described John G Taylor's work on neural bubble dynamics [48], as well as various extensions that have revealed a rich repertoire of 2D network dynamics, including breathers [15,20,39], multiple bumps [32], labyrinthine structures [16,39], rotating waves [21], spiral waves [26,29,31,43], and travelling spots [15,16]. Traditionally, in vitro methods for studying the generation and propagation of electrical activity in networks of neurons have involved removing a thin vertical slice of brain tissue [13,23,40,42], which makes it difficult to observe the various forms of 2D dynamics predicted by neural field theory. However, recent studies of tangential slices have demonstrated that various phenomena predicted by neural field theory are also observed in living tissue [26,57]. Interestingly the type of neural field models discussed here (and their localised solutions) has been adopted as a formal framework for thinking about embodied cognitive dynamics by Gregor Schöner and colleagues, see for example [28]. Initially, this work focused on developing a theory of how eye movements are planned, though has since been expanded to a variety of topics including motor planning and visuospatial cognition. Moreover, the same style of conceptual modelling has proven fruitful in developing a theory of cognitive robotics [18]. There is clearly still much to be gained by the mathematical study of neural field models along the lines developed by John G Taylor and others. where J 0 is the Bessel function of the first kind, we express w in terms of its Hankel transform of order zero, wðrÞ ¼ which, when substituted into Eq. (4), gives UðrÞ ¼ where 0 ¼ 1 and n ¼ 2 for n C 1, we thus obtain (10), using the identity J 1 ðqDÞD ¼ q R D 0 J 0 ðqr 0 Þr 0 dr 0 .