We investigate two-dimensional neural fields as a model of the dynamics of macroscopic activations in a cortex-like neural system. While the one-dimensional case was treated comprehensively by Amari 30 years ago, two-dimensional neural fields are much less understood. We derive conditions for the stability for the main classes of localized solutions of the neural field equation and study their behavior beyond parameter-controlled destabilization. We show that a slight modification of the original model yields an equation whose stationary states are guaranteed to satisfy the original problem and numerically demonstrate that it admits localized noncircular solutions. Typically, however, only periodic spatial tessellations emerge on destabilization of rotationally invariant solutions.
Neural fields (Amari, 1977; Wilson & Cowan, 1973) describe the dynamics of distributions of activity on a layer of neurons. Neural fields have been suggested as models of internal representations in natural agents (Takeuchi & Amari, 1979; Gross, Stephan, & Krabbes, 1998) as well as in robots (Steinhage, 2000; Iossifidis & Steinhage, 2001; Erlhagen & Bicho, 2006; Giese, 1998). Various modalities are covered, such as spatial localization, viewing direction, attentional spotlight, the dynamics of decision making, elementary behaviors, and positions of other agents in the environment. More extensive studies in theoretical neuroscience (Suder, Wörgötter, & Wennekers, 2001; Mayer, Herrmann, & Geisel, 2002; Bressloff, 2005) concentrate on primary visual cortex (Lieke et al., 1989), superior culliculus (Schierwagen & Werner, 1996), the representation of motoric primitives (Thelen, Schöner, Scheier, & Smith, 2001; Erlhagen & Schöner, 2002), and working memory in prefrontal cortex (Schutte, Spencer, & Schöner, 2003; Camperi & Wang, 1998). Generally neural fields serve as nonparametric representations of probability densities, and their dynamics may perform operations on the densities such as Bayesian computations (Herrmann, Pawelzik, & Geisel, 1999).
Recently neural fields have attracted attention in modeling and analysis of brain imaging data, because they are able to represent the dynamic interaction of an active medium with time-varying inputs and because the spatial and temporal scales in the data and neural field models are starting to become comparable. Especially if information about connectivity is available from the data (Jirsa, Jantzen, Fuchs, & Kelso, 2002), neural fields are of high explanatory power. Moreover the analysis has reached a level where applications directly benefit from the theoretical progress, and at the same time, computational power has become available that allows us to perform online simulations of two-dimensional neural fields.
One-dimensional problems were comprehensively studied in the 1970s (Amari, 1977; Kishimoto & Amari, 1979; Takeuchi & Amari, 1979). While for spatially extended (e.g., periodic) patterns, the transition to the more relevant two-dimensional case is nontrivial but fairly well understood (Ermentrout & Cowan, 1979; Ermentrout, 1998), localized activities in dimensions larger than one have not yet been treated with the same rigor. The situation, however, does not seem to be complex: a large body of numerical studies, together with theoretical considerations (Laing, Troy, Gutkin, & Ermentrout, 2002; Laing & Troy, 2003), imply a general instability of multibump solutions if the interactions are excitatory at small and inhibitory at large distances (see also Laing & Chow, 2001, for stability analysis in one-dimensional models of spiking neurons). Further, there has been numerical evidence that single-bump solutions in two dimensions for radially symmetric interactions are essentially circular, which was exploited as an assumption in Taylor's early attack on the two-dimensional case (Taylor, 1999). Werner and Richter (2001) provided evidence for the existence of ring-shaped solutions that are possible for certain types of neural interactions. Along these lines, one may conjecture that finite mesh-like structures of higher genus exist as well.
The situation became more spirited recently when Herrmann, Schrobsdorff, and Geisel (2004), Bressloff (2005), and Doubrovinski (2005) tackled the stability problem of localized activations in two-dimensional fields. Although the generality that has been achieved in the one-dimensional neural fields cannot be expected in two dimensions, a number of interesting variants of the circular activity configurations have been analyzed.
Here we present a concise and reproducible scheme for the analysis of the stability of localized activity distributions in neural fields. We show the applicability of the scheme to the special case of circular solutions and also study ring-shaped and bar-shaped solutions, which present the complete set of known localized solutions for simple kernels (Werner & Richter, 2001). In addition to the stability proofs based on the classical scheme (Pismen, 1999; Bressloff, 2005; Doubrovinski, 2005), we are interested in behavior beyond the phase transitions toward unstable regions. The destabilization of a circular solution is known to lead to a transient elongation of the activity bubble (Bressloff, 2005; Doubrovinski, 2005), which ultimately causes the localized solution to split or form meandering bands. Either case is unstable in a strict sense: the splitting into two continues toward a plane-filling hexagonal pattern, while the banded patterns develop a global stripe pattern or a quasiperiodic arrangement.
The destabilization is thus fundamental since for typical Mexican hat interaction kernels, there is no nearby stable state that is approached after the bifurcation, while the spatially extended patterns are not approached in finite time (unless a general criterion for convergence is drawn into consideration). Yet the destabilization, at least in neurobiological applications, is the most interesting part of the theory. Divergences are usually very slow and may halt completely due to reasonable boundary conditions, such that an activity-based correlational learning scheme may organize anisotropies in the connections, which stabilizes the anisotropic activities as assumed in Bressloff (2005) and exploited in Schierwagen and Werner (1996). A theoretical account of the interaction of activity dynamics and learning was studied in Dong and Hopfield (1992) in relation to the activity effect on feature maps (cf. Mayer et al., 2002).
2. The Neural Field Equation
The neural field model describes the activations of a layer of neurons when the geometry of interactions rather than the specific connectivity among the neurons is relevant. We assume positions for neurons with continuous-valued activations u(r, t). The synaptic weights between neurons at the positions r and r′ are expressed by isotropic interaction kernel w(r, r′) = w(|r − r′|) of Mexican hat shape. Neurons are activated if their total input is greater than zero. We will study only equilibrium solutions without external input, and we neglect slow learning effects, so the synaptic weights are constant over time.
If nonrotationally symmetric solutions are excluded from the beginning from the consideration of equation 2.2, the situation simplifies dramatically, and it can be shown (Taylor, 1999) that one-bump solutions u(‖r‖) of certain radii are stationary states of the dynamics 2.1. Analogously, a ring-shaped solution (Werner & Richter, 2001) or a stripe-like solution (i.e., a degenerate ring of infinite radius) can occur. However, when an arbitrarily small perturbation of the solution is being considered, the symmetry might be broken and new phenomena can appear.
It has previously been shown that equation 2.1 admits rotationally invariant stationary solutions with disc-shaped activated region (one-bump solutions). Generically these arise in the course of a “blue sky bifurcation” (Strogatz, 1994): no solution is present in the subcritical parameter region, whereas two solution branches bifurcate as the control parameter exceeds the critical value. The two solutions are rotationally invariant one-bumps. Stability analysis of these states with respect to rotationally invariant perturbation is essentially equivalent to stability analysis of the one-bump solution of the one-dimensional model. It reveals that the unstable branch generically corresponds to the bump of smaller radius. Upon destabilization, the region of activation expands as the solution approaches the stable branch, corresponding to the circular bump of larger radius (see Figure 1).
Seeking a stationary solution of the two-dimensional model, equation 2.1, assuming rotational invariance of the field (i.e., with r = ‖r‖), leads to a problem that is essentially equivalent to that of finding stationary states in the one-dimensional model (see the appendix). Two types of solutions besides the circular one-bumps are readily constructed: solutions with annular activated regions and solutions with a stripe-shaped region of activation (see Figure 1). The possibility of the existence of the former has been pointed out previously, while the latter can be seen as degenerate annuli of infinite inner radius.
Previous work (Werner & Richter, 2001) conjectured that the bifurcation branch corresponding to the stable one-bump remained stable for all values of the control parameter. Using equation 3.6, this is readily proved false: higher and higher frequency eigenmodes progressively turn unstable as the radius of the stationary solution is increased (see Figure 3).
Strictly speaking, the assertion that the linear stability analysis as outlined above correctly determines the stability properties of stationary solutions relies on additional assumptions on operators appearing on the right-hand side of equation 2.1 (Golubitsky, Schaeffer, & Stewart, 1985). For infinitely dimensional nonsmooth fields, these generically need not hold. In order to check whether stability analysis is adequate, correctly determining the stability properties of the stationary solutions, we performed a number of numerical simulations. Figure 4 depicts a simulation of unstable one-bump whose corresponding stability spectrum reveals that the maximal (positive) eigenvalue is that of the D2-symmetric eigenmode. As time elapses, the initially circular activated region keeps deforming, forming bloblike protrusions. These subsequently bud off from the middle bump. The newly formed activated domains keep splitting, progressively tiling the plane in a hexagonal pattern. We stress that the pattern that is formed immediately after destabilization of the stationary state is D2-symmetric, as would be expected from the properties of the eigenvalue spectrum. That is, stability properties of the stationary solutions, as well as qualitative aspects of the pattern, forming upon destabilization of the steady state, are correctly predicted from the linear stability analysis.
Symmetry breaking accompanying destabilization of a stationary bump was examined for a broad range of parameters (e.g., parameters that yielded solutions with a maximal eigenvalue of the respective spectrum corresponding to D3-, D4-, and D8-symmetric perturbations; see Figure 4). In all of the cases, the course of the symmetry breaking was appropriately determined by linear stability.
Stability analysis of annular solutions proceeds along the same lines as that of one-bumps (see the appendix). The essential difference is that instead of a single equation, 3.3, a system of two equations results, meaning that to every nonnegative whole number corresponds a pair of real eigenvalues. (In general, to every boundary of an activated domain, there corresponds an equation in the corresponding eigenvalue problem.) We find that in the case of the annular solution shown in Figure 1c, the largest eigenvalue corresponds to a D3-symmetric perturbation. Simulations demonstrate that initially a rotationally invariant annulus splits into three adjacent blobs that gradually drift apart (not shown). The symmetry of the resulting state is the same as that of the largest eigenvalue. Again, the stability of the stationary solution, as well as qualitative aspects of the emerging pattern, are readily predicted from the analytically computed spectrum.
Finally, let us consider the stripe-shaped solutions. Their stability is governed by an eigenvalue problem, similar to that governing the stability of the annuli. However, the corresponding linear operator is no longer compact (this is a consequence of activated region being unbounded), and the spectrum no longer needs to remain discrete. Actually, in this case, the spectrum is continuous: to every real corresponds a pair of (real) eigenvalues. Figure 3b shows results of the linear stability analysis of the stripe solution, depicted in Figure 1d. In the corresponding simulation, the stripe is seen to split into a row of separate bumps—a “chain-of-pearls” configuration. The interbump separation is the same as the wavelength of the eigenmode, corresponding to the largest eigenvalue. Again, stability and semiquantitative properties of patterns resulting from the stationary state destabilization are readily predicted from the respective eigenvalue spectrum.
In summary, the preceding section describes all of the stationary nonhomogeneous solutions of the two-dimensional Amari equation known to date and exhaustively examines their stability properties. Quite strikingly, stability analysis of the nonhomogeneous steady state is possible since the eigenvalue problem, equation 3.5, is effectively one-dimensional, although a two-dimensional system is being considered.
4. Modified Equation
A long-standing question regarding the Amari model is the existence of nonrotationally invariant stationary solutions with bounded and connected regions of activation. Such states are likely to bifurcate from circular solutions on destabilization at a nonrotationally invariant eigenvalue. In fact, the stability loss of one branch is always accompanied by the emergence of another in its vicinity, provided that the mappings defining the dynamical system under consideration are sufficiently smooth (Crandall & Rabinowitz, 1971). However, in the simulations, exclusively periodic patterns resulted on destabilization of circular solutions.
In order to address the existence of the nonrotationally invariant solutions of equation 2.1 with a bounded and connected activated region, we modify the original Amari model, equation 2.1, so as to obtain a related (modified) equation fulfilling the following three conditions:
According to condition 3, rotationally invariant stationary solutions behave like those of the unmodified Amari model, admitting symmetry breaking at a nonrotationally invariant eigenvalue. However, the symmetry breaking cannot result in a spatially extended periodic pattern according to condition 2. Consequently, a nonrotationally invariant localized state is likely to emerge. It will be a solution of the original (unmodified) Amari model with desired properties, provided that its region of activation remains connected in the course of destabilization. Note, however, that conditions 1 through 3 do not suffice to ensure that the activated domain will not start splitting into a separate disconnected region upon symmetry breaking.
Equation 4.1 does not admit stationary solutions with an unbounded region of activation. Indeed, assuming that such a solution exists, substituting into the equation, we would conclude that the integral term of the right-hand side of equation 4.1 is infinite, which forms a contradiction. Therefore, condition 2 above is satisfied.
Finally, we show that modification 4.1 preserves stability properties of the stationary solution (which is a stationary solution of both the modified and the unmodified problems, 2.1 and 4.1, respectively, by construction). Recall that in deriving equation 3.5, we did not make use of any particular assumptions on the form of the integral kernel w(|r − r′|). Consequently, this expression for the eigenvalue spectrum is equally valid for the modified problem as well as for the unmodified one. Note, however, that when deriving a stability spectrum in the case of the modified problem, equation 4.1, we shall exchange g(θ − θ′) by (see the derivations in the appendix). Using , it now follows from equation 3.5 that all eigenvalues of the stability spectrum of (except for λ0 corresponding to perturbations of the radius of the bump) remain unaltered by the modification. Consequently, if is unstable with respect to some nonrotationally invariant perturbation in the original problem, 2.1, it is unstable with respect to such a perturbation in the modified problem, 4.1, whereby condition 3 holds.
These arguments imply that a rotationally invariant solution of equation 2.1 that is unstable at a noncircular eigenmode also solves the modified equation, 4.1, and is unstable with respect to the same eigenmode of the dynamics of equation 4.1. Furthermore, contrary to the case of equation 2.1, the dynamics of equation 4.1 can never result in a periodic pattern with unbounded activated region if q < 0.
As stated above, conditions 1 to 3 do not suffice to guarantee that the activated region will remain connected as follows the dynamics of equation 4.1. Nevertheless, by tuning the parameter q in equation 4.1, one is able to trap the dynamics in the vicinity of instability in a state with connected activated region.
Assume that for q = 0, the dynamics in the vicinity of the bifurcation tends to increase the area of the activated region. Note that q could be understood as a Lagrange multiplier, which ensures the constancy of the area for one special value, which we certainly exclude. Suppose now that q is chosen to be negative. The additional term in the integrand of equation 4.1 will tend to counterbalance the area increase. We can now choose q such that two effects counterbalance, and a nonrotationally invariant steady state with bounded and connected region of activation will result. The calculation of effective values of q requires the consideration of higher orders of the dynamics. Here, we instead resort to numerical simulations. These confirm the emergence of nonrotationally invariant steady states with bounded and connected regions of activation. For example, Figure 5 shows the destabilization of a circular one-bump that is unstable at a D2-symmetric eigenmode, developing into a nonrotationally invariant steady state with an ellipse-shaped activated region. Only one-quarter of the domain was simulated (the field in the other three quadrants is determined by that on the simulated quadrant due to the Euclidean symmetry of dynamic equations and the D2-symmetry of initial conditions) on a grid of 300 × 300 pixels in order to increase accuracy. Symmetries corresponding to higher eigenvalues are irrelevant because the eigenvalues λn are stable for n ⩾ 3 for the given parameters. The length difference of the half-axes of the activated domain of the resulting D2-symmetric stationary state was much larger than the spatial discretization step (some 30 pixels), allowing the conclusion that the stationary solution obtained is not a discretization artifact.
In conclusion, we stress that the particular choice of modification according to equation 4.1 is very restrictive. In fact, many other equations whose stationary states satisfy equation 2.1 and do not admit solutions with unbounded regions of activation are readily constructed along the same lines. For instance, ∂tu = −u + ∫R[u]w(|r − r′|) dr′ − ∫R[u]|r − r′|dr′dr + h can be shown to have these properties, arguing essentially as when proving that conditions 1 and 2 are satisfied by stationary states of equation 4.1. We believe that further investigation of such modifications will provide insight into the properties of unstable solutions of the Amari equation.
Bifurcation theory describes the time course to critical behavior in low-dimensional systems. Under certain conditions, the parametric destabilization of an activity distribution does not lead to a nearby stable state but initiates a cascade of symmetry-breaking events that eventually approaches a new distant stable configuration. A motivation for this study is to demonstrate the complex evolution of the state of the field after the loss of stability. It is these configurations that bear the greatest computational potentials.
For the Amari equation with a simple kernel, only the existence of rotationally invariant bumps has previously been proven. This work demonstrates the existence of two other nonperiodic stationary states: stripe shaped and annular solutions. Stability analysis of these would be expected to be very involved. Surprisingly, the special form of the Amari equation makes this stability problem amenable to analytical treatment. The total synaptic input to a given neuron is dependent only on the shape of the boundary of the activated domain, making the eigenvalue problem effectively one-dimensional, thereby allowing fairly straightforward calculation of the spectrum for each of these cases.
Strictly speaking, spectral properties of the linearized operator do not guarantee the stability of the stationary state unless additional assumptions are satisfied. However, our extensive numerical investigation shows that stability is indeed correctly predicted by eigenvalue analysis. Furthermore, the spectrum allows predicting certain semiquantitative features of solutions approached after the onset of instability.
Our numerical experiments showed that exclusively spatially extended solutions (those with unbounded activated domain) appeared in the unstable parameter region. Yet a slight modification of the interaction kernel, introducing long-range interactions, circumvents this, yielding nonrotationally invariant stationary states with bounded and connected regions of activation that at the same time are stationary states of the original unmodified equation. This settles a long-standing question regarding the existence of solutions of this type. Note, however, that stability analysis of these is nontrivial, requiring the knowledge of analytical parameterization of the boundary of the activated domain. The existence of such states is also relevant for biological systems. Assuming a certain degree of shift-twist symmetry, the existence of elongated blobs suggests mechanisms for the emergence of orientation selectivity in neurons of the primary visual cortex (Tanaka, Ribot, Imamura, & Tani, 2006). Although the degree of asymmetry of nonrotationally invariant solutions was moderate, these effects could in principle be enhanced by (Hebbian) learning, which we disregarded in our treatment.
Appendix: Linear Stability of Rotationally Invariant Solutions
A.1. Circular One-bump Solution.
A.2. Annular Solutions.
A.3. Stripe-Shaped Solutions.
We thank S. Kaijser, T. Geisel, and S. Amari for their kind support of this work and for helpful remarks. We also thank Hecke Schrobsdorff for stimulating discussions. This work was partially supported by the BMBF in the framework of the Bernstein Centers for Computational Neuroscience, grant number 01GQ0432.