Robust working memory in a two-dimensional continuous attractor network

Abstract

Continuous bump attractor networks (CANs) have been widely used in the past to explain the phenomenology of working memory (WM) tasks in which continuous-valued information has to be maintained to guide future behavior. Standard CAN models suffer from two major limitations: the stereotyped shape of the bump attractor does not reflect differences in the representational quality of WM items and the recurrent connections within the network require a biologically unrealistic level of fine tuning. We address both challenges in a two-dimensional (2D) network model formalized by two coupled neural field equations of Amari type. It combines the lateral-inhibition-type connectivity of classical CANs with a locally balanced excitatory and inhibitory feedback loop. We first use a radially symmetric connectivity to analyze the existence, stability and bifurcation structure of 2D bumps representing the conjunctive WM of two input dimensions. To address the quality of WM content, we show in model simulations that the bump amplitude reflects the temporal integration of bottom-up and top-down evidence for a specific combination of input features. This includes the network capacity to transform a stable subthreshold memory trace of a weak input into a high fidelity memory representation by an unspecific cue given retrospectively during WM maintenance. To address the fine-tuning problem, we test numerically different perturbations of the assumed radial symmetry of the connectivity function including random spatial fluctuations in the connection strength. Different to the behavior of standard CAN models, the bump does not drift in representational space but remains stationary at the input position.

Appendices

Appendix A

The double integral in (5) can be calculated using the Fourier transforms and Bessel function identities (Bressloff 2012). We start with expressing w(r) as a 2D Fourier transform using polar coordinates

$$\begin{aligned}{} & {} w(r) = \frac{1}{2\pi } \int _{{\mathbb {R}}^{2}} \textrm{e}^{i{\textbf {r}}\cdot {\textbf {k}}}\widehat{w} ({\textbf {k}}) \textrm{d} {\textbf {k}}\nonumber \\{} & {} \quad = \frac{1}{2\pi } \int _{0}^{\infty } \left( \int _{0}^{2\pi }\textrm{e}^{ir\rho \cos \phi } \widehat{w} (\rho ) \textrm{d} \phi \right) \rho \textrm{d} \rho , \end{aligned}$$

(14)

where $\widehat{w}$ denotes the Fourier transform of w and ${\textbf {k}} = (\rho ,\phi )$. Using the integral representation

$$\begin{aligned} \frac{1}{2\pi } \int _{0}^{2\pi }\textrm{e}^{ir \rho \cos \phi } \textrm{d} \phi = J_{0}(\rho r), \end{aligned}$$

(15)

where $J_{0}$ is the Bessel function of the first kind, we express w in terms of its Hankel transform of order zero

$$\begin{aligned} w(r) = \int _{0}^{\infty } \widehat{w} (\rho ) J_{0}(\rho r) \rho \textrm{d} \rho , \end{aligned}$$

(16)

which, when substituted into (5), gives

$$U(r) = V(r) + \int_{0}^{{2\pi }} {\int_{0}^{R} {\left( {\int_{0}^{\infty } {\hat{w}} (\rho )J_{0} (\rho |{\mathbf{r}} - {\mathbf{r}}^{\prime } |)\rho {\text{d}}\rho } \right)r^{\prime}{\text{d}}r^{\prime}{\text{d}}\zeta ^{\prime},} } $$

(17a)

$$V(r) = U(r) - \int_{0}^{{2\pi }} {\int_{0}^{R} {\left( {\int_{0}^{\infty } {\hat{w}} (\rho )J_{0} (\rho |{\mathbf{r}} - {\mathbf{r}}^{\prime } |)\rho {\text{d}}\rho } \right)r^{\prime}{\text{d}}r^{\prime}{\text{d}}\zeta ^{\prime}.} } $$

(17b)

We reverse the order of integration and use the addition theorem

