Stability analysis of a neural field self-organizing map

We provide theoretical conditions guaranteeing that a self-organizing map efficiently develops representations of the input space. The study relies on a neural field model of spatiotemporal activity in area 3b of the primary somatosensory cortex. We rely on Lyapunov’s theory for neural fields to derive theoretical conditions for stability. We verify the theoretical conditions by numerical experiments. The analysis highlights the key role played by the balance between excitation and inhibition of lateral synaptic coupling and the strength of synaptic gains in the formation and maintenance of self-organizing maps.


Introduction
Self-organizing maps (SOMs) are neural networks mapping a high-dimensional space to a low-dimensional one through unsupervised learning. They were first introduced by Grossberg [14] and later by Kohonen [19]. SOMs are widely used in computer science and data analysis for quantization and visualization of high-dimensional data [25,38]. They also constitute a suitable tool in computational neuroscience to study the formation and maintenance of topographic maps in primary sensory cortices such as the visual cortex [24,31] and the somatosensory cortex [13,34]. Many variations and applications of Kohonen's SOM algorithm can be found in [16] and [27].
A type of self-organizing map based on neural fields theory has been introduced in [8], where neural fields are used to drive the self-organizing process. Neural fields are integrodifferential equations that describe the spatiotemporal dynamics of a cortical sheet [3][4][5]. The SOM proposed in [8] describes the topographic organization of area 3b of the primary somatosensory cortex of monkeys [21,28]. The model relies on an earlier work [29] known as the dynamic SOM (DSOM) algorithm. DSOM provides an online SOM learning algorithm, where the Kohonen's SOM time-dependent learning rate and neighborhood function have been replaced by time-invariant ones. The DSOM neighborhood function and learning rate solely depend on the distance of the winner unit (i.e., the most active neuron) from the input. The model proposed in [8,9] combines the DSOM timeinvariant learning rate and neighborhood function with Oja's learning rule [26]. As thor-oughly described in [8,9], the model is compatible with anatomical evidence of how area 3b in monkeys develops, maintains, and reorganizes topographic representations of a skin patch of the index finger.
In this work, we provide theoretical insights on the stability and convergence of the neural field SOM algorithm proposed in [8,9] by studying a more general class of systems than that originally proposed in [8]. We use Lyapunov's stability theory adapted to neural field dynamics [11]. Since typical activation functions employed in the model (such as absolute values or rectification functions) are not necessarily differentiable, we do not rely on linearization techniques but rather directly assess the stability of the original nonlinear dynamics. Yet, the obtained results are local, meaning that they are valid only for initial conditions in the vicinity of the considered equilibrium. Nonetheless, we show that they agree with numerical simulations. The stability conditions derived in this work can be used toward the direction of tuning neural field models such that they achieve the best possible results in developing self-organizing maps and thus more generalized representations. Moreover, the conditions we propose indicate that the balance between lateral excitation and inhibition keeps the system stable, thus ruling out possible configurations in which learning does not take place properly. These findings are in line with both experimental observations [18,30] and computational modeling [35][36][37].
The paper is organized as follows. In Sect. 2, we recall the SOM model under concern and its basic mechanisms. In Sect. 3, we present our main theoretical results, which we confront to numerical simulations in Sect. 4. A discussion on the obtained results is provided in Sect. 5. Mathematical proofs are given in Sect. 6.

