remarkRemark hypothesisHypothesis claimClaim Distinct excitatory/inhibitory
bump wanderingH.L. Cihak, T.L. Eissa, and Z.P. Kilpatrick

# Distinct excitatory and inhibitory bump wandering in a stochastic neural
field ††thanks: Submitted to the editors DATE. This work was funded by
National Institutes of Health grants R01MH115557-01 and R01-EB029847-01 as
well as an NSF grant DMS-1853630.

Heather L Cihak Department of Applied Mathematics, University of Colorado
Boulder, Boulder, CO<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Tahra L Eissa22footnotemark: 2 Zachary P
Kilpatrick22footnotemark: 2

###### Abstract

Localized persistent cortical neural activity is a validated neural substrate
of parametric working memory. Such activity ‘bumps’ represent the continuous
location of a cue over several seconds. Pyramidal (excitatory) and
interneuronal (inhibitory) subpopulations exhibit tuned bumps of activity,
linking neural dynamics to behavioral inaccuracies observed in memory recall.
However, many bump attractor models collapse these subpopulations into a
single joint excitatory/inhibitory (lateral inhibitory) population, and do not
consider the role of interpopulation neural architecture and noise
correlations. Both factors have a high potential to impinge upon the
stochastic dynamics of these bumps, ultimately shaping behavioral response
variance. In our study, we consider a neural field model with separate
excitatory/inhibitory (E/I) populations and leverage asymptotic analysis to
derive a nonlinear Langevin system describing E/I bump interactions. While the
E bump attracts the I bump, the I bump stabilizes but can also repel the E
bump, which can result in prolonged relaxation dynamics when both bumps are
perturbed. Furthermore, the structure of noise correlations within and between
subpopulations strongly shapes the variance in bump position. Surprisingly,
higher interpopulation correlations reduce variance.

###### keywords:

neural fields, bump attractor, wave stability, stochastic differential
equations, perturbation theory

68Q25, 68R10, 68U05

##  1 Introduction

Storing information ‘in mind’ for short periods of time is essential for
performance of daily tasks [24]. Parametric working memory (as used in tasks
requiring a delayed estimate of a continuous quantity) uses spatially
localized persistent neural activity in the prefrontal cortex, parietal
cortex, and frontal eye fields [16, 24] generated in neural circuits with
strong local recurrent excitation and broad inhibition [9, 31, 44].
Neurophysiological recordings from non-human primate subjects performing
visuospatial working memory tasks have shown that localized bumps of
persistent activity encode remembered parametric quantities for a few seconds
[9, 31, 21, 24]. For example, in the oculomotor delayed response task, a
location cue is momentarily presented on a circle displayed on a monitor, then
a delay period occurs during which the video is blank, and finally the subject
is prompted to report their memory of the cued location. During the delay
period, neural recordings identify cells tuned to specific cue locations
around the circle and reveal that the strongest (peak) firing neurons
represent the encoded location [24, 21, 22]. Peaked and localized activity
wanders stochastically during the delay period, representing a time-dependent
degradation of cue location memory consistent with subsequent response errors
[8, 48].

Neural field and spiking models are capable of producing peaked, localized
activity that wanders via spatially structured connectivity that weakens as
the distance between neural cue location preference increases with excitation
having a narrower spatial footprint [47, 27, 44, 31]. Bumps can be defined as
standing pulse solutions usually with two key dynamical properties important
for encoding memories of continuum variables: (1) they are self-sustained in
the absence of stimulus (representing memory over a delay period) and (2) they
exhibit marginal stability, such that they can occur at any location in the
network and integrate translational perturbations [27]. Such solutions have
garnered interest due to the rich dynamical features that arise when varying
the form of connectivity, introducing propagation delays, adding slow negative
feedback, or considering stochasticity [4, 27, 31, 12, 19]. Response
statistics and neural activity during oculomotor delayed response tasks are
also well characterized by the output of bump attractor models with some form
of stochasticity incorporated [48, 9, 3].

Previous psychophysical studies of delayed estimation of continuous quantities
show subject errors scale linearly with the delay period [43, 46]. Such
behavior can be accounted for by models whose low-dimensional dynamics evolve
as Brownian motion, like a particle subject to diffusion. This is consistent
with a bump attractor stochastically perturbed by noise, wandering along a
marginally stable ring attractor [31, 8]. Extended models have also considered
neural circuit mechanisms that help stabilize the movement of bumps to noise
perturbations by breaking the marginal stability of the ring attractor. For
example, spatially heterogeneous recurrent excitation leads to low-dimensional
dynamics akin to a particle stochastically perturbed along a washboard
potential, slowing the rate of diffusion [32]. Likewise, short-term
facilitation locally potentiates synaptic excitation where the bump is
instantiated, akin to a slowly moving local potential well that traps the
particle near the true memorandum location [44, 30]. In addition to
stabilizing bumps within trials, short-term facilitation also transfers memory
of the previous trial stimulus to the next, creating systematic errors
referred to as serial dependence [35, 7].

Despite advancements in understanding how these modifications to bump
attractor models affect their stochastic dynamics, a detailed examination of
the role of separate excitatory and inhibitory (E and I) population dynamics
has been overlooked, and I-I interactions are often ignored entirely.
Inhibition is often assumed to be flat and global in stochastic bump attractor
models [9], but we know prefrontal cortical synaptic inhibition exhibits
nontrivial preference-dependent interactions [24], which could have yet
unidentified effects on the dynamics of bumps and the fidelity of memory
systems that rely on them.

Here, we use a stochastic neural field model to investigate the wandering
dynamics and interactions of separate E and I activity bumps. We focus on how
a separate and spatially-extended I population shapes the variance of a bump
attractor encoding a remembered location along a continuum. First, we
introduce the neural field equations, notation conventions, and parameters
whose effects on the variance of the bumps’ final position will be analyzed
(Section 2). We then review the existence and stability of both E and I bump
solutions [47, 27, 4, 14, 37, 20]; paving the way for examining how
deterministic and stochastic perturbations deform and shift the attractor
solution. We examine the impact of modifying both the amplitude and spatial
extent of neural architecture, as well as the firing rate thresholds on the
stability of solutions (Section 3). Stable and unstable branches of bumps
annihilate at a discontinuous saddle-node bifurcation given a sufficiently
high firing rate thresholds.

Our stochastic analysis involves two different methods of identifying the
impact of noise on the bump attractor, (1) a linear perturbative analysis
(Strongly Coupled Limit between the E and I bumps) in the limit of weak noise
and (2) an Interface Based analysis where we track the threshold crossing
points of the E/I bumps (Section 4). Thus, we develop two corresponding
theoretical predictions of bump position variance, the second of which
captures nonlinear interactions between the E/I bumps in a Langevin system
that is a better match to full model simulations. This higher-order analysis
starts by obtaining distinct and nonlinearly coupled equations for the bump
interfaces. Bump position variance changes non-monotonically as a function of
most model parameters, so intermediate tuning of network architecture
maximizes memory degradation. Finally, we identify the impact of noise
correlations between the E and I populations; finding they actually serve to
reduce bump position variance.

##  2 The Model

To determine the effects of separately evolving excitatory (E) and inhibitory
(I) populations on bump attractor dynamics, we analyze a stochastic neural
field model: a coupled system of nonlinear integro-differential equations
describing interactions of spatially-extended E and I neural sub-populations
(Figure 1A and [47]). Recurrent connectivity targeting position $x$ from $y$
at time $t$ is described by convolving
($a(x)*b(x)=\int_{-\infty}^{\infty}a(x-y)b(y)dy$) the synaptic kernel and the
firing rate output of the neurons from which synapses originate. Along with
additive noise terms, we obtain the stochastic neural field equations:

(1a) $\displaystyle
du(x,t)=\left[-u(x,t)+{w_{ee}(x)*f(u(x,t))}-{w_{ei}(x)*f(v(x,t))}\right]dt+\epsilon^{\frac{1}{2}}dW_{e}$
(1b) $\displaystyle\tau
dv(x,t)=\left[-v(x,t)+{w_{ie}(x)*f(u(x,t))}-{w_{ii}(x)*f(v(x,t))}\right]dt+\epsilon^{\frac{1}{2}}dW_{i}.$

where $u(x,t)$ and $v(x,t)$ denote the E and I synaptic input profiles at
location $x$ at time $t$. We drop the time constant from the E population and
take time units to be 10ms, on the order of typical E membrane and synaptic
time constants [25]. The sigmoidal firing rate function
$f(u)=\frac{1}{1+e^{-\eta(u-\theta)}}$ with threshold $\theta$ and gain $\eta$
determines the strength of the firing rate output based on the input $u$.
Analytical results can be obtained in the high gain limit ($\eta\to\infty$)
such that

(2) $\displaystyle f(u)=H(u-\theta)=\begin{cases}1&u-\theta\geq 0,\\\
0&u-\theta<0,\end{cases}$

is the Heaviside function. The E and I population firing rate thresholds,
$\theta_{u}$ and $\theta_{v}$, can differ. Given a Heaviside firing rate
nonlinearity, it is possible to exactly characterize the dynamics of the model
using interface equations that track the evolution of level sets
$u=\theta_{u}$ and $v=\theta_{v}$ in space and time [14, 18]. The strength of
synaptic connectivity weakens with the spatial distance between positions $x$
and $y$ and is given by the distance-dependent synaptic weight profile
functions $w_{ab}(x-y)$ ($a,b\in\\{e,i\\}$). For explicit calculations, we
take these to be symmetric exponentials

(3) $\displaystyle w_{ab}(x)=A_{ab}e^{-\frac{|x|}{\sigma_{ab}}},$

where $a,b\in\\{e,i\\}$ and $A_{ab},\sigma_{ab}\in\mathbb{R}^{\geq 0}$ (non-
negative constants). Commonly used values throughout for the weight profiles
are presented in Table 1. Parameters for $w_{ee}$ were chosen to non-
dimensionalize the E to E connectivity (setting $A_{ee}=0.5$ and
$\sigma_{ee}=1$). Other parameters for the I to I and cross population
synaptic profiles were generally chosen to be broader than the synaptic
footprint of E to E connections, and we generally assume the commonly
identified $80\%$ E and $20\%$ I neuron ratio [1] in turn approximately
determines the effective synaptic amplitude from those populations.

Table 1: Model parameters for Eq. Eq. 1 Parameter |  $A_{ee}$ |  $A_{ei},A_{ie}$ |  $A_{ii}$ |  $\sigma_{ee}$ |  $\sigma_{ei},\sigma_{ie},\sigma_{ii}$ |  $\tau$
---|---|---|---|---|---|---
Definition |  E-E strength |  I-E strength |  I-I strength |  E-E spatial scale |  other spatial scales |  I time constant
Value |  0.5 |  0.15 |  0 or 0.01 |  1 |  2 |  1

The spatially-structured multiplicative noise terms
$dW_{e}=\sqrt{|u(x,t)|}dW_{u}(x,t)$ and $dW_{i}=\sqrt{|v(x,t)|}dW_{v}(x,t)$
are are introduced with small amplitude $0<\epsilon\ll 1$, allowing asymptotic
approximations of their stochastic effects. Multiplicative noise is predicted
in neural fields by linearly approximating noise arising in a system-size
expansion of a fully stochastic model [5]. Increments of spatially-extended
Wiener processes have zero mean $\langle dW_{a}(x,t)\rangle=0$
($a\in\\{u,v\\}$) and local spatial correlations $\langle
dW_{a}(x,t)dW_{a}(y,s)\rangle=C_{aa}(x-y)\delta(t-s)dtds$. For simplicity, we
start by assuming there are no correlations between the noise to the E and I
populations $\langle dW_{u}(x,t)dW_{v}(y,s)\rangle\equiv 0$.

