# Self-organization of nonlinearly coupled neural fluctuations into
synergistic population codes

Hengyuan Ma1,$\ast$, Yang Qi1,2, , Pulin Gong4, Jie Zhang1,2 Wenlian
Lu1,2,$\dagger$, Jianfeng Feng1,2,3,
1 Institute of Science and Technology for Brain-inspired Intelligence, Fudan
University, Shanghai 200433, China
2 Key Laboratory of Computational Neuroscience and
Brain-Inspired Intelligence (Fudan University), Ministry of Education, China
3 Department of Computer Science, University of Warwick, Coventry, CV4 7AL, UK
4 School of Physics, University of Sydney, Sydney, NSW 2006, Australia
https://github.com/AwakerMhy/moment_nn/tree/main These authors contributed
equally.Correspondence should be addressed to<EMAIL_ADDRESS>
###### Abstract

Neural activity in the brain exhibits correlated fluctuations that may
strongly influence the properties of neural population coding. However, how
such correlated neural fluctuations may arise from the intrinsic neural
circuit dynamics and subsequently affect the computational properties of
neural population activity remains poorly understood. The main difficulty lies
in resolving the nonlinear coupling between correlated fluctuations with the
overall dynamics of the system. In this study, we investigate the emergence of
synergistic neural population codes from the intrinsic dynamics of correlated
neural fluctuations in a neural circuit model capturing realistic nonlinear
noise coupling of spiking neurons. We show that a rich repertoire of spatial
correlation patterns naturally emerges in a bump attractor network and further
reveals the dynamical regime under which the interplay between differential
and noise correlations leads to synergistic codes. Moreover, we find that
negative correlations may induce stable bound states between two bumps, a
phenomenon previously unobserved in firing rate models. These noise-induced
effects of bump attractors lead to a number of computational advantages
including enhanced working memory capacity and efficient spatiotemporal
multiplexing and can account for a range of cognitive and behavioral phenomena
related to working memory. This study offers a dynamical approach to
investigating realistic correlated neural fluctuations and insights to their
roles in cortical computations.

##  1 Introduction

The firing activity of cortical neurons in the brain exhibits large
fluctuations over time and across trials that are characterized by rich
spatiotemporal correlation structures as revealed by large-scale simultaneous
cortical recordings [52, 63]. It has been suggested that correlated neural
variability may play synergistic or destructive roles in neural probabilistic
representations such as neural population coding [30, 44] and working memory
[46, 11]. However, existing theoretical analyses of the stochastic dynamics of
neural systems primarily rely on firing rate models with additive gaussian
noise or simplified Poisson neuron models without correlations [50, 11, 46,
31, 68]. As a result, how correlated neural fluctuations may arise from the
intrinsic neural circuit dynamics and subsequently affect the computational
properties of neural population remains less well understood [28, 18]. The
main difficulty lies in resolving the nonlinear coupling between correlated
fluctuations with the overall dynamics of the system.

It is well known that simple structures in a recurrent circuitry can give rise
to many emergent properties to neural population activity including symmetry
breaking and spontaneous pattern formation [25, 14]. One such example is found
in bump attractors, a type of self-sustained, spatially localized patterns
arising from neural circuit models with short-range excitation and long-range
inhibition [3, 14]. Due to their translation symmetry, bump attractors can be
used to encode and store continuous features such as spatial location and
orientation [73], and have been suggested as a possible neural mechanism for a
range of perceptual and cognitive processes including spatial navigation [74,
40, 42], spatial working memory [13, 69, 71, 75], sensory evidence
accumulation [22], and other cortical functions [8, 56]. Modeling studies have
shown thatmultiple bump states can undergo complex interactions including
repulsion, annihilation,merging, and splitting [48, 33]. It has been proposed
that neural fluctuations may influence the coding properties of bump
attractors in two ways. At short timescales, noisy spike count directly
contributes to the coding error of a bump attractor, whereas at long
timescales, noise further induces a random drift of the bump attractor along
the attractor manifold, resulting in degradation of working memory over time
[11]. Thus, it is critical to understand how correlated fluctuations may arise
intrinsically from recurrent neural circuitry and how they affect the dynamics
of bump attractors and subsequently the relevant cognitive functions.

In this study, we investigate how synergistic neural population codes
originate from the intrinsic dynamics of correlated neural fluctuations by
using a biologically realistic model known as the moment neural network (MNN)
[23, 38], which faithfully captures the nonlinear noise coupling of spiking
neurons up to second-order statistical moments [51]. As suggested by
experimental and theoretical studies [53, 18], second-order statistical
moments can provide a minimalistic yet adequate way for capturing the
statistical properties of fluctuating neural activity.We construct a bump
attractor network based on this model and show that complex spatial
correlations naturally emerge from the intrinsic dynamics of recurrent
circuits. Analysis of the model reveals the interplay between differential and
noise correlations as a key factor determining the computational properties of
neural population codes in bump attractors. We also find that negatively
correlated neural fluctuations can induce stable bound states between adjacent
bumps, a phenomenon previously unobserved in firing rate models. These novel
dynamical properties of bump attractors can account for a diversity of
experimental observations of cognitive and behavioral responses related to
working memory [65, 60], and provide a circuit mechanism for efficient neural
coding based on spatiotemporal multiplexing [1, 12] as well as enhanced
working memory capacity.

This study offers new insights into how nonlinear noise coupling shapes
correlated neural fluctuations into synergistic neural population codes and
opens up a dynamical approach to investigating correlated neural fluctuations
and their computational properties in neural systems.

Figure 1:  Nonlinear noise coupling in a bump attractor network. (a)
Comparison of neural circuit models with and without nonlinear noise coupling.
In models without nonlinear noise coupling (left panel), the recurrent
dynamics affect the spike count mean (solid arrow) with spike count covariance
injected additively to the system (dashed arrow). In models with nonlinear
noise coupling (right panel), both spike count mean and spike count covariance
are coupled recurrently with each other through nonlinear interactions (purple
and green arrows). Here, $\mathbf{n}$ represents spike count. (b) A ring
network with short-range excitation (red curves) and long-range inhibition
(blue curves) supporting bump activity pattern. (c) Synaptic strengths $w$ as
a function of the distance $\Delta x$ between two neurons. (d) The mean firing
rate $\mu_{i}$ of a bump attractor forms a ring manifold (black line)
parameterized by its center ofmass, additionally equipped with a firing
covariability $C_{ij}$(orange ellipsoid; scaled for visibility) at each point
on the attractor manifold. The plot shows the attractor manifold in the space
spanned by the activity (spike count) of three neurons.

##  2 A dynamical model for correlated neural variability

A major challenge faced by firing rate models is difficulties in capturing
realistic neural fluctuations that are nonlinearly coupled to the recurrent
circuit dynamics. To illustrate this point, consider the schematic diagram
shown in Figure. 1(a). In simplified neural circuit models, neural fluctuation
is often modeled as either an external additive input or independent Poisson
variability. In a more realistic setting, correlated neural fluctuation may
arise intrinsically through the recurrent neural circuit dynamics, as
represented by the strong nonlinear coupling between the spike count mean and
spike count covariance [38, 28, 18]. Therefore, to understand the coding and
computational properties of bump attractors, it is imperative to develop a
model that can accurately describe the statistics and nonlinear coupling of
fluctuating neural activity. To overcome this challenge, we employ a model
known as the moment neural network (MNN), which faithfully captures the
nonlinear coupling of correlated fluctuations in a recurrent population of
spiking neurons [38]. The model is governed by the following closed system of
equations describing the dynamics of the second-order statistical moments of
neural activity,

$\displaystyle\tau\frac{\partial\mu_{i}}{\partial
t}=-\mu_{i}+\phi_{\mu}(\bar{\mu}_{i},\bar{\sigma}^{2}_{i}),$ (1)
$\displaystyle\tau\frac{\partial C_{ij}}{\partial
t}=-C_{ij}+\begin{cases}\phi_{\sigma}(\bar{\mu}_{i},\bar{\sigma}^{2}_{i}),&\text{for
}i=j,\\\
\psi(\bar{\mu}_{i},\bar{\sigma}^{2}_{i})\psi(\bar{\mu}_{j},\bar{\sigma}^{2}_{j})\bar{C}_{ij},&\text{for
}i\neq j,\end{cases}$ (2)

where the mean firing rate $\mu_{i}$ and firing covariability $C_{ij}$
represent the mean and covariance of the spike count per unit time,
respectively, and $\tau$ is the membrane time constant. The functions
$\phi_{\mu}$, $\phi_{\sigma}$ and $\psi$ are pointwise activation functions
describing the relationship between the input current statistics and the
output spike train statistics. For the leaky integrate-and fire (LIF) spiking
neuron model, the functional form of the activation functions are given by
equations (8) and (9) in Methods. Unlike previous models with heuristic
activation functions [3, 14, 15, 73], these activations are derived through a
mathematical technique known as the diffusion approximation [5, 4]which
faithfully captures the nonlinear coupling of mean firing rate and firing
variability across populations of neurons. The quantities $\bar{\mu}_{i}$ and
$\bar{C}_{ij}$ correspond to the mean and the covariance of the total synaptic
current, respectively, and are defined as

