Synaptic mechanisms of interference in working memory

TL;DR

This paper presents a probabilistic inference model and a recurrent network mechanism explaining how interference from previous trials biases working memory responses, with implications for understanding neural dynamics and memory retention.

Contribution

It introduces a novel probabilistic model linked to a recurrent network with short-term facilitation, explaining trial-to-trial bias dynamics in working memory tasks.

Findings

Model accurately predicts bias in response locations.

Dynamic connectivity improves memory retention for correlated targets.

Low-dimensional description captures trial-to-trial bias dynamics.

Abstract

Information from preceding trials of cognitive tasks can bias performance in the current trial, a phenomenon referred to as interference. Subjects performing visual working memory tasks exhibit interference in their trial-to-trial response correlations: the recalled target location in the current trial is biased in the direction of the target presented on the previous trial. We present modeling work that (a) develops a probabilistic inference model of this history-dependent bias, and (b) links our probabilistic model to computations of a recurrent network wherein short-term facilitation accounts for the dynamics of the observed bias. Network connectivity is reshaped dynamically during each trial, providing a mechanism for generating predictions from prior trial observations. Applying timescale separation methods, we can obtain a low-dimensional description of the trial-to-trial bias based on the history of target locations. The model has response statistics whose mean is centered at the true target location across many trials, typical of such visual working memory tasks. Furthermore, we demonstrate task protocols for which the plastic model performs better than a model with static connectivity: repetitively presented targets are better retained in working memory than targets drawn from uncorrelated sequences.

Figures

Tables

Table 1. Table 1: Variables and parameters of the probabilistic inference model.

symbol	description
${\theta }_n$	target observed on trial $n$
$s_{n+1}\left(\theta \right)$	unknown distribution of targets observer is attempting to infer before trial $n+1$
$f_{{\theta }_n}\left(\theta \right)$	probability of observing target $\theta$ in the present trial
$ϵ$	assumed change rate of the environment observer uses to discount evidence
$L_{n+1,\theta }$	predictive distribution: observer’s expected probability of seeing target $\theta$ in trial $n+1$,
	given their past observations ${\theta }_{1:n}$ and an assumed change rate $ϵ$ of the environment
$r_n$	run-length: number of trials preceding trial $n+1$ with the same target distribution
${\overline{P}}_0$	uniform prior: distribution for which all targets $\theta \in [-{180}^{\circ },{180}^{\circ })$ have same probability

Table 2. Table 2: Variables and parameters of the particle evolving on a potential model.

symbol	description
$\theta \left(t\right)$	observer’s estimate of the remembered target during the delay-period
${𝒰}_{n+1}\left(\theta \right)$	potential determining the systematic dynamics of $\theta$ during the delay-period
${\sigma }_{\theta }$	standard deviation of the dynamic fluctuations
$p_{n+1}(\theta ,t)$	probability density function (pdf) for finding the observer’s estimate at $\theta$ at time $t$
${\overline{p}}_{n+1}\left(\theta \right)$	stationary density: pdf of observer’s estimate $\theta$ in the long time limit, $t\to \infty$

Table 3. Table 3: Variables and parameters of the recurrent network model.

symbol	description
$u(x,t)$	synaptic input to location $x$ at time $t$
$q(x,t)$	short term facilitation (STF) variable at location $x$ and time $t$: increases the strength of
	connectivity from neurons at location $x$ in response to their activation
${\tau }_u$	timescale of synaptic excitation, set to 10ms in simulations [60]
$w\left(x-y\right)$	baseline strength and polarity of effective connectivity from neurons at location $y$ to $x$
$F\left(u\right)$	nonlinearity that converts synaptic input to output firing rate of local neural population
$J(x,t)$	spatiotemporal input: external input due to visual stimuli plus dynamic fluctuations
$\tau$	timescale of STF, set to 1ms in simulations [17]
$\beta$	onset rate of STF, set to 0.01 in simulations [18]
$q_{+}$	maximal increase in synaptic utilization, set to 2 in simulations [18]

Equations

$$L_{n+1,\theta }={\bar{\mathrm{P}}}_0\cdot \left[\frac{(1-ϵ{)}^n}{\mathrm{P}({\theta }_{1:n})}\prod _{j=1}^n{f}_{{\theta }_j}(\theta )+ϵ\sum _{r=0}^{n-1}\frac{(1-ϵ{)}^r}{\mathrm{P}({\theta }_{n-r+1:n})}\prod _{j=n-r+1}^n{f}_{{\theta }_j}(\theta )\right],$$

$$L_{n+1,\theta }=\frac{{\bar{\mathrm{P}}}_0}{\mathrm{P}({\theta }_{1:n})}\prod _{j=1}^n{f}_{{\theta }_j}(\theta ).$$

$$L_{n+1,\theta }=\frac{\mathrm{P}({\theta }_{1:n-1})}{\mathrm{P}({\theta }_{1:n})}{f}_{{\theta }_n}(\theta )L_{n,\theta },$$

$${\tilde{L}}_{n+1,\theta }=ϵ{\bar{\mathrm{P}}}_0+(1-ϵ)f_{{\theta }_n}(\theta ).$$

$$\mathrm{d}\theta (t)=-\frac{\mathrm{d}\mathcal{U}(\theta (t))}{\mathrm{d}\theta }\mathrm{d}t+{\sigma }_{\theta }\mathrm{d}\xi (t).$$

$${\mathcal{U}}_{n+1}(\theta )\propto -f_{{\theta }_n}(\theta ).$$

$$\mathrm{d}\theta (t)$$

$$\tau {\dot{\theta }}_q(t)$$

$$s_{n+1}(\theta )={\mathcal{N}}_s\exp \left[\sum _{j=1}^N{a}_j\cos ({\omega }_j\cdot (\theta -{\psi }_j))\right],$$

$$s_{n+1}(\theta )={\mathcal{N}}_s\exp \left[a_1\cos ({\omega }_1\cdot (\theta -{\psi }_1))\right],$$

$$f_{{\theta }^{\prime }}(\theta )={\mathcal{N}}_{\theta }\exp \left[\sum _{j=1}^N{a}_j\cos ({\omega }_j\cdot (\theta -{\theta }^{\prime }))\right],$$

$$L_{n+1,\theta }=\sum _{r=0}^n\mathrm{P}({\theta }_{n+1}|r_n=r,{\theta }_{1:n}^r)\mathrm{P}(r_n=r|{\theta }_{1:n},ϵ),$$

$$\mathrm{P}({\theta }_{n+1}|r_n=r,{\theta }_{1:n}^r)=\frac{{\bar{\mathrm{P}}}_0}{\mathrm{P}({\theta }_{n-r+1:n})}\prod _{j=n-r+1}^n{f}_{{\theta }_j}(\theta ),$$