Based on prior studies of deterministic E/I population models [42, 4, 8] bump
(standing pulse) solutions emerge in parameter regions we can identify using
an existence/stability analysis. Recurrent excitation sustains both the E and
I populations and inhibition prevents the spread of the E population activity.
In the absence of noise, we obtain bump profiles that are even symmetric and
translation symmetric. Using threshold-crossing conditions, we develop
implicit formulas for the half-width variable $a_{u}$ ($a_{v}$) for the E (I)
bump, which ensure existence, defined as the distance from the bump’s center
of mass to the interface (Fig. 1B). The I bump’s half-width $a_{v}$ (Fig. 1C)
can vary non-monotonically along certain parameter axes, such as increasing
firing rate thresholds, $\theta_{u}$ and $\theta_{v}$, simultaneously. When
noise is applied to stable E/I activity bumps, they ‘wander’ due to their
neutral stability. In addition to small deformations of the bump profile
itself, a bump’s center of mass wanders stochastically (Fig. 1D and see Fig.
A.1 for an additional example). While past asymptotic analyses assumed this
wandering was roughly Brownian motion [8, 31], we will show that by relaxing
this assumption we can better characterize nonlinear interactions between the
bumps in the distinct E and I populations, and how this shapes wandering
dynamics.

Each bump has a region over which neural activity ($u$ or $v$) is
superthreshold (above $\theta_{u}$ or $\theta_{v}$). We define the
corresponding E/I population active regions to be $[x_{1}(t),x_{2}(t)]$ and
$[x_{3}(t),x_{4}(t)]$ respectively [27]. These active regions are bounded by
the interfaces (threshold crossings), since $u(x_{1,2}(t),t)=\theta_{u}$ and
$v(x_{3,4}(t),t)=\theta_{v}$. Note, large deviations could lead to multiple
disjoint active region segments in each layers, but we assume each bump’s
active region is fully connected here (See [14, 18, 36] for elaborations on
this problem). In the absence of noise, the interfaces relate directly to the
half-widths in that $a_{u}=(x_{2}-x_{1})/2$ and $a_{v}=(x_{4}-x_{3})/2$.

Stochastic motion of the bumps will be tracked by estimating the center of
mass $\Delta_{u}(t)$ ($\Delta_{v}(t)$) of the active region of the E (I) bump

(4) $\displaystyle\Delta_{u}(t)=\frac{x_{1}(t)+x_{2}(t)}{2},\hskip
56.9055pt\Delta_{v}(t)=\frac{x_{3}(t)+x_{4}(t)}{2}.$

Overall, we are interested both in how network parameters shape the form and
stability of bumps as well as how this translates into bump’s stochastic
dynamics in the presence of noise. Our ensuing analysis will initially rely on
techniques in local stability as well as weak perturbations, but a major
advancement will be the use of interface techniques to provide higher-order
corrections to the effective nonlinear equations describing stochastic bump
motion.

Figure 1: Independence of Excitatory and Inhibitory bumps. (A) Model schematic
with separate E and I populations and the synaptic connections between them.
(B) An example pair of E and I bump profiles with firing rate thresholds
$\theta_{u}=0.35$ and $\theta_{v}=0.25$ and $A_{ii}=0$. (C) Example plot of
solution half-widths as thresholds are varied, obtained by numerically solving
the threshold conditions, Eq. (8), for $a_{u}$ and $a_{v}$. The particular
half-width solution from panel B is identified by the corresponding dots. (D)
Example bumps from panel B wandering over time with noise amplitude
$\epsilon=0.002$. Traces represent the center of mass (green) and the E/I
interfaces (red/blue, respectively). Colorscales represent the profile
amplitude of $u(x,t)$ and $v(x,t)$. Other parameters are as in Table 1.

##  3 Deterministic Analysis

Our estimates for expected bump wandering (variance of center of mass)
utilizes linearization about stable stationary solutions. This necessitates an
analysis of stationary bump solutions’ (static bump solutions in the absence
of perturbations) existence and stability in varying parameter regimes for the
deterministic system. Using Eq. Eq. 1 without noise [4, 27], we demonstrate
the existence of stationary solutions through direct construction. For proofs
of uniqueness of solutions with a Heaviside firing rate see [4, 10].

###  3.1 Stationary Solutions

Assuming solutions are stationary ($u(x,t)=U(x)$ and $v(x,t)=V(x)$) and using
the Heaviside firing rate function Eq. Eq. 2, Eq. Eq. 1 becomes

(5a) $\displaystyle
U(x)=\int_{\mathbb{R}}w_{ee}(x-y)H(U(y)-\theta_{u})dy-\int_{\mathbb{R}}w_{ei}(x-y)H(V(y)-\theta_{v})dy$
(5b) $\displaystyle
V(x)=\int_{\mathbb{R}}w_{ie}(x-y)H(U(y)-\theta_{u})dy-\int_{\mathbb{R}}w_{ii}(x-y)H(V(y)-\theta_{v})dy$

We seek bumps with simply connected active regions ($U>\theta_{u}$ and
$V>\theta_{v}$) $[-a_{u},a_{u}]$ and $[-a_{v},a_{v}]$. Exploiting the
solutions’ translational invariance along $\mathbb{R}$ (opting for a center of
mass at $x=0$) and expected evenness (i.e. $U(-x)=U(x)$ and $V(-x)=V(x)$) we
obtain the system

(6a) $\displaystyle
U(x)=\int_{-a_{u}}^{a_{u}}w_{ee}(x-y)dy-\int_{-a_{v}}^{a_{v}}w_{ei}(x-y)dy$
(6b) $\displaystyle
V(x)=\int_{-a_{u}}^{a_{u}}w_{ie}(x-y)dy-\int_{-a_{v}}^{a_{v}}w_{ii}(x-y)dy$

Integrating our chosen exponential synaptic weight functions
($w_{ab}(x)=A_{ab}e^{\frac{-|x|}{\sigma_{ab}}}$, $a,b\in\\{e,i\\}$), we find
the explicit formulas for $c\in\\{a_{u},a_{v}\\}$:

(7)
$\displaystyle\int^{c}_{-c}w_{ab}(x-y)dy=\begin{cases}2A_{ab}\sigma_{ab}e^{\frac{-x}{\sigma_{ab}}}\sinh\left(\frac{c}{\sigma_{ab}}\right)&x>c,\\\
2A_{ab}\sigma_{ab}\left[1-e^{\frac{-c}{\sigma_{ab}}}\cosh\left(\frac{x}{\sigma_{ab}}\right)\right]&|x|<c,\\\
2A_{ab}\sigma_{ab}e^{\frac{x}{\sigma_{ab}}}\sinh\left(\frac{c}{\sigma_{ab}}\right)&x<-c.\end{cases}$

Substituting Eq. Eq. 7 into Eq. Eq. 6 and utilizing the threshold crossing
conditions $\theta_{u}=U(\pm a_{u})$ and $\theta_{v}=V(\pm a_{v})$, we obtain
an implicit set of equations for the half-widths, which depends on whether the
E or I bump is wider:

(8a)
$\displaystyle\theta_{u}=2A_{ee}\sigma_{ee}e^{\frac{-a_{u}}{\sigma_{ee}}}\sinh\left(\frac{a_{u}}{\sigma_{ee}}\right)-\begin{cases}2A_{ei}\sigma_{ei}e^{\frac{-a_{u}}{\sigma_{ei}}}\sinh\left(\frac{a_{v}}{\sigma_{ei}}\right)&a_{u}\geq
a_{v},\\\
2A_{ei}\sigma_{ei}\left[1-e^{\frac{-a_{v}}{\sigma_{ei}}}\cosh\left(\frac{a_{u}}{\sigma_{ei}}\right)\right]&a_{u}<a_{v},\end{cases}$
(8b)
$\displaystyle\theta_{v}=-2A_{ii}\sigma_{ii}e^{\frac{-a_{v}}{\sigma_{ii}}}\sinh\left(\frac{a_{v}}{\sigma_{ii}}\right)+\begin{cases}2A_{ie}\sigma_{ie}e^{\frac{-a_{v}}{\sigma_{ie}}}\sinh\left(\frac{a_{u}}{\sigma_{ie}}\right)&a_{v}\geq
a_{u},\\\
2A_{ie}\sigma_{ie}\left[1-e^{\frac{-a_{u}}{\sigma_{ie}}}\cosh\left(\frac{a_{v}}{\sigma_{ie}}\right)\right]&a_{v}<a_{u}.\end{cases}$

Eq. Eq. 8 is thus defined piecewise continuously (Fig. 2A, gold circle
indicates the change in cases occurs) and, as we will show, cusps appear at
these case switches in plots of eigenvalues and variance estimates.

Note, the above set of threshold conditions assume that both the E and I
population have nontrivial active regions. There is a second class of bump
solutions where there is no active region in the I population. As we will
show, these bumps always tend to be linearly unstable as there is no active
inhibition to prevent the spread of excitation upon perturbation. Applying
this assumption to Eq. Eq. 6, we obtain the simplified system,

(9a) $\displaystyle U(x)=\int_{-a_{u}}^{a_{u}}w_{ee}(x-y)dy,$ (9b)
$\displaystyle V(x)=\int_{-a_{u}}^{a_{u}}w_{ie}(x-y)dy,$

with E bump half-width given by

$\theta_{u}=U(\pm a_{u})=\int_{-a_{u}}^{a_{u}}w_{ee}(\pm a_{u}-y)dy$

or (for our particular weight functions) the formula

(10)
$\theta_{u}=2A_{ee}\sigma_{ee}e^{\frac{-a_{u}}{\sigma_{ee}}}\sinh\left(\frac{a_{u}}{\sigma_{ee}}\right).$

For such a solution to be self-consistent, we also must ensure that
$V<\theta_{v}$ everywhere. Assuming the peak of the subthreshold I bump occurs
at $x=0$, this is ensured by the condition

$\displaystyle
V(0)=\int_{-a_{u}}^{a_{u}}w_{ie}(y)dy=2A_{ie}\sigma_{ie}\left[1-e^{\frac{-a_{u}}{\sigma_{ie}}}\right]<\theta_{v}.$

There are thus, two branches of bump solutions. One branch of ‘Broad’ bumps
has superthreshold active regions in both the E and I populations ($a_{u}>0$
and $a_{v}>0$). For sufficiently high thresholds, we typically obtain
marginally stable solutions, but these destabilize through a Hopf bifurcation
for sufficiently low firing rate thresholds (See Fig. 2A and subsequent
stability analysis). In contrast, when the I bump is subthreshold
($V(0)<\theta_{v}$), we obtain unstable ‘Narrow’ solutions that create a
separatrix between the Broad solutions and the rest state (Fig. 2A). Under
each set of assumptions we obtain two distinct systems that approach the same
solution defining a discontinuous saddle node bifurcation. The discontinuity
arises from there being either two or one threshold condition, and the peak of
the I Narrow bump grazing the threshold $\theta_{v}$. We derive these results
in detail in the following subsection.

###  3.2 Eigenvalues and Noiseless Perturbations

Our stability analysis utilizes linearization about stationary solutions and
localizes the perturbation evolution problem to the bump edges, as in several
previous analyses of bump dynamics in neural fields [4, 12, 19, 42, 27]. Since
there are four threshold crossing points, analysis of the Broad solution’s
stability results in four corresponding eigenvalues and equations associated
with the degrees of freedom in the stability problem (Fig. 2B). To derive our
linearized system whose spectrum defines the stability of bumps, we start by
perturbing with small smooth functions, $u(x,t)\approx U(x)+\epsilon\psi(x,t)$
and $v(x,t)\approx V(x)+\epsilon\phi(x,t)$. Substituting into Eq. Eq. 5,
expanding about the stationary solution, and simplifying to first order we
obtain the system:

(11a) $\displaystyle\psi_{t}+\psi$
$\displaystyle=w_{ee}(x)*[H^{\prime}(U(x))\psi(x,t)]-w_{ei}(x)*[H^{\prime}(V(x))\phi(x,t)],$
(11b) $\displaystyle\tau\phi_{t}+\phi$
$\displaystyle=w_{ie}(x)*[H^{\prime}(U(x))\psi(x,t)]-w_{ii}(x)*[H^{\prime}(V(x))\phi(x,t)],$

where the distributional derivatives are given

$\displaystyle H^{\prime}(U(x))$
$\displaystyle=\frac{1}{|U^{\prime}(a_{u})|}(\delta(x-a_{u})+\delta(x+a_{u})),$
$\displaystyle H^{\prime}(V(x))$
$\displaystyle=\frac{1}{|V^{\prime}(a_{v})|}(\delta(x-a_{v})+\delta(x+a_{v})).$

Assuming separability of our perturbations, $\psi(x,t)=\psi(x)e^{\lambda t}$
and $\phi(x,t)=\phi(x)e^{\lambda t}$, and substituting our exponential weight
functions and $H^{\prime}(U(x))$ and $H^{\prime}(V(x))$, we obtain the
corresponding eigenvalue problem:

(12a)
$\begin{split}(\lambda+1)\psi(x)=\frac{A_{ee}}{|U^{\prime}(a_{u})|}\left[e^{\frac{-|x-a_{u}|}{\sigma_{ee}}}\psi(a_{u})+e^{\frac{-|x+a_{u}|}{\sigma_{ee}}}\psi(-a_{u})\right]\\\
-\frac{A_{ei}}{|V^{\prime}(a_{v})|}\left[e^{\frac{-|x-a_{v}|}{\sigma_{ei}}}\phi(a_{v})+e^{\frac{-|x+a_{v}|}{\sigma_{ei}}}\phi(-a_{v})\right],\end{split}$
(12b)
$\begin{split}(\tau\lambda+1)\phi(x)=\frac{A_{ie}}{|U^{\prime}(a_{u})|}\left[e^{\frac{-|x-a_{u}|}{\sigma_{ie}}}\psi(a_{u})+e^{\frac{-|x+a_{u}|}{\sigma_{ie}}}\psi(-a_{u})\right]\\\
-\frac{A_{ii}}{|V^{\prime}(a_{v})|}\left[e^{\frac{-|x-a_{v}|}{\sigma_{ii}}}\phi(a_{v})+e^{\frac{-|x+a_{v}|}{\sigma_{ii}}}\phi(-a_{v})\right].\end{split}$

Note, in the above, we have assumed $\psi(\pm a_{u})\neq 0$ and $\phi(\pm
a_{v})\neq 0$ to obtain an equation for the point spectrum, but taking these
values to vanish would give us an equation for the essential spectrum which
will not contribute to instabilities [13, 45]. We form a system of equations
by evaluating the perturbation functions at the bump edges $x=\pm a_{u}$ and
$x=\pm a_{v}$ into Eq. (12). Values $\lambda$ for which the resulting $4\times
4$ system is singular correspond to the four distinct eigenvalues in the point
spectrum. Forming this system, taking the determinant, and factoring, we find
the eigenvalues are the roots of the following pair of quadratics:

(13a) $\displaystyle\tau\lambda^{2}-(I+J+\tau B+\tau
C)\lambda+(I+J)(B+C)-(E+D)(F+G)$ (13b) $\displaystyle\tau\lambda^{2}-(I-J+\tau
B-\tau C)\lambda+(I-J)(B-C)+(E-D)(F-G)$

where

$\displaystyle B=\frac{A_{ee}}{|U^{\prime}(a_{u})|}-1$ $\displaystyle\hskip
72.26999ptC=\frac{A_{ee}}{|U^{\prime}(a_{u})|}e^{-\frac{2a_{u}}{\sigma_{ee}}}$
$\displaystyle
D=-\frac{A_{ei}}{|V^{\prime}(a_{v})|}e^{\frac{-|a_{v}-a_{u}|}{\sigma_{ei}}}$
$\displaystyle\hskip
72.26999ptE=-\frac{A_{ei}}{|V^{\prime}(a_{v})|}e^{\frac{-|a_{v}+a_{u}|}{\sigma_{ei}}}$
$\displaystyle
F=\frac{A_{ie}}{|U^{\prime}(a_{u})|}e^{\frac{-|a_{v}-a_{u}|}{\sigma_{ie}}}$
$\displaystyle\hskip
72.26999ptG=\frac{A_{ie}}{|U^{\prime}(a_{u})|}e^{\frac{-|a_{v}+a_{u}|}{\sigma_{ie}}}$
$\displaystyle I=-\frac{A_{ii}}{|V^{\prime}(a_{v})|}-1$ $\displaystyle\hskip
72.26999ptJ=-\frac{A_{ii}}{|V^{\prime}(a_{v})|}e^{-\frac{2a_{v}}{\sigma_{ii}}}$

Each quadratic could also be obtained by restricting the form of perturbation
(scaling and shifting) at the interfaces (See Fig. 2D for examples). Scaling
perturbations expand or contract the bump ($\psi(a_{u})=\psi(-a_{u})$ and
$\phi(a_{v})=\phi(-a_{v})$), and we obtain Eq. Eq. 13a resulting in the
associated ‘scale’ eigenvalues (Fig. 2B):

(14)
$\lambda_{\text{scale}}=\frac{I+J+\tau(B+C)}{2\tau}\pm\frac{\sqrt{(I+J-\tau(B+C))^{2}+4\tau(E+D)(F+G)}}{2\tau},$

which are both nonzero in general. The sign of their real part determines the
linear stability of the bumps. For shifting perturbations
($\psi(a_{u})=-\psi(-a_{u})$ and $\phi(a_{v})=-\phi(-a_{v})$) we obtain Eq.
Eq. 13b. A straightforward calculation shows $(I-J)(B-C)+(E-D)(F-G)=0$,
resulting in the zero eigenvalue (green) corresponding to the system’s
translational invariance (Fig. 2B). Correlated shifts in the E/I bumps’ center
of mass simply translate the solution along the real line. Such perturbations
are integrated, and neither grow nor decay. The nonzero shift eigenvalue is
given by the simple formula

(15) $\lambda=\frac{I-J+\tau(B-C)}{\tau}$

which is real and negative. Hence the bumps are stable with respect to shifts
of the E and I bump in the opposite direction.

Figure 2: Half-widths of bump solutions and linear stability. (A) Broad/Narrow
half-width solutions $a_{u}$ and $a_{v}$ as a function of firing rate
threshold $\theta_{u}=\theta_{v}=\theta$. There is a Hopf bifurcation (HB,
green dots) in the Broad solutions and a semi-stable discontinuous saddle node
bifurcation (SN, white dots) where the Narrow and Broad solutions meet. The
gold circle identifies where the E/I bump half-widths exchange size ordering,
Eq. Eq. 9. Eigenvalues are plotted in association with stability of the (B)
Broad and (C) Narrow solutions. (D) Examples of four perturbations types
related to the scale and shift Broad solution eigenvalues. (E) An example of
an unstable Broad solution destabilized via the oscillatory instability
emerging from the Hopf bifurcation, $\theta_{u}=\theta_{v}=0.1$. (F) A bump
collapsing for threshold parameters near the discontinuous saddle node,
$\theta_{u}=\theta_{v}=0.4996$. We take the bump-like solution initially
narrower than the Broad solution. The I-I connectivity $A_{ii}=0$ in all.
Other parameters are as in Table 1.

Eigenvalues of the Narrow solutions are obtained by following the same
process, which is greatly simplified since the I population is subthreshold
and perturbations of this part do not contribute to instabilities or enter
into the point spectrum equations. We therefore investigate the effect of
small smooth perturbations to the E population, with nonzero part along the
bump edges, $u(x,t)\approx U(x)+\epsilon\psi(x,t)$. The resulting two
eigenvalues are $\lambda_{1}=0$ (corresponding to the marginal stability of
shifts) and
$\lambda_{2}={2e^{\frac{-2a_{u}}{\sigma_{ee}}}}/{(1-e^{\frac{-2a_{u}}{\sigma_{ee}}})}$,
which is clearly positive, confirming that Narrow solutions are always
unstable (Fig. 2C).

We identify bifurcations in the bump solutions by checking the signs of the
real parts of each solution branch’s eigenvalues. We find two types of
bifurcations. The first is a Hopf bifurcation (green dot, HB in Fig. 2A,B)
occurring as the firing threshold is reduced, destabilizing bumps in a
pattern-destroying oscillation (Fig. 2E) [42]. This Hopf bifurcation boundary
occurs when $Re(\lambda_{\text{scale}})=0$, which is given by the implicit
formula

(16) $\frac{\tau
A_{ee}}{|U^{\prime}(a_{u})|}(1-e^{\frac{-2a_{u}}{\sigma_{ee}}})=\tau+1+\frac{A_{ii}}{|V^{\prime}(a_{v})|}(1+e^{\frac{-2a_{v}}{\sigma_{ii}}}).$

The second bifurcation observed is a discontinuous saddle node bifurcation (SN
in Fig. 2A) [15, 39] where the unstable Narrow and marginally stable Broad
solutions meet. The peak of the Narrow I population bump rises to meet the
firing threshold and the Broad I bump active region shrinks to a single point
at threshold. The meeting of two linearizations of different discrete
dimension (4 for Broad, 2 for Narrow) results in a discontinuous change in the
corresponding Jacobian matrices yielding a semistable nonsmooth bifurcation
where the solution branches annihilate one another (labelled “contraction
instability” in Fig. 2F).

###  3.3 Bounding stability regions in parameter space

Figure 3: Bump stability/instability across parameter space. Main panels take
$A_{ii}=0$ and insets take $A_{ii}=0.01$. In all plots, blue regions indicate
parameters for which a stable bump exists, red regions are where all bumps are
unstable. Black boundaries identify where Hopf bifurcations occur. White
regions are where no bump solutions exist. (A) Stability as the firing rate
thresholds $\theta$ are varied. The silver line identifies where
$\theta_{u}=\theta_{v}$. (B) Stability regions with varied E/I interaction
amplitude $A_{ei}$ and $A_{ie}$. $\theta_{u}=\theta_{v}=0.25$. (C) Stability
regions with varied inhibition timescale $\tau$ and firing thresholds
$\theta_{u}=\theta_{v}=\theta$. (D) Stability regions with varied E/I
interaction spatial extent $\sigma_{ei}=\sigma_{ii}$ and $\sigma_{ie}$.
$\theta_{u}=\theta_{v}=0.15$. Silver line identifies where
$\sigma_{ei}=\sigma_{ii}=\sigma_{ie}$. Other parameters are as in Table 1.

To identify how noise degrades stable representations of parametric working
memory, we focus on the effect of stochastic perturbations to bumps along the
Broad branch. Given that varying certain parameters can stabilize the
deterministic system Eq. Eq. 1, we identify how the bifurcation boundaries
change and stable solution regions expand/contract as parameters are varied
(Fig. 3). As discussed above, stable solutions only exist if Broad bumps have
non-positive eigenvalues associated with their linear stability.

When varying the firing rate thresholds, we find larger I thresholds,
$\theta_{v}$, expand the range of stable solutions (Fig. 3A). We speculate
that setting the I population firing rate threshold too low makes it strongly
responsive to network activity perturbations, as shown in oscillatory
instability simulations (Fig. 2E). Varying the strengths of synapses from E to
I and I to E populations, we also observed bumps destabilize if E to I
($A_{ie}$) is sufficiently strong (Fig. 3B). Such changes increase the
sensitivity of the E drive to the I population as bumps are expanded, leading
to potential I bump overcompensation of as E bumps widen, eventually silencing
both bumps.