$$\begin{aligned}{} & {} J_{0}(\rho \sqrt{r^{2}+r'^{2}-2rr'\cos \zeta '})\nonumber \\{} & {} \quad =\sum ^{\infty }_{m=0} \epsilon _{m}J_{m}(\rho r)J_{m}(\rho r')\cos m \zeta ', \end{aligned}$$

(18)

where $\epsilon _{0} = 1$ and $\epsilon _{n} = 2$ for $n \ge 1$. Then using the identity $J_{1}(\rho R)R = \rho \int _{0}^{R}J_{0}(\rho r')r' \textrm{d}r'$, we obtain (6). Note that the Fourier transform of (4) is easily calculated using the result that the Fourier transform of $K_0 \left( \dfrac{r}{\sigma } \right) = \dfrac{ 2\pi }{r^2 + \sigma ^2}$.

Appendix B

Using polar coordinates we can rewrite system (8) as

$$\begin{aligned} \begin{aligned}&\lambda \psi (r,\phi )= - \psi (r,\phi )+ \zeta (r,\phi ) \\&\quad + \int _{0}^{2\pi }\textrm{d}\phi ' \int _{0}^{\infty } r'\textrm{d}r' \\&\quad w(\sqrt{r^{2}+r'^{2}-2rr'\cos \phi })\delta (U(r')-\theta ) \psi (r',\phi -\phi '), \end{aligned} \end{aligned}$$

(19a)

$$\begin{aligned} \begin{aligned}&\lambda \zeta (r,\phi )= - \zeta (r,\phi )+ \psi (r,\phi ) \\ {}&\quad - \int _{0}^{2\pi }\textrm{d}\phi ' \int _{0}^{\infty } r'\textrm{d}r'\\&\quad w(\sqrt{r^{2}+r'^{2}-2rr'\cos \phi })\delta (U(r')-\theta ) \psi (r',\phi -\phi '). \end{aligned} \end{aligned}$$

(19b)

We look for solutions of the form

$$\begin{aligned} (\psi (r,\phi ), \zeta (r,\phi ) ) = \textrm{e}^{in\phi } (\psi (r), \zeta (r)), \end{aligned}$$

(20)

where n is the number of modes of the boundary perturbation. System (19) then takes the form