Self-organizing neural fields 2.1 Neural population dynamics
We consider the following neural fields equation: τ ∂u ∂t (r, t) = -u(r, t) + w l rr rect u r , t dr + I, where is a connected compact subset of R q (q = 1, 2, 3). For q = 2, the integral of a function g : = 1 × 2 → R is to be understood as g(r) dr = 1 2 g(r 1 , r 2 ) dr 2 dr 1 with r = (r 1 , r 2 ), and similarly for q = 3; u(r, t) represents the mean membrane potential at position r ∈ and time t ≥ 0, τ is a positive decay time constant, I denotes an external input, and w l is a function that represents the strength of lateral synaptic coupling. It is given by where the excitation and inhibition synaptic weights are typically given by and with K e , K i , σ e , σ i > 0. In [8,9] the input is provided through a two-dimensional skin model. The skin model is composed of a two-dimensional grid and receptors. The receptors are points distributed on the surface of the grid (uniformly). When a stimulus is applied on the grid, the receptors sample the input signal and convey the information to the cortical model. The skin stimulus is a noisy Gaussian-like function, and the input to the neural fields model is provided by the following function: where | · | 1 denotes the 1-norm: |x| 1 = m i=1 |x i |, and s : R 2 → [0, 1] m is a function that maps the raw input from the two-dimensional skin space to [0, 1] m . For instance, for a tactile stimulus at position p ∈ R 2 on the skin, s(p) ∈ R m could be defined as the normalized distance from p to the location of each receptor, thus potentially of much higher dimension than 2. For a more detailed description of receptor model, see [9]. The function w f : × R ≥0 → R m represents feed-forward synaptic weights with value updated according to where γ is a positive constant that represents the learning rate, and rect(x) = max{x, 0}. It is worth observing that since s(p) ∈ [0, 1] m , w f (r, t) ∈ [0, 1] m for all r ∈ and t ≥ 0 given any initial conditions satisfying w f (r, 0) ∈ [0, 1] m for all r ∈ (this can be seen by observing that the entries of ∂w f ∂t (r, t) are negative as soon as the corresponding entries of w f (r, t) become greater than 1; similarly, they are positive when the corresponding entries of w f (r, t) get below 0: see (5)). Hence interpreted as a high input when the feedforward weights are close to s(p) and as a lower input when these are more distant. The overall model (1), (4), (5) reflects the dynamics of a cortical neural population in combination with a learning rule of the feed-forward connections w f , which convey information from receptors to the cortical sheet. As described in [8,9], this model can express a variety of different behaviors, depending on the lateral connectivity kernels w e and w i .
The main advantage of the learning rule given by Eq. (5) is that it is a biologically plausible modification of the DSOM learning rule [29]. In DSOM the learning rate and neighborhood function are time-invariant and can adapt to the input according to one single parameter, called elasticity. This particular modification leads to the following behavior: if the winner neuron (i.e., the neuron that has the shortest distance from the input stimulus to its corresponding codebook-weight) is close to the stimulus, then the neighborhood function shrinks around it. This results in making the weights of neurons within the dynamic neighborhood stronger and the weights of the other units weaker. However, when the winning unit is very far from the current input, the neighborhood function exhibits a broad activity pattern, promoting learning of every unit in the network. Therefore in [8] the neighborhood function has been replaced by the term w e (|rr |) rect(u(r , t)) dr , providing a more realistic and biological plausible learning algorithm for self-organizing maps in the context of neuroscience.

Self-organizing maps
We start by briefly describing how the SOM model introduced in [8] and [9] works. The algorithm starts by initializing the feed-forward weights randomly (usually uniformly), and the neural field activity u(r, 0) is set to zero. The second step is sampling the input space by randomly drawn samples of dimension m from an input distribution. At every epoch, one sample is given to the neural field (1) and (5) through Eq. (4). This first step is depicted in Fig. 1(A), where a two-dimensional point p = (p 1 , p 2 ) is sampled from a uniform distribution, p 1 , p 2 ∼ U(0, 1). The samples are mapped to the neural space through the function s and then are passed to Eq. (4). At this point, we should point out that there are two ways of presenting stimuli while training a self-organizing map. The first is predetermining an amount of input samples and present one at each epoch (online learning) and the second is collecting all the input samples into a batch and giving all of them at once to the network (batch learning). In this work, we use the former (online learning) since it is biologically plausible.
Then the algorithm proceeds with computing the numerical solution of Eqs. (1) and (5). To that aim, Eqs. (1) and (5) are discretized and solved numerically using Euler's forward method. The numerical solution of Eq. (1) is typically a bell-shaped curve (bump) centered on the neuron that is the closest unit to the input sample and therefore is called the winner neuron or best matching unit (BMU). In Fig. 1(B), this is depicted as a black disc on a discrete lattice. The lattice represents a discretization of the field where each tile corresponds to a neuron. Neurons that lie within the vicinity (within the black disc in Fig. 1(B)) defined by the solution of Eq. (1) update their weights based on Eq. (5). The rest of the neurons feed-forward weights remain in their previous state. Once the temporal integration of Eqs. (1) and (5) is complete, the activity of the field is reset to its baseline activity. Then another input sample is drawn, and the whole process repeats itself. Once the number of epochs is exhausted, the learning stops, and the mapping process is completed.
To make the aforementioned algorithm directly comparable to Kohonen SOM [19], we provide some insights. First, in Kohonen's SOM, we compute the distance between the input and the codebooks. Here we do the same using Eq. (4). The neighborhood function that Kohonen's SOM uses to update the feed-forward weights is replaced here by the numerical solution of the neural field (Eq. (1)) and more precisely by the term w e (|rr |) rect(u(r , t)) dr . Both the learning rate and the width of the neighborhood function are time-independent in our case, as opposed to Kohonen's SOM, where they are both time-dependent. Our learning rule is different since we use a modified Oja rule [26], which is based on Hebbian learning [15], and it is therefore biologically plausible [1]. The dimensionality reduction in both models, the Kohonen and ours, takes place at the level of the learning rule. This means that Eq. (5) is responsible for learning the representations and mapping the input distribution (of dimensions m) on a manifold of lower dimension q ∈ {1, 2, 3}.