Setting firing rate thresholds equal ($\theta_{u}=\theta_{v}=\theta$) and
varying them (silver line in Fig. 3A) along with changing the I timescale,
$\tau$, results in a narrowing of the range of stable bumps as firing rate
thresholds are decreased (Fig. 3C). Specifically, as the timescale of the I
population is increased, the system reacts more slowly to being out of
equilibrium, so the E and I bump do not restabilize once perturbed. For
instantaneous inhibition ($\tau\to 0$), such oscillatory (Hopf) instabilities
never occur.

We also quantified the parameter ranges of stable bumps when varying the
spatial scale of interpopulation synaptic footprints (Fig. 3D). Stability is
again largely dependent on the profile of E to I population synaptic
connectivity; wider connectivity (high $\sigma_{ie}$) yields an unstable
region. Similar to our previous findings, broadening the E to I profile leads
to overcompensation of the I population in response to dynamical increases in
the E population width.

The main panels of Fig. 3 assume $A_{ii}=0$, common in studies of Eq. Eq. 1
due to the small amplitude of such connections. We also examined small,
nonzero amplitude I to I connections ($A_{ii}=0.01$, See insets in Fig. 3) and
found stable regions expand likely due to dampening of I population reactions.

##  4 Analysis of stochastic bump motion

Considering stochasticity emerging from neural and synaptic variability [17],
multiplicative noise in the neural field model Eq. (1) (taking $\epsilon>0$)
causes stable bumps to wander like a particle subject to Brownian motion. The
extent of this wandering has been linked to subject response errors on delayed
estimation tasks [48, 8]. We are primarily interested in identifying how
architectural features of the E/I neural circuit impact wandering [31, 30, 44,
32] as described by the time-dependent variance of bumps’ center of mass. We
can estimate this variance analytically using two different approaches. The
first we refer to as the ‘Strongly coupled limit’ approximation, motivated by
work from [31], which assumes the E/I bump has a single position (the E and I
bump do not stray too far from one another). To first order, this
approximation treats the bump as wandering by pure diffusion, and so its
stochastic motion is only characterized by a diffusion coefficient. However,
the resulting formulas have limitations for certain parameter regimes on the
timescale of interest, since they do not consider that the bumps can drift
apart. Our second approach, the ‘Interface-based approximation,’ tracks bump
interfaces (following threshold crossing points) [36, 14, 27]. Pairing this
together with a weak coupling assumption as in [28] allows us to estimate the
time-dependent changes in the distinct E and I bump centers of mass.

###  4.1 Strongly coupled limit approximation

In the Strongly coupled limit, we assume solutions to Eq. Eq. 1 take the form
of stable bumps both with position weakly perturbed by the same amount
($U(x-\Delta_{u})$ and $V(x-\Delta_{v})$) with $\Delta_{u}=\Delta_{v}=\Delta$.
Thus, the E and I bumps are assumed to move together, and we have

(17a) $\displaystyle u(x,t)$
$\displaystyle=U(x-\Delta(t))+\epsilon^{1/2}\Phi(x-\Delta(t),t)+\epsilon\Phi_{1}(x-\Delta(t),t)\dots,$
(17b) $\displaystyle v(x,t)$
$\displaystyle=V(x-\Delta(t))+\epsilon^{1/2}\Psi(x-\Delta(t),t)+\epsilon\Psi_{1}(x-\Delta(t),t)\dots,$

where $\Delta(t)$ is the $\mathcal{O}(\epsilon^{1/2})$ position perturbation
(with $d\Delta=\mathcal{O}(\epsilon^{1/2})$) and $\Phi$ and $\Psi$ are
respectively the leading order profile perturbations to the E and I bump.
Substituting Eq. Eq. 17 and truncating Eq. Eq. 1 to first order and taking
averages, we find that the stationary solutions remain the same form as
before. Moving to $\mathcal{O}(\epsilon^{1/2})$, we find

(26) $\displaystyle\left(\begin{array}[]{c}d\Phi\\\ \tau\cdot
d\Psi\end{array}\right)=\mathcal{L}\left(\begin{array}[]{c}\Phi\\\
\Psi\end{array}\right)dt+\epsilon^{-1/2}d\Delta(t)\left(\begin{array}[]{c}U^{\prime}(x)\\\
\tau\cdot
V^{\prime}(x)\end{array}\right)+\left(\begin{array}[]{c}\sqrt{|U(x)|}dW_{u}\\\
\sqrt{|V(x)|}dW_{v}\end{array}\right),$

where

$\mathcal{L}\left(\begin{array}[]{c}p\\\
q\end{array}\right)=\left(\begin{array}[]{c}-p+w_{ee}*[f^{\prime}(U)p]-w_{ei}*[f^{\prime}(V)q]\\\
-q+w_{ie}*[f^{\prime}(U)p]-w_{ii}*[f^{\prime}(V)q]\end{array}\right).$

Note that due to translation symmetry of the noise-free system,
$\mathcal{N}(\mathcal{L})$ is spanned by $\\{U^{\prime},V^{\prime}\\}$. To
briefly show this, recall the stationary solutions take the form

$\displaystyle U(x)$
$\displaystyle=\int_{\mathbb{R}}w_{ee}(x-y)f(U(y))dy-\int_{\mathbb{R}}w_{ei}(x-y)f(V(y))dy$
$\displaystyle V(x)$
$\displaystyle=\int_{\mathbb{R}}w_{ie}(x-y)f(U(y))dy-\int_{\mathbb{R}}w_{ii}(x-y)f(V(y))dy.$

Differentiating with respect to $x$ and applying integration by parts yields

$\displaystyle U^{\prime}(x)$
$\displaystyle=\int_{\mathbb{R}}w_{ee}(x-y)f^{\prime}(U(y))U^{\prime}(y)dy-\int_{\mathbb{R}}w_{ei}(x-y)f^{\prime}(V(y))V^{\prime}(y)dy$
$\displaystyle V^{\prime}(x)$
$\displaystyle=\int_{\mathbb{R}}w_{ie}(x-y)f^{\prime}(U(y))U^{\prime}(y)dy-\int_{\mathbb{R}}w_{ii}(x-y)f^{\prime}(V(y))V^{\prime}(y)dy,$

exactly the form of the functions $p$ and $q$ spanning the
$\mathcal{N}(\mathcal{L})$.

We ensure a bounded solution to Eq. Eq. 26 by requiring the inhomogeneous part
to be orthogonal to $\mathcal{N}(\mathcal{L^{*}})$. There exists a single
vector spanning the ${\mathcal{N}}(\mathcal{L^{*}})$, denoted
$(\varphi_{1}(x),\varphi_{2}(x))$. Enforcing our conditions for bounded
solutions via requiring the inner product of the nullspace
$(\varphi_{1}(x),\varphi_{2}(x))$ and inhomogeneity $(h_{1}(x),h_{2}(x))$
vanishes $\int_{\mathbb{R}}\varphi_{1}(x)h_{1}(x)+\varphi_{2}(x)h_{2}(x)dx=0$,
we isolate the bump position increment:

$d\Delta(t)=-\epsilon^{1/2}\frac{\int_{\mathbb{R}}[\varphi_{1}(x)\sqrt{|U(x)|}dW_{u}(x,t)+\varphi_{2}(x)\sqrt{|V(x)|}dW_{v}(x,t)]dx}{\int_{\mathbb{R}}[\varphi_{1}(x)U^{\prime}(x)+\tau
V^{\prime}(x)\varphi_{2}(x)]dx}.$

Since the above is simply a weighted integral over the spatiotemporal noises,
to first order, we have a Brownian stochastic differential equation (SDE) for
the bump position. Given noise is white in time, we find that the mean over
realizations is $\langle\Delta(t)\rangle=0$, and the variance is:

(27)
$\displaystyle\langle\Delta(t)^{2}\rangle=\epsilon\frac{\int_{\mathbb{R}}\varphi_{1}(x)\cdot\left[C_{u}(x)*\varphi_{1}(x)+2C_{c}(x)*\varphi_{2}(x)\right]+\varphi_{2}(x)\cdot\left[C_{v}(x)*\varphi_{2}(x)\right]dx}{\left[\int_{\mathbb{R}}\varphi_{1}(x)U^{\prime}(x)+\tau
V^{\prime}(x)\varphi_{2}(x)dx\right]^{2}}t$

where the spatial correlation functions are defined: $\langle
W_{u}(x,t)W_{u}(y,t)\rangle=C_{u}(x-y)t$, $\langle
W_{v}(x,t)W_{v}(y,t)\rangle=C_{v}(x-y)t$, and $\langle
W_{u}(x,t)W_{v}(y,t)\rangle=C_{c}(x-y)t$. Note, there is a weighted
contribution to the linear scaling variance from both the noise of the E and I
populations.

To calculate the nullspace of the adjoint linearized operator, we rearrange
the equation
$\mathcal{L}^{*}\left(\varphi_{1}(x),\varphi_{2}(x)\right)^{T}=0$, so

(28a) $\displaystyle\varphi_{1}(x)$
$\displaystyle=f^{\prime}(U)\int_{\mathbb{R}}[w_{ee}(x-y)\varphi_{1}(y)+w_{ie}(x-y)\varphi_{2}(y)]dy,$
(28b) $\displaystyle\varphi_{2}(x)$
$\displaystyle=-f^{\prime}(V)\int_{\mathbb{R}}[w_{ei}(x-y)\varphi_{1}(y)+w_{ii}(x-y)\varphi_{2}(y)]dy.$

In the case of the Heaviside firing rate function we again can determine the
distributional derivative of the Heavisde acting on the bump solution and
obtain

(29) $\displaystyle\begin{pmatrix}\varphi_{1}(x)\\\
\varphi_{2}(x)\end{pmatrix}=\begin{pmatrix}\delta(x+a_{u})+\mathcal{A}\delta(x-a_{u})\\\
-\mathcal{B}\delta(x+a_{v})+\mathcal{C}\delta(x-a_{v})\end{pmatrix}.$

Plugging Eq. Eq. 29 into Eq. Eq. 28 and solving for the constants, we find
$\mathcal{A}=-1$, $\mathcal{C}=\mathcal{B}$, and
$\mathcal{B}=\frac{w_{ei}(a_{u}-a_{v})-w_{ei}(a_{u}+a_{v})}{w_{ie}(a_{u}-a_{v})-w_{ie}(a_{u}+a_{v})}.$
Thus we find that

(30) $\displaystyle\begin{pmatrix}\varphi_{1}(x)\\\
\varphi_{2}(x)\end{pmatrix}=\begin{pmatrix}\delta(x+a_{u})-\delta(x-a_{u})\\\
-\mathcal{B}[\delta(x+a_{v})-\delta(x-a_{v})]\end{pmatrix}.$

Finally, substituting Eq. Eq. 30 into Eq. Eq. 27, we have

(31) $\displaystyle\langle\Delta(t)^{2}\rangle$
$\displaystyle=\epsilon\frac{D_{1}-D_{2}+D_{3}}{2\left[|U^{\prime}(a_{u})|+\mathcal{B}\tau|V^{\prime}(a_{v})|\right]^{2}}t$

where

$\displaystyle D_{1}$ $\displaystyle=\theta_{u}[C_{u}(0)-C_{u}(2a_{u})],$
$\displaystyle D_{2}$
$\displaystyle=2\mathcal{B}\sqrt{\theta_{u}\theta_{v}}[C_{c}(a_{u}-a_{v})-C_{c}(a_{u}+a_{v})],$
$\displaystyle D_{3}$
$\displaystyle=\theta_{v}\mathcal{B}^{2}[C_{v}(0)-C_{v}(2a_{v})].$

