Mathematical Modeling of Working Memory in the Presence of Random Disturbance using Neural Field Equations

In this paper, we describe a neural field model which explains how a population of cortical neurons may encode in its firing pattern simultaneously the nature and time of sequential stimulus events. Moreover, we investigate how noise-induced perturbations may affect the coding process. This is obtained by means of a two-dimensional neural field equation, where one dimension represents the nature of the event (for example, the color of a light signal) and the other represents the moment when the signal has occurred. The additive noise is represented by a Q-Wiener process. Some numerical experiments reported are carried out using a computational algorithm for two-dimensional stochastic neural field equations.


Introduction
In recent years, significant progress has been made in understanding the brain electrodynamics using mathematical techniques. Neural field models represent the large-scale dynamics of spatially structured networks of neurons in terms of nonlinear integro-differential equations. Such models are becoming increasingly important for the interpretation of experimental data, including those obtained from EEG, fMRI and optical imaging [2].
In Robotics, NFE are a powerful tool for modeling and analyzing the dynamic behavior of neural populations in cognitive tasks [4]. They explain the existence of stimulus-specific, persistent population activity which bridges gaps between sensation and action in order to mediate higher-order cognitive functions like working memory, planning and decision making. The neural activity is maintained in the absence of external stimuli due to strong excitatory and inhibitory feedback loops within the population. Since NFE are defined over continuous perceptual or motor variables such as orientation, color, or movement direction, the neural interactions can be characterized by their distance in feature space. Typically, neurons tuned to similar values of the continuous variable excite each other and those tuned to dissimilar values inhibit each other. The steady states that the field dynamics develops in response to a transiently presented cue are often called bump attractors since they represent bell-shaped activity patterns in metric space. Due to the assumed translation-invariant connection structure, the network can hold a continuous family of bumps, with each of the attractor states encoding a specific stimulus value. However, since bumps are neutrally stable, their positions in the continuous attractor manifold can be shifted by random noise, potentially leading to a deterioration of the encoded information over time (see, for example, [10]).
In numerical simulations, we test the robustness of a dynamic neural field model of serial order in the presence of additive noise in an experiment in which a sequence of color cues is presented. The twodimensional field model implements the idea that the serial order of stimulus events can be represented in terms of their location along the feature dimension and the temporal dimension. The position of each individual bump of a multi-bump pattern encodes the memory of a specific color cue presented at a specific time after the sequence onset. Importantly, the joint cue-timing representation allows us to address the problem of cue repetitions which has been identified as a major challenge for models of serial order [7]. A deterministic model describing this process is presented in [16], where numerical simulations are also reported. Here, we investigate how noise-induced perturbations may affect the coding process.
Some numerical experiments reported are carried out using a computational algorithm for two-dimensional stochastic neural field equations. This numerical algorithm is MATLAB-based and presented in [14] in detail. The numerical results are discussed and their physical interpretation is explained.

Deterministic Neural Field Equation
The neural field equation was first introduced in the seminal work of Wilson and Cowan [23], as a model of spatiotemporal dependence of neural activity. Here, we apply the NFE proposed by Amari [1], as a biophysically mechanistic model of pattern formation in neural tissue: 1 When describing the working memory problems the coordinates are given by the vector x = (X, Y).
(1) where • U(x, t) (the unknown function) denotes the membrane potential in point X at time t; • I(x, t) represents the external sources of excitation; • α is the potential decay rate; • S is the dependence between the firing rate of the neurons and their membrane potentials (sigmoidal or Heaviside function); • F(x−y2) defines the coupling strength between any two neurons at the positions x and y which is assumed to depend on their distance, only. To support the existence of stable multi-bump solutions, the connectivity function is chosen to be of oscillatory type with several zero crossings. They define the distances where excitation or inhibition dominates [6,15].