$\displaystyle\bar{\mu}_{i}=\sum_{k}w_{ik}\mu_{k}+\mu_{\mathrm{ext}}$ (3)
$\displaystyle\bar{C}_{ij}=\sum_{k,l}w_{ik}C_{kl}w_{jl}+\sigma_{\mathrm{ext}}^{2}\delta_{ij},$
(4)

where $w_{ij}$ is the synaptic weight from the $i$th neuron to the $j$th
neuron, $\bar{\sigma}^{2}_{i}=\bar{C}_{ii}$ denotes the current variability,
and $\mu_{\mathrm{ext}}$ and $\sigma^{2}_{\mathrm{ext}}$ are the mean and
variance of a spatially uniform external input, respectively. We fix
$\sigma^{2}_{\mathrm{ext}}$ throughout this study and systematically vary
$\mu_{\mathrm{ext}}$. See Methods for details of the parameter settings. We
note a couple of important features of this model. First, unlike models with
linearly coupled covariance as considered in [10, 62], equations (1) and (4)
capture the nonlinear interactions between mean firing rate and firing
covariability derived from realistic spiking neuron models. Second, any
covariance structure that emerges from the model is due to intrinsic dynamics
of the recurrent circuit, not determined by the external input, which is
uncorrelated.

To implement the bump attractor, each neuron $i$ is assigned to spatial
coordinates $x_{i}\in[-\pi,\pi)$ over a feature space with a periodic boundary
condition and is connected to each other via synaptic connections with short-
range excitation and long-range inhibition according to equation (10) in
Methods. A schematic diagram of this neural circuit structure is shown in
Figure 1b, and the synaptic connections are shown in Figure 1c.

##  3 Correlated variability in bump attractors

In firing rate-based models without correlation, distance-dependent synaptic
coupling such as that shown in Figure 1b can give rise to stable activity
state known as bump attractors [3]. Due to the translation invariance of the
synaptic coupling, the mean firing rate $\mu_{i}$of a bump attractor forms a
ring manifold parameterized by its center of mass. Experimentally, bump or
ring attractors have been found in various biological neural systems such as
the fruit fly ellipsoid body [29], primate V1 [52], and primate prefrontal
cortex [71]. However, the firing rate–based model is inadequate, since it
omits the fact that neurons in the ring attractors are correlated (e.g., Fig.7
in [71]). The MNN overcomes this shortcoming by integrating the ring manifold
with the intrinsic covariance structure of neural activity that is nonlinearly
coupled with the mean firing rate. As a result, the topology of the attractor
manifold in our model is no longer represented by points on a ring embedded in
the Euclidean space $\mathbb{R}^{n}$ spanned by the mean firing rates $\mu$,
but also extends to the positive semidefinite cone $S^{n}_{+}(\mathbb{R})$
spanned by the firing covariability $C$. This topological structure can be
visualized as a family of ellipsoids with their centers located on the ring
manifold and with their sizes and orientations representing the covariance
structure of neural activity (Figure 1d).

Figure 2: Correlated variability emerging from a bump attractor network. (a)
Mean firing rate $\mu(x)$, firing variability $\sigma^{2}(x)$ and the
correlation coefficients $\rho(x,{x}^{\prime})$ of neural activity states for
varying external input mean $\mu_{\mathrm{ext}}$, including spatially uniform
resting state (state I, $\mu_{\mathrm{ext}}=0.920$), bump state with positive
correlations (state II, $\mu_{\mathrm{ext}}=0.929$), bump state with both
positive and negative correlations (state III, $\mu_{\mathrm{ext}}=0.948$),
and globally active state with spatially periodic correlations (state IV,
$\mu_{\mathrm{ext}}=0.986$). (b) The height $h$ (top panel) and width $b$
(lower panel) for varying external input mean $\mu_{\mathrm{ext}}$, revealing
distinct phase transitions between different activity states shown in panel a.
(c) The correlation coefficient $\rho$ between two neurons with preferred
features $x_{1}$ and $x_{2}$ as the bump location $s$ varies in state II and
state III.

In the following, we reveal the emergence of intrinsic spatial correlations in
the bump attractor network through numerical simulations of the system defined
by equations (1)-(4). Similar to conventional firing rate–based models [3,
14], when the external input mean $\mu_{\mathrm{ext}}$ is weak, the system has
only one stable solution, the resting state (see state I in Figure 2a). As
$\mu_{\mathrm{ext}}$ increases, stable bump states emerge (see State II in
Figure 2a). Consistent with previous findings [3, 14], the mean firing rate
$\mu_{i}=\mu(x_{i}|s)$ of the bump state exhibits a unimodal spatial profile
peaked at its center $s$. In contrast, the firing variability
$\sigma^{2}_{i}=\sigma^{2}(x_{i}|s)$ is bimodal, with two peaks near the edges
of the bump and a local minimum at the center. This is in marked contrast to
the commonly assumed independent Poisson activity, in which the firing
variability is proportional to the mean firing rate or uncorrelated gaussian
noise with constant magnitude [39, 11]. We find that neural activity is
positively correlated within the bump and uncorrelated outside it. Note that
the stable bump state coexists with the stable resting state, allowing the
system to switch between them when subjected to appropriate external stimuli.

Furthermore, unlike firing rate–based models, additional phase transitions
occur in the correlation structure of the bump attractor. As
$\mu_{\mathrm{ext}}$further increases, negative correlations emerge between
neurons on the opposite sides of the bump, while the shape of the mean firing
rate remains qualitatively the same (see state III in Figure 2a). Note that
both the mean and variance of the external input are uniform, suggesting that
the spatial correlation structure is an emergent property of the system.
Similar to firing rate–based models, further increasing
$\mu_{\mathrm{ext}}$eventually leads to the saturation of mean firing rate
into a spatially uniform activity. However, we find that the same case is not
true for spatial correlations, which transform into a periodic wave pattern
(see state IV in Figure 2a). Such a distance-dependent correlation structure
may potentially arise from synchronized firing with the same and the opposite
phases.

To quantify the transition between these different activity patterns, we
calculate the height $h$ and the full width at half maximum $b$ of the bump
solution for varying the external input mean $\mu_{\mathrm{ext}}$. As shown in
Figure 2b, the phase diagrams constructed from these measurements reveal four
distinct regimes occupied by the neural activity patterns I-IV. We find that
both $h$ and $b$ vary smoothly as a function of $\mu_{\mathrm{ext}}$ within
each regime and that they exhibit discontinuous jumps at phase transitions.
The discontinuous jumps are qualitatively consistent with saddle-node
bifurcations found in firing rate bump attractor models [3, 14]. It is worth
noting that both the mean firing rate and the correlation coefficients
exhibited by our model are within a biologically plausible range, consistent
with experimental observations [71].

To further illustrate the internal correlation structure of the bump
attractor, we focus on two neurons and inspect how their mean firing rate and
correlation depend on the bump center $s$. Specifically, given two neurons
located at $x_{1}=0.2\pi$ and $x_{2}=0.4\pi$, their mean firing rates $\mu$
follow bell-shaped tuning functions peaked at $x_{1}$ and $x_{2}$
respectively. This holds for both states II and III, as shown in the top
panels in Figure 2(c). For each state, a distinct pattern is found in the
correlation coefficient between these two neurons. For state II, the
correlation coefficient $\rho(x_{1},x_{2})$ reaches the maximum when s is
located in the middle between $x_{1}$ and $x_{2}$, and monotonically decays to
zero when s move away from the middle point in both directions. For State III,
when $x_{1}<s<x_{2}$ (in-flank cases), we have $\rho(x_{1},x_{2})<0$; when
$s=x_{1}$ or $s=x_{2}$ (peak cases), we have $\rho(x_{1},x_{2})=0$; and when
$s<x_{1}$ or $s>x_{2}$ (out-flank cases), we have $\rho(x_{1},x_{2})>0$. We
find that the correlation pattern predicted by state III of our model is
consistent with experimentally observed pairwise correlation structures in the
primate prefrontal cortex [71]. This study analyzed neuronal recordings in
monkeys during a visuospatial delayed response task, in which they were
required to memorize the spatial locations of transient stimuli and found that
the tuning curves and pairwise correlation structures in the prefrontal cortex
corroborated the bump attractor hypothesis, that is, the neural population
held a bump activity that encoded the cue location $s\in[-\pi,\pi)$ by the
bump center during the delay period. As shown in Figure 7f of their paper, in
the out-flank, in-flank and peak cases, the sign of the correlations of two
neurons with two preferred stimuli $x_{1}$ and $x_{2}$ are the same as our
prediction in State III of Figure 2c. Additionally, our model predicts the
existence of a pairwise correlation structure as in state II of Figure 2c,
which has not been experimentally observed before.