Explicit conditions for stability
The most important question when one trains a self-organizing map is: Will the learning process converge and properly map the input space to the neural one? In most of the cases, it is not possible to predict this. However, in the specific case of the self-organizing algorithm provided by [8], here we show that it is possible to obtain an analytical condition that guarantees the stability of the equilibrium point of system (1)- (5). Stability during learning is a prerequisite to generate a meaningful mapping and thus a proper topographic map. Moreover, a byproduct of deriving such a stability condition is providing some insights on how to properly tune model parameters.
To this end, we now proceed to the mathematical analysis of the model. For generality, the adopted mathematical framework is slightly wider than merely Eqs. (1), (4), (5) and encompasses more general classes of activation functions and synaptic kernels. We start by introducing the considered class of systems and then provide sufficient conditions for its stability and convergence.

Model under study
The self-organizing neural field (1), (4), (5) is a particular case of the more general dynamics where τ , γ > 0, w l , w e ∈ L 2 ( 2 , R), the set of all square-integrable functions from 2 to R, and f e , f l and f s are Lipschitz continuous functions.

Stability analysis of Eq. (6a) and (6b)
We recall that an equilibrium x * of a systemẋ(t) = f (x(t)), where x(t) : → R n for each fixed t ≥ 0, is called globally exponentially stable if there exist k, ε > 0 such that, for all admissible initial conditions, where · denotes the spatial L 2 -norm. This property ensures that all solutions go to the equilibrium configuration x * in the L 2 sense (global convergence) and that the transient overshoot is proportional to the L 2 -norm of the distance between the initial configuration and the equilibrium (stability). The equilibrium pattern x * is said to be locally exponentially stable if (8) holds only for solutions starting sufficiently near from it (in the L 2 sense). We refer the reader to [11] for a deeper discussion on the stability analysis of neural fields. Our main result proposes a sufficient condition for the local exponential stability of Eq. (6a) and (6b). Its proof is given in Sect. 6.1.

Theorem 1 Let be a compact connected set of
the equilibrium pattern (u * , w * f ) is locally exponentially stable for Eq. (6a) and (6b).
Condition (9) imposes that the synaptic weights of the lateral coupling w l are sufficiently small: stronger lateral synaptic weights can be tolerated if the maximum slope l of the activation function f l is low enough, meaning that the system given by Eq. (6a) is less selfexcitable. Recall that if f l is a differentiable function, then l can be picked as the maximum value of its derivative. Nonetheless, Theorem 1 does not impose such a differentiability requirement, thus allowing us to consider nonsmooth functions such as absolute values, saturations, or rectification functions. Note that it was shown in [33] that condition (9) ensures that the system owns a single equilibrium pattern. It is also worth stressing that the slopes of the functions f s and f e do not intervene in the stability conditions. Condition (10) requires a sufficient excitation in the vicinity of the equilibrium u * . Roughly speaking, it imposes that the considered equilibrium pattern u * does not lie in a region where f e is zero.

Stability analysis of the SOM neural fields
Theorem 1 provides a stability condition for the model described by Eq. (6a) and (6b). We next apply it to the model given in [8] to derive more explicit and testable stability conditions. More precisely, the self-organizing neural fields (1), (4), (5) can be put in the form of Eq. (6a) and (6b) by letting f e (x) = f l (x) = rect(x), f s (x) = 1 -|x| 1 m , and w e r, r = K e e -|r-r | 2 /2σ 2 e , (11a) In view of (7a) and (7b), the equilibrium patterns of Eqs. (1), (4), (5) are given by The Lipschitz constant of f l is l = 1. Based on this, we can also derive the following corollary, whose proof is provided in Sect. 6.2.