Working Memory Application
Working memory (WM) defined as the ability to actively retain stimulus information over short periods of time is crucial for cognitive control of behavior. Persistent, stimulus-selective activity observed in many brain areas is commonly believed to represent a neural substrate of WM. In neural field models, WM is implemented in a bistable regime of the field dynamics in which a brief localized input may switch from a homogeneous resting state below threshold to a self-stabilized bump. Fig.1 shows an example of a 1D field representing the continuous dimension color. When two or more transient inputs are presented simultaneously or sequentially, a multi-item bump pattern stabilizes. Fig.2 illustrates the formation of a twobump pattern in response to two simultaneously presented inputs. Notice that the population activity decays slightly after cessation of the inputs but the recurrent interactions are strong enough to maintain the pattern above the threshold.
The 2D-field model of serial order of sequential events adds a temporal dimension (X) to the color dimension (Y) 1 . This choice is motivated by the experimental observation of neurons that are tuned to the temporal intervals and/or the ordinal structure of sequential tasks [20]. We use a traveling wave (assumed to be triggered by a sequence onset signal) as a second input to the field to integrate the information about the event timing relative to the sequence onset. Traveling waves which are observed in many cortical areas provide a subthreshold depolarization to individual neurons and increase their firing probability to the presentation of external stimuli [21]. Consisting with this view, we assume that only at field positions receiving the traveling wave and the color input simultaneously, the combined input is strong enough to trigger the evolution of a bump. A projection of the bump on the X and Y axes reveals bell-shaped activity profiles which are consistent with the notion of broadly tuned neurons in the feature (color) dimension and the time dimension [3,20]. Figs. 3 and 4 illustrate, respectively, the localized input and the wave input to the 2D-field. If a signal of color Y occurs at a certain time, we have a ridge-like input which is localized in the feature space Y but extends in the time dimension X. The wave has the form of a ridge which extends in the Y direction and propagates in the direction of X with elapsed time t since the sequence onset occurred at t = 0. One bump evolves in response to the combined input (see Fig.5), which persists after all inputs have been switched off (see Fig. 6).     May a stationary solution be transformed into another one just as a result of noise? What is the probability of this event?
In [12], a stochastic version of the neural field equation with additive noise is introduced as follows: is supposed to satisfy the periodic boundary conditions in space x ∈ R 2

Numerical Algorithm
Since NFE cannot in general be solved analytically, numerical methods for their approximation have been developed by many researchers. In [5], a computational method is introduced which applies quadrature rule in space to reduce the problem to a system of delay differential equations, which is then solved by a standard algorithm for this kind of equations. A more efficient approach is proposed in [8], [9], where the authors introduce a new algorithm to deal with the convolution kernel of the equation and use the Fast Fourier Transforms to reduce significantly the computational effort required by numerical integration. Recently, a new numerical method for the approximation of twodimensional neural fields was introduced, with use of an implicit second order scheme for the integration in time and applying the Chebyshev interpolation to reduce the dimensions of the matrices [17]. Some applications of this algorithm to Neuroscience problems are discussed in [18].
In the stochastic case, the numerical approximation of NFE becomes a harder challenge.
Numerical algorithms for the stochastic NFE are proposed in [19] and [13].
The numerical integration scheme proposed in the present work is rooted in the Karhunen-Loève presentation in the following form:  (2) under examination. Here and below, these are chosen to satisfy the formula (4) We point out that functions (4) establish an orthonormal basis in the Hilbert space associated with the problem at hand.
Next, the uncorrelated random processes u kl (t) defined on the time interval [0, T] obey the following SDE: (5) whose nonlinear term refers to the double integral, that is (6) the inner product in SDE (5) enjoys the standard definition (7) at any time t, and the scalars λ kl and the stochastic process increments dβ kl (t) come from the infinite series representation of the differential dW(x, t) of the traceclass space-valued Q-Wiener process in 2D-SNFE (2). It is known that the mentioned representation has the following fashion: (8) with the stochastic processes β kl (t) being the mutually independent Brownian motions with zero mean and unit variance and λ kl standing for the eigenvalues of the covariance operator Q. The correlation in this noise obeys the formula for any space points x, y ∈ Ω and time instants t, s > 0.
Here and below, E{·} refers to the expectation operator, and the fixed parameter ξ determines the spatial correlation length. Following [22], we impose the condition ξ << 2L and, then, yield explicitly values of the eigenvalues λ kl of the covariance operator Q in series (8). These are found to be: (9) In what follows, the Galerkin-type approximation underlying our numerical scheme implies that the infinite series (3) is truncated to a finite summation formula. We choose a sufficiently large positive integer K and replace solution (3) to 2D-SNFE (2) with its approximation of the fashion (10) where the deterministic functions v kl (x) are given by rule (4) and the stochastic processes u kl (t) satisfy SDE (5).
Taking into account the 2D-square-form of the domain Ω, we introduce then an equidistant 2D-spacemesh as follows: with its step size h x = 2L/N, where N is a user-supplied quantity of the subdivision steps fulfilled in each direction of the square [−L, L] × [−L, L] underlying 2D-SNFE (2). More formally, this mesh X consists of 2D-vectors in the domain Ω, that is X is the square matrix of size N+1 whose (i, j)th entry is the vector with the coordinates x i 1 = x 0 1 +ih x x j 2 = x 0 2 +jh x . We recall that the number of summation terms in the approximate solution (10) and that of the discretization steps used in the domain Ω must obey the condition N >> K, as explained in [14].
First, we substitute the Galerkin-type solution (10) into SDE (5). Second, we replace the continuous spacedependent functions in equations (5) and (6) with their values calculated at nodes of the mesh X and approximate the integrals arisen by means of the Trapezoidal Rule summations, as follows: Further, we explain our effective MATLAB-based solution technique for attacking the space-discretized version of SDE (5)-(9) rooted in approximations (10)- (12). At first we need to vectorize all values and functions defined at entries of the space-discterizedmeshmatrix X. Here, we employ the built-in MATLAB command reshape for converting any matrix into the requested vectorized fashion. For instance, the mesh X is translated into its column-wise vector representation by the command X(:). We note that formulas (11) and (12) imply that the last row and column of the mesh array X are not used and must be excluded from the calculation, as that explained in [14]. That is why we further remove them from the mesh. This entails that the (i, j)th entry of such a truncated matrix X and the (j × N + i + 1)th entry in its vectorized representation X(:) coincide. Then, we vectorize the scalar stochastic processes u kl (t) in SDE (5) to the form (13) where each random variable u kl (t) is presented by the (k × K + l + 1)th entry of vector (13) for any k, l = 0, 1, . . . , K.
Similarly, we reshape and vectorize all the Fourierbasis-functions (4) evaluated at nodes of our truncated space mesh X. For that, we assemble the rectangular matrix V of size (K + 1) 2 × N 2 whose rows are sorted at the same manner as the stochastic processes u kl (t) in vector (13). More formally, every eigenfunction v kl (x) is calculated first at space points (x i , x j ), i, j=0, 1, . . . , N−1, in use. Second, the resulting matrix derived for the function v kl (X) evaluated at all nodes of our truncated mesh X is further reshaped to the requested vectorized fashion of size N 2 by the MATLAB command v kl (X(:)).
Finally, we utilize the Euler-Maruyama scheme [11, section 10.2] for attacking SDE (15) on some equidistant mesh predefined by the user in the time interval [0, T] of 2D-SNFE (2). The time-integration step size ht is assumed to be sufficiently small. The effective MATLAB-based implementation of the abovepresented-SNFE-solution-method is based on the computation of the term FS(U K (X(:), t)) in SDE (15) by the formula (16) where the time-variant column-vector s is of size N 2 and, for any i, j = 0, 1, . . . , N − 1, its (j × N + i + 1)th entry is given by the value S(u т (t)V(:, j × N + i + 1)) at each time instant t. We recall that the matrix V represents the values of Fourier-basis-functions (4) evaluated at the vectorized-truncated-space-mesh X(:), but the MATLAB command V(:, j × N + i + 1) returns the (j × N + i + 1)th column in the latter matrix. In other words, the mentioned column contains the values of all eigenfunctions v kl (x), k, l = 0, 1, . . . , K, computed at the space position (x i , x j ). Lastly, the notation F refers to the symmetric matrix of size N 2 whose (j × N + i + 1, q × N + p + 1)-entry means the following factor in calculation (12):