Importantly, in both [71] and our model,the correlations of each neuronal pair
are dependent on the location of the bump center $s$, and this is more general
than previous models, where the correlations are assumed to be translation
invariant [58, 6, 20]. Although both the rate-based model [71]and theMNN
display a bell-shaped firing rate profile during the delay period of working
memory (compare Figure 8a in their paper with State III in Figure 2), the MNN
offers a number of unique advantages. First, in the rate-based model, the
spiking activity of each neuron is modeled as an inhomogeneous Poisson process
whose instantaneous firing rate is determined by the dynamics of the rate-
based model. In contrast, the MNN directly captures the dynamics of neural
spiking activity in terms of the nonlinear coupling between mean firing rate
and firing covariability, without imposing such simplifying abstractions.
Moreover, in the rate-based model, the diffusion of bump state along the
attractor manifold due to external noise is required to explain several
experimental observations such as negative noise correlation and stimulus-
dependent Fano factor. However, by considering the intrinsic covariance
structure in the MNN, many of these experimental observations can be explained
without invoking the diffusion of bump state. See the supplementary
information (Supplementary Information) for detailed comparisons between the
stochastic rate-based model and the MNN. These analyses demonstrate that by
considering the nonlinear noise coupling in recurrent neural circuits, the MNN
can account for the correlation structures of neural fluctuations in
experimental data and lay the theoretical foundation for exploring the role of
correlated neural fluctuations in brain functions such as neural coding and
working memory.

##  4 Noise covariance influences coding accuracy

Due to the translation symmetry of bump attractor, its center s can be used to
encode continuous variables. We analyze how the coding accuracy of bump
attractor is influenced by correlated activity state. The coding accuracy can
be quantified using the linear Fisher information$I_{\mathrm{linear}}(\Delta
t)=\mathbf{u}(s)^{T}C(s)^{-1}\mathbf{u}(s)\Delta t$, where
$\mathbf{u}_{i}(s)=\frac{d\mu_{i}(s)}{ds}$corresponds to the tangent vector
along the attractor manifold and $\Delta t$ is the spike count time window.
The linear Fisher information $I_{\mathrm{linear}}(\Delta t)$ is approximately
independent of $s$ for large population size due to translation invariance.
Linear Fisher information provides an upper bound on the certainty of the
bump’s location as can be determined by an optimal linear decoder from random
spike count [26, 9, 24]. Figure 3a shows the linear Fisher information rate
($I:=\frac{I_{\text{linear}}(\Delta t)}{\Delta t}$) of a bump attractor for
encoding $s$ under varying the external input mean $\mu_{\mathrm{ext}}$. We
observe that the phase transition from state II to state III caused by
increasing $\mu_{\mathrm{ext}}$ coincides with an abrupt decrease in $I$.
Since cognitive performance at the behavioral level is limited by the coding
accuracy in upstream neural populations [35], these two dynamical regimes of
bump activity state may correspond to different levels of cognitive
performance. In particular, our model predicts that variations in global
circuit parameters such as task-independent external input mean
$\mu_{\mathrm{ext}}$are sufficient to induce shifts in cognitive state, which
may manifest as a sudden deterioration of task performance.

Figure 3: Noise and differential covariance together influence coding
accuracy. (a) Linear Fisher information rate $I$ under varying external input
mean $\mu_{\mathrm{ext}}$. (b) Differential covariance and noise covariance of
state II ($\mu_{\mathrm{ext}}=0.925$) and state III
($\mu_{\mathrm{ext}}=0.931$). (c) Projecting the system state onto the space
spanned by the activity of two neurons whose preferred $s$ differs by
$0.2\pi$. The black line shows the ring manifold represented by the the mean
firing rates of two neurons. Orange ellipses represent the firing
covariability when the bump center is located at different values represented
by the red dots. Both state II (left panel) and state III (right panel) show
similar structures in mean firing rate, but not the firing covariability. The
major axis (blue lines) of the firing covariability in state II is roughly
orthogonal to the ring manifold, whereas that in state III is roughly parallel
to the ring manifold. (d) Linear Fisher information rate $I$ grows linearly
with the neural population size $N$.

To understand what causes the drastic difference in the coding accuracy
between bump states II and III, we turn to analyze the contribution of
correlated neural variability to the linear Fisher information. We find that
the linear Fisher information is jointly influenced by two types of firing
covariability, differential covariance and noise covariance [6]. Differential
covariance $\mathbf{u}(s)\mathbf{u}(s)^{T}$ measures the covariance due to
changes in the mean firing rate $\mu$ under infinitesimal shifts in $s$ [30],
whereas noise covariance $C(s)$ in our model measures the spike count
covariance per unit time between two neurons for fixed $s$. Figure 3b compares
the spatial structures of differential covariance and noise covariance in
states II and III. Both states exhibit similarly shaped differential
covariance with positively correlated neurons at the same side of the bump and
negatively correlated neurons at the different sides of the bump, consistent
with previous firing rate–based models of bump attractor [73]. However, our
model reveals distinct shapes in the noise covariance that could not be
captured by previous models. For state II, the noise covariance is
nonnegative, whereas for state III, the noise covariance contains both
positive and negative values. It turns out that unlike state II, both noise
and differential covariances have the same sign in state III, which could
cause the sudden decrease in linear Fisher information. This result is
consistent with previous theoretical analysis of neural population coding
[41], which attributes the information decrease to the overlap between these
two kinds of covariances.

The interplay between the differential covariance and noise covariance can be
further understood geometrically by considering two-dimensional projections of
the bump attractor state onto the space spanned by the spike count of two
neurons. As shown in Figure 3c, the black lines represent the mean firing rate
of the bump attractor in the space spanned by the activity of two neurons
whose preferred $s$ differ by $0.2\pi$. Orange ellipses represent the firing
covariability (scaled for visibility) for when the mean firing rates take
specific values on the ring manifold as marked by the red dots. In this
figure, the differential covariance corresponds to the tangent space of the
attractor manifold represented by the mean firing rate, whereas noise
covariance corresponds to the orange ellipses. Remarkably, the major axis
(blue line) of the noise covariance is roughly orthogonal to the attractor
manifold in state II but is roughly parallel to the attractor manifold in
state III. As a result, the coding error, which is dominated by the projection
of the noise covariance onto the attractor manifold, is smaller in state II
than in state III.

Previous theoretical studies on neural population coding suggest that Fisher
information may saturate with increasing neural population size due to
information-limiting correlations [6, 41, 44]. It is shown analytically that
in a neural population with homogeneous tuning function and distance-dependent
noise correlation $\rho(x_{1},x_{2})=\rho(x_{1}-x_{2})$, Fisher information
grows linearly or superlinearly with population size when the noise
correlation is zero or negative, respectively, but saturates to a finite value
when the noise correlation is positive [59]. However, this saturation may
disappear when the neural population has heterogeneous tuning [20]. It is
pointed out that such information-limiting correlation can be largely
attributed to the overlap between noise correlation C and differential
correlation $\mathbf{u}\mathbf{u}^{T}$ [41].

We run simulations of the neural circuit model with varying numbers of neurons
to investigate the effect of population size on linear Fisher information rate
$I$ for states II and III, respectively (see Figure 3d). We find that in both
cases, linear Fisher information grows approximately linearly with the
population size, though this growth is severely reduced in the case of state
III. This outcome can be explained by observing that the noise covariance in
state III significantly overlaps with differential covariance but not in state
II (see Figure 3b). While information is reduced in state III, it does not
saturate to a finite limit as the noise covariance is similar in sign but not
completely parallel to the differential covariance. This is confirmed by
calculating the ratio between the information and the population size, which
converges to a nonvanishing value in the limit of large population size
(Supplementary Information, including Figure S1),and is further elucidated by
theoretical analysis using simplified covariance structures (see equations S1
and S30 in Supplementary Information).

These results are thus consistent with the theory that differential
correlation serves as an information-limiting correlation [41]. Furthermore,
our model captures the emergence of mean firing rate $\mu(x)$ and noise
covariance $C(x,{x}^{\prime})$ through the intrinsic dynamics of the recurrent
circuit, thereby offering an explanation about the dynamical origin of
information-limiting correlations. Importantly, our model demonstrates that
such information-limiting correlation can be switched on and off by changing a
global parameter (the external input mean $\mu_{\mathrm{ext}}$ in this case),
which determines the dynamical regime of the bump attractor.

Previously, it has also been found that the information content in a bump
attractor is more concentrated around the edges of the bump state with
independent Poisson activity [57]. We also investigate the contribution of
linear Fisher information from each neuron as shown in Supplementary
Information, Figure S2. We find that in state II, the neurons located at the
boundaries of the bump are more important for coding than other neurons, as
consistent with previous models [57]. Unexpectedly, in state III, the neurons
near the center of the bump provide negative contributions to the linear
Fisher information, suggesting that some neurons may be harmful for coding.
Since the correlated fluctuation is coupled with the neural dynamics in our
model, our findings suggest a way for the neural population to directly
regulate its coding efficiency through manipulating its intrinsic correlation
structure.

The analysis above primarily focuses on the coding property of stationary bump
activity. In the context of working memory, this can be thought as
representing the neural activity state during the delay epoch of working
memory. In addition to such time-independent mnemonic coding subspace, [43]
have also identified the presence of a dynamic coding subspace by analyzing
neural spike data of the primate prefrontal cortex during working memory
tasks. Following [43], we investigate the dynamic coding in our MNN by
analyzing the temporal dynamics of the bump state during and after cue
presentation for each of states II and III (see Figures S5-S6 in Supplementary
Information). We observe that the population correlation (equation S1 of [43])
between the population activity states and a reference population activity at
the cue epoch reaches its peak during the cue epoch and gradually decays
during the delay epoch (see Figure S5c, Supplementary Information). This
observation suggests the existence of a dynamic subspace that undergoes decay
during the delay epoch. Furthermore, we find that the linear Fisher
information of the dynamic coding in state II (see Figure 2)is significantly
higher than that of the mnemonic coding, and the linear Fisher information of
both types of coding becomes similar during the delay epoch (see Figure S5d,
Supplementary Information). This is consistent with that study, which
demonstrated that the decoding accuracy of dynamic coding is substantially
higher than that of mnemonic coding during the cue epoch, but both accuracies
become similar during the delay epoch. Interestingly, we find the opposite
result in state III (see Figure 2), where the linear Fisher information of
mnemonic coding is significantly higher than that of the dynamic coding during
the cue period. This suggests that the brain might adopt a correlation
structure in state II during the cue epoch for better coding efficiency. (See
the Supplementary Information for details.)