Corollary 1
Assume that is a compact connected set of R q , and let w e , w i , and w l be as in (11a)-(11c). Then, under the condition that the equilibrium (u * , w * f ), as defined in Eq. (12a) and (12b), is locally exponentially stable for Eqs. (1)-(5).
A particular case for which local exponential stability holds is when the excitation and inhibition weight functions are sufficiently balanced. Indeed, it appears clearly that Eq. (13) is fulfilled if K e K i and σ e σ i . See the discussion in Sect. 5 for further physiological insights on this condition.
The integral involved in (13) can be solved explicitly. For instance, in the twodimensional case (q = 2) the condition boils down to the following: and let w e , w i , and w l be as in (11a)-(11c). Define where Erf : R → (-1, 1) denotes the Gauss error function. Then, under the condition the equilibrium (u * , w * f ), as defined in Eq. (12a) and (12b), is locally exponentially stable for Eq. (1)- (5).
Plenty of approximations are available for the Erf function in the literature. For instance, the following expression approximates it with a 5.10 -4 error: with a 1 = 0.278393, a 2 = 0.230389, a 3 = 0.000972, and a 4 = 0.078108; see, for instance, [2]. The Erf function is also commonly implemented in mathematical software, thus making Eq. (15) easily testable in practice.

Numerical assessment on a two-dimensional map
To numerically assess whether the above stability condition correctly predicts the performance of the learning process, we focus on a simple example of a two-dimensional map (q = 2) and a two-dimensional input space (n = 2). Furthermore, we choose s(p) to be the identity function since we do not consider any receptors: the position of the tactile stimuli is assumed to be directly available. This choice is motivated by the fact that the presence or absence of a receptors grid does not affect the theoretical results of the current work. We refer to [8,9] for a more complex application of the neural field self-organizing algorithm.
We sample two-dimensional inputs from a uniform distribution. Therefore we have s i (p) = (p 1 , p 2 ), where i indicates the ith sample, and p 1 , p 2 ∼ U(0, 1). In all our simulations, we use 7000 sample points and train the self-organizing map over each of them (7000 epochs). It is worth stressing the difference between the training time (epochs) and the simulation time. The former refers to the iterations over all the input samples (stimuli): one such input is presented to the model at each epoch. The latter is attributed to the numerical temporal integration of Eqs. (1)- (5). Thus each epoch corresponds to a predefined number of simulation steps. At the end of each epoch the activity of the neural field is reset to baseline activity before proceeding to the next epoch.