The contributions arising due to noise correlations within the E and I
populations, $D_{1}$ and $D_{3}$, are positive. In contrast, the diffusion
contribution from the cross-population correlation term $C_{c}$ is negative
assuming the correlation function is monotone decreasing. Thus, noise
correlations between the E and I populations reduce wandering. Here, we assume
$C_{c}\equiv 0$. Later, we discuss the effects of nonzero correlated noise
components across the E and I populations and compare these theoretical
results to numerical simulations.

###  4.2 Interface Based Approximation Theory

Complementing our Strongly coupled limit approximation of bump wandering, we
also derive interface equations that approximate the stochastic dynamics of
the bump in response to multiplicative noise in Eq. Eq. 1. While the Strongly
Coupled Limit approximation assumes the E and I bump move as one, numerical
simulations reveal that this is not the case (See Fig. 4A and C). For
instance, the I bump tends to be more susceptible to profile deformations and
may wander more, making the Strongly coupled limit less valid. Higher-order
corrections of the variance estimate can be obtained using interface theory to
develop nonlinearly coupled Langevin equations for estimating the E and I
bumps’ coupled centers of mass.

Interface theory for bumps in neural fields tracks the level sets of neural
activity variables where firing rate thresholds are crossed. This approach was
originally pioneered in the case of noiseless single-bump neural field models
[27] and then subsequently for traveling fronts in inhomogeneous neural fields
[11] as well as solutions in planar neural fields [14]. More recently, these
approaches were adapted to obtain higher-order approximations for the
timescale of front initiation [18] as well as the stochastic motion of
multiple interacting bumps [36]. Here, we further adapt this work [36] to
account for the motion of separate bumps in the E and I neural populations of
Eq. Eq. 1.

As defined previously, the active regions of the E and I bump are
$[x_{1}(t),x_{2}(t)]$ and $[x_{3}(t),x_{4}(t)]$ respectively. Without noise
perturbations, we would expect these interfaces to remain constant for the
equilibrium bump solution, but noise perturbs these values so they wander over
time. However, unlike in the Strongly coupled limit, we do not expect this
wandering to be pure Brownian motion, but rather the result of a nonlinearly
coupled system of SDEs. To obtain these SDEs, we start by writing the level
set condition for the interfaces in each neural population (E and I):

(32) $\displaystyle u(x_{1}(t),t)=u(x_{2}(t),t)=\theta_{u}\hskip
28.45274ptv(x_{3}(t),t)=v(x_{4}(t),t)=\theta_{v},$

where $\theta_{u}$ and $\theta_{v}$ are the firing rate thresholds of the E
and I populations. For a Heaviside firing rate, we substitute Eq. Eq. 32 into
Eq. Eq. 1 and obtain

(33a) $\displaystyle du$
$\displaystyle=\left[-u+\int^{x_{2}(t)}_{x_{1}(t)}w_{ee}(x-y)dy-\int^{x_{4}(t)}_{x_{3}(t)}w_{ei}(x-y)dy\right]dt+\epsilon^{\frac{1}{2}}dW_{e}$
(33b) $\displaystyle\tau dv$
$\displaystyle=\left[-v+\int^{x_{2}(t)}_{x_{1}(t)}w_{ie}(x-y)dy-\int^{x_{4}(t)}_{x_{3}(t)}w_{ii}(x-y)dy\right]dt+\epsilon^{\frac{1}{2}}dW_{i}$

Differentiating Eq. Eq. 32, we obtain consistency equations for the interfaces
$x_{j}(t)$:

$\displaystyle\partial_{x}u(x_{j}(t),t)dx_{j}+du(x_{j}(t),t)$
$\displaystyle=0,\ \ j=1,2;$
$\displaystyle\partial_{x}v(x_{j}(t),t)dx_{j}+dv(x_{j}(t),t)$
$\displaystyle=0,\ \ j=3,4.$

We then approximate the above exact evolution equations by assuming the
spatial gradients at the interfaces remain constant and odd symmetric
throughout the evolution motivated by the form of the Strongly coupled limit
expansion:

$\displaystyle|U^{\prime}(a_{u})|$
$\displaystyle\equiv\alpha_{u}(t)\approx\frac{\partial u(x_{1}(t),t)}{\partial
x}=-\frac{\partial u(x_{2}(t),t)}{\partial x},$
$\displaystyle|V^{\prime}(a_{v})|$
$\displaystyle\equiv\alpha_{v}(t)\approx\frac{\partial v(x_{3}(t),t)}{\partial
x}=-\frac{\partial v(x_{4}(t),t)}{\partial x}.$

Subsequently, we drop $\mathcal{O}(\epsilon)$ terms to obtain the relations

(34a) $\displaystyle du(x_{1}(t),t)$ $\displaystyle=-\alpha_{u}dx_{1}(t),$
(34b) $\displaystyle du(x_{2}(t),t)$ $\displaystyle=\alpha_{u}dx_{2}(t),$
(34c) $\displaystyle dv(x_{3}(t),t)$ $\displaystyle=-\alpha_{v}dx_{3}(t),$
(34d) $\displaystyle dv(x_{4}(t),t)$ $\displaystyle=\alpha_{v}dx_{4}(t).$

Plugging the formulas in Eq. Eq. 34 into Eq. 33, we obtain the following
system of nonlinear Langevin equations describing the stochastic evolution of
the interfaces:

(35a) $\displaystyle dx_{1}$
$\displaystyle=-\frac{1}{\alpha_{u}}\left(\left[-\theta_{u}+\mathcal{W}_{ee}(x_{1};x_{1},x_{2})-\mathcal{W}_{ei}(x_{1};x_{3},x_{4})\right]dt+\epsilon^{\frac{1}{2}}dW_{e}(x_{1},t)\right)$
(35b) $\displaystyle dx_{2}$
$\displaystyle=\frac{1}{\alpha_{u}}\left(\left[-\theta_{u}+\mathcal{W}_{ee}(x_{2};x_{1},x_{2})-\mathcal{W}_{ei}(x_{2};x_{3},x_{4})\right]dt+\epsilon^{\frac{1}{2}}dW_{e}(x_{2},t)\right)$
(35c) $\displaystyle\tau dx_{3}$
$\displaystyle=-\frac{1}{\alpha_{v}}\left(\left[-\theta_{v}+\mathcal{W}_{ie}(x_{3};x_{1},x_{2})-\mathcal{W}_{ii}(x_{3};x_{3},x_{4})\right]dt+\epsilon^{\frac{1}{2}}dW_{i}(x_{3},t)\right)$
(35d) $\displaystyle\tau dx_{4}$
$\displaystyle=\frac{1}{\alpha_{v}}\left(\left[-\theta_{v}+\mathcal{W}_{ie}(x_{4};x_{1},x_{2})-\mathcal{W}_{ii}(x_{4};x_{3},x_{4})\right]dt+\epsilon^{\frac{1}{2}}dW_{i}(x_{4},t)\right),$

where the coupling functions are given by the integrals over the active
regions

$\displaystyle\mathcal{W}_{jk}(x;x_{a},x_{b})=\int^{x_{b}}_{x_{a}}w_{jk}(x-y)dy.$

To derive estimates for the evolution of the centers of mass of the E and I
bump, we apply the definitions $\Delta_{u}=(x_{1}+x_{2})/2$ and
$\Delta_{v}=(x_{3}+x_{4})/2$ from Eq. Eq. 4 to Eq. Eq. 35 and combine
equations to obtain

(36a)
$\displaystyle\begin{split}d\Delta_{u}=\frac{1}{2\alpha_{u}}\left(\left[-2\mathcal{W}_{ee}(x_{1};x_{1},x_{2})+\mathcal{W}_{ei}(x_{1};x_{3},x_{4})-\mathcal{W}_{ei}(x_{2};x_{3},x_{4})\right]dt\right.\\\
\left.+\sqrt{\epsilon\theta_{u}}[dW_{u}(x_{2},t)-dW_{u}(x_{1},t)]\right),\end{split}$
(36b)
$\displaystyle\begin{split}d\Delta_{v}=\frac{1}{2\tau\alpha_{u}}\left(\left[2\mathcal{W}_{ii}(x_{3};x_{3},x_{4})+\mathcal{W}_{ie}(x_{4};x_{1},x_{2})-\mathcal{W}_{ei}(x_{3};x_{1},x_{2})\right]dt\right.\\\
\left.+\sqrt{\epsilon\theta_{v}}[dW_{v}(x_{4},t)-dW_{v}(x_{3},t)]\right).\end{split}$

Assuming position perturbations remain small we approximate
$x_{1}=\Delta_{u}-a_{u}$, $x_{2}=\Delta_{u}+a_{u}$, $x_{3}=\Delta_{v}-a_{v}$,
and $x_{4}=\Delta_{v}+a_{v}$ with
$\Delta_{u},\Delta_{v}=\mathcal{O}(\epsilon^{1/2})$. Plugging in these
approximations, linearizing about the stationary solutions, and simplifying
yields the multivariate Ornstein-Uhlenbeck (OU) process:

(37a)
$\displaystyle\begin{split}d\Delta_{u}(t)=\frac{1}{\alpha_{u}}\left((\Delta_{u}-\Delta_{v})\cdot\left[w_{ei}(a_{u}-a_{v})-w_{ei}(a_{u}+a_{v})\right]dt\right.\\\
\left.+\sqrt{\epsilon\theta_{u}}[dW_{u}(\Delta_{u}+a_{u},t)-dW_{u}(\Delta_{u}-a_{u},t)]\right)\end{split}$
(37b) $\displaystyle\begin{split}\tau
d\Delta_{v}(t)=\frac{1}{\alpha_{v}}\left((\Delta_{v}-\Delta_{u})\cdot\left[w_{ie}(a_{u}+a_{v})-w_{ie}(a_{u}-a_{v})\right]dt\right.\\\
\left.+\sqrt{\epsilon\theta_{v}}[dW_{v}(\Delta_{v}+a_{v},t)-dW_{v}(\Delta_{v}-a_{v},t)]\right).\end{split}$

Note, $w_{ei}$ and $w_{ie}$ are even and monotonic in distance, so we define
the positive quantaties

$\displaystyle M_{u}$
$\displaystyle=\alpha_{u}^{-1}\cdot\left[w_{ei}(a_{u}-a_{v})-w_{ei}(a_{u}+a_{v})\right]>0,$
$\displaystyle M_{v}$
$\displaystyle=\tau^{-1}\alpha_{v}^{-1}\cdot\left[w_{ie}(a_{u}-a_{v})-w_{ie}(a_{u}+a_{v})\right]>0.$

Then we have the matrix-vector representation of the multivariate OU process,
$d\boldsymbol{\Delta}=\boldsymbol{\mathcal{K}}\boldsymbol{\Delta}dt+\boldsymbol{d\mathcal{W}}$,
where $\boldsymbol{\Delta}=\begin{pmatrix}\Delta_{u}\\\
\Delta_{v}\end{pmatrix}$,
$\boldsymbol{\mathcal{K}}=\begin{pmatrix}M_{u}&-M_{u}\\\
M_{v}&-M_{v}\end{pmatrix}$ is the coupling matrix, and

$\boldsymbol{d\mathcal{W}}=\begin{pmatrix}\frac{\sqrt{\epsilon\theta_{u}}}{2\alpha_{u}}[dW_{u}(\Delta_{u}+a_{u},t)-dW_{u}(\Delta_{u}-a_{u},t)]\\\
\frac{\sqrt{\epsilon\theta_{v}}}{2\tau\alpha_{v}}[dW_{v}(\Delta_{v}+a_{v},t)-dW_{v}(\Delta_{v}-a_{v},t)]\end{pmatrix}$

is the correlated Wiener process noise. Following methods for solving linear
SDEs [23], we can determine the mean and variance of $\Delta_{u}$ and
$\Delta_{v}$ to obtain estimates of the diffusion of each population. First we
diagonalize $\boldsymbol{\mathcal{K}}$:

$\boldsymbol{\mathcal{K}}=\frac{1}{M_{v}-M_{u}}\begin{pmatrix}1&M_{u}\\\
1&M_{v}\end{pmatrix}\begin{pmatrix}0&0\\\
0&-(M_{v}-M_{u})\end{pmatrix}\begin{pmatrix}-M_{v}&M_{u}\\\ -1&1\end{pmatrix}$

Thus our eigenvalues are $\lambda_{1,2}=0,-(M_{v}-M_{u})$ and the eigenvectors
are $\mathbf{v}_{1}=(1,1)^{T}$ and $\mathbf{v}_{2}=(M_{u},M_{v})^{T}$
respectively. Note, this implies that there is a marginally stable direction
along perturbations of the bump that move both the E and I bump the same
amount, and there is an attractive (stable) direction for perturbations that
move the E and I bump differently. Similar results have been found for coupled
lateral I layers for which each layer individually supports a self-sustaining
bump in the absence of cross-population coupling [20, 28, 6]. The mean is
given by
$\langle\boldsymbol{\Delta}\rangle=e^{\boldsymbol{\mathcal{K}}t}\boldsymbol{\Delta}(0)$,
and so

$\displaystyle\begin{pmatrix}\langle\Delta_{u}\rangle\\\
\langle\Delta_{v}\rangle\end{pmatrix}=\frac{1}{M_{v}-M_{u}}\begin{pmatrix}M_{v}\Delta_{v}(0)-M_{u}\Delta_{u}(0)-M_{u}[\Delta_{u}(0)-\Delta_{v}(0)]e^{-(M_{v}-M_{u})t}\\\
M_{v}\Delta_{v}(0)-M_{u}\Delta_{u}(0)+M_{v}[\Delta_{v}(0)-\Delta_{u}(0)]e^{-(M_{v}-M_{u})t}\end{pmatrix}$

Next we seek the covariance
$\langle\boldsymbol{\Delta}(t)\boldsymbol{\Delta}^{T}(t)\rangle=\int^{t}_{0}e^{\boldsymbol{\mathcal{K}}(t-s)}\boldsymbol{D}e^{\boldsymbol{\mathcal{K}^{T}}(t-s)}ds$
where the covariance matrix of the noise term, $\boldsymbol{D}$ is found to
be:

(38) $\displaystyle\boldsymbol{D}=\begin{pmatrix}D_{u}&D_{c}\\\
D_{c}&D_{v}\end{pmatrix}$

with

$\displaystyle
D_{u}=\frac{\epsilon\theta_{u}}{2\alpha_{u}^{2}}[C_{u}(0)-C_{u}(2a_{u})]$
$\displaystyle
D_{v}=\frac{\epsilon\theta_{v}}{2\tau^{2}\alpha_{v}^{2}}[C_{v}(0)-C_{v}(2a_{v})]$
$\displaystyle\begin{split}D_{c}=\frac{\epsilon\sqrt{\theta_{u}\theta_{v}}}{4\tau\alpha_{u}\alpha_{v}}[C_{c}(\Delta_{u}-\Delta_{v}+a_{u}-a_{v})-C_{c}(\Delta_{u}-\Delta_{v}+a_{u}+a_{v})\\\
-C_{c}(\Delta_{u}-\Delta_{v}-a_{u}-a_{v})+C_{c}(\Delta_{u}-\Delta_{v}-a_{u}+a_{v})].\end{split}$

Multiplying out
$e^{\boldsymbol{\mathcal{K}}(t-s)}\boldsymbol{D}e^{\boldsymbol{\mathcal{K}^{T}}(t-s)}$
and then integrating yields our predictions of E and I bump center of mass
variance:

(39a)
$\displaystyle\begin{split}\langle\Delta_{u}(t)^{2}\rangle&=\frac{D_{v}M_{u}^{2}-2D_{c}M_{u}M_{v}+D_{u}M_{v}^{2}}{(M_{v}-M_{u})^{2}}t\\\
&-2\frac{e^{-(M_{v}-M_{u})t}-1}{(M_{v}-M_{u})^{3}}[D_{c}M_{u}^{2}-D_{v}M_{u}^{2}+D_{c}M_{u}M_{v}-D_{u}M_{u}M_{v}]\\\
&-M_{u}^{2}\frac{e^{-2(M_{v}-M_{u})t}-1}{2(M_{v}-M_{u})^{3}}[D_{u}+D_{v}-2D_{c}],\end{split}$
(39b)
$\displaystyle\begin{split}\langle\Delta_{v}(t)^{2}\rangle&=\frac{D_{v}M_{u}^{2}-2D_{c}M_{u}M_{v}+D_{u}M_{v}^{2}}{(M_{v}-M_{u})^{2}}t\\\
&-2\frac{e^{-(M_{v}-M_{u})t}-1}{(M_{v}-M_{u})^{3}}[D_{c}M_{v}^{2}-D_{u}M_{v}^{2}+D_{c}M_{u}M_{v}-D_{v}M_{u}M_{v}]\\\
&-M_{v}^{2}\frac{e^{-2(M_{v}-M_{u})t}-1}{2(M_{v}-M_{u})^{3}}[D_{u}+D_{v}-2D_{c}].\end{split}$

In the limit as $t\rightarrow\infty$ we find that both variances are dominated
by the term:

$\displaystyle\langle\Delta_{u}(t)^{2}\rangle=\langle\Delta_{v}(t)^{2}\rangle=\frac{D_{v}M_{u}^{2}-2D_{c}M_{u}M_{v}+D_{u}M_{v}^{2}}{(M_{v}-M_{u})^{2}}t,$

which is essentially an estimate of the variance derived from assuming the E
and I bump are co-located as in the Strongly coupled limit approximation.
Relating these two routes to one another, we find

$\mathcal{B}=\frac{w_{ei}(a_{u}-a_{v})-w_{ei}(a_{u}+a_{v})}{w_{ie}(a_{u}-a_{v})-w_{ie}(a_{u}+a_{v})}=\frac{\alpha_{u}M_{u}}{\tau\alpha_{v}M_{v}}.$

Using this relation and plugging in the expressions for $D_{u}$, $D_{v}$, and
$D_{c}$ we find that we obtain the diffusion coefficient expression from the
Strongly coupled limit prediction, Eq. Eq. 31. At short times, the Interface-
based approximation, Eq. Eq. 39, has contributions based on the interactions
of the E and I bumps, which decay over time.

The main difference between the Strongly coupled limit and Interface-based
methods is that the bumps are allowed to drift apart in the Interface-based
method, which allows us to separately estimate the I bump’s variance. The most
general form of the approximation can further describe stochastic widening and
contraction of bump, described by a full four-dimensional and nonlinear
approximation. Even collapsing to centers of mass approximations, using the
Interface-based interactions as a starting point is more accurate than the
Strongly coupled limit, as we will subsequently show. On the other hand, the
Strongly coupled limit is more straightforward to obtain, and there are fewer
arbitrary assumptions we must make (e.g., the static gradient approximation).
Yet, the Interface-based method allows us to obtain greater precision by using
the fully nonlinear approximation over all the interfaces, $x_{j}$, making it
a better metric for variance predictions.

Figure 4: Variance predictions and simulations. Numerical simulations of Eq.
Eq. 1 were run using an approximate version of the interface equations derived
from Eq. Eq. 33. Euler-Maruyama was used for time-stepping with noise
amplitude $\epsilon=0.001$, the spatial interval was truncated to
$[-3\pi,3\pi]$ with steps $dx=\frac{3\pi}{1000}$, timesteps are $dt=1$ms, and
variances were calculated by marginalizing over $10^{4}$ realizations per
point. (A) The Strongly coupled limit prediction and corresponding simulations
over 1 second. Although this prediction works reasonably well for small
thresholds $\theta$, it breaks down for higher thresholds where the I bump
center of mass differs considerably, e.g., $\theta_{u}+0.05=\theta_{v}=0.45$.
Note the kink and change in the variance trend when $\theta$ passes through a
value at which the half-widths $a_{u}$, $a_{v}$ exchange order (gold circle).
Insets below show single E and I bump simulations at indicated threshold
values. (B) When comparing to estimates of variance made using the Interface-
based approach, theory more closely tracks the simulation results at higher
firing rate threshold $\theta$. (C) The Strongly coupled limit predicts pure
diffusion and linearly scaling variance, which underestimates variance
calculated from simulations at $\theta_{u}+0.05=\theta_{v}=0.45$. (D) The
Interface-based estimate tracks the drifting apart of the E and I bump,
leading to more accurate variance predictions when
$\theta_{u}+0.05=\theta_{v}=0.45$. All other parameters are as in Table 1.

###  4.3 Variance Predictions and Simulations Comparison

To validate our Strongly coupled limit and Interface-based predictions of bump
variance, we ran stochastic simulations of Eq. Eq. 1 using spatiotemporal
noises to the E and I populations that are not correlated between populations;
$C_{c}(x)\equiv 0$. As firing rate thresholds are increased (Fig. 4A-B), both
our theoretical predictions and the averaged numerical simulations suggest
that variance changes non-monotonically. This stands in stark contrast to
results from previous studies, which found that the effective diffusion of
bumps generally tends to increase monotonically with the firing threshold for
single-population lateral I networks [31]. Interestingly, we find that the
variance peaks at the precise point in parameter space where the noise-free E
and I bump have the same width (gold circles, Fig. 4A,B).

Moreover, we see that the Strongly coupled limit approximation adequately
captures the effective motion of the E bump, but not the I bump, which tends
to stray further (Fig. 4A). In contrast, the Interface-based approximation is
able to capture both the distinct I and E bumps’ variance by accounting for
the stochastic dynamics of the I bump being perturbed away from the E bump
center of mass (Fig. 4B).

The dynamics of the I bump vary with firing rate thresholds and timescales. At
higher firing rate thresholds, the I bump wanders more than the E bump, and
the two bumps tend to be more weakly coupled. As the timescale increases, the
I bump center of mass relaxes to wander slightly away from that of the E bump,
while the E bump variance scales linearly in time like pure diffusion (Fig.
4C,D). Since the I bump is sustained by input from the E population, the I
bump appears to weakly track the input(see Fig. A.1 for single simulation of
such activity), which can be estimated by OU processes [31, 6]. Hence, the I
bump variance can be higher than that of the E population.

To further analyze how the variance in bump position depends on changes in
parameters, we determined how it changes along several different parameter
axes (Fig. 5). For simplicity, we studied the variations in the E bump’s
position variance and compared to the Interface-based approximation. We
started by varying the amplitude of inter-population connectivity, either
$E\to I$ or $I\to E$. Weakening either of these projection amplitudes tended
to increase bump variance (Fig. 5A,B). We speculate that weakening cross-
population connectivity leads to less input between the bumps. Thus, the bumps
are more weakly stable to noise perturbations and re-equilibrate more slowly,
ultimately leading to more wandering.

Figure 5: Excitatory bump position variances as a function of network
connectivity, spatial extent, and firing rate threshold. Noise amplitude
$\epsilon=0.001$ throughout. (A) Bump position variance mostly decreases as
connectivity amplitude $A_{ie}$ is increased. Other parameters are
$A_{ii}=0,\tau=1$, and $\theta_{u}=\theta_{v}=\theta$. (B) Bump position
variance mostly decreases as $A_{ei}$ is increased. Other parameters are
$A_{ii}=0,\tau=1$. $\theta_{u}=\theta_{v}=\theta$. (C) Bump position variance
changes non-monotonically as the spatial extent of the I projections
($\sigma_{ei}=\sigma_{ii}$) are increased. Other parameters are
$A_{ii}=0,\theta_{u}=\theta_{v}=0.3$ and $\tau=1$. $\sigma_{ie}$ is varied
vary. (D) Bump positive variance primarily increases with I population
threshold $\theta_{v}$. Other parameters are $A_{ii}=0$ and $\tau=1$. We also
vary firing threshold $\theta_{u}$. Other parameters are as in Table 1.