Figure 4: Interactions of bump attractors mediated by negatively correlated
fluctuations.  (a) Mean firing rate $\mu$, firing variability $\sigma^{2}$ and
correlation coefficient $\rho$ of two bumps in stable bound state. (b) Two
bumps separated by a distance closer to or farther from the stable equilibrium
at $d^{\ast}\approx 0.63\pi$ experience repulsion or attraction, respectively.
(c) The effective interaction energy $E$ as a function of interbump distance
$d$ (solid line) exhibits a local minimum at $d\approx 0.65\pi$ (dot), which
disappears when the negative correlation is removed (dashed line). We set
$\mu_{\mathrm{ext}}=0.922$ in panels a-c. (d) Stable equilibrium distance
$d^{\ast}$ between two bumps for varying external input mean
$\mu_{\mathrm{ext}}$ (black lines), overlaid on the effective interaction
energy $E$ (heat map).

##  5 Stable bound states of bump attractors

In the previous section, we revealed the influence of correlated neural
variability on neural coding by a single bump attractor state. In this
section, we further investigate its impact on the dynamical interaction of
multiple bump attractor states. In bump attractor networks, multiple bump
states elicited by different stimuli may coexist in the same neural circuit
and interact through distance-dependent synaptic connections. Conventional
firing rate–based models predict that such interactions typically lead to
either repulsion due to long-range synaptic inhibition or merging due to
short-range excitation [33, 72].

In our model, however, we find that bump attractors may become dynamically
coupled through negatively correlated fluctuations, forming stable bound
states previously unobserved in firing rate–based models. Figure 4a shows two
bumps (see State II in Figure 2a) forming such a stable bound state. Both the
mean firing rate $\mu$ and firing variability $\sigma^{2}$ of each of them
largely retain the similar shapes as that of a single isolated bump. However,
the system now also exhibits negative correlations between two bumps, in
addition to the positive correlations within each bump. Remarkably, the
interbump distance $d$ of such a two-bump state forms a stable equilibrium,
that is, if the system is initialized with the interbump distance $d$ smaller
or larger than the stable equilibrium, the bumps will repulse or attract each
other to eventually reach the stable distance, as shown in Figure 4b. Similar
phenomena of stable bound states have been found in other noise-coupled
nonlinear systems such as interacting solitons in mode-locked lasers [27, 70]
and electron-electron attraction due to lattice vibration in superconductors
[7].

To quantify the interaction between two bumps, we employ a projection method
to calculate the effective interaction energy. In this method, we first clamp
one of the two bumps in Figure 4b and consider its influence on the other bump
(the free-moving bump) as if it acts like an external input. We next assume
that the shape of the free-moving bump does not vary significantly as it is
shifted around the stable equilibrium $s^{\mathrm{f}}$ by $z$. Under these
simplifications, we can project the entire system state on to the neutrally
stable left eigenspace along the attractor manifold represented by the mean
firing rate of the free-moving bump to arrive at an approximate equation of
motion for the bump’s spatial deviation $z$ [11]. Due to the complexity of our
model, explicit analytical solution is prohibitive, and we instead resort to
numerical analysis as follows. We approximate the mean firing rate and firing
covariability of the free bump as
$\mu^{\mathrm{f}}_{i}(z)={\mu}^{\mathrm{f}}(x_{i}|s^{\mathrm{f}}+z)$ and
$C^{\mathrm{f}}_{ij}(z)={C}^{\mathrm{ff}}(x_{i},x_{j}|s^{\mathrm{f}}+z)$,
respectively, where ${\mu}^{\mathrm{f}}(x|s^{\mathrm{f}})$ and
${C}^{\mathrm{ff}}(x,{x}^{\prime}|s^{\mathrm{f}})$ represent the shape of the
free bump’s mean and covariance at the stable equilibrium $s^{\mathrm{f}}$.
Following [62], we denote $\langle
f,g\rangle=\tfrac{2\pi}{N}\sum_{i}f(x_{i})g(x_{i})$ as the inner product, and
calculate the left eigenvector of the system for the free bump near the
equilibrium as $v=a\frac{d\bar{\mu}^{\mathrm{f}}}{dz}$, where
$\bar{\mu}^{\mathrm{f}}$ is the synaptic current received by the free bump
(equation S37), and
$a=1/\langle\frac{d\bar{\mu}^{\mathrm{f}}}{dz},\frac{d\mu^{\mathrm{f}}}{dz}\rangle$
is a normalization coefficient. Projecting both sides of equation (1) for the
full system to the left-eigenspace spanned by $v$, we obtain

$\displaystyle\frac{dz}{dt}=\left\langle
v,-\mu(z)+\phi\big{(}\bar{\mu}(z),\bar{C}(z)\big{)}\right\rangle:=b(d),$ (5)

where $\mu(z)$ and $C(z)$ denote the system state after shifting the free bump
by $z$, $d=s^{\mathrm{c}}-s^{\mathrm{f}}-z$ is the interbump distance, and
$s^{\mathrm{c}}$is the peak location of the clamped bump. Based on this
equation of motion, given the distance $d$ between two bumps, we can define
the corresponding effective interaction energy as
$E(d)=-\int_{\pi}^{d}b(x)dx$. (See equations S32 to S46 for details of the
derivations of$b(d)$.) As shown in Figure 4c, the stable bound state of bump
attractors corresponds to the local minimum of the effective interaction
energy with an interbump distance of
$d^{\ast}=s^{\mathrm{c}}-s^{\mathrm{f}}=0.65\pi$. Interestingly, when the
negative correlation between two bumps is clamped to zero, the bound state
loses its stability as indicated by the disappearance of the local minimum in
the effective interaction energy (see Figure 4c), resulting in repulsive
interaction. We find that as we increase the strength of the external input
mean $\mu_{\mathrm{ext}}$, the stable distance $d^{\ast}$ between two bump
attactors increases, as shown by the black curve in Figure 4d. Interestingly,
the stable distance $d^{\ast}$ exhibits a steep increase near
$\mu_{\mathrm{ext}}=0.93$, which coincides with the transition point between
states II and III of a single bump in Figure 2. However, the internal
covariance structures of these two bumps do not change until
$\mu_{\mathrm{ext}}$is much higher, presumably due to the interaction between
the bumps (see Figure S3). The stable distance curve in Figure 4d is laced
over the effective interaction energy, whose local minimum for each fixed
$\mu_{\mathrm{ext}}$ coincides with $d^{\ast}$. These results indicate that
negative interbump correlation is responsible for generating the stable-bound
states of bump attractors in our model, preventing strong repulsion between
two bumps as typically observed in firing rate-based models.

To gain further theoretical understanding of the origin of negative interbump
correlations, we consider the overall interactions between the neural
populations of these two bumps, which can be roughly characterized by two
neural masses with self-excitation and mutual inhibition. Denote the mean
firing rate and firing covariability of these two neural masses at the steady
state as

$\displaystyle\mu=\begin{pmatrix}\mu_{1}\\\ \mu_{2}\end{pmatrix},\quad
C=\begin{pmatrix}c_{11}&c_{12}\\\ c_{21}&c_{22}\end{pmatrix},$ (6)

and the average synaptic weights between the two populations as
$w=\bigl{(}\begin{smallmatrix}w_{11}&w_{12}\\\
w_{21}&w_{22}\end{smallmatrix}\bigr{)}$, such that $w_{11},~{}w_{22}>0$ and
$w_{12},~{}w_{21}<0$, as required by the self-excitation and the mutual
inhibition. Under mild conditions over the neural activation, which are
satisfied in our model, the activity of these two populations of neurons
becomes negatively correlated: $c_{12}<0$ (see equations S47 and S55 for
proof). This implies that the negative correlation is a direct consequence of
the inhibition between two neural populations. This analysis leads to the
surprising conclusion that mutual inhibition between bumps increases working
memory capacity by reducing the minimal distance between adjacent bumps,
contrary to conventional wisdom suggesting that inhibitions should limit
working memory capacity through repulsive interactions between bumps [21, 33].
We turn to this point in the next section.

##  6 Intrinsic negative correlations regulate dynamical interactions between
memorized items