$$\begin{aligned} \begin{aligned}&\lambda \psi (r) \textrm{e}^{in\phi } = - \psi (r)\textrm{e}^{in\phi }+ \zeta (r)\textrm{e}^{in\phi } \\ {}&\quad + \int _{0}^{2\pi }\textrm{d}\phi ' \int _{0}^{\infty } r'\textrm{d}r'\\&\quad w(\sqrt{r^{2}+r'^{2}-2rr'\cos (\phi -\phi ')})\delta (U(r')-\theta ) \psi (r')\textrm{e}^{in(\phi -\phi ')}, \end{aligned} \end{aligned}$$

(21a)

$$\begin{aligned} \begin{aligned}&\lambda \zeta (r) \textrm{e}^{in\phi } = - \zeta (r)\textrm{e}^{in\phi }+ \psi (r)\textrm{e}^{in\phi } \\ {}&\quad - \int _{0}^{2\pi }\textrm{d}\phi ' \int _{0}^{\infty } r'\textrm{d}r'\\&\quad w(\sqrt{r^{2}+r'^{2}-2rr'\cos (\phi -\phi ')})\delta (U(r')-\theta ) \psi (r')\textrm{e}^{in(\phi -\phi ')}. \end{aligned} \end{aligned}$$

(21b)

We set $r=R$ and after dividing both sides by $\textrm{e}^{in\phi }$ we get

$$\begin{aligned} \begin{aligned} \lambda \psi (R)&= - \psi (R)+ \zeta (R) \\ {}&+ \int _{0}^{2\pi }\textrm{d}\phi R w(R\sqrt{2-2\cos \phi }))\frac{\psi (R)\textrm{e}^{-in\phi }}{|U'(R) |}, \end{aligned} \end{aligned}$$

(22a)

$$\begin{aligned} \begin{aligned} \lambda \zeta (R)&= - \zeta (R)+ \psi (R) \\ {}&- \int _{0}^{2\pi }\textrm{d}\phi R w(R\sqrt{2-2\cos \phi }))\frac{\psi (R)\textrm{e}^{-in\phi }}{|U'(R) |}. \end{aligned} \end{aligned}$$

(22b)

The system (22) can be written as

$$\begin{aligned} A \begin{bmatrix} \psi (R) \\ \zeta (R) \\ \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ \end{bmatrix}, \end{aligned}$$

where the matrix A is given by

$$\begin{aligned} A = \begin{bmatrix} \lambda +1-S_{n} &{} -1 \\ -1+S_{n} &{} \lambda +1 \\ \end{bmatrix}, \end{aligned}$$

with

$$\begin{aligned} S_{n} = \frac{R}{|U'(R) |}\int _{0}^{2\pi } w(R\sqrt{2-2\cos \phi })\textrm{e}^{-in\phi }\textrm{d}\phi . \end{aligned}$$

(23)

Then, we find that

$$\begin{aligned} (\lambda +1+S_{n})(\lambda +1)-(-1+S_{n})(-1+S_{n})=0. \end{aligned}$$

(24)

Hence the eigenvalues of A are

$$\begin{aligned}{} & {} \lambda _{-1}= 0, \end{aligned}$$

(25)

$$\begin{aligned}{} & {} \lambda _{n}= -2+S_{n}. \end{aligned}$$

(26)

Note that $\lambda _{n}$ is real, since after setting $\sqrt{2-2\cos \phi } = 2\sin \left( \frac{\phi }{2}\right)$ and rescaling $\phi$ we have

$$\begin{aligned} \text {Im}\{\lambda _{n}\}= -\frac{2R}{|U'(R) |}\int _{0}^{\pi } w(2R\sin (\phi ))\sin (2n\phi )\textrm{d}\phi = 0, \end{aligned}$$

(27)

i.e., the integrand is odd-symmetric about $\frac{\pi }{2}$. Hence,

$$\begin{aligned}{} & {} \lambda _{n}= \text {Re}\{\lambda _{n}\} = -2+ \frac{R}{|U'(R) |} \nonumber \\{} & {} \quad \int _{0}^{2\pi } w(2R\sin (\phi /2))\cos (n\phi )\textrm{d}\phi , \end{aligned}$$

(28)

with the integrand even-symmetric about $\frac{\pi }{2}$.

We then evaluate the integral in (28) using Bessel functions

$$\begin{aligned} & \int_{0}^{{2\pi }} w (2R\sin (\phi ^{\prime}/2))\cos (n\phi ^{\prime}){\text{d}}\phi ^{\prime} \\ & = \int_{0}^{{2\pi }} {\left( {\int_{0}^{\infty } {\hat{w}} (\rho )J_{0} (\rho (2R\sin (\phi ^{\prime}/2)))\rho {\text{d}}\rho } \right)} \cos \phi ^{\prime}{\text{d}}\phi ^{\prime} \\ & = 2\pi \int_{0}^{\infty } {\hat{w}} (\rho )J_{n} (\rho R)J_{n} (\rho R)\rho {\text{d}}\rho . \\ \end{aligned}$$

(29)

We differentiate (6a) with respect to r, and, knowing that $U(r)+V(r)=K$ we have

$$\begin{aligned} U'(R) = - \pi R \int _{0}^{\infty } \widehat{w} (\rho ) J_{1}(\rho R) J_{1}(\rho R)\rho \textrm{d} \rho . \end{aligned}$$

(30)

We can now write the eigenvalues of A as (9) and (10).

Appendix C

Numerical simulations of the model were done in MATLAB using a forward Euler method with uniform spatial mesh with $\textrm{dx}=0.05$ and time step $\textrm{dt}=0.01$. To compute the two-dimensional spatial convolution of w and f we employ a two-dimensional fast Fourier transform (2D FFT), using MATLAB’s in-built functions fft2 and ifft2 to perform the Fourier transform and the inverse Fourier transform, respectively. Periodic boundary conditions are used. By choosing a sufficiently large domain size, we make sure that the localized patterns evolve sufficiently far from the boundaries.

For performing numerical continuation, we use the method described in (Rankin et al. 2014) and adapt MATLAB code available in (Avitabile 2016). The main advantage of this method is that it can be applied directly to the full integral model. This is possible due to the usage of Newton-GMRES solvers combined with a fast Fourier transform (FFT) employed for computing the convolution term (Rankin et al. 2014).