Aside from non-monotonicity arising as a function of the firing rate
threshold, we also observe that narrowing the I bump synaptic profile (varying
the spatial extent of the I projections) yields a peak in variance (Fig. 5C)
which may be due to a greater susceptibility to noise perturbations when the
width of I to E projections $\sigma_{ei}=\sigma_{ii}$ is decreased, resulting
in a peak of variance. For smaller $\sigma_{ei}=\sigma_{ii}$, we speculate
that the equilibrium E bump width is narrower and can be more easily
stabilized by I feedback. Note that the simulations do not match the peaks as
well as in other panels likely due to approximations made in the Interface-
based approach, though the trends are still largely captured.

We also observe that lower I firing rate thresholds lead to bumps that are
more stabilized to noise perturbations (Fig. 5D). However, the largest peaks
in bump position variance occur where E and I bump half-widths are most
similar, i.e., where E (I) threshold crossing gradients are the lowest
(highest) (Fig. A.2).

These results assumed I to I connections are non-existent $A_{ii}=0$. However,
we found that even low amounts of connectivity, $I\rightarrow I$ connectivity
$A_{ii}=0.01$, decreased bump variance (Fig. A.3).

Figure 6: Effects of inter-population noise correlations. (A) Bump position
variance as a function of firing rate threshold $\theta_{u}=\theta_{v}=\theta$
for a network with completely uncorrelated ($c=0$) or completely correlated
($c=1$) inter-population noise correlations. Interface-based theory (solid)
agrees well with numerical simulations. (B) Bump position variance for fixed
firing rate threshold $\theta_{u}=\theta_{v}=0.45$ as a function of inter-
population noise correlations $c$. (C) Bump profile evolution in absence of
noise for firing rate threshold $\theta_{u}=\theta_{v}=0.25$ given a
correlated center shift perturbation at $t=0$. (D) Bump profile evolution in
absence of noise over 1 second for firing rate threshold
$\theta_{u}=\theta_{v}=0.25$ given an uncorrelated center shift perturbation
at $t=0$. (E) Bump profile evolution for firing rather threshold
$\theta_{u}=\theta_{v}=0.25$ given two correlated ‘kicks’ (0.1 amplitude
Gaussian bumps shifted by random amount) applied to each population. (F) Bump
profile evolution for firing rate threshold $\theta_{u}=\theta_{v}=0.25$ with
two uncorrelated ‘kicks’ (0.1 amplitude Gaussian bumps shifted by random
amount) applied to each population. Other parameters are $A_{ii}=0$,
$\epsilon=0.001$ and otherwise as in Table 1.

###  4.4 Correlated versus uncorrelated noise

The impact of noise correlations on neural circuit codes for delayed estimates
are varied. Correlated noise within a neural subpopulation can improve working
memory coding [38], but cross-population noise correlations can lead to an
increase in bump attractor wandering that degrades memory [28]. We find in the
subsequent investigation that cross-population noise correlations between E
and I populations lead to less bump wandering than uncorrelated noise.

Our aim is to explore the system when there is cross-population noise
(corresponding to the case where $C_{c}\neq 0$). To obtain greater control
over the extent of correlated noise we opt to express our spatiotemporal noise
sources in the E and I populations as

(40) $dW_{u,v}(x,t)\rightarrow\sqrt{1-c^{2}}dW_{u,v}(x,t)+cdW_{c}(x,t),$

where the independent spatially correlated and temporally white noise process
$dW_{c}$ represents a correlated stochastic component in system Eq. Eq. 1. We
define these three noise terms the same as before though now we have nonzero
cross population spatial correlation $C_{c}(x-y)t=\langle
W_{c}(x,t)W_{c}(y,t)\rangle$. The correlation parameter $c\in[0,1]$ such that
$c=0$ implies uncorrelated noise, and $c=1$ fully correlated noise.
Correlating noise across the E and I populations drastically alters to our
variance predictions (Fig. 6A and Fig. A.4), with increased correlated noise
decreasing the predicted and simulated bump variances (Fig. 6B).

We analyze the response of the E/I bump structure to correlated as opposed to
uncorrelated shifts. Correlated shifts translated both bumps centers, due to
translational invariance of the system (Fig. 6C). Uncorrelated shifts move the
E and I bumps in opposite directions, but then show an attraction of the I
bump to the E bump, while the E bump is repulsed by the I bump leading the I
bump to ‘catch up’ to the shifted E bump and the E bump moves further away
from its original position (Fig. 6D). Thus, uncorrelated shifts lead to
additional drift of the bump and higher variance. This stands in stark
contrast to the uncorrelated/correlated perturbation analysis carried out for
two coupled identical lateral I layers [28], in which opposing perturbations
of two weakly connected bumps are effectively canceled by the attractive force
between the bumps.

We also considered the effects of two small additive Gaussian inputs, more
akin to the noisy kicks arising in stochastic simulations. When the position
of these kicks was strongly correlated (Fig. 6E; i.e. the same application to
each population), the E and I bumps stayed together and relaxed back to their
original position. Uncorrelated kicks led to different perturbations in the
position and profile of the E and I bumps (e.g., Fig. 6F), with a relaxation
where the I bump was attracted the E bump, but the E bump was repelled by the
I bump.

Overall, we found that the separate E/I network model shows parametrically-
dependent bump wandering. First, bump position variance depends strongly on
the relative widths of the E and I bumps, and relaxation dynamics from noise
perturbations can cause the I bump to stray from the E bump. Second, non-
monotonic variance trends arise with respect to network connectivity amplitude
and spatial scale parameters. Lastly, increases in inter-population noise
correlations reduce bump wandering by eliminating the relaxation effects that
would otherwise extend the effects of stochastic perturbations on the E and I
bumps.

##  5 Discussion

Stochastic bump attractor models have become a useful tool for characterizing
variability of systems that code for memory-based estimates of continuous
variables [8, 48, 3, 34]. As we show here, it is important to appreciate the
nuances of E/I architecture and noise-correlations [24] when making
predictions about continuous variable estimates. Our study has provided a
suite of new predictions concerning how inter-population connectivity
amplitude and spatial profiles impact bump wandering. To obtain tractable
expressions for bump position variance predictions, one form of our
approximations relied upon the marginal stability of solutions to the noise-
free system. Prior to performing our variance estimates, we observed that
along several parameter axes, we obtain non-monotonic changes in half-widths
and two types of instabilities: an oscillatory (Hopf) instability located in
the red unstable regions, and a contraction instability located at semistable
discontinuous saddle nodes (Fig. 3). Partial versions of the results have been
observed previously in E/I population models [4]. We also found that stability
of solutions was significantly affected by nonzero I$\rightarrow$I synaptic
strength (compare Fig. 3 and insets). Thus, even ‘weak’ I$\rightarrow$I
connectivity can strongly affect the linearized dynamics of stationary bump
solutions.

Our asymptotic analysis aimed at predicting the stochastic motion of bumps
subject to noise moved beyond the standard single center-of-mass approximation
often used to estimate the stochastic motion of bumps [37, 8, 31]. While a
single center-of-mass approximation works reasonably well across some parts of
parameter space, at high firing rate threshold, close to the discontinuous
saddle-node bifurcation, we find this breaks down and is better captured by a
stochastic interface based approximation previously developed in [36]. In
particular, the I bump diffuses more than the E bump, which is well captured
by the nonlinear Langevin approximation derived by approximating the motion of
bump interfaces. The amplitude of bump wandering changes non-monotonically in
most network parameters, due to the change in bump half-width amplitudes, well
captured by our Interface-based approximation. Notably, inter-population noise
correlations reduced bump wandering. Uncoordinated E and I bump motion,
arising in networks with uncorrelated inter-population noise, lead to
additional bump drift during relaxation periods while the I bump ‘chased’ the
E bump.

There are several possible extensions of our analysis of bump stochastic
motion in the E/I network. It is important to note that we made multiple
approximations to collapse our stochastically evolving bump interface
equations to a pair of SDEs describing the coupling between the E and I bump.
Alternatively, we could have retained a higher order approximation of the bump
interface gradients in order to obtain a more accurate approximation [18].
Recall, the Interface Based approximation begins as a fully nonlinear
description and can be used to describe dynamics near the oscillatory (Hopf)
bifurcation or the discontinuous saddle node. Alternatively, such near-
bifurcation approximations could also be determined by choosing a scaling for
a bifurcation parameter similar to the weak noise amplitude as in [33, 29].

Our analysis of bump position variance across multiple parametric axes helped
identify a number of ways to reduce bump wandering via network architectural
tuning. We largely chose parameters roughly assuming $80\%$ E and $20\%$ I
neurons present in the prefrontal cortex [1], but we could certainly explore
broader ranges of parameter space beyond this typical fraction. Another
natural extension for this work would be to consider more complex mechanisms
for synaptic tuning, such as short term plasticity in the E and I populations,
to reduce the propensity for bumps to wander. Recent studies in mean field
reductions of spiking networks have demonstrated that short term facilitation
(depression) on the E population with global inhibition tends to decrease
(increase) bump drift and diffusion [44], ultimately improving parametric
working memory. Extensions of our model and present analysis could be used to
further investigate how the introduction of different forms of short term
plasticity into either the E or I population would impact bump wandering. Not
only could short term plasticity reduce the effect of stochastic perturbations
on bumps, but also make them more robust to distraction inputs [40].

Finally, our investigation into the role of cross-population noise
correlations raises questions regarding the coding advantages brought about by
noise correlations. Noise correlations can increase, decrease, or not affect
the amount of information encoded by a neural circuit [2]. Most such results
have been derived in network models devoid of spatial structure. However,
recently work has demonstrated how disruptive broadly correlated
spatiotemporal noise can be to information transmission in spatially-organized
neural circuits [26]. Our work adds to this ongoing line of inquiry by
demonstrating variability-reducing mechanisms possible via increased
correlation in noise between E and I populations. Our bump position variance
predictions and model are sensitive to changes in the structure of noise. An
improved understanding of the precise form and structure of noise in
prefrontal cortex and other areas [41] could help further constrain neural
circuit models of memory-encoding persistent activity. Mechanistic models that
connect synaptic architecture, psychophysical performance, stochastic and
spatiotemporal dynamics, as well as the rich structure of internal and
external noise can help us further understand the dynamical principles
underlying information coding in the brain.

## Acknowledgments

We thank Sage Shaw for assistance with python code development.