Previous studies using a firing rate–based model found that the interactions
between two bumps are attractive when they are close, and they become
repulsive past an unstable equilibrium [72], similar to the energy landscape
indicated by the dashed line in Figure 4c. However, such attractive or
repulsive interactions will either cause the two bumps to merge, preventing
the system from discriminating two distinct items [2],or to drift away from
each other, thereby limiting the coding resolution. First, our model
demonstrates that through negative correlations, two bumps can maintain a
closer distance than the case without negative correlations. Moreover, when
two bumps are too close, the repulsive interaction will restore their distance
to the stable equilibrium. As a result, our findings provide a mechanism
enabling discrimination of similar items encoded by bump states [36]. Second,
our results suggest that nonlinearly coupled, correlated noise can have
unexpected effects of stabilizing the spatial location of bump states, in
contrast to previous findings that bump states should drift away along the
attractor manifold under the influence of additive uncorrelated noise [11]. In
particular, when the position of one bump is maintained through external
mechanisms such as external stimulus or top-down attention, another adjacent
bump state can remain stabilized through the attraction toward the equilibrium
distance. Our model thus provides a potential way to improve the coding
accuracy of a bump attractor against this drifting effect.

The preceding analysis of bump interactions can be applied to understand a
range of experimental observations of cognitive and behavioral processes
involved in working memory [71, 73] and eye movement presumably reflecting
higher cognitive states [65, 60, 72]. It has been widely observed that the
trajectory of saccadic eye movement (typically from an initial fixation to a
target location) can deviate either toward [66] or away from [61]an element
other than the target in visual field. In one study [66],human participants
were instructed to perform two saccades from fixation toward different targets
in consecutive trials. When the delay between two trials was short and the
distance between two targets was large, the trajectory of the second saccade
deviated toward the first target instead of directly going to the second
target through a straight line. In another study [61], on the contrary, the
trajectory of eye movement deviated away from a nearby visual stimulus whose
location was required to be memorized. Existing theoretical proposals for
explaining these phenomena involve top-down modulation whose ability to
control rapid eye movements, however, appears to be limited both temporally
and spatially [64], and an exact mechanism is still unclear [60].

The stable bound states found in our model provide an alternative neural
mechanism for the above phenomena without resorting to top-down modulation. We
first establish the connection between our model and the above two
experiments. In our bump attractor model, the ring attractor manifold
represents the angular direction in visual space, whereas the previously
memorized and the current targets of eye movement are encoded by the centers
of bump activities. The clamped bump in our model as in Figure 4correspond to
the direction of the previously shown target in [66] or that of the stimulus
stored in working memory in [61]. In both cases, the clamped bump induces an
energy landscape like that in Figure 4(c),which interacts with the free-moving
bump encoding the angular direction of the affected eye movement. It turns out
that the influence of the energy landscape induced by the clamped bump is
consistent with these experiments: the interaction of the two bumps is
attractive when their distance is large, which can explain the deviation of
eye movement toward the previously shown target in [66],and is repulsive when
their distance is small, which can explain the deviation of the eye movement
away from the stimulus location in [61]. The nontrivial dependence of bump
interaction on their distance as predicted by our model can be tested by
future experiments through systematically investigating the dependence of
behavioral errors during working memory tasks on the difference between
memorized items or between one memorized item and a distractor stimulus. These
results demonstrate that negatively correlated neural fluctuations modulate
the interaction between bump states, which in turn influences cognitive
functions and behavior and highlights the importance of considering correlated
neural fluctuations in modeling neural circuit dynamics. In the next section,
we analyze how these noise-mediated dynamical effects can facilitate neural
representation in working memory.

##  7 Intrinsic negative correlations lead to efficient coding strategies for
working memory

Due to their translation invariance and bistability, bump attractors have been
proposed as the neural substrate for spatial working memory [71, 75, 73, 46,
21]. The number of bumps the system can simultaneously maintain is interpreted
as working memory capacity. To investigate how working memory capacity is
influenced by correlated neural fluctuations, we elicit multiple bump states
in our model with bell-shaped transient inputs.

Figure 5: Multiple bump attractors with intricate correlation structures. The
mean firing rate, firing variability, and correlation coefficient of the
neural circuit with three (left panels) or four (right panels) bumps. In the
three-bump case, correlations between different bumps are mainly negative. In
the four-bump case, adjacent bumps are negatively correlated, while every
other bump is positively correlated. (b) The circuit fails to support more
than two bumps when the negative correlation is clamped to zero.

Neural activity states holding multiple (three and four) bumps (see state II,
$\mu_{\mathrm{ext}}=0.935$) are shown in Figure 5a. The shape of mean firing
rate $\mu$ and firing variability $\sigma^{2}$ of individual bumps are
qualitatively consistent with that of an isolated bump (Figure 2), but with
slight deformations, lower mean firing rate, and higher firing variability.
When the number of bumps is three, correlations with different bumps are
mainly negative, similar to that in a two bumps’ case (see Figure 4). When the
number of bumps is four, adjacent bumps are negatively correlated (spaced by
$\pi/2$), while every other bump is positively correlated (spaced by $\pi$).
As we have shown with the two-bump case, negative correlation can reverse the
strong repulsion between adjacent bumps, suggesting that negative correlation
could be a key factor for enabling the multibump configuration shown above. To
test this idea, we clamp the negative correlations in the four-bump case in
Figure 5a to zero and find that the activity of some bumps is quenched by the
mutual inhibition from adjacent bumps, resulting in the failure for the system
to hold more than two bumps, as shown in Figure 5b. This observation is
consistent with that the removal of negative correlation turns the interaction
between two bumps from a balance of attraction and repulsion into pure
repulsion, as described in Figure 4. These results suggest that negative
correlations can enable multiple bumps to become more closely packed, thereby
enhancing the memory capacity of the system.

Experimental evidence suggests that working memory capacity is limited [16,
17]. There are two theoretical mechanisms through which such limit may occur.
The first is a spatial mechanism that mutual inhibition limits the number of
stable bump states that can be simultaneously maintained. The second of them
is a temporal mechanism based on the theta-gamma coupling of neural
oscillation, in which a limited number of gamma cycles, representing distinct
memorized items, can be fitted in each theta cycle [37]. Interestingly,
although our model does not explicitly represent temporal oscillations, it is
compatible with the interpretation that the observed correlations may reflect
the synchrony or antisynchrony between different bump states. Thus, by taking
into account correlated neural variability, our model hints toward a unifying
perspective for the limit of working memory capacity based on spatiotemporal
working memory representation.

These results also suggest an efficient computational strategy for memory
storage by fluctuating neural activity, in which limited resources (i.e.,
neural spikes) are distributed across both space and time. Spatially, negative
correlations increase memory capacity by reducing the distance between
adjacent bumps (see Figure 4). Temporally, negative correlations can be
interpreted as a form of antisynchrony that corresponds to temporally less
concentrated bump activity states for each stored item. Conceptually, this is
similar to time-division multiplexing recently proposed in a spiking neural
network model [1, 12] and can also be understood as a form of sparse coding.

We also find that negative correlations can improve the robustness of working
memory by preventing the forgetting of a previously stored bump when a new
stimulus arrives. To demonstrate this, we sequentially inject transient
stimuli with durations of $T=0.1$ s and various spatial separations and
investigate the resulting dynamics of bump attractors. As shown in Figure 6,
if the initial distance between two bumps is sufficiently large ($0.6\pi$),
both remain persistent after the removal of the transient stimuli and reach a
stable distance under either attractive (top panel) or repulsive (middle
panel) interactions. If the initial distance between two bumps is sufficiently
small ($0.5\pi$), the two bumps will merge into one bump (bottom panel). Next,
we clamp the negative correlations to zero and observe the effect on the
dynamics of sequentially memorized bump states (see the right column in Figure
6). If the initial distance between two bumps is sufficiently large
($d=0.8\pi$, top panel), the second bump is maintained, while the previously
stored bump is quenched. If the initial distance between two bumps is
sufficiently small ($d=0.6\pi$, middle panel), the two bumps will annihilate
each other. If the initial distance between two bumps is too small
($d=0.5\pi$, bottom panel), the two bumps will merge into one bump, which is
the same outcome as the case with negative correlations. These results suggest
that the negative correlations help to prevent forgetting due to competition
from similar memorized items

Figure 6: Negative correlations prevent forgetting of sequentially memorized
items. The left column shows the cases in which we keep the negative
correlations. When the initial distance between two bumps is $0.8\pi$ (top
panel) or $0.6\pi$ (middle panel), both of them remain persistent after the
removal of the transient stimuli and reach a stable distance. If the initial
distance between two bumps is sufficiently small ($0.5\pi$),the two bumps will
merge into one bump (bottom panel). The right column shows the cases in which
we clamp the negative correlations to zero. If the initial distance between
two bumps is $d=0.8\pi$, the second bump is maintained, while the previously
stored bump is quenched (top panel). If the initial distance between two bumps
is $d=0.6\pi$, the two bumps will annihilate each other (middle panel). If the
initial distance between two bumps is $d=0.5\pi$, the two bumps will merge
into one bump (bottom panel). The system is initialized with zero mean and
covariance, and two transient stimuli of duration $T=0.1$ s (marked by the red
rectangles) are exerted at two locations sequentially. The grayscale indicates
mean firing rate (sp/s).

##  8 Discussion