$$\displaystyle{\rm P}(r_{n}=r|\theta_{1:n},\epsilon)={\rm P}(r_{n}=r|\epsilon)=\left\{\begin{array}[]{cl}\epsilon(1-\epsilon)^{r},&r<n,\\ (1-\epsilon)^{n},&r=n.\end{array}\right.$$

$$L_{n+1,\theta }={\bar{\mathrm{P}}}_0\cdot \left[\frac{(1-ϵ{)}^n}{\mathrm{P}({\theta }_{1:n})}\prod _{j=1}^n{f}_{{\theta }_j}(\theta )+ϵ\sum _{r=0}^{n-1}\frac{(1-ϵ{)}^r}{\mathrm{P}({\theta }_{n-r+1:n})}\prod _{j=n-r+1}^n{f}_{{\theta }_j}(\theta )\right].$$

$${\tilde{L}}_{n+1,\theta }={\mathcal{N}}_s\cdot \left[(1-nϵ)\prod _{j=1}^n\frac{f_{{\theta }_j}(\theta )}{{\bar{\mathrm{P}}}_0}+ϵ\sum _{r=0}^{n-1}\prod _{j=n-r+1}^n\frac{f_{{\theta }_j}(\theta )}{{\bar{\mathrm{P}}}_0}\right],$$

$$L_{n+1,\theta }\approx ϵ{\bar{\mathrm{P}}}_0\left[1+\frac{1-ϵ}{\mathrm{P}({\theta }_n)}{f}_{{\theta }_n}(\theta )\right].$$

$${\tilde{L}}_{n+1,\theta }=\frac{{\bar{\mathrm{P}}}_0+(1-ϵ)f_{{\theta }_n}(\theta )}{2-ϵ},$$

$${\tilde{L}}_{n+1,\theta }=ϵ{\bar{\mathrm{P}}}_0+(1-ϵ)f_{{\theta }_n}(\theta ),$$

$$L_{n+1,\theta }$$

$$\frac{\partial p_{n+1}(\theta ,t)}{\partial t}=\frac{\partial }{\partial \theta }\left[\frac{\mathrm{d}{\mathcal{U}}_{n+1}(\theta )}{\mathrm{d}\theta }{p}_{n+1}(\theta ,t)\right]+\frac{{\sigma }_{\theta }^2}{2}\frac{{\partial }^2{p}_{n+1}(\theta ,t)}{\partial {\theta }^2},$$

$$\underset{t\to \infty }{\lim }{p}_{n+1}(\theta ,t)={\bar{p}}_{n+1}(\theta )=L_{n+1,\theta }.$$

$${\bar{p}}_{n+1}(\theta )={\chi }_{n+1}\exp \left[-\frac{2{\mathcal{U}}_{n+1}(\theta )}{{\sigma }_{\theta }^2}\right],$$

$${\mathcal{U}}_{n+1}(\theta )=\frac{{\sigma }_{\theta }^2}{2}\ln \frac{{\chi }_{n+1}}{L_{n+1,\theta }}.$$

$${\mathcal{U}}_{n+1}(\theta )$$

$$\approx \frac{{\sigma }_{\theta }^2}{2}\left[\ln \frac{{\chi }_{n+1}}{{\bar{\mathrm{P}}}_0}-(1-ϵ)\frac{f_{{\theta }_n}(\theta )-{\bar{\mathrm{P}}}_0}{{\bar{\mathrm{P}}}_0}\right],$$

$${\mathcal{U}}_{n+1}(\theta )\propto -f_{{\theta }_n}(\theta ).$$

$${\tau }_u\mathrm{d}u(x,t)$$

$$\tau \dot{q}(x,t)$$

$$\displaystyle I(x,t)=\left\{\begin{array}[]{lll}I_{0}\exp\left[I_{1}(\cos(x-\theta_{n})-1)\right],&t\in[t_{n},t_{n}+T_{C}),\\ 0,&t\in[t_{n}+T_{C},t_{n}+T_{C}+T_{D}^{n}),\\ -I_{R},&t\in[t_{n}+T_{C}+T_{D}^{n},t_{n}+T_{C}+T_{D}^{n}+T_{A}),\\ 0,&t\in[t_{n}+T_{C}+T_{D}^{n}+T_{A},t_{n+1}),\end{array}\right.$$

$$U(x)={\int }_{-180}^{180}w(x-y)F(U(y))\mathrm{d}y.$$

$$\mathrm{d}\Phi (x,t)-\mathcal{L}\Phi (x,t)\mathrm{d}t$$

Full text

**Synaptic mechanisms of interference in working memory **

Zachary P Kilpatrick1,2

1 Department of Applied Mathematics, University of Colorado, Boulder CO, USA

2 Department of Physiology and Biophysics, University of Colorado School of Medicine, Aurora CO, USA

***** zpkilpat@colorado.edu

Summary

Information from preceding trials of cognitive tasks can bias performance in the current trial, a phenomenon referred to as interference. Subjects performing visual working memory tasks exhibit interference in their trial-to-trial response correlations: the recalled target location in the current trial is biased in the direction of the target presented on the previous trial. We present modeling work that (a) develops a probabilistic inference model of this history-dependent bias, and (b) links our probabilistic model to computations of a recurrent network wherein short-term facilitation accounts for the dynamics of the observed bias. Network connectivity is reshaped dynamically during each trial, providing a mechanism for generating predictions from prior trial observations. Applying timescale separation methods, we can obtain a low-dimensional description of the trial-to-trial bias based on the history of target locations. The model has response statistics whose mean is centered at the true target location across many trials, typical of such visual working memory tasks. Furthermore, we demonstrate task protocols for which the plastic model performs better than a model with static connectivity: repetitively presented targets are better retained in working memory than targets drawn from uncorrelated sequences.

Keywords: working memory; interference; short-term facilitation; probabilistic inference; bump attractor

Introduction

Parametric working memory experiments are a testbed for behavioral biases and errors, and help identify neural mechanisms that underlie them [1, 2]. In visuospatial working memory, subjects identify, store, and recall target locations in trials lasting a few seconds. Response errors are normally distributed [3, 4], and tend to accumulate during the delay-period, while subjects retain the target location in memory [1, 4]. Complementary neural recordings suggest these working memories are implemented in circuits comprised of stimulus-tuned neurons with slow excitation and broad inhibition [5, 6]. Persistent activity emerges as a tuned pattern of activity called a bump state, whose peak encodes the remembered target position [7, 4].

Neuronal studies of visual working memory typically focus on population activity within a single trial, ignoring serial correlations across trials [8]. Several authors have identified behavioral biases that cause the previous trial’s visual target to interfere with the subject’s response on the subsequent trial [9]. For instance, in delayed match-to-sample tests, false alarms occur more often when comparison stimuli match samples from previous trials [10]. Originally, interference was observed in verbal working memory tasks [11], and evidence suggests the effect impacts working memory capacity [12]. One consistent observation is that interference is reduced by increasing the time interval between trials [13, 9], suggesting the effect persists for a few seconds. Investigations of interference in visuospatial working memory reveal other effects: Increasing the delay-period of working memory trials increases the bias strength, and responses are biased in the direction of the stimulus from the previous trial [9].

Our study focuses on why and how interference biases arise visuospatial working memory. First, what evidence accumulation strategy accounts for the bias introduced by the previous trial’s target? We will show these biases emerge in observers using sequential Bayesian updating to predict the location of the next target. Such models are obtained by iteratively applying Bayes’ rule to a stream of noisy measurements, updating an observer’s belief of the most likely choice. In changing environments, older measurements are discounted at a rate that increases with the assumed change rate of the environment [14, 15]. In our model, the sequence of targets observed on each trial is used to predict the next target. When subjects assume the environment changes rapidly, only the most recent target is used to make their prediction, leading to suboptimal inference of the subsequent target [16].

What neurophysiological processes could account for intertrial biases? Both the decay and activation timescales of the bias appear to be on the order of seconds. We propose short-term facilitation (STF), which acts on the timescale of seconds [17], can account for the dynamics of the bias. In a recurrent network that sustains persistent activity during a delay-period in the form of an activity bump, facilitated synapses from neurons tuned to the previous target attract the activity bump in the subsequent trial. Previous models identified STF as a possible mechanism for lengthening the timescale of working memory [18, 19, 20]. Our study proposes interference arises as a result of an irrelevant working memory remaining from the previous trial.

Our neurocomputational model accounts for recent observations of interference from visual working memory experiments [9], and makes novel predictions linking behavioral responses to corresponding neural and synaptic mechanisms. The separation in timescales between the neural activity dynamics and STF variable allows us to derive a low-dimensional model describing the bump’s interaction with the network’s evolving synaptic weights. We find that protocols with a uniform distribution of possible target angles lead to response error distributions that are normally distributed about zero as found previously [1, 4, 8]. Conversely, target protocol sequences with strong serial correlations can lead to a biased distributions in recalled target positions. Such biases may be advantageous in more complex tasks, where information from previous trials provides information about the target location in subsequent trials, as we show. Finally, we demonstrate that a recurrent network with STF supports bump attractors whose diffusion time course possesses two distinct phases, a prediction we propose to validate our model.

Results

Our study presents two frameworks for generating interference in a sequence of visual working memory trials. Both models use information about the target location on the previous trial to bias the response on the current trial. First, we develop a probabilistic inference model that predicts the distribution of possible target angles on the current trial based on observations of past trials. When the observer assumes the environment changes rapidly, the predictive distribution is primarily shaped by the previous trial’s target. Second, we analyze a recurrent network model with STF wherein a localized bump of activity represents the observer’s belief on the current trial and the spatial profile of STF represents the observer’s evolving predictive distribution for the subsequent target. We show the attractor structure of the network model can be directly related to the predictive distribution of the inference model.

Interference in a visual working memory task

We focus specifically on an oculomotor delayed-response task with a single target presented in each trial [1, 4]. On each trial, the subject views a target $\theta_{n}$ during a short cue period (Fig. 1A). They must remember the target location during a delay-period and saccade to the remembered location at the end. Response errors depend on the previous trial in three distinct ways, originally reported in Papadimitriou et al. (2015) [9]: (1) responses are attracted to the location of the previous target, graded with the difference between the current and previous target (Fig. 1B); (2) the bias decreases as the interval between trials is increased (Fig. 1C); and (3) the bias increases as the delay-period increases (Fig. 1D). As we will show, these biases are captured by a model of an observer that predicts the current target based on the previous target. These effects also emerge in a recurrent network model with slow excitation, subject to STF, and broad inhibition (Fig. 1E). This network represents the memory of the presented target as a bump of neural activity, which drifts in the direction of the target presented on the previous trial (Fig. 1F). Before analyzing the mechanics of this network model in more detail, we derive a probabilistic inference model that accounts for these tendencies.

Inference model for updating target predictions

Our goal now is to show that sequential Bayesian updating can account for interference observed in working memory, given specific constraints on a probabilistic updating algorithm. The observer attempts to predict the probability of observing target angle $\theta_{n+1}=\theta$ in trial $n+1$, given the targets $\theta_{1:n}=\{\theta_{1},\theta_{2},...,\theta_{n}\}$ observed in the previous $n$ trials (Fig. 2A). However, the target $\theta_{j}$ on the $j^{\text{th}}$ trial will only help predict the target $\theta_{n+1}$ on the $n+1^{\text{th}}$ trial if the distribution $s_{n+1}(\theta)$ from which targets are drawn remains the same between trial $j$ and trial $n+1$ [21]. The observer assume the distribution from which presented targets are drawn changes stochastically at a fixed rate $\epsilon:={\rm P}(s_{n+1}(\theta)\not\equiv s_{n}(\theta))$. Most visual working memory protocols fix the distribution of target angles throughout the task ($\epsilon=0$) [1, 22, 4, 9], as we do for most of our study, so the observer employs a potentially incorrect model to estimate this distribution ($\epsilon>0$). Subjects in psychophysical tasks can have a strong bias toward assuming environments change on a timescale of several seconds [14], and this bias is not easily trained away [16, 23]. Combining these features of the model, the observer updates their predictive distribution for the target during the $(n+1)^{\text{th}}$ trial.

Our algorithm is based on models that compute a predictive distribution for a stochastically moving target, given a sequence of noisy observations [24, 21]. The predictive distribution is computed using sequential analysis [25, 15]: Prior to trial $n+1$, the observer has seen $n$ targets $\theta_{1:n}=\{\theta_{1},\theta_{2},...,\theta_{n}\}$. The observer computes $f_{\theta^{\prime}}(\theta):={\rm P}(\theta_{n+1}=\theta|\theta_{j}=\theta^{\prime},s_{n+1}(\theta)\equiv s_{j}(\theta))$ (Fig. 2B), the probability of observing the target $\theta_{n+1}$ in the $(n+1)^{\text{th}}$ trial assuming the underlying probability distribution from which targets are sampled does not change from trial $j$ to $n+1$ ($s_{n+1}(\theta)\equiv s_{j}(\theta)$), for each trial $j=1,...,n$. The true distribution of target angles $\theta$ remains uniform throughout most our study, so the observer is applying suboptimal inference. Further details of our Bayesian nonparametric model are given in Methods.

The observer thus computes a predictive distribution $L_{n+1,\theta}={\rm P}(\theta_{n+1}=\theta|\theta_{1:n},\epsilon)$, using the previous targets $\theta_{1:n}$ (Fig. 2A) to predict the subsequent target $\theta_{n+1}$. If the observer assumes the distribution $s_{n+1}(\theta)$ from which targets are drawn in trial $n+1$ changes stochastically with a rate $\epsilon\in(0,1)$, recent observations will be weighted more in determining $L_{n+1,\theta}$ [21, 14, 15]. Each observation $\theta_{j}$ contributes to the current estimate of $L_{n+1,\theta}$ via the probability $f_{\theta_{j}}(\theta)$ (Fig. 2B). Observations are weighted by assuming the observer has a fixed belief about the value $\epsilon$, specifying the average number of trials they expect the distribution $s_{n}(\theta)$ to remain the same. Leveraging techniques in probabilistic inference (See Methods), we find

[TABLE]

where $\bar{{\rm P}}_{0}:=1/360$ is the uniform density for $-180^{\circ}\leq\theta<180^{\circ}$. To understand Eq. (1), it is instructive to examine limits of the parameter $\epsilon$ that admit approximations or exact updates.

Static environments ($\epsilon\to 0$). In the limit $\epsilon\to 0$, the observer assumes the environment is static, so the predictive distribution is comprised of equal weightings of each observation (See Fig. 2C) [26, 27]:

[TABLE]

As has been shown previously, Eq. (2) can be written iteratively [28]:

[TABLE]

suggesting such a computation could be implemented and represented by neural circuits. Temporal integration of tuned inputs has been demonstrated in both neural recordings [29, 30, 31] and circuit models [32, 27, 28] of decision-making tasks. Most oculomotor delayed-response tasks use a distribution of targets $s(\theta)$ that is constant across trials [1, 22, 4, 9]. Therefore, Eq. (2) is the optimal strategy for obtaining an estimate of $s(\theta)$, assuming the observer has a correct representation of the probability $f_{\theta_{j}}(\theta)$. For instance, if the distribution $s(\theta)$ were peaked, repeated observations $\theta_{1:n}$ would gradually improve the observer’s estimate of that peak in Eq. (2). In changing environments ($\epsilon>0$), recently observed targets are weighted more strongly than older targets, and the predictive distribution should down-weight the influence of past targets at a rate that increases with $\epsilon$ [15].

Rapidly-changing environment ($\epsilon\approx 1$). Our work focuses on the limit where the environment changes rapidly, $\epsilon\approx 1$ ($0<(1-\epsilon)\ll 1$), to account for biases that depend on the previous trial’s target $\theta_{n}$ (See Methods for other cases). In this case, the predictive distribution for trial $n+1$ is a single peaked function centered at $\theta_{n}$ (Fig. 2C). The observer assumes the environment changes fast enough that each subsequent target is likely drawn from a new distribution ($s_{n+1}(\theta)\not\equiv s_{n}(\theta)$). This is a suboptimal strategy, but matches the typical trends of interference in working memory. Applying this assumption to Eq. (1), the formula for $L_{n+1,\theta}$ is dominated by terms of order $(1-\epsilon)$ and larger. Truncating to ${\mathcal{O}}(1-\epsilon)$ and normalizing the update equation (See Methods) then yields

[TABLE]

Thus, the dominant contribution from $\theta_{1:n}$ to $\tilde{L}_{n+1,\theta}$ is the target $\theta_{n}$ observed during the previous trial $n$ (Fig. 2C), similar to the behavioral findings of Papadimitriou et al. (2015) [9].

Note, sequential computations are trivial in the limit of a constantly-changing environment $\epsilon\to 1$, since the observer assumes the environment is reset after each trial. Prior observations provide no information about the present distribution, so the predictive distribution is always uniform: $L_{n+1,\theta}\equiv\bar{{\rm P}}_{0}$.

In summary, a probabilistic inference model that assumes the distribution of targets is predictable over short timescales leads to response biases that depend mostly on the previous trial. We now demonstrate that this predictive distribution can be incorporated into a low-dimensional attractor model which describes the degradation of target memory during the delay-period of visual working memory tasks [33, 7, 34].

Incorporating suboptimal predictions into working memory

We model the loading, storage, and recall of a target angle $\theta$ using a low-dimensional attractor model spanning the space of possible target angles $\theta\in[-180,180)^{\circ}$. These dynamics can be implemented in recurrent neuronal networks with slow excitation and broad inhibition [35, 6, 4]. Before examining the effects of neural architecture, we discuss how to incorporate the predictive distribution update, Eq. (3), into an associated low-dimensional model. Our analysis links the update of the predictive distribution to the spatial organization of attractors in a network. Importantly, working memory is degraded by dynamic fluctuations, so the stored target angle wanders diffusively during the delay-period [6, 34, 4].

During the delay-period of a single trial, the stored target angle $\theta(t)$ evolves according to a stochastic differential equation [7]:

[TABLE]

Here $\theta(t)$ is restricted to the periodic domain $\theta\in[-180,180)^{\circ}$ and ${\rm d}\xi$ is a standard white noise process. Eq. (4) can be derived as the low-dimensional projection for the location of a bump attractor in a recurrent network. The potential gradient $-{\mathcal{U}}^{\prime}(\theta)$ models spatial heterogeneity in neural architecture that shapes attractor dynamics (Fig. 3A). During trial $n+1$, we label the potential ${\mathcal{U}}(\theta)\equiv{\mathcal{U}}_{n+1}(\theta)$. Classic models of bump attractors on a ring assume distance-dependent connectivity [35, 6]. The case ${\mathcal{U}}_{n+1}^{\prime}(\theta)\not\equiv 0$ accounts for spatial heterogeneity in connectivity that may arise from a combination of training and synaptic plasticity [7], or random components of synaptic architecture [36]. The potential landscape can be updated during each trial, so at the beginning of trial $n+1$ it has the form ${\mathcal{U}}_{n+1}(\theta)$. When ${\mathcal{U}}_{n+1}(\theta)\equiv 0$, the potential is flat, so $\theta(t)$ evolves along a line attractor [33]. On the other hand, when the potential is heterogeneous, ${\mathcal{U}}_{n+1}(\theta)\not\equiv 0$, $\theta(t)$ tends to drift toward one of a finite number of discrete attractors [7, 34]. We will incorporate a process whereby previous targets are used to update the potential, so ${\mathcal{U}}_{n+1}(\theta)$ is typically heterogeneous. The observer sees the target at the beginning of trial $n+1$, $\theta(0)=\theta_{n+1}$ (Fig. 3A), and the angle $\theta(t)$ evolves according to Eq. (4) during the delay-period, lasting $T_{D}$ time units. After the delay-period, $\theta(T_{D})$ is the recalled angle. Depending on the underlying potential ${\mathcal{U}}_{n}(\theta)$, there will be a strong bias to a subset of possible targets.

We derive a correspondence between the probabilistic inference model and attractor model by assuming stationarity of ${\mathcal{U}}_{n+1}(\theta)$ within each trial (See Methods). In the recurrent network model (Fig. 1E), we take these within-trial dynamics into account. Freezing ${\mathcal{U}}_{n+1}(\theta)$ during a trial allows us to relate the statistics of the position $\theta(t)$ to the shape of the potential. Specifically, we relate the stationary density of Eq. (4) to the desired predictive distribution $L_{n+1,\theta}$ (See Methods). Thus, if information about the current trial’s target $\theta_{n+1}$ is degraded, the probability distribution associated with the recalled target angle $\theta$ is $L_{n+1,\theta}$. Focusing on interference trends in Fig. 1, we aim to have the attractor structure of Eq. (4) represent the predictive distribution in Eq. (3). Our calculations relate the potential function in trial $n+1$ to the probability generated by the trial $n$ target (Fig. 3A) as

[TABLE]

The potential ${\mathcal{U}}_{n+1}(\theta)$ can be implemented by a decaying plasticity process that facilitates synapses from neurons tuned to the previous target $\theta_{n}$. As we will show, this can be accomplished via STF (Fig. 3B).

Short-term facilitation generates interference in working memory

We now show a neuronal network model describing neural activity $u(x,t)$ subject to STF $q(x,t)$ can incorporate predictive distribution updates derived above. Predictions are stored in the dynamically changing synaptic weights of a recurrent neuronal network. The recurrent network model spatially labels neurons according to their target orientation preference, determining the distance-dependent structure of inputs to the network. This is captured by a network with local excitation and effective inhibition that is fast and broad. Connectivity is shaped dynamically by STF (Fig. 1E). See Methods for more details.

A sequence of delayed-response protocols is implemented in the recurrent network by specifying a spatiotemporal input $I(x,t)$ across trials (top of Fig. 1F). Each trial has a cue period of time length $T_{C}$; a delay-period of time length $T_{D}^{n+1}$; and an intertrial period of time length $T_{I}^{n+1}$ before the next target is presented. The network receives a peaked current centered at the neurons preferring the presented target angle $\theta_{n+1}$ during the cue period of trial $n+1$; no external input during the delay-period; and a strong inactivating current after the delay-period [6, 4]. The resulting bump attractor drifts in the direction of the bump from trial $n$, due to the STF at the location of the trial $n$ bump (Figs. 1F and 3B).

The mechanism underlying intertrial bias is determined by projecting our recurrent network model to a low-dimensional system that extends Eq (4) to account for STF. To reduce the recurrent network, we project the fast dynamics of bump solutions to an equation for the bump’s position $\theta(t)$ in trial $n+1$ [19, 34]. The STF variable $q(x,t)$ determines an evolving potential function ${\mathcal{U}}(\theta,t)$ that shapes the bump’s position (Fig. 4). We use timescale separation methods (See Methods) to derive a set of stochastic differential equations, which approximates the motion of the bump’s position $\theta(t)$ and the location of STF $\theta_{q}(t)$:

[TABLE]

during trial $n+1$ ($t_{n}<t<t_{n+1}$). The slowly-evolving potential gradient $-\frac{\partial}{\partial\theta}{\mathcal{U}}(\theta,t)$ shaping the dynamics of $\theta(t)$ is a mixture of STF contributions from trial $n$ (decaying ${\mathcal{A}}_{n}(t)$) and trial $n+1$ (increasing ${\mathcal{A}}_{n+1}(t)$). The bump position $\theta(t)$ moves towards the minimum of this dynamic potential, ${\rm argmin}_{\theta}\left[{\mathcal{U}}(\theta,t)\right]$ (Fig. 4). The variable $\theta_{q}(t)$ is the location of STF originating in trial $n+1$, and its position slowly moves toward the bump location $\theta(t)$. The parametrized timescale $\tau$ of STF is inversely related to the observer’s perceived environmental change rate $\epsilon$ in Eq. (3), since increasing $\epsilon$ corresponds to decreasing $\tau$.

The presence of STF provides two contributions to the slow dynamics of the bump position $\theta(t)$. The memory of the previous trial’s target $\theta_{n}$ is reflected by the potential $\bar{\mathcal{U}}(\theta-\theta_{n})$, whose effect decays slowly during trial $n+1$. This attracts $\theta(t)$, but the movement of $\theta(t)$ towards $\theta_{n}$ is slowed by the onset of the STF variable initially centered at $\theta_{n+1}$. The STF variable’s center-of-mass $\theta_{q}(t)$ slowly drifts towards $\theta_{n}$, which allows $\theta(t)$ to drift there as well, $\bar{\mathcal{U}}(\theta-\theta_{q}(t))$. This accounts for the slow build-up of the bias that increases with the length of the delay-period [9].

Target- and time-dependent trends match experimental observations

We now demonstrate that the bias observed in the visual working memory experiments of Papadimitriou et al. (2015) can be accounted for by our recurrent network model (Fig. 1E) and our low-dimensional description of bump motion dynamics (Fig. 4). To represent a sequence of working memory trials, targets $(\theta_{1},\theta_{2},\theta_{3},...)$ were presented to the recurrent network, and the center-of-mass of the bump was recorded at the end of each delay-period, representing the response $(r_{1},r_{2},r_{3},...)$ (See Methods). The bias of responses was determined by computing the difference between the response and the presented target, $r_{n}-\theta_{n}$. Means and variances of the bias were determined under each condition.

Our results are summarized in Fig. 5, focusing on three conditions considered by Papadimitriou et al. (2015) [9]. First, we calculated the bias when conditioning on the angle between the trial $n$ and trial $n+1$ targets, $\theta_{n}-\theta_{n+1}$ (Fig. 5A). Positive (negative) angles lead to positive (negative) bias; i.e. the bump drifts in the direction of the previous target $\theta_{n}$. The bias is graded with the difference, $\theta_{n}-\theta_{n+1}$. To expose this effect, we averaged across trials, since the recurrent network incorporates dynamic input fluctuations, as in bump attractor models of visuospatial working memory [6, 4]. We also calculated the peak bias as a function of the intertrial interval (ITI), the time between the trial $n$ response ($r_{n}$) and the trial $n+1$ target presentation $\theta_{n+1}$. Consistent with Papadimitriou et al. (2015) [9], the peak bias decreased with the ITI (Fig. 5B). The mechanism for this decrease is the slow decay in the STF of synapses utilized in the previous trial. Finally, the peak bias increased with the delay within a trial, since persistent activity was slowly attracted to the location of the previous target (Fig. 5C). This slow saturation arises due to the slow kinetics of STF within a trial.

Not only did our recurrent network model recapitulate the main findings of Papadimitriou et al. (2015) [9], we also found our low-dimensional description of the bump and STF variable dynamics had these properties (blue curves in Fig. 5). The mechanics underlying the bias can be described with a model of a particle evolving in a slowly changing potential (Fig. 4), shaped by the dynamics of STF. Having established a mechanism for the bias, we consider how different protocols determine the statistics of responses, not conditioned with sequential trial information.

Task protocol shapes ensemble statistics

Visual working memory tasks are often designed such that sequential target locations are uncorrelated [6, 4]. In such protocols, there is no advantage in using previous trial information to predict targets within the current trial. Nonetheless, these biases persist in the intertrial response correlations discussed in Papadimitriou et al. (2015) [9] and Fig. 5. On the other hand, interference might be advantageous for tasks with correlations between successive target angles, $\theta_{n}$ and $\theta_{n+1}$. Consider object motion tracking tasks with transiently occluded objects [37], so the object’s location prior to occlusion predicts its subsequent location when it reappears. Memory of object location that persists beyond a single trial can therefore be useful for naturally-inspired tasks.

We demonstrate this idea by comparing the network’s performance in working memory tasks whose targets are drawn from distributions with different intertrial correlation profiles (Fig. 6). As a control, we consider the case of uniform correlations between target $\theta_{n}$ and target $\theta_{n+1}$ (Fig. 6A). Responses are normally distributed about the true target angle. The dynamics of the bump encoding the target are shaped by both input fluctuations and a bias in the direction of the previous target on individual trials. However, the directional bias is not apparent in the entire distribution of response angles, since it samples from all possible pairs $(\theta_{n+1},\theta_{n})$. An ensemble-wide measure of performance is given by the standard deviation of the response distribution ($\sigma\approx 4.42$). When target angles are correlated between trials, the relative response distribution narrows (Fig. 6B). Memory of the previous trial’s target $\theta_{n}$ stabilizes the memory of the current trial’s target $\theta_{n+1}$, decreasing the standard deviation of responses ($\sigma\approx 3.20$). There is a high probability the current target $\theta_{n+1}$ will be close to the previous target $\theta_{n}$, so the timescale of the network’s underlying inference process is reasonably well matched to the environment. However, such effects can be deleterious when the previous angle $\theta_{n}$ is skewed in comparison to the current angle $\theta_{n+1}$. Protocols with this correlation structure lead to a systematic bias in the relative response distribution, so its peak is shifted away from zero (Fig. 6C).

Our neuronal network model predicts that, if an intertrial bias is present in the computations of a neural circuit, it should be detectable by varying the intertrial correlation structure of target angles $\theta_{n}$. Furthermore, when there are strong local correlations between adjacent trials (${\rm corr}(\theta_{n+1},\theta_{n})$ is large for $|\theta_{n+1}-\theta_{n}|$ small), responses are more accurate than for protocols with uncorrelated adjacent trial angles. Since the strength of the bias increases as the intertrial interval is decreased, due to the decay of STF, we expect the effect to be more pronounced for trials taken closer together.

Two timescales of memory degradation

Wimmer et al. (2014) have shown that both the normal distribution of saccade endpoints and observed changes in neural firing rates during the delay-period can be replicated by a diffusing bump attractor model [4]. We have shown that a recurrent network with STF (Fig. 1E) still leads to a normal distribution of predicted response angles (Fig. 6A). Our model also provides new predictions for the dynamics of memory degradation, which we now compare with the standard diffusing bump attractor model [6] (Fig. 7). In a network with STF (Fig. 7A), bump trajectories evolve in a history-dependent fashion (Fig. 7B). Initially, bumps diffuse freely, but are eventually drawn to their starting location by facilitated synapses (See also Fig. 4). This leads to two distinct phases of diffusion, as shown in plots of the bump variance (Fig. 7C). Rapid diffusion occurs initially, as the bump equilibrates to the quasistationary density determined by the slowly evolving potential (Fig. 4). Slower diffusion occurs subsequently, as spatial heterogeneity in synaptic architecture gradually responds to changes in bump position via STF. Stabilizing effects of STF on bump attractors have been analyzed previously [19], but our identification of these multiple timescale dynamics is novel. This feature of the bump dynamics is not present in networks with static synapses (Fig. 7D). Here, bumps evolve as a noise-driven particle over a flat potential landscape (Fig. 7E), described by Brownian motion: a memoryless stochastic process [33]. Variance in bump position scales purely linearly with time (Fig. 7F), and the diffusion coefficient can be computed using a low-dimensional approximation [34].

The qualitative differences between the bump attractor with and without dynamic synapses should be detectable in both behavioral and neurophysiological recordings [4]. Moreover, the observed intertrial bias identified in recent analyses of behavioral data requires some mechanism for transferring information between trials that is distinct from neural activity [9], as dynamic synapses are in our model. In total, our model provides both an intuition for the behavioral motivation as well as neurophysiological mechanisms that produce such interference.

Discussion

Neural circuit models of visual working memory tend to use neural activity variables as the encoders of target locations. Our computational models account for interference in visual working memory using both suboptimal Bayesian inference and STF acting on a recurrent network model of delay-period activity. The timescale and prior target dependence of attractive biases we observe correspond to psychophysical observations of behavioral experiments in monkeys [9]. STF evolves dynamically over seconds [36], matching the kinetics of interference in visual working memory. The link we have drawn between our two models suggests neural circuits can implement probabilistic inference using short-term plasticity.

Experimental predictions

More complete descriptions of the neural mechanics of visual working memory could be accomplished by analyzing the effects of correlations in sequential target presentations. Since responses in subsequent trials are shaped by the previous trial’s target [9], computational models can be validated by determining how well their response distributions reflect trial-to-trial target correlations (Fig. 6). It is also possible that introducing target sequences whose distributions change in time could impact quantitative features of interference. For instance, implementing tasks with target sequences that have multiple trial correlations may extend the timescale of interference beyond a single trial. Furthermore, our model predicts that multiple timescales emerge in the statistics of delay-period activity during a working memory task (Fig. 7). Variance of recall error increases sublinearly in our model, consistent with a recent reanalysis of psychophysical data of saccades to remembered visual targets [3, 38]. The dynamics of our model are thus inconsistent with the purely linear diffusion of recall error common in bump attractor models with static synapses [6, 4].

The idea that STF may play a role in working memory is not new [39, 18], and there is evidence that prefrontal cortex neurons exhibit dynamic patterns of activity during the delay-period, suggestive of an underlying modulatory process [40]. However, it remains unclear how the presence of STF may shape the encoding of working memories. Our model suggests STF can transfer attractive biases between trials. Recent findings on the biophysics of STF could be harnessed to examine how blocking STF shapes behavioral biases in monkey experiments [41]. We predict that reducing the effects of the STF would both decrease the systematic bias in responses and increase the amplitude of errors, since the stabilizing effect of STF on the persistent activity will be diminished [19].

Alternative neurophysiological mechanisms for intertrial bias

Our study accounts for biases observed by Papadimitriou et al. (2015) [9], who identified an attraction between the previous target and current response. Strengthening synapses that originate from recently active neurons can attract neural activity states in subsequent trials. This is consistent with recent experiments showing latent and “activity-silent” working memories can be reactivated in humans using transcranial magnetic stimulation [42], suggesting working memory is maintained by mechanisms other than target-tuned persistent neural activity [18, 40]. The principle of using short-term plasticity to store memories of visual working memory targets could be extended to account for longer timescales and more intricate statistical structures. Short-term depression (STD) could effect a repulsive bias on subsequent responses, since neural activity would be less likely to persist in recently-activated depressed regions of the network. In this way, STD could encode a predictive distribution for targets that are anti-correlated to the previously present target.

Other physiological mechanisms could also shape network responses to encode past observations in a predictive distribution. Long-term plasticity is a more viable candidate for encoding predictive distributions that accumulate observations over long timescales. Consider a protocol that uses the same distribution of target angles throughout an entire experiment, but this distribution is biased towards a discrete set of possible angles [34]. For a recurrent network to represent this distribution, it must retain information about past target presentations over a long timescale. Many biophysical processes underlying plasticity are slow enough to encode information from such lengthy sequences [43]. Furthermore, the distributed nature of working memory suggests that there may be brain regions whose task-relevant neural activity partially persists from one trial to the next [8]. Such activity could shape low-level sensory interpretations of targets in subsequent trials.

Synaptic plasticity can stabilize working memory

Several modeling studies of working memory have focused on the computational capability of synaptic dynamics [44]. For instance, STF can prolong the lifetime of working memories in spatially heterogeneous networks, since facilitated synapses slow the systematic drift of bump attractor states [19]. This is related to our finding that STF reduces the diffusion of bumps in response to dynamic fluctuations (Fig. 7B), generating two timescales of memory degradation, corresponding to the bump variance (Fig. 7C). This scaling may be detectable in neural recordings or behavioral data, since recall errors may saturate if stabilized by STF. Facilitation can also account for experimentally observed increases in spike train irregularity during the delay-period in circuit models that support tuned persistent activity [45]. Alternatively, homeostatic synaptic scaling can compensate for spatial heterogeneity, which would otherwise cause persistent states to drift [7]. However, the short homeostatic timescales often suggested in models do not often match experimental observations [46].

Models of working memory have also replaced persistent neural firing with stimulus-selective STF, so neuronal spiking is only required for recall at the end of the delay-period [18]. One advantage of this model is that multiple items can be stored in the dynamic efficacy of synapses, and the item capacity can be regulated by external excitation for different task load demands [20]. Our model proposes that STF plays a supporting rather than a primary role, and there is extensive neurophysiological evidence corroborating persistent neural activity as a primary working memory encoder [4].

Robust working memory via excitatory/inhibitory balance

Computational modeling studies have demonstrated that a balance of fast inhibition and slow excitation can stabilize networks, so they accurately integrate inputs [32]. Drift in the representation of a continuous parameter can be reduced by incorporating negative-derivative feedback into the firing rate dynamics of a network, similar to introducing strong friction into the mechanics of particle motion on a sloped landscape [47]. Fast inhibition balanced by slower excitation produces negative feedback that is proportional to the time-derivative of population activity. A related mechanism can be implemented in spiking networks wherein fast inhibition rapidly prevents runaway excitation, and the resulting network still elicits highly irregular activity characteristic of cortical population discharges [48]. Mutually inhibiting balanced networks are similarly capable of representing working memory of continuous parameters [49], and extending our framework by incorporating STF into this paradigm would be a fruitful direction of future study.

Extensions to multi-item working memory

Working memory can store multiple items at once, and the neural mechanisms of interference between simultaneously stored items are the focus of ongoing work [50]. While there is consensus that working memory is a limited resource allocated across stored items, controversy remains over whether resource allocation is quantized (e.g., slots) [51] or continuous (e.g., fluid) [50]. Spatially-organized neural circuit models can recapitulate inter-item biases observed in multi-item working memory experiments, and provide a theory for how network interactions produce such errors [52]. In these models, each remembered item corresponds to an activity bump, and the spatial scale of lateral inhibition determines the relationship between recall error and item number [53]. The model provides a theory for attractive bias and forgetting of items since nearby activity bumps merge with one another. This is related to the mechanism of attractive bias in our model, but a significant difference is that previous models only required localized excitation whereas we use STF. It would be interesting identify the temporal dynamics of biases in multi-item working memory to see if they require slower timescale processes like short-term plasticity.

Tuning short-term plasticity to the environmental timescale

We have not identified a mechanism whereby our network model’s timescale of inference could be tuned to learn the inherent timescale of the environment. There is recent evidence from decision-making experiments that humans can learn the timescale on which their environment changes and use this information to weight their observations toward a decision [14]. Our model suggests that the trial-history inference utilized by subjects in Papadimitriou et al. (2015) is significantly suboptimal [9], perhaps because it is difficult to infer the timescale of relevant past-trial information. The complexity, sensitivity, and resource expense of optimal inference in most contexts likely makes it impossible to implement exactly in neural circuits [54]. This may explain why humans often use suboptimal methods for accumulating evidence [16]. Plasticity processes that determine the timescale of evidence accumulation may be shaped across generations by evolution, or across a lifetime of development. Nonetheless, metaplasticity processes can internally tune the dynamics of plasticity responses in networks without changing synaptic efficacy itself, and these changes could occur in a reward-dependent way [55]. Recently, a model of reward-based metaplasticity was proposed to account for adaptive learning observed in a probabilistic reversal learning task [56]. Such a process could modify the timescale and other features of short-term plasticity in ways that improve task performance in working memory as well.

Conclusions

Our results suggest that interference observed in visual working memory tasks can be accounted for by a persistently active neural circuit with STF. Importantly, interference is graded by the time between trials and during a trial. The interplay of synaptic and neural processes involved in interference may have arisen as a robust system for processing visual information that changes on the timescale of seconds. More work is need to determine how information about the environment stretches across multiple timescales to shape responses in cognitive tasks. We expect that identifying the neural origin of such biases will improve our understanding of how working memory fits into the brain’s information-processing hierarchy.

Online Methods

Assumptions of the inference model

Our model performs nonparametric density estimation to approximate the distribution $s_{n+1}(\theta)$ from which a target $\theta$ will be drawn before trial $n+1$. The observer assumes the possible distributions $s(\theta)$ are drawn from a function space $s\in S$ according to the prior $p(s)$. We assume that marginalizing over all such distributions yields the uniform density $\bar{{\rm P}}_{0}=\int_{S}s(\theta)p(s){\rm d}s=1/360$. One possibility is that the distribution $s_{n+1}(\theta)$ is constructed by drawing $N$-tuples ${\bf a}$ and $\boldsymbol{\psi}$ from a uniform distribution over the hypercubes $[0,a_{max}]^{N}$ and $[-180^{\circ},180^{\circ})^{N}$ and using the entries to construct an exponential distribution of a sum of cosines:

[TABLE]

where $\omega_{j}=j\pi/180$ and ${\mathcal{N}}_{s}$ is a normalization constant. For instance, when $N=1$,

[TABLE]

peaked at $\psi_{1}$. For the main instantiation and reduction of our model, knowing the specific family of distributions is unnecessary.

The probability $f_{\theta^{\prime}}(\theta):={\rm P}(\theta_{n+1}=\theta|\theta_{n}=\theta^{\prime},s_{n+1}(\theta)\equiv s_{n}(\theta))$ is defined under static conditions ($s_{n+1}(\theta)\equiv s_{n}(\theta)$) to separate the dynamic effects of sampling distribution $s_{n}(\theta)$ changes. We are performing nonparametric Bayesian estimation of the distribution, and the probability $f_{\theta^{\prime}}(\theta)$ is already marginalized over the space of distributions $s(\theta)$. Thus, we do not model the intermediate step of inferring the probability of each distribution $s(\theta)$ and marginalizing, but it could be computed by integrating over the prior on the function space, $f_{\theta^{\prime}}(\theta)=\int_{S}s(\theta)f(s|\theta^{\prime}){\rm d}s$. Each observation $\theta^{\prime}$ would give the probability $f(s|\theta^{\prime})$ that the current distribution is $s(\theta)$. Integrating over the space of all distributions $s\in S$ provides the probability the next target will be $\theta$, based on the previous observation $\theta^{\prime}$ alone and the assumption that the distribution remains the same from trial $n$ to $n+1$. Details on the difference between parametric and nonparametric Bayesian estimation of densities can be found in Orbanz and Teh (2011) [57]. Note, we assume self-conjugacy of $f_{\theta^{\prime}}(\theta)=f_{\theta}(\theta^{\prime})$, which follows since the order of observations does not matter while the environment remains fixed. This relationship will also make the predictiveness of our model more apparent. It is important to note that the observer assumes the form of $f_{\theta^{\prime}}(\theta)$, but this is not necessarily the distribution an ideal observer should use. For illustration, we consider a family of distributions given by an exponential of cosines:

[TABLE]

which is self-conjugate: $f_{\theta^{\prime}}(\theta)\equiv f_{\theta}(\theta^{\prime})$ [58]. A distribution like Eq. (7) would emerge from a generative model with distance-dependent spatial correlations in the ensemble of produced targets. The example $f_{\theta^{\prime}}(\theta)$ we use for comparison with our recurrent network with STF is close to the case of Eq. (7) with $N=1$. A description of the parameters and variables in our model is provided in Table 1.

Derivation of the probabilistic inference model

The observer’s predictive distribution $L_{n+1,\theta}={\rm P}(\theta_{n+1}|\theta_{1:n},\epsilon)$ is derived by computing the probability of observing $\theta_{n+1}$ given each prior observation $\theta_{j}$, $j=1,...,n$. Importantly, we must compute the probability of each run length $r_{n}=r$, $r=0,...,n$, corresponding to the number of trials the assumed underlying distribution $s_{n}(\theta)$ has remained the same [24, 21]. Knowing the probability of each run length will inform us of how much to weight each observation $\theta_{j}$, $j=1,...,n$. In particular, $r_{n}=n$ indicates the environment has remained the same since the first trial, and $r_{n}=0$ indicates the environment changes between trial $n$ and $n+1$. Summing over all possible run lengths, the marginal predictive distribution is

[TABLE]

where ${\rm P}(\theta_{n+1}|r_{n}=r,\theta_{1:n}^{r})$ is the conditional predictive distribution assuming run length $r_{n}=r$ and ${\rm P}(r_{n}=r|\theta_{1:n})$ is the conditional probability of the run length $r_{n}=r$ given the series of target angles $\theta_{1:n}$. We further simplify Eq. (8) as follows: First, utilizing sequential analysis, we find that if the present run length is $r_{n}=r$, the conditional predictive distribution is given by the product of probabilties from the last $r$ observations [15]:

[TABLE]

where again $\bar{{\rm P}}_{0}$ is the uniform distribution and we have utilized our self-conjugacy assumption for $f_{\theta^{\prime}}(\theta)\equiv f_{\theta}(\theta^{\prime})$. Next, we assume that observations provide no information about the present run length $r_{n}$, which would be a consequence of the observer making no a priori assumptions on the overall distribution from which targets $\theta_{1:n}$ are drawn. Thus, the observer only uses their knowledge of the change rate of the environment $\epsilon$ to determine the probability of a given run length $r_{n}=r$, and the conditional probability can be computed

[TABLE]

Plugging Eqs. (9–12) into the update Eq. (8), we find the probability of the next target being at angle $\theta_{n+1}=\theta$, given that the previous $n$ targets were $\theta_{1:n}$, is:

[TABLE]

Limit of slowly-changing environment (small $\epsilon$)

Here, we examine the case $0<\epsilon\ll 1$, where the environment changes very slowly. Assuming independence of the target angles selected on each trial $\theta_{1:n}$ [27], ${\rm P}(\theta_{n-r:n})={\rm P}(\theta_{n-r:n-1}){\rm P}(\theta_{n})$, we can split the probabilities over the target sequences $\theta_{n-r:n}$ into products: ${\rm P}(\theta_{n-r:n})=\prod_{j=n-r}^{n}{\rm P}(\theta_{j})=\bar{{\rm P}}_{0}^{r+1}$. The last equality holds since the family of possible distributions $s_{n}(\theta)$ averages to a constant $\bar{{\rm P}}_{0}$, the uniform density. Applying this assumption to Eq. (1) and truncating to ${\mathcal{O}}(\epsilon)$, we have

[TABLE]

noting we must choose ${\mathcal{N}}_{s}$ so $\int_{-180}^{180}L_{n+1,\theta}{\rm d}\theta=1$, normalized at each step.

Limit of rapidly-changing environment ($\epsilon\approx 1$)

Here, we examine the case $\epsilon\approx 1$ ($0<(1-\epsilon)\ll 1$), a rapidly-changing environment. Applying this assumption to Eq. (1), we find $L_{n+1,\theta}$ is dominated by terms of order $(1-\epsilon)$ and larger. Terms of order $(1-\epsilon)^{2}$ are much smaller. For instance, we can approximate to linear order, dropping terms of ${\mathcal{O}}((1-\epsilon)^{2})$, to reduce Eq. (1) to

[TABLE]

Furthermore, we ensure the expression in Eq. (13) is normalized by writing

[TABLE]

since $\int_{-180}^{180}\left[\bar{{\rm P}}_{0}+(1-\epsilon)f_{\theta_{n}}(\theta)\right]{\rm d}\theta=2-\epsilon$. Alternatively, we can truncate by multiplying through by $[1-(1-\epsilon)]/[1-(1-\epsilon)]$, truncating to ${\mathcal{O}}(1-\epsilon)$ and normalizing to yield

[TABLE]

the key update equation in our Results (Figs. 2 and 3A). Higher order approximations are obtained by keeping more terms from Eq. (1); e.g., a second order approximation yields

[TABLE]

successively downweighting the influence of previous observations ($\theta_{n-1}$).

Relating the predictive distribution to the potential of an attractor model

A predictive distribution can be represented by an attractor model by first determining the formula of the stationary distribution of Eq. (4), given an arbitrary potential function ${\mathcal{U}}_{n+1}(\theta)$. Eq. (4) can be reformulated as an equivalent Fokker-Planck equation for the represented angle $\theta$ during trial $n+1$ assuming the present potential function is ${\mathcal{U}}_{n+1}(\theta)$ [59],

[TABLE]

where $p_{n+1}(\theta,t)$ is the probability density corresponding to the target angle estimate $\theta$ at time $t$. The initial estimate of the target is exact, $\theta(0)=\theta_{n+1}$, so $p_{n+1}(\theta,0)=\delta(\theta-\theta_{n+1})$ is the initial condition. We summarize the constituent variables and model parameters in Table 2.

We now derive the form of ${\mathcal{U}}_{n+1}(\theta)$ that leads to a stationary density corresponding the predictive distribution $L_{n+1,\theta}$ in the limit $t\to\infty$ in Eq. (14). The stationary density $\bar{p}_{n+1}(\theta)$ is analogous to a predictive distribution represented by Eq. (4) since it is the probability the system represents when no information about the current trial’s target $\theta_{n+1}$ remains. Thus, we build a rule to update ${\mathcal{U}}_{n+1}(\theta)$ to mirror the update of $L_{n+1,\theta}$ in Eq. (3). To obtain this result, we match the stationary density for Eq. (14) to the updated predictive distribution:

[TABLE]

Solving Eq. (14) for its stationary solution, we find that during trial $n+1$:

[TABLE]

where $\chi_{n+1}$ is a normalization factor chosen so that $\int_{-180}^{180}\bar{p}_{n+1}(\theta){\rm d}\theta=1$. Plugging Eq. (16) into Eq. (15) and solving for ${\mathcal{U}}_{n+1}(\theta)$, we obtain

[TABLE]

For a rapidly changing environment $0<(1-\epsilon)\ll 1$, we approximate $L_{n+1,\theta}$ using Eq. (3) so that

[TABLE]

where we have linearized in $(1-\epsilon)$. However, for Eq. (4), only the derivative of ${\mathcal{U}}_{n+1}(\theta)$ impacts the dynamics, so we drop the additive constants and examine the proportionality

[TABLE]

In the limit of weak interactions between trials, the potential ${\mathcal{U}}_{n+1}(\theta)$ should be shaped like the negative of the probability $f_{\theta_{n}}(\theta)$ based on the previous trial’s target $\theta_{n}$.

Bump attractor model with short-term facilitation

Our neuronal network model is comprised of two variables evolving in space $x\in[-180,180)^{\circ}$, corresponding to the stimulus preference of neurons at that location, and time $t>0$. Variables and parameters are summarized in Table 3, and the evolution equations consist of one stochastic integrodifferential equation and one auxiliary differential equation:

[TABLE]

where $u(x,t)$ describes the evolution of the normalized synaptic input at location $x$. The model Eq. (17) can be derived as the large system size limit of a population of synaptically coupled spiking neurons [61], and similar dynamics have been validated in spiking networks with lateral inhibitory connectivity  [6, 4]. We fix the timescale of dynamics by setting $\tau_{u}=10$ms, so time evolves according to units of a typical excitatory synaptic time constant [60]. This population rate model can be explicitly analyzed to link the architecture of the network to a low-dimensional description of the dynamics of a bump attractor as described by Eq. (4).

Each location $x$ in the network receives recurrent coupling defined by the weight function $w(x-y)$. We take this function to be peaked when $x=y$ and decreasing as the distance $|x-y|$ grows, in line with anatomical studies of delay-period neurons in prefrontal cortex [5]. We do not separately model excitatory and inhibitory populations, but Eq. (17) can be derived from a model with distinct excitatory and inhibitory populations in the limit of fast inhibitory synapses [35, 62]. Thus, we have combined excitatory and inhibitory populations, so $w(x-y)$ takes on both positive and negative values. Our analysis can be applied to a general class of distance-dependent connectivity functions, given by an arbitrary sum of cosines $w(x-y)=\sum_{n=0}^{\infty}\alpha_{n}\cos(\omega_{n}(x-y))$ where $\omega_{n}=n\pi/180$, and we will use a single cosine to illustrate in examples: $w(x-y)=\cos(\omega_{1}(x-y))$. The nonlinearity $F(u)$ converts the normalized synaptic input $u(x,t)$ into a normalized firing rate, $F(u)\in[0,1]$. We take this to be sigmoidal $F(u)=1/\left[1+{\rm e}^{-\gamma(u-\kappa)}\right]$ [63], with a gain of $\gamma=20$ and a threshold of $\kappa=0.1$ in numerical simulations. In the high-gain limit ($\gamma\to\infty$), a Heaviside step function $F(u)=H(u-\kappa)$ allows for explicit calculations [35, 61].

Recurrent coupling is shaped by STF in active regions of the network ($F(u)>0$), as described by the variable $q(x,t)\in[0,q_{+}]$; $q_{+}>0$ and $\beta$ determine the maximal increase in synaptic utilization and the rate at which facilitation occurs [17, 64]. For our numerical simulations, we consider the parameter values $q_{+}=2$ and $\beta=0.01$, consistent with previous models employing facilitation in working memory circuits [19, 18, 20] and experimental findings for facilitation responses in prefrontal cortex [36]. The timescale of plasticity is slow, $\tau=1000{\rm ms}\gg 10$ms, consistent with experimental measurements [17]. Our qualitative results are robust to parameter changes. Information from the previous trial is maintained by the slow-decaying kinetics of the facilitation variable $q(x,t)$, even in the absence of neural activity [18, 20].

Effects of the target and the response are described by the deterministic spatiotemporal input $I(x,t)$, which we discuss more in detail below. The noise process $W(x,t)$ is white in time and has an increment with mean $\langle{\rm d}W(x,t)\rangle\equiv 0$ and spatial correlation function $\langle{\rm d}W(x,t){\rm d}W(y,s)\rangle=C(x-y)\delta(t-s){\rm d}t{\rm d}s$. In numerical simulations, we take our correlation function to be $C(x-y)=\sigma_{W}^{2}\cos(x-y)$ with $\sigma_{W}=0.005$, so the model recapitulates the typical 1-5% standard deviation in saccade endpoints observed in oculomotor delayed-response tasks with delay-periods from 1-10s [1, 3, 4].

Implementing sequential delayed-response task protocol

A series of oculomotor delayed-response tasks is executed by the network Eq. (17) by specifying a schedule of peaked inputs occurring during the cue periods of length $T_{C}$, no input during trial $n$’s delay-period of length $T_{D}^{n}$, and brief and strong inhibitory input of length $T_{A}$ after the response has been recorded, and then no input until the next trial. This is described by the spatiotemporal function

[TABLE]

for all $n=1,2,3,...$, where $t_{n}$ is the starting time of the $n^{\text{th}}$ trial which has cue period $T_{C}$, delay-period $T_{D}^{n}$, inactivation period $T_{A}$, and subsequent intertrial interval $T_{I}^{n}$. Note that the delay and intertrial interval times may vary trial-to-trial, but the cue is always presented for the same period of time as in Papadimitriou et al. (2015) [9]. The amplitude of the cue-related stimulus is controlled by $I_{0}$, and $I_{1}$ controls is sharpness. Activity from trial $n$ is ceased by the global inactivating stimulus of amplitude $I_{R}$.

In numerical simulations, we fix the parameters $T_{C}=500$ms; $T_{A}=500$ms; $I_{0}=1$; $I_{1}=1$; and $I_{R}=2$. Target locations $\theta_{n}$ are drawn from a uniform probability mass function (pmf) for the discrete set of angles $\theta_{n}\in\{-180^{\circ},-162^{\circ},...,162^{\circ}\}$ to generate statistics in Fig. 5A, which adequately resolves the bias effect curves for comparison with the results in Papadimitriou et al. (2015) [9]. Intertrial intervals are varied to produce Fig. 5B by drawing $T_{I}^{n}:=t_{n+1}-(T_{C}+T_{D}^{n}+T_{A})$ randomly from a uniform pmf for the discrete set of times $T_{I}^{n}\in\{1000,1200,...,5000\}$ms and $\theta_{n}$ randomly as in Fig. 5A and identifying the $\theta_{n}$ that produces the maximal bias for each value of $T_{I}^{n}$. Delay-periods are varied to produce Fig. 5C by drawing $T_{D}^{n}$ randomly from a uniform pmf for the discrete set of times $T_{I}^{n}\in\{0,200,...,5000\}$ms and following a similar procedure to Fig. 5B. Draws from a uniform density function ${\rm P}(\theta_{n})\equiv\bar{{\rm P}}_{0}$, defined on $\theta_{n}\in[-180,180)^{\circ}$ are used to generate the distribution in Fig. 6A and plots in Fig. 7. Nontrivial correlation structure in target selection is defined by the sum of a von Mises distribution and uniform distribution ${\rm corr}(\theta_{n+1},\theta_{n})=(1-\epsilon){\mathcal{N}}_{v}{\rm e}^{25\cos(\theta_{n}-\theta_{n+1}-\mu)}+\epsilon\bar{{\rm P}}_{0}$ for fixed $\theta_{n}$ with $\epsilon=0.5$; $\mu=0$ for local correlations (Fig. 6B) and $\mu=90$ for skewed correlations (Fig. 6C).

The recurrent network, Eq. (17), is assumed to encode the initial target $\theta_{n}$ during trial $n$ via the center-of-mass $\theta(t)$ of the corresponding bump attractor. Representation of the cue at the end of the trial is determined by performing a readout on the neural activity $u(x,t)$ at the end of the delay time for trial $n$: $t=t_{n}+T_{C}+T_{D}^{n}$. One way of doing this would be to compute a circular mean over $x$ weighted by $u(x,t)$, but since $u(x,t)$ is a roughly symmetric and peaked function in $x$, computing $\theta(t):={\rm argmax}_{x}u(x,t)$ (when $t\in[t_{n},t_{n}+T_{C}+T_{D}^{n})$) is an accurate and efficient approximation [34, 4]. The bias and relative saccade endpoint on each trial $n$ are then determined by computing the difference $\theta(t)-\theta_{n}$ (Figs. 5, 6, and 7).

Deriving the low-dimensional description of bump motion

We analyze the mechanisms by which STF shapes the bias on subsequent trials by deriving a low-dimensional description for the motion of the bump position $\theta(t)$. To begin, note that in the absence of facilitation ($\beta\equiv 0$), the variable $q(x,t)\equiv 0$. In the absence of noise ($W(x,t)\equiv 0$), the resulting deterministic Eq. (17) has stationary bump solutions that are well studied and defined by the implicit equation [35, 61, 65]:

[TABLE]

Assuming the stimulus $I(x,t)$ presented during the cue period of trial $n$ ($t\in[t_{n},t_{n}+T_{C})$) is strong enough to form a stationary bump solution, the impact of the facilitation variable $q(x,t)$ and noise $W(x,t)$ on $u(x,t)$ during the delay-period ($t\in[t_{n}+T_{C},t_{n}+T_{C}+T_{D}^{n})$) can be determined perturbatively, assuming $|q|\ll 1$ and $|{\rm d}W|\ll 1$. Since $\tau\gg\tau_{u}$, $u(x,t)$ will rapidly equilibrate to a quasi-steady-state determined by the profile of $q(x,t)$. We thus approximate the neural activity dynamics as $u(x,t)\approx U(x-\theta(t))+\Phi(x,t)$, where $\theta(t)$ describes the dynamics of the bump center-of-mass during the delay-period ($|\theta|\ll 1$ and $|{\rm d}\theta|\ll 1$), and $\Phi(x,t)$ describes perturbations to the bump’s shape ($|\Phi|\ll 1$). Plugging this approximation into Eq. (17) and truncating to linear order yields

[TABLE]

where ${\mathcal{L}}u=-u+\int_{-180}^{180}w(x-y)F^{\prime}(U(y))u(y){\rm d}y$ is a linear operator and $q(x,t_{s})$ is the facilitation variable evolving on the slow timescale $t_{s}=\tau_{u}t/\tau\ll t$, quasi-stationary on the fast timescale of $u(x,t)$. We ensure a bounded solution by requiring the right hand side of Eq. (18) is orthogonal to the nullspace $V(x)$ of the adjoint linear operator ${\mathcal{L}}^{*}v=-v+F^{\prime}(U)\int_{-180}^{180}w(x-y)v(y){\rm d}y$. Orthogonality is enforced by requiring the inner product $\langle u,v\rangle=\int_{-180}^{180}u(x)v(x){\rm d}x$ of the nullspace $V(x)$ with the inhomogeneous portion of Eq. (18) is zero. It can be shown $V(x)=F^{\prime}(U(x))U^{\prime}(x)$ spans the nullspace of ${\mathcal{L}}^{*}$ [65]. This yields the following equation for the evolution of the bump position:

[TABLE]

where the slowly evolving nonlinearity

[TABLE]

is shaped by the form of $q(x,t_{s})$ and the noise $\xi(t)$ is a standard Wiener process that comes from filtering the full spatiotemporal noise process ${\rm d}W(x,t)$, so the diffusion coefficient

[TABLE]

Eq. (19) has the same form as Eq. (4). Thus, if the facilitation variable $q(x,t_{s})$ evolves trial-to-trial such that $K(\theta,t_{s})$ has similar shape to $\displaystyle-\frac{{\rm d}{\mathcal{U}}_{n+1}}{{\rm d}\theta}(\theta)$ at the beginning of the $(n+1)^{\text{th}}$ trial ($t=t_{n+1}$), the dynamics of the network Eq. (17) can reflect a prior distribution based on the previous target(s). Given the approximation we derived in Eq. (5), we enforce proportionality $K(\theta,t_{n+1})\propto-\frac{{\rm d}{\mathcal{U}}_{n+1}}{{\rm d}\theta}(\theta)$:

[TABLE]

where $\alpha$ is a scaling constant and $t_{n+1}$ is the starting time of trial $n+1$ in the original time units $t=\tau t_{s}/\tau_{u}$. The form of the probability $f_{\theta^{\prime}}(\theta)$ that can be represented is therefore restricted by the dynamics of the facilitation variable $q(x,t)$. We can perform a direct calculation to identify how $q(x,t)$ relates to the predictive distribution it represents in the following special case.

Explicit solutions for high-gain firing rate nonlinearities

To explicitly calculate solutions, we take the limit of high-gain, so that $F(u)\to H(u-\kappa)$ and $w(x)=\cos(\omega_{1}x)$, note $\omega_{1}=180/\pi$. In this case, the bump solution $U(x-x_{0})=(2\sin(a)/\omega_{1})\cos(\omega_{1}(x-x_{0}))$ for $U(\pm a)=\kappa$ and null vector $V(x-x_{0})=\delta(x-x_{0}-a)-\delta(x-x_{0}+a)$ (without loss of generality we take $x_{0}\equiv 0$) [65]. Furthermore, we can determine the form of the evolution of $q(x,t)$ by studying the stationary solutions to Eq. (17) in the absence of noise ($W\equiv 0$). For a bump $U(x)$ centered at $x_{0}=0$, the associated stationary form for $Q(x)$ assuming $H(U(x)-\kappa)=1$ for $x\in(-a,a)$ and zero otherwise is $Q(x)=\beta q_{+}/(1+\beta)$ for $x\in(-a,a)$ and zero otherwise. Thus, if the previous target was at $\theta_{n}$, we expect $q(x,t)$ to have a shape resembling $Q(x-\theta_{n})$ after trial $n$. Assuming the cue plus delay time during trial $n$ was $T_{C}+T_{D}^{n}$ and the intertrial interval is $T_{I}^{n}$, slow dynamics will reshape the amplitude of $q(x,t)$ so ${\mathcal{A}}_{n}(T^{n})=(1-{\rm e}^{-(T_{C}+T_{D}^{n})/\tau}){\rm e}^{-T_{I}^{n}/\tau}$ ($T^{n}=T_{C}+T_{D}^{n}+T_{I}^{n}$ is the total time block of each trial) and so $q(x,t)\approx{\mathcal{A}}_{n}(T^{n})\cdot Q(x-\theta_{n})$ at the beginning of trial $n+1$. A lengthy calculation of Eq. (20) combined with the relation Eq. (21) yields:

[TABLE]

for $|\theta-\theta_{n}|<2a$, and $\frac{{\rm d}f_{\theta_{n}}(\theta)}{{\rm d}\theta}\equiv 0$ otherwise. Integrating, we find this implies

[TABLE]

for $|\theta-\theta_{n}|<2a$, and $f_{\theta_{n}}(\theta)$ constant otherwise. Thus, the STF dynamics allows the network architecture to represent a predictive distribution that is peaked at the previous target location (Fig. 3). The amplitude of the $\theta$-dependent portion of the predictive distribution during trial $n+1$ is then controlled by cue, delay, and intertrial times ($T_{C},T_{D}^{n+1},T_{I}^{n+1}$) and the facilitation parameters ($\beta$, $q_{+}$, $\tau$).

To derive a coupled pair of equations (Fig. 4) describing the dynamics of the bump location $\theta(t)$ and the slow evolution of the nonlinearity $K(\theta,t)$, we focus on the limit $F(u)\to H(u-\kappa)$. We approximate $q(x,t)$ by summing the contributions from each of the $n+1$ trials. This yields

[TABLE]

where the slowly evolving function ${\mathcal{A}}_{n}(t)$ defines the rising and falling kinetics of the facilitation variable originating in trial $n$:

[TABLE]

increasing towards saturation (${\mathcal{A}}_{n}\to 1$) during the cue and delay-period $[t_{n},t_{n}+T_{C}+T_{D}^{n})$ and decaying afterward (${\mathcal{A}}_{n}\to 0$). The variable $\theta_{q}(t)$ describes the slow movement of the center-of-mass of the saturating portion of the facilitation variable $q(x,t)$ due to the drift of the neural activity $u(x,t)$ described by $\theta(t)$. However, since ${\mathcal{A}}_{1}(t)\ll{\mathcal{A}}_{2}(t)\ll\cdots\ll{\mathcal{A}}_{n}(t)$, we only keep the terms ${\mathcal{A}}_{n}(t)$ and ${\mathcal{A}}_{n+1}(t)$ in Eq. (22). Furthermore, since ${\mathcal{A}}_{n}(t)$ becomes much smaller than ${\mathcal{A}}_{n+1}(t)$ for most times $t>t_{n+1}$ in trial $n+1$, we approximate $\theta_{q}(t_{n}+T_{C}+T_{D}^{n})\approx\theta_{n}$. This provides intuition as to why it is sufficient to only consider the previous target rather than the response in trial $n$ as the variable influencing the bias in Papadimitriou et al. (2015) [9]. Therefore, we start with the following ansatz for the evolution of the facilitation variable during trial $n+1$:

[TABLE]

A bump centered at $\theta(t)$, $U(x-\theta(t))$, attracts the STF variable to the same location $q\to Q(x-\theta(t))$, but the dynamics of $q$ are much slower ($\tau\gg 1$). Thus, we model the evolution of $\theta_{q}(t)$ by linearizing the slow dynamics of Eq. (17b) about $(u,q)=(U(x-\theta(t)),Q(x-\theta(t)))+(0,\phi(x,t))$ (with $|\phi|\ll 1$) to find

[TABLE]

The perturbation $\phi(x,t)$ describes the displacement of the variable $q$ away from its equilibrium position. We thus introduce the field $\Phi(x,t)=\int_{-180}^{180}w(x-y)\phi(y,t)F(U(y-\theta(t))){\rm d}y$, which reduces Eq. (24) to

[TABLE]

so separating variables $\Phi(x,t)=\bar{\Phi}(x){\rm e}^{\lambda t}$ we see that perturbations of the facilitation variable’s center-of-mass $\theta_{q}(t)$ away from $\theta(t)$ should relax at rate $\lambda_{\tau}=-(1+\beta)/\tau$.

Therefore, the slow evolution of the potential gradient function $K(\theta,t_{s})$ in Eq. (19) can be described by integrating Eq. (20) using the ansatz Eq. (23) for $q(x,t)$. Our low-dimensional system for the dynamics of the bump location $\theta(t)$ and leading order facilitation bump $\theta_{q}(t)$ during the delay-period of trial $n+1$ ($t\in[t_{n+1}+T_{C},t_{n+1}+T_{C}+T_{D}^{n+1})$) is given by the set of non-autonomous stochastic differential equations:

[TABLE]

where we have defined a parametrized time-invariant potential gradient $\frac{{\rm d}\bar{\mathcal{U}}(\theta-\theta^{\prime})}{{\rm d}\theta}$ corresponding to the stationary profile of the facilitation variable centered at $\theta^{\prime}$: $Q(x-\theta_{n})$. For our specific choices of weight function and firing rate nonlinearity, we find the potential gradient is:

[TABLE]

and

[TABLE]

calculates the shorter difference on the periodic domain. As in our recurrent network, we use the parameters $\kappa=0.1$; $q_{+}=2$; $\beta=0.01$; and $\tau/\tau_{u}=100$ to compare with network simulations in Fig. 5.

Numerical simulations of the neuronal network model

Numerical simulations of the recurrent network Eq. (17) were done in MATLAB using an Euler-Maruyama method with timestep ${\rm d}t=0.1$ms and spatial step ${\rm d}x=0.18^{\circ}$ with initial conditions generated randomly by starting $u(x,0)\equiv q(x,0)\equiv 0$ and allowing the system to evolve in response to the dynamic fluctuations for $t=2$s prior to applying the sequence of stimuli $I(x,t)$ described for each numerical experiment in Figs. 5, 6, and 7. Numerical simulations of Eq. (6) were also performed using an Euler-Maruyama method with timestep ${\rm d}t=0.1$ms. The effects of the target $\theta_{n}$ on each trial $n$ were incorporated by holding $\theta(t)=\theta_{n}$ during the cue period $t\in[t_{n},t_{n}+T_{C})$. Otherwise, the dynamics were allowed to evolve as described.

Data Analysis

MATLAB was used for statistical analysis of all numerical simulations. The bias effects in Fig. 5 were determined by identifying the centroid of the bump at the end of the delay-period. Means were computed across $10^{5}$ simulations each, and standard deviations were determined by taking the square root of the var command applied to the vector of endpoints. Histograms in Fig. 6 were computed for $10^{5}$ simulations using the hist and bar commands applied to the vector of endpoints for each correlation condition. Bump positions were computed in Fig. 7 by determining the centroid of the bump at each timepoint, and $10^{5}$ simulations were then used to determine the standard deviation and variance plots (using var again).

Author Contributions

ZPK conceived the study, performed simulations and analysis, and wrote the paper.

Acknowledgments

This work was funded by National Science Foundation grants NSF-DMS-1615737 and NSF-DMS-1517629. We thank Krešimir Josić and Brent Doiron for helpful conversations and comments on the manuscript.