## References

  * [1] M. Abeles, Corticonics: Neural circuits of the cerebral cortex, Cambridge University Press, 1991.
  * [2] B. B. Averbeck, P. E. Latham, and A. Pouget, Neural correlations, population coding and computation, Nature reviews neuroscience, 7 (2006), pp. 358–366.
  * [3] J. Barbosa, H. Stein, R. L. Martinez, A. Galan-Gadea, S. Li, J. Dalmau, K. Adam, J. Valls-Solé, C. Constantinidis, and A. Compte, Interplay between persistent activity and activity-silent dynamics in the prefrontal cortex underlies serial biases in working memory, Nature neuroscience, 23 (2020), pp. 1016–1024.
  * [4] P. Blomquist, J. Wyller, and G. T. Einevoll, Localized activity patterns in two-population neuronal networks, Physica D, 206 (2005), pp. 180–212, https://doi.org/DOI:lookup.
  * [5] P. C. Bressloff, Stochastic neural field theory and the system-size expansion, SIAM Journal on Applied Mathematics, 70 (2010), pp. 1488–1521.
  * [6] P. C. Bressloff and Z. P. Kilpatrick, Nonlinear langevin equations for wandering patterns in stochastic neural fields, SIAM Journal on Applied Dynamical Systems, 14 (2015), pp. 305–334.
  * [7] G. M. Cicchini, K. Mikellidou, and D. C. Burr, The functional role of serial dependence, Proceedings of the Royal Society B, 285 (2018), p. 20181722.
  * [8] A. Compte, N. Brunel, P. S. Goldman-Rakic, and X.-J. Wang, Synaptic mechanisms and network dynamics underlying spatial working memory in a cortical network model, Cerebral Cortex, 10 (2000), pp. 910–923, https://doi.org/DOI:lookup.
  * [9] C. Constantinidis and T. Klingberg, The neuroscience of working memory capacity and training, Nature Reviews, 206 (2016), pp. 1–12, https://doi.org/DOI:lookup.
  * [10] S. Coombes, P. beim Graben, R. Potthast, and J. Wright, Neural Fields: Theory and Applications, Springer, Berlin Heidelberg, 2014.
  * [11] S. Coombes and C. Laing, Pulsating fronts in periodically modulated neural field models, Physical Review E, 83 (2011), p. 011912.
  * [12] S. Coombes, G. J. Lord, and M. R. Owen, Waves and bumps in neuronal networks with axo-dendritic synaptic interactions, Physica D: Nonlinear Phenomena, 178 (2003), pp. 219–241.
  * [13] S. Coombes and M. R. Owen, Evans functions for integral neural field equations with heaviside firing rate function, SIAM Journal on Applied Dynamical Systems, 3 (2004), pp. 574–600.
  * [14] S. Coombes, H. Schmidt, and I. Bojak, Interface dynamics in planar neural field models, Journal of Mathematical Neuroscience, 2 (2012), pp. 1–27, https://doi.org/DOI:10.1186/2190-8567-2-9.
  * [15] S. Coombes, R. Thul, and K. Wedgwood, Nonsmooth dynamics in spiking neuron models, Physica D, 241 (2012), pp. 2042–2057, https://doi.org/DOI:10.1016/j.physd.2011.05.012.
  * [16] C. E. Curtis, Prefrontal and parietal contributions to spatial working memory, Neuroscience, 139 (2006), pp. 173–180, https://doi.org/DOI:10.1016/j.neuroscience.2005.04.070.
  * [17] A. A. Faisal, L. P. Selen, and D. M. Wolpert, Noise in the nervous system, Nature reviews neuroscience, 9 (2008), pp. 292–303.
  * [18] G. Faye and Z. P. Kilpatrick, Threshold of front propagation in neural fields: An interface dynamics approach, SIAM Journal on Applied Mathematics, 78 (2018), pp. 2575–2596.
  * [19] S. E. Folias and P. C. Bressloff, Breathing pulses in an excitatory neural network, SIAM Journal on Applied Dynamical Systems, 3 (2004), pp. 378–407.
  * [20] S. E. Folias and G. B. Ermentrout, New patterns of activity in a pair of interacting excitatory-inhibitory neural fields, Phys. Rev. Lett., 107 (2011), pp. 1–4.
  * [21] S. Funahashi, C. Bruce, and P. Goldman-Rakic, Mnemonic coding of visual space in the monkey’s dorsolateral prefrontal cortex,, J Neurophysiol, 61 (1989), pp. 331–349, https://doi.org/DOI:.
  * [22] J. M. Fuster, Unit activity in prefrontal cortex during delayed-response performance: Neuronal correlates of transient memory, J Neurophysiol, 36 (1973), pp. 61–78, https://doi.org/DOI:.
  * [23] C. W. Gardiner, Handbook of stochastic methods: for physics, chemistry, and the natural sciences, Springer complexity, Berlin Heidelberg, 2009\.
  * [24] P. S. Goldman-Rakic, Cellular basis of working memory, Neuron, 14 (1995), pp. 477–485, https://doi.org/DOI:lookup.
  * [25] M. Häusser and A. Roth, Estimating the time course of the excitatory synaptic conductance in neocortical pyramidal cells using a novel voltage jump method, Journal of Neuroscience, 17 (1997), pp. 7606–7625.
  * [26] C. Huang, A. Pouget, and B. Doiron, Internally generated population activity in cortical networks hinders information transmission, bioRxiv, (2020).
  * [27] S. ichi Amari, Dynamics of pattern formation in lateral-inhibition type neural fields, Biological Cybernetics, 27 (1977), pp. 77–87, https://doi.org/DOI:lookup.
  * [28] Z. P. Kilpatrick, Interareal coupling reduces encoding variability in multi-area models of spatial working memory, Frontiers in Computational Neuroscience, 7 (2013), pp. 1–14, https://doi.org/DOI:lookup.
  * [29] Z. P. Kilpatrick, Ghosts of bump attractors in stochastic neural fields: Bottlenecks and extinction, Discrete & Continuous Dynamical Systems-B, 21 (2016), p. 2211.
  * [30] Z. P. Kilpatrick, Synaptic mechanisms of interference in working memory, Nature: Scientific Reports, 8 (2018), pp. 1–20, https://doi.org/DOI:lookup.
  * [31] Z. P. Kilpatrick and B. Ermentrout, Wandering bumps in stochastic neural fields, SIAM J. Applied Dynamical Systems, 12 (2013), pp. 61–94, https://doi.org/DOI:lookup.
  * [32] Z. P. Kilpatrick, B. Ermentrout, and B. Doiron, Optimizing working memory with heterogeneity of recurrent cortical excitation, Journal of neuroscience, 33 (2013), pp. 18999–19011.
  * [33] Z. P. Kilpatrick and G. Faye, Pulse bifurcations in stochastic neural fields, SIAM Journal on Applied Dynamical Systems, 13 (2014), pp. 830–860.
  * [34] S. S. Kim, H. Rouault, S. Druckmann, and V. Jayaraman, Ring attractor dynamics in the drosophila central brain, Science, 356 (2017), pp. 849–853.
  * [35] A. Kiyonaga, J. M. Scimeca, D. P. Bliss, and D. Whitney, Serial dependence across perception, attention, and memory, Trends in Cognitive Sciences, 21 (2017), pp. 493–497.
  * [36] N. Krishnan, D. B. Poll, and Z. P. Kilpatrick, Synaptic efficacy shapes resource limitations in working memory, Journal of Computational Neuroscience, (2018), pp. 1–23, https://doi.org/DOI:lookup.
  * [37] C. R. Laing and C. C. Chow, Stationary bumps in networks of spiking neurons, Neural Computation, 13 (2001), pp. 1473–1494.
  * [38] M. L. Leavitt, F. Pieper, A. J. Sachs, and J. C. Martinez-Trujillo, Correlated variability modifies working memory fidelity in primate prefrontal neuronal ensembles, Proceedings of the National Academy of Sciences, 114 (2017), pp. E2494–E2503.
  * [39] R. I. LEINE, D. H. V. CAMPEN, and B. L. V. D. VRANDE, Bifurcations in nonlinear discontinuous systems, Nonlinear Dynamics, 23 (2000), pp. 105–164.
  * [40] E. S. Lorenc, R. Mallett, and J. A. Lewis-Peacock, Distraction in visual working memory: Resistance is not futile, Trends in Cognitive Sciences, 25 (2021), pp. 228–239.
  * [41] M. D. McDonnell and L. M. Ward, The benefits of noise in neural systems: bridging theory and experiment, Nature Reviews Neuroscience, 12 (2011), pp. 415–425.
  * [42] D. J. Pinto and G. B. Ermentrout, Spatially structured activity in synaptically coupled neuronal networks: Ii. lateral inhibition and standing pulses, SIAM Journal on Applied Mathematics, 62 (2001), pp. 226–243.
  * [43] C. J. Ploner, B. Gaymard, S. Rivaud, Y. Agid, and C. Pierrot-Deseilligny, Temporal limits of spatial working memory in humans, European Journal of Neuroscience, 10 (1998), pp. 794–797.
  * [44] A. Seeholzer, M. Deger, and W. Gerstner, Stability of working memory in continuous attractor networks under the control of short- term plasticity, PLOS Computational Biology, 15 (2019), pp. 1–48, https://doi.org/DOI:lookup.
  * [45] R. Veltz and O. Faugeras, Local/global analysis of the stationary solutions of some neural field equations, SIAM Journal on Applied Dynamical Systems, 9 (2010), pp. 954–998.
  * [46] J. M. White, D. L. Sparks, and T. R. Stanford, Saccades to remembered target locations: an analysis of systematic and variable errors, Vision Res, 34 (1994), pp. 79–92.
  * [47] H. R. Wilson and J. D. Cowan, A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue, Biol. Cybernet., 13 (1973), pp. 55–80, https://doi.org/DOI:lookup.
  * [48] K. Wimmer, D. Q. Nykamp, C. Constantinidis, and A. Compte, Bump attractor dynamics in prefrontal cortex explains behavioral precision in spatial working memory, Nature Neuroscience, 17 (2014), pp. 431–439, https://doi.org/DOI:10.1038/nn.3645.

##  Appendix A Additional Figures

Figure A.1: Single simulation of E and I bump wandering. The noise amplitude
is $\epsilon=0.001$. The I bump wanders much further in certain parameter
regimes. Synaptic weight strength $A_{ii}=0$; and firing rate thresholds
$\theta_{u}=0.42$, and $\theta_{v}=0.47$. (A) Initial and final profiles of E
and I bumps. (B) E and I center of mass positions evolving over time. The I
bump generally follows E wandering, but also additionally wanders more wildly
about E. (C) Bumps wandering over time. Traces represent the center of mass
(green) and the E/I interfaces (red/blue, respectively). Colorscales represent
the profile amplitude of $u(x,t)$ and $v(x,t)$. Other parameters are as in
Table 1.  Figure A.2: Isolating factors affecting variance of bumps’
wandering. Synaptic strength $A_{ii}=0.01$. (A) Colorscale represents the
difference in relative E and I halfwidths for varying $\theta_{u}$ and
$\theta_{v}$. Bump position variances are maximized when $a_{u}=a_{v}$. (B)
Colorscale of E bump position variance for $\epsilon=0.001$ amplitude noise.
Black line bounds unstable region. (C,D) E,I bump threshold crossing gradient
($|U^{\prime}(a_{u})|$ and $|V^{\prime}(a_{v})|$) as a function of firing rate
thresholds. Note, the E variance is largest when the threshold crossing
gradient for E (I) is low (high). Slices further taken through these plots
correspond to predictions shown in Fig. 5D.  Figure A.3: Predicted and
simulated center of mass variances. Here, we take I to I connectivity
$A_{ii}=0.01$ and noise amplitude $\epsilon=0.001$. (A) Bump position variance
as a function of $A_{ei}=A$ vary and $\theta_{u}=\theta_{v}=\theta$ are
varied. (B) Variance as a function of $A_{ie}=A$ as
$\theta_{u}=\theta_{v}=\theta$ are also varied. (C) Variance as a function of
spatial scale $\sigma_{ei}=\sigma_{ii}$ and varying $\sigma_{ie}$. We fix
$\theta_{u}=\theta_{v}=0.3$. (D) Variance as a function of $\theta_{v}$ as
$\theta_{u}$ is varied. Parameters not mentioned are as defined in Table 1.
Figure A.4: Predicted and simulated bump position variances for correlated and
uncorrelated noise. Noise amplitude is set to $\epsilon=0.001$ and we take
$A_{ii}=0$. (A)Bump position variance as a function of
$\theta_{u}=\theta_{v}=\theta$. (B) Bump position variance as a function
$A_{ei}=A$ with $\theta_{u}=\theta_{v}=0.25$. (C) $\theta_{u}=\theta_{v}=0.2$
and $\sigma_{ie}=3$. Bump position variance as a function of
$\sigma_{ei}=\sigma_{ii}$. (D) Bump position variance as a function of
$\theta_{v}$. We take $\theta_{u}=0.2$. Other parameters are as in Table 1.