In this work, we have revealed the emergence of synergistic neural population
codes and their mechanistic origin from nonlinearly coupled, correlated neural
fluctuations in a bump attractor network. Our model demonstrates that complex
spatial correlation structures can emerge from a bump attractor network model,
as shown in Figure 2, suggesting the need for closer experimental
investigations of the pairwise correlations of neural population activity in
the brain. As shown in Figure 3a, the model additionally shows that global
neural circuit parameters (such as the external input mean) can drastically
change the dynamical regimes of bump attractor states with distinct noise
correlation patterns, which can be interpreted as corresponding to different
cognitive states. Specifically, our model predicts that behavioral error
should be larger when the noise covariance is parallel to the tangent space of
the ring manifold than when the noise covariance is orthogonal to the tangent
space of the ring manifold. Future experimental studies could investigate this
relationship between the neural population codes and behavioral performance
using simultaneous recordings of cortical neurons during cognitive tasks.
Furthermore, we reveal correlation induced dynamical interactions between two
bump states, which lead to a stable bound state of two bumps as characterized
by a local minimum in the effective interaction energy between them as shown
in Figure 4. As a result, the interaction between two bumps, which represent
the memorized items, can be attractive or repulsive depending on their
distance (i.e., difference in the encoded features). To test our theory,
future experimental studies could systematically investigate how behavioral
errors during working memory tasks depend on the difference between memorized
items or between one memorized item and a distractor stimulus. We also show
that the negative correlation between the neural activity of two bumps are
originated from mutual inhibitions in the recurrent network, a property often
overlooked in previous firing rate models of bump attractors. These results
thus highlight the importance of incorporating correlated neural variability
in modeling studies of cortical dynamics.

The validity of analyzing stochastic neuronal dynamics through approximations
up to second-order moments is supported by both empirical evidence and
theoretical analysis. Through examining multielectrode array recordings in the
vertebrate retina [53],it is found that the strong collective behavior
exhibited by a neural population is well accounted for by a model that
captures the observed pairwise correlations but assumes no higher-order
interactions, even in the case where the pairwise correlations are weak. A
theoretical study [18] investigates the effect of higher-order cumulants to
the dynamics in a stochastic binary neural network model. It is found that in
a weakly correlated state, the contributions from higher-order cumulants to
the network dynamics are significantly smaller than pairwise interactions.
This suggests that the first two moments are fairly accurate at capturing the
stochastic dynamics of the system. In summary, it can be argued that our MNN
offers a minimalistic yet effective representation of fluctuating neural
activity and is suitable for the purpose of modeling working memory in this
work.

The MNN model presented in this study has a couple of limitations. First, the
current version of MNN lacks the depiction of temporal statistics of neural
activity, which may encode computationally useful information such as phase
and time-delayed synchrony. Future studies are required to develop moment
mappings that account for temporal statistics such as auto- and cross-
covariances of neural spike trains [18]. Second, in order to simulate the
dynamics of the MNN, it is necessary to store the covariance matrix of the
system, which entails a quadratic space complexity. Furthermore, multiplying
the covariance matrix with the connection weights results in a cubic time
complexity. As a consequence, simulating large-scale neural systems could
require significant computational resources. However, this trade-off is
worthwhile considering that high dimensional joint probability distribution of
neural activity is computationally intractable.

In the context of ring attractor networks, the MNN based on the gaussian field
approximation can be generalized in a number of ways. First, we can consider a
two-dimensional spatial structure instead of just one to incorporate more
complex two-dimensional patterns [25, 47]. A straightforward 2D extension of
the model considered in this paper is included in Supplementary Information,
where the firing covariability manifests as more complex patterns compared to
the one-dimensional case. Second, this study considers stationary bump
activity with correlated variability by using a model derived from the
current-based leaky integrate-and-fire neurons without adaptation. Previous
theoretical studies of bump attractors using firing rate models suggest that
neural adaptations can lead to the formation of complex dynamical patterns
such as propagating and rotating waves [25, 48, 14, 49, 55] and that short-
term synaptic facilitation and depression can have a significant impact on the
dynamics of bump diffusion [54, 49]. It is thus of great theoretical interest
to further extend our model to incorporate neural adaptation, such as spike
frequency adaptation, or synaptic depression, and investigate the influence of
correlated neural fluctuations on the dynamics of the resulting wave patterns
as well as its functional implications. Third, neural fluctuations may also
cause bump states to drift away along the attractor manifold over time,
resulting in further degradation of working memory. It has been shown that
uncorrelated noise imposes a bound on the coding capability of bump attractor
[11]. The model presented in this paper may be extended to analyze the
influence of correlated noise on this diffusive effect and its impact on
working memory performance.

Although we have primarily focused on bump attractor networks in this paper,
the general approach to modeling correlated neural variability using the MNN
can be easily extended to other types of models. For instance, [43] proposed a
parsimonious linear model to interpret the coexistence of dynamical coding and
mnemonic coding in working memory. The network structure of this model can be
applied to our MNN to gain a deeper understanding of the specific role played
by firing covaribility in dynamical coding and mnemonic coding in working
memory. Another extension is to divide the neural population into excitatory
and inhibitory ones, as consistent with Dale’s principle [19], for analyzing
the dynamics and computational roles of firing covaribility in balanced
network [67] and in neural oscillations [45]. Moreover, the MNN considered in
this paper is derived from the simple leaky integrate-and-fire neuron model.
The general framework of MNN can further incorporate more complex spiking
neuron models with slow and fast synaptic dynamics, which are important in
generating complex dynamics such as bursting activity [32]. Finally, the
modeling framework based on the MNN also opens up new opportunities for
investigating the role of neural correlations in synaptic plasticity,
learning, and attention [34].

## Appendix

###  8.1 The moment neural network

To capture the nonlinear coupling of correlated neural activity, we employ a
model known as the moment neural network derived from the following leaky
integrate-and-fire neuron model [23, 38]

$\displaystyle\dfrac{dV_{i}(t)}{dt}=-LV_{i}(t)+I_{i}(t),\quad
I_{i}(t)=\sum_{j=1}^{N}w_{ij}S_{j}(t),$ (7)

where $V_{i}$ is the membrane potential of the $i$th neuron, $L$ is the leak
conductance, $I_{i}(t)$ is the total synaptic current, $w_{ij}$ is the
synaptic weight from neuron $j$ to neuron $i$, and $S_{j}(t)$ is the spike
train from pre-synaptic neurons. When $V_{i}$ reaches the firing threshold
$V_{\text{th}}$ the neuron will release a spike which is transmitted to other
connected neurons, and then $V_{i}$ is reset to the resting potential
$V_{\text{res}}$ and enters a refractory period $T_{\text{ref}}$. The
functions $\phi_{\mu}$ and $\phi_{\sigma}$ together map the mean $\bar{\mu}$
and variance $\bar{\sigma}^{2}$ of the input current $I_{i}(t)$ to that of the
output spikes according to [23, 38]

$\displaystyle\phi_{\mu}:(\bar{\mu},\bar{\sigma}^{2})\mapsto\mu,\quad\mu=\dfrac{1}{T_{\rm
ref}+\tfrac{2}{L}\int_{I_{\rm lb}}^{I_{\rm ub}}g(x)dx},$ (8)
$\displaystyle\phi_{\sigma}:(\bar{\mu},\bar{\sigma}^{2})\mapsto\sigma^{2},\quad\sigma^{2}=\tfrac{8}{L^{2}}\mu^{3}\textstyle\int_{I_{\rm
lb}}^{I_{\rm ub}}h(x)dx,$ (9)

where $T_{\rm ref}$ is the refractory period with integration bounds $I_{\rm
ub}(\bar{\mu},\bar{\sigma})=\tfrac{V_{\rm
th}L-\bar{\mu}}{\sqrt{L}\bar{\sigma}}$ and $I_{\rm
lb}(\bar{\mu},\bar{\sigma})=\tfrac{V_{\rm
res}L-\bar{\mu}}{\sqrt{L}\bar{\sigma}}$. The constant parameters $L$, $V_{\rm
res}$, and $V_{\rm th}$ are identical to those in equation (7). The pair of
Dawson-like functions $g(x)$ and $h(x)$ appearing in equations (8) and (9) are
$g(x)=e^{x^{2}}\int_{-\infty}^{x}e^{-u^{2}}du$ and
$h(x)=e^{x^{2}}\int_{-\infty}^{x}e^{-u^{2}}[g(u)]^{2}du$. The mapping $\psi$,
which we refer to as the linear response coefficient, is equal to
$\psi(\bar{\mu},\bar{\sigma}^{2})=\tfrac{\partial\mu}{\partial\bar{\mu}}$ and
it is derived using a linear perturbation analysis around correlation
coefficient $\bar{\rho}_{ij}=0$ [51, 38]. This approximation is justified as
pairwise correlations between neurons in the brain are typically weak.
Throughout this paper, we set neuron parameters to be $V_{\rm th}=20$ mV,
$V_{\rm res}=0$ mV, $T_{\rm ref}=5$ ms, and $\tau=1/L=20$ ms. An efficient
numerical algorithm is used for implementing the moment mappings [47].

###  8.2 Network parameters

We construct a ring network by setting distance-dependent synaptic weights
with short-range excitation and long-range inhibition as follows,

$\displaystyle
w_{ij}=\tfrac{2\pi}{N}[w_{E}\kappa(x_{i}-x_{j};d_{E})-w_{E}\kappa(x_{i}-x_{j};d_{I})],$
(10)