Numerical Results
The numerical results presented in [16] support the conjecture that if the external input has appropriate intensity and duration, and if the connectivity kernel is of a certain type (see equation (17)) the neural activity can generate stable multi-bump solutions which contain the information carried by the external signals.
Here, we begin by validating those results with the new numerical algorithm described in Section 2 and, then, we introduce some noise and observe how this noise affects the behavior of the solutions, in particular, how the structure of the multi-bump solutions can be changed.
The connectivity kernel F is of the oscillating type (17) where A, a 1 , l are certain positive constants (see [16]), The activation function S is of the Heaviside type, that is, where b > 0 is a certain threshold.
The function I represents the sum of the external inputs with several components of the form where I 0 is a traveling wave of the form The remaining inputs are localized (have the centers at C i ) and can be written as: (19) where C i ∈ [−L, L], α i , γ i , are given positive numbers, i = 1, ..., n.
As the initial condition, we take We carry out a set of numerical experiments. In each of them, we have a certain sequence of external stimuli, and simulate three different cases, with E = 0 (deterministic case), E = 0.04 (weak noise) and E = 0.08 (strong noise). In all the cases where the stochastic equation is simulated, we perform ns = 100 Monte Carlo runs and display the average value (as an approximation of the mathematical expectation of solution).
The snapshots of the solution at the time t = 0.5 are displayed in Fig. 7. In the deterministic case, we see that this solution reproduces the inputs I 0 (traveling wave) and I 1 . Moreover, at the place where the graphics of the two inputs intersect, we observe suprathreshold activity (U > b = 0.1). In the case of the weak noise (E = 0.04), we observe some perturbation of the signals, but the main features of the solution are preserved. Finally, in the case of the strong noise (E = 0.08), the perturbation of the activity is even higher, but the bump is still well visible at the intersection of the stimuli. Here, in the deterministic case, the effect of the input I 0 (traveling wave) is no more visible; the image of the input I 1 is considerably reduced (as expected, since the corresponding signals have been turned off at the time t = 1.5). However, in the place where the graphics of the two inputs intersect, the onebump solution still exists and carries the information about the nature of the signal and the time at which it has occurred. In the case of the weak noise (E = 0.04), we observe some perturbation of this one-bump pattern, but the coordinates of the peak still reproduce the characteristics of the input. Finally, in the case of the strong noise (E = 0.08), we still observe suprathreshold activity in the place of the intersection, but the peak position deviates slightly from the location of the external input. In this case, we have again the traveling wave I 0 and two localized inputs I 1 and I 2 , with α 1 = 0.12, γ 1 = 1, C = −10 and α 2 = 0.12, γ 2 = 1, C 2 = 10, which correspond to the experiment with two different colors. As in the previous example, the inputs I 0 , I 1 and I 2 are constant in the time interval 0 < t < 1, and, then, are removed. The other parameters have the same values as in the previous example.
The snapshots of the solution at the time t = 1 are displayed in Fig. 9 In the deterministic case, we see that the solution reproduces the inputs I 0 (traveling wave), I 1 and I 2 . Moreover, at the places where the graphic of I 0 intersects I 1 and I 2 , we observe suprathreshold activity (U > b = 0.1). In the case of the weak noise (E = 0.04), we observe some perturbation of the patterns, but the main features of this solution are preserved. Finally, in the case of the strong noise (E = 0.08), the perturbation of the activity pattern is even higher, but the bump is still well visible at the intersections of the stimuli. Here, in the deterministic case, the effect of inputs I 0 (traveling wave), I 1 and I 2 are no more visible (since these inputs have been turned off at time t = 1.) However, similarly to what happens in the previous example, the suprathreshold activity remains at the places of the intersections. In this case, this suprathreshold activity takes the form of a two-bump solution, where the Y-coordinate of each peak corresponds to one color of the signal. Note that these two peaks have the same X-coordinate, which means that the two inputs occurred simultaneously. The effects of the weak noise (E = 0.04) and strong noise (E = 0.08) are less visible in this case in comparison with those observed in the previous example, because we consider the later time instant, that is, t = 5. we employ I 0 , I 1 and I 2 (with the same parameters as in Example 2); then, in [1,5], we implement only the traveling wave I 0 ; finally, in [5,6], we have I 0 , I 3 and I 4 , where I 3 = 1.2I 1 and I 4 = 1.2I 2 . After t = 6, all the inputs are removed.
In Fig. 11, we observe the formation of two bumps, which correspond to the signals that occurred in the time interval [0, 1]. The graph in Fig. 12 shows the formation of a new set of two bumps. These relate to the signals in the time interval [5,6]. In Fig. 13, in the deterministic case, we observe a solution with four bumps, as can be expected: two of them correspond to the first input signal, and two of them to the second one. This solution contains the whole information about the sequence of the stimuli in use. In the case of the strong noise (E = 0.08), Fig.13 shows an interesting feature that is not observed in the previous examples. Instead of 4 bumps, the field activity displays 6 ones, that is, two additional bumps have been formed just as the effect of the implemented noise. Note that this does not happen in the deterministic case, nor in the case of the weak noise.
In order to observe how this picture changes in time, we have carried out computations over a longer time interval, that is, we show the results at the time t = 12 and t = 15. Such outputs are exhibited in Fig.14 and Fig.15. In these graphs, we observe that the two additional bumps, which appear as the result of the noise in use, do not disappear or even decrease in time. This is an example of false memories that can be induced by random noise. Differently from the inputinduced bumps, their activity level is subthreshold (less than b = 0.1). This can be explained as follows: only part of the paths of the stochastic solution contain these two additional bumps (whose height should be above the suprathreshold value), while the remaining ones are similar to the deterministic solution. Since what we show in the picture is the mean solution, the heights of these bumps are averaged on all the paths computed and, hence, these are subthreshold).
It is also worth to mention that the distance between a 'false bump' and the closest real memory representation is approximately the coordinate of one of the local maxima of the connectivity function, that is, the location of false memory receives excitation from existing bumps.

Conclusion
We have described a two-dimensional neural field model which explains how a population of cortical neurons may encode in its self-sustained firing pattern simultaneously the nature and time of sequential stimulus events.
The postulated wave mechanism explains how a nervous system lacking specific sensors for temporal perception may develop neurons that respond to specific interval durations.
The numerical results presented support the conjecture that if the external input has appropriate intensity and duration, and if the connection kernel is of the oscillatory type described here, the neural activity can generate stable multi-bump solutions which reliably contain the information carried by the external signal.
In the case of the weak noise, the structure of the stationary solutions is preserved in the sense that the mean solution of the stochastic equation has the same number of bumps as the solution of the deterministic one has.
In the case of the strong noise, the structure of the stationary solutions can be significantly changed, affecting the precision of the memorized information or even generating false memories.