Parameters and simulation details
The neural fields equations are discretized using k = 40 × 40 units. Accordingly, the twodimensional model (1)-(5) is simulated over a spatial uniform discretization d of the spatial domain = [0, 1] × [0, 1], namely d = 40 i,j=1 ( i 40 , j 40 ). The considered input space, over which the stimuli are uniformly distributed, is also [0, 1] × [0, 1] (two-dimensional input vectors). The temporal integration is performed using the forward Euler method, whereas the spatial convolution in Eqs. (1)-(5) is computed via the fast Fourier transform (FFT). The learning process runs for 7000 epochs. The components of the feed-forward weights are initialized from a uniform distribution U(0, 0.01), and the neural field activity is set to zero. At each epoch, we feed a stimulus to Eqs. (1)-(5), and the system evolves according to its dynamics, whereas the feed-forward weights are being updated. Then we reset the neural fields activity to zero. We run each experiment ten times using a different pseudorandom number generator (PRNG) seed each time (the PRNG seeds are given in Appendix 7: the same initial conditions and the set of PRNG seeds were used in each experimental condition).
The source code is written in Python (Numpy-, Numba-, Sklearn, and Matplotlibdependent) and are freely distributed under the GPL 3-Clause License (https://github.com/gdetor/som_stability). All the parameters used in numerical simulations are summarized in Table 1. All simulations ran on an Intel NUC machine equipped Table 1 Simulation parameters. K e and K i are the amplitudes of excitatory and inhibitory lateral connections, respectively; σ e and σ i are the variances of excitatory and inhibitory lateral connections, respectively; τ is the decay time constant, dt is the integration time step in ms, t is the simulation time in seconds, and γ is the learning rate. In each epoch, one stimulus is presented to the model

SOM quality measures
We measure the quality of the self-organizing maps using two performance indicators, the distortion D [6] and the δxδy representation [7]. We recall here that d is the spatial uniform discretization of = [0, 1] × [0, 1] and k = 40 × 40 is the number of nodes (neurons). Furthermore, for each j ∈ {1, . . . , k}, w j f (t * ) denotes the steady-state value of the feed-forward weights at the jth node of the spatial discretization, and t * corresponds to the time at the end of an epoch.
The distortion assesses the quality of a self-organizing map. It measures the loss of information over the learning process. In other words, it indicates how good a reconstruction of an input will be after the mapping of all inputs to a lower-dimensional neural map. In a sense, distortion measures how well a SOM algorithm "compresses" the input data with respect to the neighborhood structure. Mathematically, the distortion is computed according to its discrete approximation: where n is the number of samples we use during the training of the self-organizing map. Distortion is essentially an indicator of the map convergence, but it is not a reliable tool for assessing its quality. To gauge the quality of the map, we use the δxδy representation [7]. It shows when a map preserves the topology of the input space and hence how well a topographic map is formed. To estimate the δxδy, we compute all the pairwise distances between the feed-forward weights, δx = δx(i, j) = |w i f (t * )w j f (t * )|, and all the distances between the nodes of the uniform discretization of the input space [0, 1] 2 , δy(i, j) = |y iy j | for i, j = 1, . . . , k, where y i are the discrete nodes of d . We plot the δxδy (i.e., δx is the ordinate, and δy the abscissa) along with a straight line, named L δx-δy , that crosses the origin and the mean of δx points. If the point cloud representation of δxδy closely follows the line L δx-δy , then the map is considered well-formed and preserves the topology of the input space.
To quantify the δxδy representation through a scalar performance index, we perform a linear regression on the point cloud of δxδy without fitting the intercept (magenta line in figures), and we get a new line named L . Then we define the measure where a i ∈ L δ x -δy and b i ∈ L . Naturally, P should approach zero as the two lines are getting closer, indicating that the self-organizing map respects the topology of the input space, and thus it is well-formed.

Stable case
We start by simulating the model described by Eqs. (1)-(5) with the parameters given in the first line of Table 1. With these parameters, condition (15) is fulfilled (0.47 < 1), and Corollary 2 predicts that the equilibrium is exponentially stable over each epoch. Accordingly, the model succeeds in building up a self-organizing map as shown in panel (A) of Fig. 2. The white discs indicate the feed-forward weights after learning, and the black dots indicate the input data points (two-dimensional rectangular uniform distribution). Panels (B) and (C) show the δxδy representation and the distortion, respectively. We observe that the δxδy representation indicates a correlation between the feed-forward weights and the rectangular grid points (aligned with the mean of δx-red line). This means that the self-organizing map is well-formed and conserves the topology of the input. Moreover, the distortion declines and converges toward 0.0025, pointing out first that the loss of information during learning is low and that the structure in the self-organizing map is preserved. However, the boundary effects (the density of points is higher at the boundary of the map in panel (A)) affect both the distortion (it does not converge to zero; see panel (C) in Fig. 2) and the δxδy representation (it is not perfectly aligned with the red line; see panel (B) in Fig. 2). In spite of these boundary effects, the obtained δxδy performance indicator is good (P = 0.01).
The evolution of the norm-2 of feed-forward weights of three randomly chosen units The oscillations around the equilibrium are due to a repeated alteration of the input stimulus, which causes a shift to the feed-forward weights values of each winner neuron (see [8] for more detail).

Unstable case
The second line of Table 1 provides parameters for which Condition (15) is violated (5.25 > 1). According to our theoretical predictions, the model might not be stable and thus may not be able to develop any self-organizing map at all. To make sure that this is the case (and not merely a transient effect), we have let the training take more epochs (20,000). Nevertheless, here we present only the 7000 first epochs for consistency with the rest of our experiments. This situation is illustrated in Fig. 3, where the self-organizing process has failed to generate a well-formed map (panel (A)). In this case, it is apparent that self-organization process has failed to generate a topographic map.
The δx -δy representation in panel (B) of Fig. 3 looks like a diffused cloud, indicating that there is no correlation between the grid points and the feed-forward weights, meaning that there is no preservation of the topology of the input space. Accordingly, the performance index reaches the value P = 0.41, thus higher than the stable case. Moreover, the distortion in panel (C) of Fig. 3 oscillates without converging to an equilibrium, pointing out that the loss of information remains high and therefore the mapping is not successful. Finally, the norm-2 of feed-forward weights of three units (r * = (0.25, 0.25), (0.1, 0.225), (0.35, 0.075)) are shown in panel (D): it is apparent that they do not converge to an equilibrium. Instead, they oscillate violently and never stabilize around an equilibrium configuration.

Numerical assessment of Corollary 2
Finally, we numerically tested condition (15) of Corollary 2 for different values of the parameters K e and K i (all other parameters remained the same as in Table 1). For each pair  (15) is fulfilled, and therefore the weights converge to an equilibrium giving rise to a well-formed topographic map (K e , K i ), we computed the left-hand side of Eq. (15), the distortion D (averaged over the last 10 epochs), and the δxδy performance index P; see Fig. 4. We observe that for high values of K e and K i , the stability condition of Corollary 2 is violated (the black solid line overpasses the black dashed line). The distortion (orange curve) closely follows the lefthand side of condition (15) (up to a scaling factor), suggesting that distortion can serve as a measure of stability of system (1)- (5). Furthermore, the distortion and the δxδy performance index P indicate that the learning process degrades for high values of (K e , K i ), in line with the fact that condition (15) is violated. Figure 5 confirms this good alignment between the theoretical stability condition and the performance of the self-organizing map: for the first five cases, it properly maps the input space to the neural one, whereas the topology of the input space is not preserved in the last two cases, and a malformed topographic map is obtained.

Conclusion
In this work, we have presented theoretical conditions for the stability of a neural field system coupled with an Oja-like learning rule [26]. Numerical assessments on a twodimensional self-organizing map indicate that the theoretical condition is closely aligned with the capacity of the network to form a coherent topographic map.
Previous works have shown through simulations that the dynamical system described by Eqs. (1)-(5) can develop topographic maps through an unsupervised self-organization process [8,9]. The model relies on the activity of a neural field to drive a learning process. This type of models are capable of developing topographic maps and reorganize them in face of several kinds of disturbances. Here we proceed to a rigorous theoretical analysis of such kind of models by employing neural field Lyapunov theory.
The obtained stability conditions are reminiscent of those obtained for general neural fields dynamics, in which the spatial L 2 -norm of synaptic weights plays an essential role [11,22,33]. In our setting, these conditions translate in a good balance between excita- Figure 4 Numerical investigation of Corollary 2. Eight different pairs of the parameters (K e , K i ) were used to investigate the conservativeness of the stability condition given by Corollary 2. We ran eight different simulations for 7000 epochs, keeping always the rest of the parameters as in Table 1 and the same PRNG seed as before (7659). The black curve indicates the numerical value of the left-hand side of (15): the stability is guaranteed if it is below the black dashed line. The green curve indicates the δxδy performance index P.
The orange curve represents the distortion D averaged of the 10 last epochs. It is apparent that as the values of (K e , K i ) increase the Corollary 2 becomes violated and the self-organizing map fails to map the input space to the neural one (see Fig. 5 for more detail) These stability conditions provide a means to identify the parameters set within which the unsupervised learning works efficiently and thus provides an indication on how to tune them in practice. In particular, they can be used to further investigate how the dynamics of an underlying system affects the learning process during an unsupervised training process and what is the effect of the parameters on the final topographic map: as Fig. 4 indicates, the parameters of the model directly affect the quality of the topographic map. However, a limitation of the present work is that it does not offer a way of choosing the parameters in an optimal way. Furthermore, although the conditions provided by Theorem 1 guarantee the stability of the neural field, they do not predict the quality of the obtained map: Stability ensures that the learning will converge to an equilibrium, but the quality of the obtained map strongly depends on the structure of this equilibrium and hence on the chosen initial values of the feed-forward weights. This is a well-known problem with selforganizing maps [19], which is generally solved using a decreasing neighborhood, starting from a very wide one. In our case the neighborhood function is directly correlated with the profile of the field activity and is fixed (stereotyped). We thus cannot always ensure the proper unfolding of the map. It is to be noted that when the neighborhood of a Kohonen is kept fixed, it suffers from similar problems. Nevertheless, the numerical assessment of the proposed theoretical stability conditions suggests that the stability condition accurately predicts the emergence of topographic maps through unsupervised learning: see

Figs. 4 and 5.
Other works have studied stability conditions for Kohonen maps and vector quantization algorithms using methods from linear systems stability theory [32] or through energy functions [10]. However, these works focus on the learning rule for the Kohonen selforganizing maps [20], and the dynamics are not explicitly given by dynamical systems.
Our work goes beyond by taking into account not only the learning dynamics, but also the neural dynamics that drives the self-organizing process.
Last but not least, it has been shown that neural adaptation is crucial in the development of the neocortex [23] and neurons tend to adapt their input/output relation according to the statistics of the input stimuli. Our theoretical results provide conditions under which this input/output adaptation successfully takes place at least at a computational level.

Proof of Theorem 1
To place the equilibrium at the origin, we employ the following change of variables: u(r, t) = u(r, t)u * (r), where u * and w * f denote the equilibrium patterns of Eq. (6a) and (6b), as defined in Eq. (7a) and (7b). Then system (6a) and (6b) can be written as where for all x ∈ R and all r ∈ , x) = f e x + u * (r) .
With this notation, we havef l (r, 0) =f s (0) = 0 for all r ∈ , meaning that (17a) and (17b) owns an equilibrium at zero. Thus the stability properties of the origin of (17a) and (17b) determine those of the equilibria of (6a) and (6b). First, observe that since w e is a bounded function and is compact, there existsω e > 0 such that w e r, r 2 dr ≤w 2 e , ∀r ∈ .
To assess the stability of (17a) and (17b), we may be tempted to rely on linearization techniques. Nevertheless, the linearized system (17a) and (17b) around the origin would necessarily involve the derivative of f s at zero, which may be undefined if f s is not differentiable at zero (which is the case for the system of interest (1)-(5), where f s involves an absolute value). Consequently, the proof we propose here relies on Lyapunov methods [17], which were extended to neural fields in [11].
Now assumption (10) ensures that inf r∈ w e (r, r )f e (r , 0) dr > 0. It follows that there exists c > 0 such that w e r, r f e r , 0 dr ≥ 2c, ∀r ∈ .
Consequently, using (24) and the Cauchy-Schwarz inequality, we get that for any v ∈ L 2 ( , R), w e r, r f e r , v r dr = w e r, r f e r , 0 dr + w e r, r f e r , v r -f e r , 0 dr ≥ 2cw e r, r f e r , v r -f e r , 0 dr ≥ 2ce w e r, r v r dr where the last bound comes from (18). Let B ε denote the ball (in L 2 -norm) of radius ε > 0, that is, B ε := {v ∈ L 2 ( , R) : v < ε}. Letting ε := c/ ewe , we conclude from the above expression that w e r, r f e r , v r dr ≥ c, ∀r ∈ , ∀v ∈ B ε .
Thus if initial conditions are picked within the L 2 -ball of radius √ α ε √ α , then ũ(·, t) + w f (·, t) < ε at all times t ≥ 0. This means that for these initial conditions, solutions never leave the ball B ε , and hence T = +∞. Thus Eq. (30) ensures the exponential stability on this set of initial conditions.
Moreover, we claim that the solution u * of the implicit Eq. (12a) and (12b) is necessarily positive on some subset of of nonzero measure. To see this, assume on the contrary that u * (r) ≤ 0 for almost all r ∈ . Then rect(u * (r)) = 0 for almost all r ∈ , which implies that w l r, r rect u * r dr = 0, ∀r ∈ .
In view of Eq. (12a) and (12b), this implies that u * (r) = 1 for all r ∈ , thus leading to a contradiction. Consequently, as claimed, u * is necessarily positive on some subset + of of nonzero measure. Recalling that here is assumed to be a compact set, it follows that inf r∈ e -|r-r | 2 /2σ 2 e rect u * (r) dr ≥ inf r∈ + e -|r-r | 2 /2σ 2 e u * (r) dr > 0, which makes Eq. (10) satisfied. The conclusion then follows from Theorem 1.

Proof of Corollary 2
The following one-dimensional relation holds: To compute its two-dimensional counterpart, let r = (r 1 , r 2 ) and r = (r 1 , r 2 ). Then, for =