where $\kappa(x;d)=\exp[\tfrac{1}{d^{2}}(\cos x-1)]$, $w_{E}$ and $w_{I}$
control the strengths of excitatory and inhibitory synapses, and $d_{E}$ and
$d_{I}$ control their spatial range, respectively. Throughout this paper, we
set the synaptic weight parameters to be $w_{E}=15$, $w_{I}=6$, $d_{E}=0.5$,
$d_{I}=1$, and $N=400$. We fix $\sigma_{\mathrm{ext}}^{2}=0.01$ throughout
this study and systematically vary $\mu_{\mathrm{ext}}$. We use Euler’s method
for simulating equations (1) and (2) with a time step of $\Delta t=0.5\tau$.

###  8.3 Linear Fisher information analysis

To quantify the amount of information content encoded by a bump attractor, we
calculate the linear Fisher information $I_{\mathrm{linear}}$ as follows [30].
Given a bump state with mean firing rate $\mu(x|s)$ and firing covariability
$C(x,{x}^{\prime}|s)$ centered at $s$, the linear Fisher information can be
calculated

$\displaystyle I_{\mathrm{linear}}(\Delta
t)=\mathbf{u}(s)^{T}C(s)^{-1}\mathbf{u}(s)\Delta t,$ (11)

where $u_{i}(s)=\frac{d\mu_{i}(s)}{ds}$ and $\Delta t$ is the spike count time
window. Note that $I_{\mathrm{linear}}(\Delta t)$ becomes approximately
independent of $s$ due to the translation invarianace when the population size
is large. For calculating $I_{\mathrm{linear}}(\Delta t)$, we apply centered
finite difference to approximate $\mathbf{u}(s)$ and use Moore-Penrose pseudo-
inverse for calculating $C(s)^{-1}$.

### Acknowledgments

Supported by STI2030-Major Projects (No. 2021ZD0200204); supported by Shanghai
Municipal Science and Technology Major Project (No. 2018SHZDZX01), ZJ Lab, and
Shanghai Center for Brain Science and Brain-Inspired Technology; supported by
the 111 Project (No. B18015).

## References

  * Akam and Kullmann, [2014] Akam, T. and Kullmann, D. M. (2014).  Oscillatory multiplexing of population codes for selective communication in the mammalian brain.  Nature Reviews Neuroscience, 15(2):111–122.
  * Almeida et al., [2015] Almeida, R., Barbosa, J., and Compte, A. (2015).  Neural circuit basis of visuo-spatial working memory precision: a computational and behavioral study.  Journal of Neurophysiology, 114(3):1806–1818.
  * Amari, [1977] Amari, S.-i. (1977).  Dynamics of pattern formation in lateral-inhibition type neural fields.  Biological Cybernetics, 27(2):77–87.
  * Amit and Brunel, [1997] Amit, D. J. and Brunel, N. (1997).  Dynamics of a recurrent network of spiking neurons before and following learning.  Network: Computation in Neural Systems, 8(4):373–404.
  * Amit and Tsodyks, [1991] Amit, D. J. and Tsodyks, M. (1991).  Quantitative study of attractor neural network retrieving at low spike rates. I. substrate-spikes, rates and neuronal gain.  Network: Computation in Neural Systems, 2(3):259.
  * Averbeck et al., [2006] Averbeck, B. B., Latham, P. E., and Pouget, A. (2006).  Neural correlations, population coding and computation.  Nature Reviews Neuroscience, 7(5):358–366.
  * Bardeen et al., [1957] Bardeen, J., Cooper, L. N., and Schrieffer, J. R. (1957).  Microscopic theory of superconductivity.  Physical Review, 106(1):162.
  * Ben-Yishai et al., [1995] Ben-Yishai, R., Bar-Or, R. L., and Sompolinsky, H. (1995).  Theory of orientation tuning in visual cortex.  Proceedings of the National Academy of Sciences, 92(9):3844–3848.
  * Berens et al., [2012] Berens, P., Ecker, A. S., Cotton, R. J., Ma, W. J., Bethge, M., and Tolias, A. S. (2012).  A fast and simple population code for orientation in primate V1.  Journal of Neuroscience, 32(31):10618–10626.
  * Buice et al., [2010] Buice, M. A., Cowan, J. D., and Chow, C. C. (2010).  Systematic fluctuation expansion for neural network activity equations.  Neural Computation, 22(2):377–426.
  * Burak and Fiete, [2012] Burak, Y. and Fiete, I. R. (2012).  Fundamental limits on persistent activity in networks of noisy neurons.  Proceedings of the National Academy of Sciences, 109(43):17645–17650.
  * Caruso et al., [2018] Caruso, V. C., Mohl, J. T., Glynn, C., Lee, J., Willett, S. M., Zaman, A., Ebihara, A. F., Estrada, R., Freiwald, W. A., Tokdar, S. T., et al. (2018).  Single neurons may encode simultaneous stimuli by switching between activity patterns.  Nature Communications, 9:2715.
  * Compte et al., [2000] Compte, A., Brunel, N., Goldman-Rakic, P. S., and Wang, X.-J. (2000).  Synaptic mechanisms and network dynamics underlying spatial working memory in a cortical network model.  Cerebral Cortex, 10(9):910–923.
  * Coombes, [2005] Coombes, S. (2005).  Waves, bumps, and patterns in neural field theories.  Biological Cybernetics, 93(2):91–108.
  * Coombes and Schmidt, [2010] Coombes, S. and Schmidt, H. (2010).  Neural fields with sigmoidal firing rates: approximate solutions.  Discrete and Continuous Dynamical Systems, 28(4):1369.
  * Cowan, [2001] Cowan, N. (2001).  The magical number 4 in short-term memory: a reconsideration of mental storage capacity.  Behavioral and Brain Sciences, 24(1):87–114.
  * Cowan, [2010] Cowan, N. (2010).  The magical mystery four: how is working memory capacity limited, and why?  Current Directions in Psychological Science, 19(1):51–57.
  * Dahmen et al., [2016] Dahmen, D., Bos, H., and Helias, M. (2016).  Correlated fluctuations in strongly coupled binary networks beyond equilibrium.  Physical Review X, 6(3):031024.
  * Dale, [1935] Dale, H. (1935).  Pharmacology and nerve-endings.
  * Ecker et al., [2011] Ecker, A., Berens, P., Tolias, A., and Bethge, M. (2011).  The effect of noise correlations in populations of diversely tuned neurons.  The Journal of Neuroscience, 31(40):14272–14283.
  * Edin et al., [2009] Edin, F., Klingberg, T., Johansson, P., McNab, F., Tegnér, J., and Compte, A. (2009).  Mechanism for top-down control of working memory capacity.  Proceedings of the National Academy of Sciences, 106(16):6802–6807.
  * Esnaola-Acebes et al., [2022] Esnaola-Acebes, J. M., Roxin, A., and Wimmer, K. (2022).  Flexible integration of continuous sensory evidence in perceptual estimation tasks.  Proceedings of the National Academy of Sciences, 119(45):e2214441119.
  * Feng et al., [2006] Feng, J., Deng, Y., and Rossoni, E. (2006).  Dynamics of moment neuronal networks.  Physical Review E, 73(4):041906.
  * Fetsch et al., [2012] Fetsch, C. R., Pouget, A., DeAngelis, G. C., and Angelaki, D. E. (2012).  Neural correlates of reliability-based cue weighting during multisensory integration.  Nature Neuroscience, 15(1):146–154.
  * Folias and Bressloff, [2005] Folias, S. E. and Bressloff, P. C. (2005).  Breathers in two-dimensional neural media.  Physical Review Letters, 95(20):208107.
  * Graf et al., [2011] Graf, A., Kohn, A., Jazayeri, M., and Movshon, J. A. (2011).  Decoding the activity of neuronal populations in macaque primary visual cortex.  Nature Neuroscience, 14(2):239–245.
  * Grelu and Akhmediev, [2012] Grelu, P. and Akhmediev, N. (2012).  Dissipative solitons for mode-locked lasers.  Nature Photonics, 6(2):84–92.
  * Helias et al., [2014] Helias, M., Tetzlaff, T., and Diesmann, M. (2014).  The correlation structure of local neuronal networks intrinsically results from recurrent dynamics.  PLoS Computational Biology, 10(1):e1003428.
  * Kim et al., [2017] Kim, S. S., Rouault, H., Druckmann, S., and Jayaraman, V. (2017).  Ring attractor dynamics in the drosophila central brain.  Science, 356(6340):849–853.
  * Kohn et al., [2016] Kohn, A., Coen-Cagli, R., Kanitscheider, I., and Pouget, A. (2016).  Correlations and neuronal population information.  Annual Review of Neuroscience, 39:237.
  * Koyluoglu et al., [2017] Koyluoglu, O. O., Pertzov, Y., Manohar, S., Husain, M., and Fiete, I. R. (2017).  Fundamental bound on the persistence and capacity of short-term memory stored as graded persistent activity.  Elife, 6:e22225.
  * Krahe and Gabbiani, [2004] Krahe, R. and Gabbiani, F. (2004).  Burst firing in sensory systems.  Nature Reviews Neuroscience, 5(1):13–23.
  * Krishnan et al., [2018] Krishnan, N., Poll, D. B., and Kilpatrick, Z. P. (2018).  Synaptic efficacy shapes resource limitations in working memory.  Journal of Computational Neuroscience, 44(3):273–295.
  * Kuśmierz et al., [2017] Kuśmierz, L., Isomura, T., and Toyoizumi, T. (2017).  Learning with three factors: modulating hebbian plasticity with errors.  Current Opinion in Neurobiology, 46:170–177.
  * Li et al., [2021] Li, H.-H., Sprague, T. C., Yoo, A. H., Ma, W. J., and Curtis, C. E. (2021).  Joint representation of working memory and uncertainty in human cortex.  Neuron, 109(22):3699–3712.
  * Lin and Luck, [2009] Lin, P.-H. and Luck, S. J. (2009).  The influence of similarity on visual working memory representations.  Visual Cognition, 17(3):356–372.
  * Lisman and Jensen, [2013] Lisman, J. E. and Jensen, O. (2013).  The theta-gamma neural code.  Neuron, 77(6):1002–1016.
  * Lu et al., [2010] Lu, W., Rossoni, E., and Feng, J. (2010).  On a gaussian neuronal field model.  NeuroImage, 52(3):913–933.
  * Ma, [2010] Ma, W. J. (2010).  Signal detection theory, uncertainty, and poisson-like population codes.  Vision Research, 50(22):2308–2319.
  * McNaughton et al., [2006] McNaughton, B. L., Battaglia, F. P., Jensen, O., Moser, E. I., and Moser, M.-B. (2006).  Path integration and the neural basis of the ‘cognitive map’.  Nature Reviews Neuroscience, 7(8):663–678.
  * Moreno-Bote et al., [2014] Moreno-Bote, R., Beck, J., Kanitscheider, I., Pitkow, X., Latham, P., and Pouget, A. (2014).  Information-limiting correlations.  Nature Neuroscience, 17(10):1410–1417.
  * Moser et al., [2014] Moser, E. I., Roudi, Y., Witter, M. P., Kentros, C., Bonhoeffer, T., and Moser, M.-B. (2014).  Grid cells and cortical representation.  Nature Reviews Neuroscience, 15(7):466–481.
  * Murray et al., [2017] Murray, J. D., Bernacchia, A., Roy, N. A., Constantinidis, C., Romo, R., and Wang, X.-J. (2017).  Stable population coding for working memory coexists with heterogeneous neural dynamics in prefrontal cortex.  Proceedings of the National Academy of Sciences, 114(2):394–399.
  * Panzeri et al., [2022] Panzeri, S., Moroni, M., Safaai, H., and Harvey, C. D. (2022).  The structures and functions of correlations in neural population codes.  Nature Reviews Neuroscience, 23:551–567.
  * Pina et al., [2018] Pina, J. E., Bodner, M., and Ermentrout, B. (2018).  Oscillations in working memory and neural binding: A mechanism for multiple memories and their interactions.  PLoS computational biology, 14(11):e1006517.
  * Polk et al., [2012] Polk, A., Litwin-Kumar, A., and Doiron, B. (2012).  Correlated neural variability in persistent state networks.  Proceedings of the National Academy of Sciences, 109(16):6295–6300.
  * Qi, [2022] Qi, Y. (2022).  An efficient numerical algorithm for the moment neural activation.  arXiv preprint.
  * Qi and Gong, [2015] Qi, Y. and Gong, P. (2015).  Dynamic patterns in a two-dimensional neural field with refractoriness.  Physical Review E, 92(2):022702.
  * Qi and Gong, [2022] Qi, Y. and Gong, P. (2022).  Fractional neural sampling as a theory of spatiotemporal probabilistic computations in neural circuits.  Nature Communications, 13(1):1–19.
  * Renart et al., [2007] Renart, A., Moreno-Bote, R., Wang, X.-J., and Parga, N. (2007).  Mean-driven and fluctuation-driven persistent activity in recurrent networks.  Neural Computation, 19(1):1–46.
  * Rocha et al., [2007] Rocha, J. d. l., Doiron, B., Shea-Brown, E., Josić, K., and Reyes, A. (2007).  Correlation between neural spike trains increases with firing rate.  Nature.
  * Rosenbaum et al., [2017] Rosenbaum, R., Smith, M. A., Kohn, A., Rubin, J. E., and Doiron, B. (2017).  The spatial structure of correlated neuronal variability.  Nature Neuroscience, 20(1):107–114.
  * Schneidman et al., [2006] Schneidman, E., Berry, M. J., Segev, R., and Bialek, W. (2006).  Weak pairwise correlations imply strongly correlated network states in a neural population.  Nature, 440(7087):1007–1012.
  * Seeholzer et al., [2019] Seeholzer, A., Deger, M., and Gerstner, W. (2019).  Stability of working memory in continuous attractor networks under the control of short-term plasticity.  PLoS Computational Biology, 15(4):e1006928.
  * Senk et al., [2020] Senk, J., Korvasová, K., Schuecker, J., Hagen, E., Tetzlaff, T., Diesmann, M., and Helias, M. (2020).  Conditions for wave trains in spiking neural networks.  Physical review research, 2(2):023174.
  * Seung et al., [2000] Seung, H. S., Lee, D. D., Reis, B. Y., and Tank, D. W. (2000).  Stability of the memory of eye position in a recurrent network of conductance-based model neurons.  Neuron, 26(1):259–271.
  * Seung and Sompolinsky, [1993] Seung, H. S. and Sompolinsky, H. (1993).  Simple models for reading neuronal population codes.  Proceedings of the National Academy of Sciences, 90(22):10749–10753.
  * Shamir and Sompolinsky, [2004] Shamir, M. and Sompolinsky, H. (2004).  Nonlinear population codes.  Neural Computation, 16(6):1105–1136.
  * Sompolinsky et al., [2001] Sompolinsky, H., Yoon, H., Kang, K., and Shamir, M. (2001).  Population coding in neuronal systems with correlated noise.  Physical Review E, 64(5):051904.
  * Theeuwes et al., [2009] Theeuwes, J., Belopolsky, A., and Olivers, C. N. (2009).  Interactions between working memory, attention and eye movements.  Acta Psychologica, 132(2):106–114.
  * Theeuwes et al., [2005] Theeuwes, J., Olivers, C. N., and Chizk, C. L. (2005).  Remembering a location makes the eyes curve away.  Psychological Science, 16(3):196–199.
  * Touboul and Ermentrout, [2011] Touboul, J. D. and Ermentrout, G. B. (2011).  Finite-size and correlation-induced effects in mean-field dynamics.  Journal of Computational Neuroscience, 31(3):453–484.
  * Urai et al., [2022] Urai, A. E., Doiron, B., Leifer, A. M., and Churchland, A. K. (2022).  Large-scale neural recordings call for new insights to link brain and behavior.  Nature Neuroscience, 25(1):11–19.
  * Van der Stigchel et al., [2009] Van der Stigchel, S., Belopolsky, A. V., Peters, J. C., Wijnen, J. G., Meeter, M., and Theeuwes, J. (2009).  The limits of top-down control of visual attention.  Acta Psychologica, 132(3):201–212.
  * Van der Stigchel et al., [2006] Van der Stigchel, S., Meeter, M., and Theeuwes, J. (2006).  Eye movement trajectories and what they tell us.  Neuroscience & Biobehavioral Reviews, 30(5):666–679.
  * Van Gisbergen et al., [1987] Van Gisbergen, J., Van Opstal, A., and Roebroek, J. (1987).  Stimulus-induced midflight modification of saccade trajectories.  In Eye Movements from Physiology to Cognition, pages 27–36. Elsevier.
  * Van Vreeswijk and Sompolinsky, [1996] Van Vreeswijk, C. and Sompolinsky, H. (1996).  Chaos in neuronal networks with balanced excitatory and inhibitory activity.  Science, 274(5293):1724–1726.
  * Wang and Kang, [2022] Wang, R. and Kang, L. (2022).  Multiple bumps can enhance robustness to noise in continuous attractor networks.  BioRxiv.
  * Wang, [2001] Wang, X.-J. (2001).  Synaptic reverberation underlying mnemonic persistent activity.  Trends in Neurosciences, 24(8):455–463.
  * Weill et al., [2016] Weill, R., Bekker, A., Smulakovsky, V., Fischer, B., and Gat, O. (2016).  Noise-mediated casimir-like pulse interaction mechanism in lasers.  Optica, 3(2):189–192.
  * Wimmer et al., [2014] Wimmer, K., Nykamp, D. Q., Constantinidis, C., and Compte, A. (2014).  Bump attractor dynamics in prefrontal cortex explains behavioral precision in spatial working memory.  Nature Neuroscience, 17(3):431–439.
  * Wojtak et al., [2021] Wojtak, W., Coombes, S., Avitabile, D., Bicho, E., and Erlhagen, W. (2021).  A dynamic neural field model of continuous input integration.  Biological Cybernetics, 115(5):451–471.
  * Wu et al., [2016] Wu, S., Wong, K. M., Fung, C. A., Mi, Y., and Zhang, W. (2016).  Continuous attractor neural networks: candidate of a canonical model for neural information representation.  F1000Research, 5.
  * Zhang, [1996] Zhang, K. (1996).  Representation of spatial orientation by the intrinsic dynamics of the head-direction cell ensemble: a theory.  Journal of Neuroscience, 16(6):2112–2126.
  * Zylberberg and Strowbridge, [2017] Zylberberg, J. and Strowbridge, B. W. (2017).  Mechanisms of persistent activity in cortical circuits: possible neural substrates for working memory.  Annual Review of Neuroscience, 40:603–627.

See pages ,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15 of Supplementary.pdf