# Forgetting leads to chaos in attractor networks

Ulises Pereira-Obilinovic  Department of Statistics, The University of
Chicago, Chicago, Illinois, United States of America  Center for Neural
Science, New York University, New York, United States of America  Johnatan
Aljadeff  Department of Neurobiology, The University of Chicago, Chicago,
Illinois, United States of America  Section of Neurobiology, Division of
Biological Sciences, University of California San Diego, La Jolla, California,
United States of America  Nicolas Brunel<EMAIL_ADDRESS>Departments
of Statistics and Neurobiology, The University of Chicago, Chicago, Illinois,
United States of America  Departments of Neurobiology and Physics, Duke
University, Durham, North Carolina, United States of America

###### Abstract

Attractor networks are an influential theory for memory storage in brain
systems [1, 2, 3, 4]. This theory has recently been challenged by the
observation of strong temporal variability in neuronal recordings during
memory tasks. In this work, we study a sparsely connected attractor network
where memories are learned according to a Hebbian synaptic plasticity rule.
After recapitulating known results for the continuous, sparsely connected
Hopfield model [5], we investigate a model in which new memories are learned
continuously and old memories are forgotten, using an online synaptic
plasticity rule. We show that for a forgetting time scale that optimizes
storage capacity, the qualitative features of the network’s memory retrieval
dynamics are age-dependent: most recent memories are retrieved as fixed-point
attractors while older memories are retrieved as chaotic attractors
characterized by strong heterogeneity and temporal fluctuations. Therefore,
fixed-point and chaotic attractors co-exist in the network phase space. The
network presents a continuum of statistically distinguishable memory states,
where chaotic fluctuations appear abruptly above a critical age and then
increase gradually until the memory disappears. We develop a dynamical mean
field theory (DMFT) to analyze the age-dependent dynamics and compare the
theory with simulations of large networks. We compute the optimal forgetting
time-scale for which the number of stored memories is maximized. We found that
the maximum age at which memories can be retrieved is given by an instability
at which old memories destabilize and the network converges instead to a more
recent one. Our numerical simulations show that a high-degree of sparsity is
necessary for the DMFT to accurately predict the network capacity. Finally,
our theory provides specific predictions for delay response tasks with aging
memoranda. In particular, it predicts a higher degree of temporal fluctuations
in retrieval states associated with older memories, and it also predicts
fluctuations should be faster in older memories. Our theory of attractor
networks that continuously learn new information at the price of forgetting
old memories can account for the observed diversity of retrieval states in the
cortex, and in particular the strong temporal fluctuations of cortical
activity.

††preprint: APS/123-QED

##  I Introduction

Attractor networks are an influential theory based on recurrent networks for
learning and memory in brain systems [1, 2, 3, 4]. In this theory, memories
correspond to stable patterns of network activity representing the stored
memoranda, which correspond to fixed-point attractor states of the network
dynamics. When each memory is learned, the external input associated with the
memory modifies the network synaptic connectivity thanks to a synaptic
plasticity rule. These synaptic modifications lead to the creation of a fixed-
point attractor of the dynamics, corresponding to a neural representation of
the learned memorandum. In the attractor state, the network activity is
correlated with, but not necessarily identical to, the original external input
to the network. Upon a brief presentation of an external cue correlated with
one of the stored memoranda, the network retrieves the corresponding
memorandum by autonomously relaxing to the associated attractor state, which
permits decoding the memorandum’s identity by downstream circuitry. The theory
parsimoniously reproduces key observations in neuronal responses during memory
tasks - in particcular in reproduces the phenomenon of selective persistent
activity, i.e., the observation that some neurons in multiple cortical areas
have persistent elevated firing rates following the presentation of specific
items to be maintained in memory during the task, but not others [6, 7, 8, 9,
10, 11, 12].

A major challenge to the attractor network scenario as the mechanism
subserving mnemonic representations is the observation of high degree of
temporal irregularity in neuronal activity during delay periods in the
prefrontal cortex [13, 14, 15, 16, 17, 18]. The observed high spike
irregularity and non-stationary of the trial-averaged firing rates challenge
the classic picture of retrieval states as discrete fixed-point attractor
states. Extensions of classic attractor network models [19, 20] successfully
account for the observed spike irregularity. However, these models fail to
reproduce the non-stationarity observed in the firing rates. Various models
have been proposed to account this non-stationarity by encoding memories using
time-varying variables [20, 21], but they cannot reproduce the stable
population coding observed in neuronal recordings during memory tasks [18]
(but see [22]).

Recently, we have proposed an alternative scenario to account for the observed
variability in an attractor network whose learning rules are inferred from in
vivo data [23]. In this scenario, memories are stored by chaotic attractors
(in contrast to fixed-point attractors in classical attractor networks [1, 2,
24]). Neural activity is correlated with the stored pattern at all times, yet
the dynamics are characterized by strong internally generated temporal
fluctuations. The network displays associative memory properties in spite of
strong chaotic fluctuations in neuronal activity. Using a dynamical mean field
theory (DMFT) [25, 26] Tirozzi and Tsodyks predicted the existence of this
chaotic associative memory phase in the sparse version of the Hopfield model
[5]. In both the Tirozzi and Tsodyks model, as well as in networks constrained
by in vivo data [23], all fixed-point attractor memory states transition to
chaos when either a parameter characterizing the overall strength of synaptic
connectivity, or the number of stored patterns, reach a critical value.

In such models, as well as in classical attractor networks [2], when the
number of stored patterns is larger than the network capacity (i.e., maximum
number of patterns the network can store) they undergo catastrophic
forgetting, and all memories are forgotten at once (i.e., no memory can be
retrieved, no matter how close the initial network state is to the stored
memory). Catastrophic forgetting occurs because of the statistical symmetry
between stored patterns. In a large network, when all patterns are identical
and independently distributed, as in classical attractor network models [1, 2,
24] forgetting one pattern is statistically equivalent to forgetting all. The
recipe for fixing catastrophic forgetting is well known - it consists in
introducing an online learning process, in which new memories are written on
top of older memories that are gradually forgotten (sometimes called
palimpsest models [27]). In this models, each memory has a different age,
which breaks the statistical symmetry between patterns, allowing successful
retrieval of the most recent patterns, at the price of a forgetting of older
ones [28, 29, 30, 31]. Online learning in attractor networks has been actively
studied using networks of binary neurons [28, 31, 32, 33, 34, 35, 36, 37, 38],
in which memory states are fixed-point attractors. However, the effect of
online learning rules on the dynamics of memory states remains an open
question in networks of continuous rate units. Building such theory is key for
contrasting theoretical predictions with firing rates recorded in
neurophysiological experiments. Additionally, since networks of rate units can
be obtained as approximations of more realistic networks of spiking neurons
[39], understating the dynamics and learning in such networks is a critical
intermediate step that combines analytically tractability and biological
realism.

Here we study an attractor network of continuous rate units, where each memory
is initially acquired and then progressively degrades according to an online
Hebbian learning rule. We provide a DMFT describing the network dynamics, in
the large $N$, sparse connectivity limit, and compare our theory with
numerical simulations. Remarkably, we show that in this network both fixed-
point and chaotic attractors memory states co-exist in a large range of
forgetting time scales. The pattern’s age determines the nature of its
retrieval state: fixed-point for newer patterns, and chaotic attractor for
older patterns, leading to a continuum of statistically distinguishable memory
states. Analysis of solutions to the DMFT equations reveals the storage
capacity (i.e., maximal age at which memories can be retrieved) and the
forgetting time-scale that maximizes capacity. We show that storage capacity
is given by an instability line, on which old memories become unstable and the
network retrieves a newer one instead. We complement analytical results by
comparison with simulations of large networks (up to 107 neurons). We find
that a good agreement between simulations and theory is obtained only for
extremely sparse networks, i.e., when the average number of connections per
neuron scales logarithmically with the network size. However, the transition
to chaos is well predicted by our theory for denser networks, i.e., when the
average number of connections scales as a power law with the network size.

##  II Model

We consider a recurrent neuronal network composed of $N$ neurons whose input
currents are described by analog variables $h_{i}$, where $i=1,2,\dots,N$. The
instantaneous firing rates of neurons are given by the the input-output single
neuron transfer function (or f-I curve) $\phi(x)=\tanh(x)$. This is a
unrealistic choice from a neurobiological point of view, but has the advantage
of simplicity. The major features presented in this theory are qualitatively
preserved in more realistic networks [23] (see Discussion). The synaptic input
currents obey the standard current-based version of the rate equations

$\frac{{\rm d}{h}_{i}(t)}{{\rm d}t}=-h_{i}(t)+\sum_{i\neq
j}^{N}c_{ij}J_{ij}\phi[h_{j}(t)]+I_{i},$ (1)

where $c_{ij}$ represents the ‘structural connectivity’ of the network, which
is generated as a sparse random Erdős-Rényi graph, i.e., all $c_{ij}$ are
i.i.d. Bernoulli random variables, $c_{ij}=1$ with probability $K/N$ where $K$
is the average in- and out-degree; $J_{ij}$ is the strength of the synapse
connecting the pre-synaptic neuron $j$ to the post-synaptic neuron $i$; and
$I_{i}$ is an external input to neuron $i$. We assume that $1\ll K\ll N$,
which means that the average number of connections is large but much smaller
that the network size. This is a relevant limit for microcircuits in the
brain, where each neuron receives a large number of inputs ($K\sim 1,000$s),
and connection probabilities are small, of order $\sim$10% in cortex[40, 41,
42, 43, 44] and 1% in hippocampus [45].

The synaptic connectivity matrix is assumed to have been structured through
past presentations of external stimuli (patterns) shown to the network. The
stimulus presentation time scale is assumed much longer than the time scale of
network dynamics, and the synaptic connectivity is assumed fixed at the scale
of retrieval of specific patterns. We use the variable $u$ to refer to the
stimulus presentation times, to avoid confusion with $t$ in Eq. (1), where
$u\in\mathbb{Z}$. External inputs to the network consist in a stream of time-
ordered binary patterns $\\{\vec{\eta}^{u}\\},u\in\mathbb{Z}$. We assume that
$\eta_{i}^{u}$ are i.i.d. binary random variables, such that $\eta_{i}^{u}=\pm
1$ with equal probability. We note that while the stored memories are binary,
the network dynamics and therefore the retrieval states are continuous. The
patterns $\eta$ can be interpreted as the population firing rates when the
dynamics in Eq. (1) is driven by a strong random external stimulus $I_{i}$
symmetrically distributed around zero. The synaptic weights are modified by
the patterns in one shot using an online version of the covariance rule [46].
According to such a rule,

$J^{u+1}_{ij}=\rho J^{u}_{ij}+\frac{A}{K}\eta^{u}_{i}\eta^{u}_{j}.$ (2)

In Eq. (2), the first term in the r.h.s. represents a decay term with a
forgetting rate $0<\rho<1$ that prevents synaptic weights from growing
unbounded, while the second term represents a ‘Hebbian’ synaptic modification
that is proportional to the covariance of pre and post-synaptic activity, with
a strength parameter $A$. On one extreme, for $\rho=0$ the learned
connectivity from previous patterns is forgotten in one time step, while in
the case $\rho=1$ there is no forgetting. The case $\rho=1$ is problematic
since it leads to an unbounded growth of the weights when the number of
presented patterns becomes very large. However, for completeness we will
describe the $\rho=1$ scenario, in which a finite number of patterns is
presented, in Section IV.

With such a learning rule, the connectivity matrix is given at any given time
by

$J^{u}_{ij}=\frac{A}{K}\sum_{u^{\prime}=-\infty}^{u-1}\rho^{u-1-u^{\prime}}\eta^{u^{\prime}}_{i}\eta^{u^{\prime}}_{j}.$
(3)

Without loss of generality, we focus on a particular time $u$, drop the index
$u$ for the connectivity matrix $J_{ij}=J_{ij}^{u}$, and rewrite $J_{ij}$ as

$J_{ij}=\frac{A}{K}\sum_{\mu=0}^{\infty}e^{-\frac{\mu}{\tau
K}}\eta^{\mu}_{i}\eta^{\mu}_{j}.$ (4)

where $\mu$ is the age of pattern $\eta^{\mu}_{i}$ (i.e., how far in the past
pattern $\mu$ was presented), and $\tau$ is a forgetting time constant defined
as $\tau=-1/(K\log\rho)$. Note that in the sum in the r.h.s. of Eq. (4), we
have changed the pattern indices for convenience, such that a pattern of age
$\mu$ corresponds to the pattern shown at time $u^{\prime}=u-1-\mu$.

##  III Dynamic mean field theory

###  III.1 Formalism

We start by deriving the DMFT for the network model defined by Eqs. (1, 4), in
the limit of infinitely many neurons $N\to\infty$, synapses per neuron
$K\to\infty$, and sparse connectivity $K/N\to 0$ [47, 5]. Although the network
in Eqs. (1,4) stores binary patterns, i.e., $\eta_{i}^{s}\in\\{-1,1\\}$, here
we present a general DMFT for arbitrary $P(\eta)$, with the constraint that
$\langle\eta\rangle=0$. This ensures that the average change in connection
strength due to learning of a single pattern is zero, and therefore the
synaptic weights average does not grow with the number of stored patterns.
This could be enforced by a homeostatic mechanism that controls the mean
changes in the incoming inputs due to learning [48, 49]. We use similar
methods as in refs. [25, 5], and anticipate that for some parameters the
dynamics of the network will be chaotic. In the above limit, the instantaneous
synaptic inputs to each neuron can then be rewritten of its mean $\nu_{i}$ and
random temporal variations around the mean $\zeta(t)$ due to chaotic dynamics,

$\sum_{i\neq j}^{N}c_{ij}J_{ij}\phi(h_{j}(t))=\nu_{i}+\zeta_{i}(t).$ (5)

Therefore, the network dynamics in Eq. (1) becomes statistically equivalent to
a stochastic differential equation (SDE) given by

$\frac{{\rm d}{h}_{i}(t)}{{\rm d}t}=-h_{i}(t)+\nu_{i}+\zeta_{i}(t).$ (6)

where the $\zeta_{i}$s are uncorrelated stochastic processes, whose
autocorrelation will be computed self-consistently in the following.

As in classical mean field theories for attractor neuronal network models [2,
24], we define the overlaps between network state and the memories stored in
the network as

$m_{\mu}=\frac{1}{N}\sum_{i=1}^{N}\eta_{i}^{\mu}\phi(h_{i})=\langle\eta^{\mu}\phi(h)\rangle_{h,\eta},$
(7)

where $\langle\dots\rangle_{h,\eta}$ represents an average over the statistics
of the stored patterns $\eta$ and the input currents $h$. The network is able
to retrieve a pattern stored in memory if it settles in a state in which there
is an overlap of order one with the corresponding pattern. This means in
particular that the stimulus identity can be retrieved from the network
activity by a linear decoder. On the other hand, if the overlap is zero in the
large $N$ limit, then the memory is forgotten and cannot be retrieved. We
assume that the overlaps do not depend on time (i.e., $m_{\mu}(t)=m_{\mu}$),
which is trivially true for fixed-point attractors, and becomes a good
approximation at large $N$ for chaotic attractor memory states. Except when
stated otherwise (see Section III.4), we will assume that in a given retrieval
state only a single memory is retrieved. With this assumption, the mean
synaptic inputs $\nu_{i}$ to neuron $i$ become

$\nu_{i}=A\eta_{i}^{\mu}e^{-\frac{\mu}{\tau K}}m_{\mu}.$ (8)

Note that in Eq. (8), the mean synaptic input depends exponentially on the age
of the pattern. For most recent memories ($\mu\ll K$) the mean synaptic
current is close to its maximum in magnitude $\nu_{\mu}\approx
A\eta_{i}^{\mu}m_{\mu}$, but it decays exponentially for older memories,
$\mu\sim O(K)$.

We next turn to the fluctuations around the mean. The variable $\zeta_{i}(t)$
is assumed to be a Gaussian random field with mean zero, representing the
variability around the mean synaptic inputs to neuron $i$. Its autocovariance
function is given by

$\text{Cov}\left(\zeta_{i}(t),\zeta_{i}(t+\tau)\right)=\gamma
A^{2}\kappa\langle\phi(h(t))\phi(h(t+\tau))\rangle_{h},$ (9)

where

$\kappa=\frac{1}{K}\sum_{\mu=1}^{\infty}e^{-\frac{2\mu}{\tau K}},$ (10)

and $\gamma=\langle\eta^{2}\rangle^{2}.$ In particular, for binary patterns
considered here, $\gamma=1$. In the large $K$ limit, Eq. (10) leads to
$\kappa=\tau/2$.

Eqs. (8,10) makes it clear that the number of patterns that can be retrieved
is of order $K$. Thus, in the large $N$ limit, we define a continuous age
index $s=\mu/K$, such that $s=0$ for recently learned patterns, while
$s\to\infty$ for patterns seen in the distant past.

For a pattern with age $s$, the dynamics of the currents can be rewritten as

$\frac{{\rm d}h_{i}}{{\rm
d}t}=-h_{i}+A\eta_{i}^{s}e^{-\frac{s}{\tau}}m_{s}+A\sqrt{\frac{\gamma\tau}{2}}y_{i}(t),$
(11)

where $y(t)$ is a Gaussian random field with auto-covariance function

$C_{s}(t^{\prime})=\langle
y(t)y(t+t^{\prime})\rangle_{y}=\langle\phi(h(t))\phi(h(t+t^{\prime})\rangle_{h},$
(12)

where the subscript $s$ in the above equation is here to remind us that the
autocovariance function depends on the age of the retreived pattern. By
defining the local currents
$u_{i}(t)=h_{i}(t)-A\eta^{s}_{i}e^{-\frac{s}{\tau}}m_{s}$, Eqs. (11,7,12) can
be re-written as

$\displaystyle\dot{u}$ $\displaystyle=$
$\displaystyle-u+A\sqrt{\gamma\alpha}y(t)$ (13) $\displaystyle m_{\mu}$
$\displaystyle=$ $\displaystyle\langle\eta\phi(u(t)+A\eta
e^{-\frac{s}{\tau}}m_{s}),\rangle_{u,\eta},$ (14)

while

$C_{s}(t^{\prime})=\langle\phi\left(u(t)+A\eta
e^{-\frac{s}{\tau}}m_{s}\right)\phi\left(u(t+t^{\prime})+A\eta
e^{-\frac{s}{\tau}}m_{s}\right)\rangle_{u,\eta}.$ (15)

As in [25, 50] we introduce the synaptic input current auto-covariance
function

$\Delta_{s}(t^{\prime})=\langle u(t)u(t+t^{\prime})\rangle_{u}.$ (16)

where again $s$ reminds us of the age dependence of this autocovariance
function. Analogously to the derivation in [26, 51], we can derive a self-
consistent equation for the local-field auto-covariance

$\frac{{\rm d}^{2}\Delta_{s}(t^{\prime})}{{\rm
d}\tau^{2}}=\Delta(t^{\prime})-\frac{A^{2}\gamma\tau}{2}C_{s}(t^{\prime}).$
(17)

Using Eq. (6), we can evaluate the age-dependent overlap. Eq. (7) becomes,

$m_{s}=\int D\eta Dx\eta\phi\left(A\left[\sqrt{\Delta_{0s}}x+\eta
e^{-\frac{s}{\tau}}m_{s}\right]\right).$ (18)

We define $\Delta_{0s}$ as the auto-covariance at zero time-lag
$t^{\prime}=0$, i.e., $\Delta_{s}(0)=\Delta_{0s}$, i.e., the variance of the
synaptic input currents $h_{i}$. The auto-covariance in Eq. (15) can be
written as

$C_{s}(t^{\prime})=\int D\eta Dz\left[\int
Dx\phi\left(A\left[\sqrt{\Delta_{0s}-|\Delta_{s}(t^{\prime})|}x+\sqrt{|\Delta_{s}(t^{\prime})|}z+\eta
e^{-\frac{s}{\tau}}m_{s}\right]\right)\right]^{2},$ (19)

where $D\eta=p_{\eta}(\eta)d\eta$, $Dx=dxe^{-x^{2}/2}/\sqrt{2\pi}$, and
$Dz=dze^{-z^{2}/2}/\sqrt{2\pi}$.

We further assume that $\Delta(t^{\prime})\geq 0$, and re-scale
$\Delta(t^{\prime})$, $\Delta(t^{\prime})\to A^{2}\Delta(t^{\prime})$. As in
[25], Eq. (17) can be re-written in terms of a potential

$\frac{{\rm d}^{2}\Delta_{s}}{{\rm d}t^{\prime 2}}=-\frac{\partial
V(\Delta_{s},\Delta_{0s})}{\partial\Delta_{s}},$ (20)

by defining

$V(\Delta_{s},\Delta_{0s})=-\frac{\Delta_{s}^{2}}{2}+\frac{\kappa\gamma}{A^{2}}\int
D\eta Dz\left[\int
Dx\Phi\left(A\left[\sqrt{\Delta_{0s}-|\Delta_{s}|}x+\sqrt{|\Delta_{s}|}z+\eta
e^{-\frac{s}{\tau}}m_{s}\right]\right)\right]^{2},$ (21)

where $\Phi(x)=\int_{0}^{x}dr\phi(r)$.

The dynamics in Eqs. (20, 21) corresponds to equation of particle in a
Newtonian potential $V(\Delta_{s},\Delta_{0s})$ that depends parametrically on
$\Delta_{0s}$. As shown by [25], the motion of the particle should be such
that $\frac{d\Delta_{s}(0)}{dt^{\prime}}=0$ and
$V(\Delta_{0s},\Delta_{0s})=V(\Delta_{s}(\infty),\Delta_{0s})$. Therefore, the
network dynamics in Eq. (11) is described by three order parameters: the
overlap $m_{s}$; the synaptic input current variance $\Delta_{0s}$; and
$\Delta_{1s}$, which is the long time limit of the auto-covariance of the
synaptic input currents
$\Delta_{1s}=\lim_{t^{\prime}\to\infty}\Delta_{s}(t^{\prime})$. These can be
computed self-consistently, as shown below.

As we will see below, there are four qualitatively different types of
solutions to the DMF equations. Solutions with $m_{s}=0$ correspond to
background states (i.e., states that are uncorrelated with all stored
memories), while solutions with $m_{s}>0$ correspond to memory (or retrieval)
states, since there exists a finite overlap between network state and the
corresponding stored pattern. Both types of solutions can be either fixed-
point attractors of the dynamics (i.e., with autocovariance functions that are
constant in time, i.e., $\Delta_{s}(t)=\Delta_{0s}$ at all times) or chaotic
attractors, with an autocovariance function that depends on time. We first
focus our attention on the transitions to chaos, and then study the capacity
of the system (largest age at which retrieval states still exist).

###  III.2 Transitions to chaos

In fixed-point attractors there are no temporal fluctuations in the input
currents. The auto-covariance of the local fields in Eq. (16) is then equal to
the variance of the local currents at all times (i.e.,
$\Delta_{s}(\tau)=\Delta_{0s}$), which leads to

$\displaystyle m_{s}$ $\displaystyle=$ $\displaystyle\int D\eta
Dx\eta\phi\left(A\left[\sqrt{\Delta_{0s}}x+\eta
e^{-\frac{s}{\tau}}m_{s}\right]\right)$ (22) $\displaystyle\Delta_{0s}$
$\displaystyle=$ $\displaystyle\gamma\kappa\int D\eta
Dx\phi^{2}\left(A\left[\sqrt{\Delta_{0s}}x+\eta
e^{-\frac{s}{\tau}}m_{s}\right]\right).$ (23)

The above equations give the overlap in a retrieval fixed point attractor
state of a memory of age $s$. However, these fixed points can destabilize and
become chaotic, leading to chaotic retrieval states. Importantly, as we will
see below in this model the dynamical properties of the attractors (i.e.,
fixed-point or chaotic) strongly depend on the age of the patterns. Eqs.
(22,23) have solutions in parameter space even beyond the transition to chaos,
when the assumption of fixed-point attractor memory states with no temporal
fluctuations in the input currents is no longer valid. However, as expected,
their predictions depart from network simulations (See Fig. 1 dashed line;
Fig. 2 A and B red lines). We refer to the solutions of these equations the
static solutions to DMFT (SMFT).

Analogous to [25], to find the transition to chaos of memory states, it is
necessary to find the point in parameter space where the static solution
$\Delta(t^{\prime})=\Delta_{0}$ becomes unstable. At this point the auto-
covariance of the local-field $\Delta(t^{\prime})$ transitions from stationary
to time-dependent. Since the dependence on time of the auto-covariance of the
local fields is ruled by the Newtonian equation for the position $\Delta$ at
“time” $t^{\prime}$ of a particle subject to a potential energy (Eq. (17)),
finding the transition point is equivalent to finding the critical point
$\Delta_{0}^{\text{chaos}}$ where the potential in Eq. (21) changes its
concavity. After this point, solutions for the auto-covariance of the local
field starting at $\Delta_{0}$ relax to
$\lim_{t^{\prime}\to\infty}\Delta(t^{\prime})\equiv\Delta_{1}$. The transition
point is given by

$A^{2}\gamma\kappa\int D\eta
Dz\left\\{\phi^{{}^{\prime}}\left(A\left[\sqrt{\Delta_{0s}}z+\eta
e^{-\frac{s}{\tau}}m_{s}\right]\right)\right\\}^{2}=1.$ (24)

Equation (24) in addition to Eqs. (22,23) describe the curve in the parameter
space that separates fixed-point from chaotic memory states.

In the background state, all the overlaps with the stored memories are zero,
i.e., $m_{s}=0$. The critical line in the space of parameters for its
transition to chaos is thus given by the following set of equations

$\displaystyle A^{2}\gamma\kappa\int
Dz\left\\{\phi^{{}^{\prime}}\left(A\sqrt{\Delta_{0}}z\right)\right\\}^{2}=1$
(25) $\displaystyle\Delta_{0}=\gamma\kappa\int
Dz\phi^{2}\left(A\sqrt{\Delta_{0}}z\right).$ (26)

###  III.3 Capacity

As there exist two types of memory states, there are two possible ways to
compute the capacity: the maximum age at which fixed-point or chaotic
retrieval states cease to exist, respectively. However, as we will see later,
in a range of parameters of the system the physically relevant capacity is the
one computed from chaotic retrieval states, since fixed point memory states
destabilize before reaching capacity for those parameters. In chaotic
attractors, the potential defined in Eq. (21) is not convex (i.e.,
$\partial^{2}V/\partial\Delta^{2}>0$), and the auto-covariance of the local
currents in Eq. (16) is time dependent. Additionally, as shown by [25, 50], a
chaotic solution is characterized by an aperiodic, decreasing potential. This
corresponds to the condition
$\lim_{t^{\prime}\to\infty}V(\Delta_{s}(t^{\prime}))=V(\Delta_{0s})$, which is
equivalent to

$-\frac{\Delta_{0s}^{2}}{2}+\frac{\kappa\gamma}{A^{2}}\int D\eta
Dx\Phi^{2}\left(A\left[\sqrt{\Delta_{0s}}x+\eta
e^{-\frac{s}{\tau}}m_{s}\right]\right)=-\frac{\Delta_{1s}^{2}}{2}+\frac{\kappa\gamma}{A^{2}}\int
D\eta Dz\left[\int
Dx\Phi\left(A\left[\sqrt{\Delta_{0s}-|\Delta_{1s}|}x+\sqrt{|\Delta_{1s}|}z+\eta
e^{-\frac{s}{\tau}}m_{s}\right]\right)\right]^{2},$ (27)

where $\Delta_{1s}$ corresponds to the large time limit of $\Delta_{s}(t)$,
$\Delta_{s}(t)\xrightarrow{t\to\infty}\Delta_{1s}$. Therefore,
$\left.\frac{\partial V}{\partial\Delta}\right|_{\Delta=\Delta_{1s}}=0$, which
is equivalent to

$\Delta_{1s}=\kappa\gamma\int D\eta Dz\left[\int
Dx\phi\left(A\left[\sqrt{\Delta_{0s}-|\Delta_{1s}|}x+s\sqrt{|\Delta_{1s}|}z+\eta
e^{-\frac{s}{\tau}}m_{s}\right]\right)\right]^{2}.$ (28)

In Eqs. (27,28), the overlap $m_{s}$ is given by Eq. (22). Thus, the three
order parameters $m_{s}$, $\Delta_{0s}$ and $\Delta_{1s}$ can be obtained by
solving numerically Eqs. (22,27,28). For fixed-point attractor memory states,
there are no temporal fluctuations in the input currents and therefore
$\Delta_{1s}=\Delta_{0s}$, recovering from Eqs. (22,27,28) the static
solutions in Eq. (22,23). Beyond the transition to chaos, the overlap curve
for chaotic attractors is given by Eqs. (22,27,28). The capacity of this
network corresponds to the older memory age $s_{c}$ beyond which no older
memory can be retrieved. At capacity, the corresponding memory state has zero
overlap $m_{s_{c}}=0$. Since the transfer function in this model
($\phi(x)=\tanh(x)$) is odd, at capacity we have that $\Delta_{1s_{c}}=0$ (see
Eq. (28)). The value of $\Delta_{0s_{c}}$ is then obtained by solving Eq.
(27), and is given by

$\tau=\frac{(A\Delta_{0s^{c}})^{2}}{\left[\int
Dx\Phi^{2}\left(A\left[\sqrt{\Delta_{0s^{c}}}x\right]\right)-\left(\int
Dx\Phi\left(A\left[\sqrt{\Delta_{0s^{c}}}x\right]\right)\right)^{2}\right]},$
(29)

with

$Ae^{-\frac{s^{c}}{\tau}}\int
Dx\phi^{{}^{\prime}}\left(A\left[\sqrt{\Delta_{0s^{c}}}x\right]\right)=1.$
(30)

Eqs. (29,30) provide the capacity curve in parameter space ($s_{c},\tau,A$)
(see the curve that separates the green and the gray region in Fig. 4). Notice
that Eq. (30) is obtained by derivating Eq. (22) with respect of the overlap
and evaluating the resulting equation at capacity (i.e., $m_{s_{c}}=0$),
considering the overlap changes continuously with the memory age (see Fig.
5A).

###  III.4 Stability of memory states

In this section, we study the stability of a memory state of a given age $s$
with respect to perturbations in neuronal activity in the direction of
another, more recent, memory states. Thus, we assume that the network state is
correlated with two memories: the currently retrieved one, of age $s$, but
also a recent memory ($s=0$). For simplicity we focus on the case of binary
patterns. The aim of this calculation is to derive dynamical equations for the
transient dynamics of the overlaps $m_{0}$ and $m_{s}$, similarly to a recent
calculation for recurrent networks storing sequences of activity [52]. Then,
we can use these equations to study numerically the stability of the memory
state corresponding to the memory $s$, i.e., $m_{s}\sim\mathcal{O}(1)$ and
$m_{0}\ll 1$.

Similarly to Eq. (11), the dynamics of the currents are approximated by a time
dependent Gaussian random field

$\frac{{\rm d}h_{i}}{{\rm
d}t}=-h_{i}+A(\eta_{i}^{0}m_{0}+\eta_{i}^{s}e^{-\frac{s}{\tau}}m_{s})+A\sqrt{\frac{\tau}{2}}y_{i}(t).$
(31)

At any time, the synaptic input current $h_{i}(t)$ can be described by a
Gaussian random variable

$h(t)=A(\eta_{i}^{0}m_{0}+\eta_{i}^{s}e^{-\frac{s}{\tau}}m_{s})+\sqrt{\frac{\tau\Delta_{0s}(t)}{2}}x,$
(32)

where $x$ is Gaussian random variable with mean 0 and variance 1, and the
variance $\Delta_{0s}$ is now time dependent. Using the above in Eq. (1) we
obtain dynamical equations for the overlaps

$\displaystyle\frac{{\rm d}m_{0}}{{\rm d}t}$ $\displaystyle=$ $\displaystyle-
m_{0}+\frac{1}{2}\int
Dx\left[\phi\left(A\left[m_{0}+e^{-\frac{s}{\tau}}m_{s}+\sqrt{\frac{\tau\Delta_{0s}}{2}}x\right]\right)+\phi\left(A\left[m_{0}-e^{-\frac{s}{\tau}}m_{s}+\sqrt{\frac{\tau\Delta_{0s}}{2}}x\right]\right)\right]$
(33) $\displaystyle\frac{{\rm d}m_{s}}{{\rm d}t}$ $\displaystyle=$
$\displaystyle-m_{s}+\frac{1}{2}\int
Dx\left[\phi\left(A\left[m_{0}+e^{-\frac{s}{\tau}}m_{s}+\sqrt{\frac{\tau\Delta_{0s}}{2}}x\right]\right)+\phi\left(A\left[e^{-\frac{s}{\tau}}m_{s}-m_{0}+\sqrt{\frac{\tau\Delta_{0s}}{2}}x\right]\right)\right].$
(34)

The dynamics of the variance $\Delta_{0s}$ and the two-time autocorrelation
function can be derived by a system of integro-differential equations (see
section 2.2 in the supplementary information of [52] for a derivation in a
closely related model). Here, we choose a simpler approach and approximate the
dynamics of $\Delta_{0s}$ and $\Delta_{1s}$ by a system of ODEs (see Appendix
A for more details) whose steady states are given by

$\begin{split}&-\Delta_{0s}^{2}+\frac{\tau}{2A^{2}}\int
Dx\left\\{\Phi^{2}\left(A\left[\sqrt{\Delta_{0s}}x+m_{0}+e^{-\frac{s}{\tau}}m_{s}\right]\right)+\Phi^{2}\left(A\left[\sqrt{\Delta_{0s}}x+e^{-\frac{s}{\tau}}m_{s}-m_{0}\right]\right)\right\\}=-\Delta_{1s}^{2}+\\\
&\frac{\tau}{2A^{2}}\int Dz\left\\{\left[\int
Dx\Phi\left(A\left[\sqrt{\Delta_{0s}-|\Delta_{1s}|}x+\sqrt{|\Delta_{1}|}z+m_{0}+e^{-\frac{s}{\tau}}m_{s}\right]\right)\right]^{2}+\right.\\\
&\left.\left[\int
Dx\Phi\left(A\left[\sqrt{\Delta_{0s}-|\Delta_{1s}|}x+\sqrt{|\Delta_{1s}|}z+e^{-\frac{s}{\tau}}m_{s}-m_{0}\right]\right)\right]^{2}\right\\}\end{split}$
(35)

and

$\begin{split}\Delta_{1s}=\frac{\tau}{4}\int Dz\left\\{\left[\int
Dx\phi\left(A\left[\sqrt{\Delta_{0s}-|\Delta_{1s}|}x+\sqrt{|\Delta_{1s}|}z+m_{0}+e^{-\frac{s}{\tau}}m_{s}\right]\right)\right]^{2}+\right.\\\
\left.\left[\int
Dx\phi\left(A\left[\sqrt{\Delta_{0s}-|\Delta_{1s}|}x+\sqrt{|\Delta_{1s}|}z+e^{-\frac{s}{\tau}}m_{s}-m_{0}\right]\right)\right]^{2}\right\\}.\end{split}$
(36)

Neglecting the dynamics of the $\Delta$ parameters, and expanding the overlap
equations Eqs. (33,34) around the memory state $m_{s}$ and $m_{0}=0$ we obtain
following instability line

$\int
Dx\phi^{{}^{\prime}}\left(A\left[e^{-\frac{s}{\tau}}m_{s}+\sqrt{\frac{\tau\Delta_{0s}}{2}}x\right]\right)=\frac{1}{A}.$
(37)

Our network simulations show that this approximation matches well the network
dynamics (see Fig. 5 and 6) suggesting that the observed memory state
instability is mainly governed by the overlap dynamics in Eqs. (33,34).

##  IV Dynamics in the absence of forgetting

We start by analyzing the dynamics of this network in the limit where there is
no forgetting ($\tau\to\infty$). In this scenario, we need to consider a
finite number of presented patterns $p$, since the variance of synaptic
weights diverges in the limit $p\rightarrow\infty$, and the synaptic
connectivity can be written as

$J_{ij}=\frac{Ac_{ij}}{Nc}\sum_{\mu=1}^{p}\eta^{\mu}_{i}\eta^{\mu}_{j}.$ (38)

Figure 1: Center: Bifurcation diagram for the TT model in the parameter space
spanned by memory load ($\alpha$) vs the learing gain ($A$). Full lines:
Boundaries of the four qualitatively different regions calculated using DMFT.
Red circles: location of the transition to chaos of retrieval states, using
simulations of networks of $N=10^{7}$ units and $K=2\log(N)$ connections per
neuron. Blue circles: storage capacity, computed using simulations. Dashed red
line: Transition to chaos of the background state. Black dashed line: Line on
which (unstable) fixed point retrieval states disappear. Surrounding panels,
labelled I-IV, show representative numerical simulations in the four
qualitatively different regions. In each simulation the network is initialized
close to one of the stored memories. Each plots show the overlaps of network
state with this memory (blue), as well as other memories (green); and the
activity of 10 randomly selected neurons, as a function of time. The networks
parameters $A$ and $\alpha$ are indicated by the location of the bottom left
corner of the corresponding roman numbers. The network parameters for the
surrounding panels are $N=10^{6}$ and $K=2\log(N)$. The parameter values for
$A$ and $\alpha$ for panels I-IV are respectively $A=0.5,2.5,2.5,2.5$ and
$\alpha=11/K,11/K,21/K,33/K$. Figure 2: (A) Overlap vs memory load. Circles:
Overlaps computed from network simulations (average over 10 realizations).
Dashed vertical line: Transition to chaos. Black lines overlaps computed from
DMFT. Red lines: Overlaps in static solutions to DMFT (SMFT) which are
unstable beyond the transition to chaos. (B) Parameters characterizing the
autocovariance function $\Delta_{0s}$ and $\Delta_{1s}$. As in A, full lines
represent solutions from MFT, while circles represent simulation results. (C)
Autocovariance function for $\alpha=0.5429$ (i.e., $\alpha=15/K$) (Black:
DMFT; Yellow: 10 realizations of network simulations). Other parameters,
$N=10^{6}$, $K=2\log(N)$, and $A=5.5$. Figure 3: Effect of connectivity
sparseness on the dynamics of the TT model. (A) Overlap vs memory load, for
different scalings of connection probability with network size $N$. Circles:
Overlap averaged time after transients and over 10 realizations. Horizontal
dashed line: transition to chaos. Black line: Overlap computed from DMFT. (B)
Capacity vs $\log(c)/\log(N)$ (blue circles) and transition to chaos vs
$\log(c)/\log(N)$ (red circles). The capacity is estimated from network
simulations as the minimal memory load for which the overlap is smaller than
0.1. Dashed lines: DMFT red line corresponds to the capacity predicted by the
DMFT. Notice that for $\lim_{N\to\infty}\log(\log(2\log(N))/N)/\log(N)=-1$.
Green square: capacity in the fully connected case [53] which is equal to the
capacity of the Hopfield model (0.14). Networks sizes corresponding to the
scalings
$c=1/{N^{0.3}},1/{N^{0.4}},1/{N^{0.5}},1/{N^{0.6}},1/{N^{0.75}},2\log(N)/N$
are $N=10^{5},2.5\cdot 10^{5},5\cdot 10^{5},10^{6},10^{6},10^{6}$
respectively. The learning gain is $A=5.5$.

Our network model therefore coincides with the model studied by Tirozzi and
Tsodyks [5] (referred to as TT model in the following) who provided an
analytical description of the network qualitative different behaviors in the
sparse coding limit. In the fully connected case, $c_{ij}=1$, it coincides
with the model introduced by Hopfield [54] and studied analytically by Kuhn et
al [53]. In this limit, the sum on the r.h.s. of Eq. (10) becomes a sum over
$p$ patterns, and $\kappa$ becomes equal to the memory load
$\alpha:=\kappa=p/K$, i.e., the number of stored patterns divided by the
average number of connections. Here we recapitulate the analytical results
provided by Tirozzi and Tsodyks and complement them with simulations of large
networks, to investigate how well the theory matches results of such
simulations at various degrees of sparsity (see Figs.1-3).

In the TT model, the network dynamics depends on two parameters, the memory
load $\alpha=p/K$, i.e., the number of stored patterns scaled by the average
number of connections, and $A$, the strength by which the patterns are
imprinted in the connectivity when learned. Figure 1 shows the network
bifurcation diagram derived from the DMFT delineating the four different
dynamical regimes in parameter space, in which the four qualitatively
different types of states described in Section III exist: I) fixed point
background state. II) fixed point memory states (i.e., states with a finite
overlap $m\sim\mathcal{O}(1)$ with the stored patterns); III) chaotic memory
states (i.e., with a finite overlap $m\sim\mathcal{O}(1)$ with stored
patterns); IV) chaotic background state. For small values of the update
strength $A$, the network is weakly coupled, and only the background state, in
which the firing rates of all neurons are equal to zero, exists. In this
regime, memory retrieval is not possible, and the overlap decays to zero when
the network is initialized with any of the stored patterns (see red region I
in Fig. 1). For larger values of $A$ and small memory load, any of the stored
patterns can be retrieved as fixed-point attractor states. When the network is
initialized close to one of the stored patterns, the network goes to a fixed-
point that is correlated with that pattern (see region II in blue in Fig. 1).
Within this parameter region, the red dashed line in Fig. 1 indicates the
transition to chaos of the background state. Below the red dashed line, for
small memory load, the background state is a fixed-point while for larger
memory loads it is a chaotic state. Larger memory loads lead a transition to
chaos of the memory states (green region III in Fig 1). In this regime, the
network dynamics are chaotic, but it retains a finite overlap with the
retrieved memory. Therefore, in spite of chaotic fluctuations, the population
activity is stably correlated with the retrieved pattern. Finally, if the
memory load further increases none of the stored patterns can be retrieved and
all memories are forgotten. This is a well known phenomenon in attractor
neuronal networks called catastrophic forgetting. The value of the memory load
when this happens is called the memory capacity. Beyond the memory capacity,
the only stable attractor is the chaotic background state (gray region in Fig
1). The dashed black line in Fig. 1 corresponds to the memory load at which
(unstable) fixed point retrieval states disappear. Interestingly, this line is
below the line at which chaotic retrieval states disappear (compare dashed
line with solid line between the green and gray region in Fig. 1). Thus,
chaotic fluctuations allow the network to increase the storage capacity of the
network.

Fig. 2 compares the results from DMFT with numerical simulations of a very
large and sparse network ($N=10^{6}$, $K=2\log(N)$). For these parameters the
DMFT is in excellent agreement with numerical simulations, as shown by
comparing the overlaps in memory states (Fig. 2A), parameters of the
autocovariance (Fig. 2B), and the full autocovariance function (Fig. 2C). We
next numerically investigate the validity of our theory for different scalings
of the sparsity with the network size. In networks in which the average number
of connections scale logarithmicaly with the network size (i.e.,
$K\sim\log(N)$ and $c\sim\log(N)/N$), the storage capacity computed from
network simulations is very close to the to the capacity predicted by the DMFT
(see Fig. 3 A and B). Denser scalings gradually decrease the capacity of the
network, and consequently the DMFT capacity predictions become less accurate
(see Fig. 3A and B). The transition to chaos is well predicted by DMFT in
networks in which the connection probability $c=1/N^{a}$ with $a\geq 0.4$. For
less sparse connectivity, retrieval states disappear before the transition to
chaos is reached, and therefore chaotic retrieval states non longer exist. For
the sake of completeness, we also plot in Fig. 3B the capacity in fully
connected networks, that was obtained analytically by Kuhn et al [53].

##  V Dynamics of the network with forgetting

Figure 4: (A) Bifurcation diagram for the network with forgetting, for $A=10$.
Solid lines: transition to chaos, and maximal age at which retrieval states
exist, calculated using the DMFT. Dashed line: Maximal age at which static
retrieval solutions exist (SMFT). Red line: Maximal age at which retrieval
states are stable. (B) Retrieval dynamics of memories of age 1, 5 and 10 for
the same realization of the connectivity matrix ($N=10^{6}$, $K=2\log(N)$,
$A=10$, $\tau=0.5$). These ages correspond to $s=1/K,5/K,10/K$ respectively,
see I-III in panel A.

We now turn to a network where connectivity is learned using an online Hebbian
synaptic plasticity rule, Eq. (2). In this rule, recent patterns are imprinted
in network connectivity using the covariance between pre and post-synaptic
activity, while older patterns are forgotten with a decay time-scale
$\tau=-1/(K\log(\rho))$. We assume that an infinite number of patterns has
been presented to the network, leading to a connectivity matrix given by Eq.
(4).

In this network, the retrieval dynamics of each memory state depends on its
age, and there is a continuum of statistically distinguishable memory states
parameterized by the memory age (see section III). In a given memory state
$s$, the mean input current depends exponentially on the age of the pattern
(see Eq. (11)). Due to the implicit dependence of the input current
autocovariance function on the mean input current (see Eq. (19)), both the
mean and the autocovariance of the synaptic input currents $h_{i}$ vary with
the age of the memory. For sufficiently large $A$, we identify using DMFT
three different regions depending on the age $s$ and the memory time scale
(Fig. 4): A region in which retrieval states are fixed point attractors,
another region in which retrieval states are chaotic attractors, and finally a
region in which no retrieval is possible. Depending on the forgetting time
scale, different scenarios as possible as memories age. For short forgetting
time scales ($\tau<0.34$ in Fig. 4A), most recent memories are retrieved as
fixed point attractors (blue region in Fig. 4A, and Fig. 4B-I). As memories
age, they become unstable (red line in Fig. 4), and the network typically
retrieves one of the most recently learned patterns instead. For larger
forgetting time scales (0.34$<\tau<$0.68 in Fig. 4), most recent memories are
also retrieved as fixed point attractors. However, unlike for shorter
forgetting time scales, fixed point attractor states become chaotic retrieval
states after a certain age (above the line separating blue and green regions
in Fig. 4 see an example in Fig. 4B-II). As memories further age, there is a
second transition line where the chaotic attractor becomes unstable and the
network retrieves one of the most recent memories instead (see Fig. 4B-III).
For even larger forgetting time scales (0.68$<\tau<$1.58), fixed point
attractor states no longer exist, and even the most recent memories are
retrieved as chaotic attractors. Finally, when $\tau>1.58$ the network is no
longer able to retrieve any memory. Note that the maximal capacity (defined as
the maximal age at which memories can be retrieved) is given by the red line
in Fig. 4. Above this line, retrieval states still exist in a finite region of
parameter space, but they are unstable and the network retrieves more recent
patterns instead. Note also that the capacity is optimized at an intermediate
value of $\tau\sim 0.6$ for which both types of retrieval states coexist in
the network phase space (fixed point attractors for recent memories, chaotic
attractors for older ones). For this value of $\tau$, the capacity is about
$s_{max}=0.18$, which is significantly lower than the capacity in the TT
model.

Figure 5: (A) Overlap vs memory age $s$ in the network with online learning.
Green line/circles: overlap with a memory of age $s$ after the network is
initialized close to the corresponding state, calculated from DMFT/ network
simulations, respectively. Blue line/circles: overlap with a recent memory
that is retrieved instead of the memory of age $s$. Red line: overlap with a
memory of age $s$ in the unstable retrieval state (DMFT) In simulations,
averages are over time and over 10 realizations. Vertical dashed gray line:
transition to chaos. Vertical dashed blue line: Age at which memories become
unstable. Red line: overlap with memory of age $s$ (DMFT) Blue line: overlap
with most recent memories, $s=0$ (DMFT). (B) Autocovariance parameters vs
memory age. Green line/circles: $\Delta_{0s}=\Delta_{1s}$ computed from
DMFT/simulations respectively. Blue line/circles: $\Delta_{00}=\Delta_{10}$
computed from DMFT/simulations respectively. Black/red line:
$\Delta_{0s}/\Delta_{1s}$ in unstable retrieval state (DMFT) (C) Dynamics of
overlaps for a memory of age $s=5/K=0.18$. The memory is retrieved
successfully. (D) Dynamics of overlaps for a memory of age $s=8/K=0.29$. The
memory is not retrieved, as the network goes instead in the attractor state
corresponding to the most recent memory. In both C and D overlaps are color-
coded by age. Network parameters: $N=10^{6}$, $K=2\log(N)$, $\tau=0.64$, and
$A=4$.  Figure 6: (A) Overlap vs memory age ($s$) for different values of the
forgetting time scale $\tau$, and $A=4$. Full lines: Overlap computed using
DMFT in stable retrieval states. Dashed lines: unstable retrieval states.
Circles: Simulations (average over time and over 10 realizations). (B)
Bifurcation diagram for the network with forgetting, for $A=4$. Black line:
maximal age at which retrieval states exist, calculated using DMFT. Blue
line/circles: transition to chaos of retrieval states (DMFT/simulations). Red
line/circles: maximal age at which retrieval states are stable
(DMFT/simulations). Simulation parameters: $N=10^{6}$, $K=2\log(N)$.

Numerical simulations in large very sparsely connected networks show a good
agreement with DMFT, as shown in Fig. 5 in a network with $A=4$ and
$\tau=0.64$. Fig. 5A,B shows that both the overlap with the retrieved pattern,
and the variance $\Delta_{0s}$ decay with age. At the age of about $s=0.18$,
retrieval states become unstable, and the network retrieves one of the most
recent memories instead. Fig. 5C,D show examples of successful retrieval (Fig.
5C), and of an unsuccessful retrieval ending in the retrieval of the most
recent memory instead (Fig. 5D). Fig. 6 shows how the overlaps with retrieved
memories depend on the forgetting time constant $\tau$. It shows that the
overlaps in retrieval states of recent memories decay with $\tau$, since
increasing $\tau$ means older patterns provide increasing interference with
the retrieval of recent patterns. Finally, Fig. 6 shows that the location of
the two instability lines of retrieval states (in red, the instability towards
more recent memories, and in blue, the instability towards chaotic retrieval
states) predicted by DMFT are in good agreement with numerical simulations.

##  VI Discussion

Figure 7: Schematics of the attractor networks phase space. (A) Extensive
transition to chaos of fixed-point attractor memory states in an attractor
network with no forgetting. When the memory load $\alpha$ surpasses a critical
point all the fixed-point attractors transition to chaos at once. For a given
value of the memory load $\alpha$, all fixed-point or chaotic memory states
are statistically equivalent. (B) In an attractor network that forgets memory
states are heterogeneous. Most recent memories are stored as fixed-point
attractors with larger basin of attraction while older memories are chaotic
with smaller basins of attraction. The chaotic temporal fluctuations become
faster for older memories. (C) As predicted by our stability analysis, after a
critical memory age, older memory states become unstable and the network jumps
to the most recent memory state (see light blue arrows).

In the present work, we have characterizez using dynamical mean field theory
the dynamics of an attractor network that learns and forgets through an online
Hebbian learning rule, and shown that in this network fixed-point and chaotic
memory states co-exist. In such a network, recent memories are strongly
imprinted in the connectivity while older memory are exponentially forgotten.
In the limit of no forgetting, the network undergoes a global transition to
chaos, where all memory states transition to chaos at once when the number of
stored patterns surpasses a critical value (see Fig. 7A). When forgetting
takes place, memory states are heterogeneous and form a continuum of
statistically different memory states. In a broad range of forgetting time
scales, recent memories are stored in the network as fixed-point attractors,
while older memories as chaotic attractors (see Fig. 7B). The magnitude of the
chaotic fluctuations increase parametrically with the memory age and the basin
of attractions shrink with age. Interactions between most recent and older
memories destabilize older memory states (see Fig. 7C). Our stability analysis
explains this effect and accurately predicts the capacity curve, i.e., the
curve in parameter space when older memory state destabilize. We contrast our
results with simulation of large networks and find that our theory matches
well the network dynamics, provided synaptic connectivity is sparse enough.
Overall, our work puts forward a network mechanism that might explain the
diversity of neuronal responses in the cortex during memory tasks, and
provides an analytical framework for dissecting the network mechanisms of such
heterogeneity.

###  VI.1 Relation with other chaotic networks

Random networks of rate units are a popular class of tractable models for
exploring the neural mechanisms underlying the temporal fluctuations and
heterogeneity of the firing rate patterns observed in the brain. Transition to
chaos has been extensively studied in randomly connected networks [25, 55, 26,
56] and the analysis has been extended to multiple scenarios such as networks
with multiple populations [57, 50, 58], networks of clusters of recurrently
connected neurons [59], networks driven by external stimuli [60, 61], and
random networks partially structured with low-rank connectivity [62, 63, 64].
The effect of chaos in the computational capabilities of recurrent networks
has also been studied extensively [65, 66, 67, 61, 68]. In the vast majority
of the above studies the focus has been in the transition to chaos of a single
fixed-point, the background (zero average activity) state. A notable exception
is the work by [5]. In this work by Tirozzi and Tsodyks, they analytically
study the highly diluted version of the Hopfield model in a network of
continuous ‘firing rate’ units. They derived a DMFT that fully characterizes
the networks dynamics. Their DMFT is a particular case of our model for
infinite forgetting time scale $\tau\to\infty$ and number of patterns
proportional to $K$, i.e., $p=\alpha K$. The DMFT for the network in [5] is
recovered in our equations for $\alpha=\kappa=\tau/2$ and $s=0$. Although in
[5] computed analytically the phase diagram of the model, no comparison with
network simulations were provided. In Section IV of the present work, we
complemented the DMFT in [5] with a systematic comparison between numerical
solutions of the DMFT and simulations of large networks. We show that a good
agreement between simulations and theory is only obtained for very sparse
networks, in which $c\sim 1/N^{a}$ and $a>0.8$. (see Section VI.4). Our work
differs from the above studies in that in our network there exist a
coexistence between a continuum of statistically distinguishable fixed-point
and chaotic memory states (see Fig. 7 B and C).

###  VI.2 Relation with other palimpsest networks

The problem of catastrophic forgetting was recognized soon after the Hopfield
model was proposed, and this issue was addressed using various online learning
rules [29, 27, 28] that lead to gradual forgetting of old memories, thereby
successfully avoiding catastrophic forgetting. This class of models are
sometimes referred to as ‘palimpsest’ models. In all these models the price to
pay for avoiding catastrophic forgetting is a substantial decrease in storage
capacity. This phenomenon is also present in our model, where the capacity is
only about $s_{max}=0.18$ for $A=10$, to be compared with a capacity of order
0.9 in the TT model.

Online learning rules has also been widely studied in networks with binary
synapses [30, 69, 31, 34, 32, 70, 36, 37, 35, 38]. In these models, synapses
have two states, one with strong, the other with weak (or zero) efficacy.
Transitions between states are induced stochastically by the activity of pre
and post-synaptic neurons. Under general conditions, these networks also
behave as palimpsests, maintaining a memory of recently shown patterns but
gradually forgetting patterns shown in the past. Most studies of such networks
have used an ideal observer approach to estimate the maximal age at which a
pattern can still be retrieved from synaptic connectivity [31, 32, 70, 37,
38]. Storage capacity has also been investigated in attractor networks with
synapses, using either binary neurons [70, 36, 35], or firing rate units [34].
We expect that the results presented here about the diversity of retrieval
states should hold in sparsely connected networks with such synapses, but the
capacity should depend drastically on the implementation of the Markov process
describing the transition between states.

###  VI.3 Relation to random networks with low-rank structure

Recently, low-rank networks have been used as a framework for modeling the
dynamics and computations in neuronal circuits in the brain [62, 63, 71, 72,
64]. In these networks, the connectivity is comprised of a low-rank matrix
plus a random matrix with i.i.d. Gaussian entries, i.e.,
$J_{ij}=\sum_{s=1}^{p}\frac{A}{N}\eta_{i}^{s}\xi_{j}^{s}+X_{ij}$ where
$X_{ij}\overset{i.i.d.}{\sim}N(0,g^{2}/N)$ and $p\sim\mathcal{O}(1)$.
Attractor networks and low-rank networks are very similar from a mean field
perspective. In the attractor network analyzed here, for a given memory state
$s$, in the large $N$ and $K$ limit, the result of our DMFT is equivalent to
the following low-rank network
$\frac{A}{N}\eta_{i}^{s}\eta_{j}^{s}e^{-s/\tau}+\frac{A\tau}{2\sqrt{N}}X_{ij}$,
where $X_{ij}\overset{i.i.d.}{\sim}N(0,1)$. Therefore, the mean input current
due to the low-rank component in low-rank networks is equivalent to the mean
input current of the retrieved memory in attractor networks, while the random
component is equivalent to the quenched noise that arises from the stored
memories (see [72]). Because these networks are similar in the mean field
limit, we hypothesize that by introducing heterogeneity in the strength of
low-rank components, the co-existence between chaotic and fixed-point
attractors displayed in our network can also be obtained in low-rank networks.

###  VI.4 Network sparsity

The DMFT derived here, as well as the ones in [73, 5], can be obtained using a
path integral formalism, expanding the generating functional in powers of
$c=K/N\ll 1$, and then neglecting the $\mathcal{O}(K^{2}/N^{2})$ corrections
[47]. However, the conditions on $c$ and $N$ for which these higher order
corrections can be neglected remained unclear. Our numerical simulations show
that our DMFT accurately predicts the network dynamics in very sparse networks
with $c\sim\log(N)/N$ (Fig. 3). For denser connectivity, the agreement between
capacity predicted by theory and simulations degrades progressively (see Fig.
3 A and B). On the other hand, the transition to chaos is well predicted by
theory up to connection probabilities $c\sim 1/N^{0.4}$, at which point the
chaotic region disappears altogether. These results are consistent with
results on sparse networks of binary neurons where the disruptive effect of
correlations between neurons in the capacity can be avoided when $K\ll\log(N)$
[74].

###  VI.5 Attractor networks with learning rules inferred from data

We have recently studied an attractor network similar to the one presented
here, but with both transfer function and synaptic plasticity rules inferred
from _in vivo_ data [23]. We showed that this network display qualitatively
similar dynamics as the TT model. For small learning gain $A$ and memory loads
$\alpha$ memory states are fixed-point attractors, while for large $A$ and
$\alpha$ memory states are chaotic (see Fig. 4 and Fig. 6 respectively in
[23]). A detailed study of this more realistic network with an online learning
rule is outside the scope of this paper. Building and analyzing more realistic
attractor network models with learning rules inferred from _in vivo_ data
should help to clarify the network mechanisms underlying the observed changes
in neural activity across days during learning novel images [75, 76].

###  VI.6 Experimental predictions

To investigate the neural mechanisms underlying the maintenance and
manipulation of information stored in memory, neuroscientists have designed a
broad class of tasks named delay response tasks. In these tasks, subjects must
keep in memory information about a previously presented stimulus during a
delay period, to be able to produce a behavioral response that depends on this
information. Neuronal recordings in delay response tasks show that a
(typically small proportion) of neurons display elevated tonic activity during
the delay period [6, 7, 8, 9]. The observed persistent elevated delay activity
is consistent with attractor network models in which fixed-point attractors
correspond the memory states [77, 78]. However, in many of these recordings,
neurons exhibit firing patterns with high degree of temporal variation. These
observations challenge attractor networks as viable models of memory storage
[79]. Interestingly, population analysis of neural activity during delay
response tasks shows that in spite of the variability observed at a single
neuron level there is a stable population encoding of the memoranda during the
delay periods [18]. A possible scenario compatible with the above observations
[23, 80] is that the single neuron temporal variability and stable population
encoding can be explained by attractor networks in the chaotic regime as
described here for the TT model in Section IV (see Fig. 1; green region and
Fig. 7A). In an attractor network with forgetting, memory states are
heterogeneous and as a memory ages, its activity ranges from fixed-point to
chaotic. Additionally, for chaotic retrieval states, the temporal fluctuations
of their activity gets faster with age (see Fig. 6B). Our theory thus predicts
that for delay response tasks with aging memoranda: 1) retrieval states
corresponding to newer memories will present a large proportion of neurons
with persistent delay activity consistent with fixed-point attractor dynamics;
2) temporal fluctuations will become more pronounced, and faster, as memories
become older. These predictions can be tested by recording neural activity
during delay response tasks performed across the learning and forgetting
process of multiple items.

##  VII Acknowledgments

UP thanks to the Swartz Fundation for its support and to the Champaign Public
Library for providing the physical space where an early version of this work
was developed. NB was supported by R01MH115555, R01NS112917 and ONR
N00014-17-1-3004.

##  Appendix A Relaxational dynamics for $\Delta_{0s}$ and $\Delta_{1s}$

For computing the order parameters $m_{0}$, $m_{s}$, $\Delta_{0s}$ and
$\Delta_{1s}$ in panels Fig. 5 A (green and blue lines) and B (red and yellow
lines), we replace the full integro-differential dynamics of the variance
$\Delta_{0s}$ and the two-time autocorrelation function (see section 2.2 in
the supplementary information of [52] for the corresponding equations in a
closely related model) by a relaxational dynamics for $\Delta_{0s}$ and
$\Delta_{1s}$ given by

$\begin{split}\frac{{\rm d}\Delta_{1s}}{{\rm
d}t}=-\Delta_{1s}+\frac{\tau}{4}\int Dz\left\\{\left[\int
Dx\phi\left(A\left[\sqrt{\Delta_{0s}-|\Delta_{1s}|}x+\sqrt{|\Delta_{1s}|}z+m_{0}+e^{-\frac{s}{\tau}}m_{s}\right]\right)\right]^{2}+\right.\\\
\left.\left[\int
Dx\phi\left(A\left[\sqrt{\Delta_{0s}-|\Delta_{1s}|}x+\sqrt{|\Delta_{1s}|}z+e^{-\frac{s}{\tau}}m_{s}-m_{0}\right]\right)\right]^{2}\right\\}\end{split}$
(39)

and

$\begin{split}&\frac{{\rm d}\Delta_{0s}}{{\rm
d}t}=-\Delta_{0s}+\sqrt{\Delta_{1s}^{2}+\psi_{1}(m_{0},m_{s},\Delta_{0s},\Delta_{1s})-\psi_{2}(m_{0},m_{s},\Delta_{0s},\Delta_{1s})},\end{split}$
(40)

with

$\begin{split}\psi_{1}(m_{0},m_{s},\Delta_{0s},\Delta_{1s})=\frac{\tau}{2A^{2}}\int
Dx\left\\{\Phi^{2}\left(A\left[\sqrt{\Delta_{0s}}x+m_{0}+e^{-\frac{s}{\tau}}m_{s}\right]\right)+\Phi^{2}\left(A\left[\sqrt{\Delta_{0s}}x+e^{-\frac{s}{\tau}}m_{s}-m_{0}\right]\right)\right\\}\end{split}$
(41)

$\begin{split}&\psi_{2}(m_{0},m_{s},\Delta_{0s},\Delta_{1s})=\frac{\tau}{2A^{2}}\int
Dz\left\\{\left[\int
Dx\Phi\left(A\left[\sqrt{\Delta_{0s}-|\Delta_{1s}|}x+\sqrt{|\Delta_{1}|}z+m_{0}+e^{-\frac{s}{\tau}}m_{s}\right]\right)\right]^{2}+\right.\\\
&\left.\left[\int
Dx\Phi\left(A\left[\sqrt{\Delta_{0s}-|\Delta_{1s}|}x+\sqrt{|\Delta_{1s}|}z+e^{-\frac{s}{\tau}}m_{s}-m_{0}\right]\right)\right]^{2}\right\\}.\end{split}$
(42)

Notice that the dynamics for the overlaps is given by Eqs. (33, 34).
Therefore, the full dynamics for the order parameters
$m_{0},m_{s},\Delta_{0s},\Delta_{1s}$ is ruled by a system of integro-
differential equations given by Eqs. (39, 40, 33, 34). We initialize the
dynamics using several initial conditions of the the form

$\displaystyle m_{0}$ $\displaystyle=$ $\displaystyle\delta_{m_{0}}$ (43)
$\displaystyle m_{s}$ $\displaystyle=$ $\displaystyle 1$ (44)
$\displaystyle\Delta_{0s}$ $\displaystyle=$ $\displaystyle d_{0}$ (45)
$\displaystyle\Delta_{1s}$ $\displaystyle=$ $\displaystyle d_{1},$ (46)

where $d_{0}$ and $d_{1}$ are fixed and $\delta_{m_{0}}$ is small and
variable. The steady state solutions are shown in Fig. 5 A and B. This method
makes an accurate prediction for the memory age $s^{*}$ above which networks
jump form an older to the most recent memory state (see Fig. 5 A and B).

##  Appendix B Network simulations and numerical solutions to the DMFT

The network simulations and numerical solutions to the DMFT were generated
using custom Python scripts. This code is available at GitHub repository
https://github.com/ulisespereira/chaos-forgetting-palimpsest.

## References

  * Hopfield [1982] J. J. Hopfield, Neural networks and physical systems with emergent collective computational abilities, Proceedings of the national academy of sciences 79, 2554 (1982).
  * Amit  _et al._ [1985] D. J. Amit, H. Gutfreund, and H. Sompolinsky, Spin-glass models of neural networks, Physical Review A 32, 1007 (1985).
  * Amit [1992] D. J. Amit,  _Modeling brain function: The world of attractor neural networks_ (Cambridge University Press, 1992).
  * Brunel [2005] N. Brunel, Network models of memory, in  _Methods and Models in Neurophysics, Volume Session LXXX: Lecture Notes of the Les Houches Summer School 2003_ , edited by C. Chow, B. Gutkin, D. Hansel, C. Meunier, and J. Dalibard (Elsevier, 2005) pp. 407–476.
  * Tirozzi and Tsodyks [1991] B. Tirozzi and M. Tsodyks, Chaos in highly diluted neural networks, EPL (Europhysics Letters) 14, 727 (1991).
  * Fuster  _et al._ [1971] J. M. Fuster, G. E. Alexander, _et al._ , Neuron activity related to short-term memory, Science 173, 652 (1971).
  * Miyashita [1988] Y. Miyashita, Neuronal correlate of visual associative long-term memory in the primate temporal cortex, Nature 335, 817 (1988).
  * Funahashi  _et al._ [1989] S. Funahashi, C. J. Bruce, and P. S. Goldman-Rakic, Mnemonic coding of visual space in the monkey’s dorsolateral prefrontal cortex, Journal of neurophysiology 61, 331 (1989).
  * Goldman-Rakic [1995] P. S. Goldman-Rakic, Cellular basis of working memory, Neuron 14, 477 (1995).
  * Liu  _et al._ [2014] D. Liu, X. Gu, J. Zhu, X. Zhang, Z. Han, W. Yan, Q. Cheng, J. Hao, H. Fan, R. Hou, _et al._ , Medial prefrontal activity during delay period contributes to learning of a working memory task, Science 346, 458 (2014).
  * Guo  _et al._ [2014] Z. V. Guo, N. Li, D. Huber, E. Ophir, D. Gutnisky, J. T. Ting, G. Feng, and K. Svoboda, Flow of cortical activity underlying a tactile decision in mice, Neuron 81, 179 (2014).
  * Inagaki  _et al._ [2017] H. K. Inagaki, L. Fontolan, S. Romani, and K. Svoboda, Discrete attractor dynamics underlying selective persistent activity in frontal cortex, Biorxiv , 203448 (2017).
  * Compte  _et al._ [2003] A. Compte, C. Constantinidis, J. Tegnér, S. Raghavachari, M. Chafee, P. S. Goldman-Rakic, and X.-J. Wang, Temporally irregular mnemonic persistent activity in prefrontal neurons of monkeys during a delayed response task, J. Neurophysiol. 90, 3441 (2003).
  * Shafi  _et al._ [2007] M. Shafi, Y. Zhou, J. Quintana, C. Chow, J. Fuster, and M. Bodner, Variability in neuronal activity in primate cortex during working memory tasks, Neuroscience 146, 1082 (2007).
  * Barak  _et al._ [2010] O. Barak, M. Tsodyks, and R. Romo, Neuronal population coding of parametric working memory, J. Neurosci. 30, 9424 (2010).
  * Barak and Tsodyks [2014] O. Barak and M. Tsodyks, Working models of working memory, Curr. Opin. Neurobiol. 25, 20 (2014).
  * Kobak  _et al._ [2016] D. Kobak, W. Brendel, C. Constantinidis, C. E. Feierstein, A. Kepecs, Z. F. Mainen, X. L. Qi, R. Romo, N. Uchida, and C. K. Machens, Demixed principal component analysis of neural population data, Elife 5 (2016).
  * Murray  _et al._ [2017] J. D. Murray, A. Bernacchia, N. A. Roy, C. Constantinidis, R. Romo, and X.-J. Wang, Stable population coding for working memory coexists with heterogeneous neural dynamics in prefrontal cortex, Proceedings of the National Academy of Sciences 114, 394 (2017).
  * Barbieri and Brunel [2007] F. Barbieri and N. Brunel, Irregular persistent activity induced by synaptic excitatory feedback, Frontiers in Computational Neuroscience 1, 5 (2007).
  * Mongillo  _et al._ [2008] G. Mongillo, O. Barak, and M. Tsodyks, Synaptic theory of working memory, Science 319, 1543 (2008).
  * Lundqvist  _et al._ [2010] M. Lundqvist, A. Compte, and A. Lansner, Bistable, irregular firing and population oscillations in a modular attractor memory network, PLoS Comput. Biol. 6, e1000803 (2010).
  * Druckmann and Chklovskii [2012] S. Druckmann and D. B. Chklovskii, Neuronal circuits underlying persistent representations despite time varying activity, Current Biology 22, 2095 (2012).
  * Pereira and Brunel [2018] U. Pereira and N. Brunel, Attractor dynamics in networks with learning rules inferred from in vivo data, Neuron 99, 227 (2018).
  * Tsodyks and Feigel’Man [1988] M. Tsodyks and M. Feigel’Man, The enhanced storage capacity in neural networks with low activity level, EPL (Europhysics Letters) 6, 101 (1988).
  * Sompolinsky  _et al._ [1988] H. Sompolinsky, A. Crisanti, and H.-J. Sommers, Chaos in random neural networks, Physical Review Letters 61, 259 (1988).
  * Crisanti and Sompolinsky [2018] A. Crisanti and H. Sompolinsky, Path integral approach to random neural networks, Physical Review E 98, 062120 (2018).
  * Nadal  _et al._ [1986] J. Nadal, G. Toulouse, J. Changeux, and S. Dehaene, Networks of formal neurons and memory palimpsests, EPL (Europhysics Letters) 1, 535 (1986).
  * Parisi [1986] G. Parisi, A memory which forgets, Journal of Physics A: Mathematical and General 19, L617 (1986).
  * Mézard  _et al._ [1986] M. Mézard, J. Nadal, and G. Toulouse, Solvable models of working memories, Journal de physique 47, 1457 (1986).
  * Tsodyks [1990] M. Tsodyks, Associative memory in neural networks with binary synapses, Modern Physics Letters B 4, 713 (1990).
  * Amit and Fusi [1994] D. J. Amit and S. Fusi, Learning in neural networks with material synapses, Neural Computation 6, 957 (1994).
  * Fusi  _et al._ [2005] S. Fusi, P. J. Drew, and L. F. Abbott, Cascade models of synaptically stored memories, Neuron 45, 599 (2005).
  * Fusi and Abbott [2007] S. Fusi and L. Abbott, Limits on the memory storage capacity of bounded synapses, Nature neuroscience 10, 485 (2007).
  * Romani  _et al._ [2008] S. Romani, D. J. Amit, and Y. Amit, Optimizing one-shot learning with binary synapses, Neural computation 20, 1928 (2008).
  * Dubreuil  _et al._ [2014] A. M. Dubreuil, Y. Amit, and N. Brunel, Memory capacity of networks with stochastic binary synapses, PLoS computational biology 10, e1003727 (2014).
  * Huang and Amit [2011] Y. Huang and Y. Amit, Capacity analysis in multi-state synaptic models: a retrieval probability perspective, Journal of computational neuroscience 30, 699 (2011).
  * Lahiri and Ganguli [2013] S. Lahiri and S. Ganguli, A memory frontier for complex synapses, in  _Advances in neural information processing systems_ (2013) pp. 1034–1042.
  * Benna and Fusi [2016] M. K. Benna and S. Fusi, Computational principles of synaptic memory consolidation, Nature neuroscience 19, 1697 (2016).
  * Ostojic and Brunel [2011] S. Ostojic and N. Brunel, From spiking neuron models to linear-nonlinear models, PLoS Comput Biol 7, e1001056 (2011).
  * Mason  _et al._ [1991] A. Mason, A. Nicoll, and K. Stratford, Synaptic transmission between individual pyramidal neurons of the rat visual cortex in vitro, Journal of Neuroscience 11, 72 (1991).
  * Markram  _et al._ [1997] H. Markram, J. Lübke, M. Frotscher, A. Roth, and B. Sakmann, Physiology and anatomy of synaptic connections between thick tufted pyramidal neurones in the developing rat neocortex., The Journal of physiology 500, 409 (1997).
  * Holmgren  _et al._ [2003] C. Holmgren, T. Harkany, B. Svennenfors, and Y. Zilberter, Pyramidal cell communication within local networks in layer 2/3 of rat neocortex, The Journal of physiology 551, 139 (2003).
  * Thomson and Lamy [2007] A. M. Thomson and C. Lamy, Functional maps of neocortical local circuitry, Frontiers in neuroscience 1, 2 (2007).
  * Lefort  _et al._ [2009] S. Lefort, C. Tomm, J.-C. F. Sarria, and C. C. Petersen, The excitatory neuronal network of the c2 barrel column in mouse primary somatosensory cortex, Neuron 61, 301 (2009).
  * Guzman  _et al._ [2016] S. J. Guzman, A. Schlögl, M. Frotscher, and P. Jonas, Synaptic mechanisms of pattern completion in the hippocampal ca3 network, Science 353, 1117 (2016).
  * Sejnowski [1977] T. J. Sejnowski, Storing covariance with nonlinearly interacting neurons, Journal of mathematical biology 4, 303 (1977).
  * Kree and Zippelius [1987] R. Kree and A. Zippelius, Continuous-time dynamics of asymmetrically diluted neural networks, Physical Review A 36, 4421 (1987).
  * Toyoizumi  _et al._ [2014] T. Toyoizumi, M. Kaneko, M. P. Stryker, and K. D. Miller, Modeling the dynamic interaction of hebbian and homeostatic plasticity, Neuron 84, 497 (2014).
  * Vogels  _et al._ [2011] T. Vogels, H. Sprekeler, F. Zenke, C. Clopath, and W. Gerstner, Inhibitory plasticity balances excitation and inhibition in sensory pathways and memory networks, Science 334, 1569 (2011).
  * Kadmon and Sompolinsky [2015] J. Kadmon and H. Sompolinsky, Transition to chaos in random neuronal networks, Physical Review X 5, 041030 (2015).
  * Schücker  _et al._ [2016] J. Schücker, S. Goedeke, D. Dahmen, and M. Helias, Functional methods for disordered neural networks, arXiv preprint arXiv:1605.06758 (2016).
  * Gillett  _et al._ [2020] M. Gillett, U. Pereira, and N. Brunel, Characteristics of sequential activity in networks with temporally asymmetric hebbian learning, Proceedings of the National Academy of Sciences 117, 29948 (2020).
  * Kühn  _et al._ [1991] R. Kühn, S. Bös, and J. L. van Hemmen, Statistical mechanics for networks of graded-response neurons, Physical Review A 43, 2084 (1991).
  * Hopfield [1984] J. J. Hopfield, Neurons with graded response have collective computational properties like those of two-state neurons, Proc. Natl. Acad. Sci. U.S.A. 81, 3088 (1984).
  * Wainrib and Touboul [2013] G. Wainrib and J. Touboul, Topological and dynamical complexity of random neural networks, Physical review letters 110, 118101 (2013).
  * Engelken  _et al._ [2020] R. Engelken, F. Wolf, and L. Abbott, Lyapunov spectra of chaotic recurrent neural networks, arXiv preprint arXiv:2006.02427 (2020).
  * Aljadeff  _et al._ [2015] J. Aljadeff, M. Stern, and T. Sharpee, Transition to chaos in random networks with cell-type-specific connectivity, Physical review letters 114, 088101 (2015).
  * Harish and Hansel [2015] O. Harish and D. Hansel, Asynchronous rate chaos in spiking neuronal circuits, PLoS Comput Biol 11, e1004266 (2015).
  * Stern  _et al._ [2014] M. Stern, H. Sompolinsky, and L. Abbott, Dynamics of random neural networks with bistable units, Physical Review E 90, 062710 (2014).
  * Rajan  _et al._ [2010] K. Rajan, L. Abbott, and H. Sompolinsky, Stimulus-dependent suppression of chaos in recurrent neural networks, Physical Review E 82, 011903 (2010).
  * Schuecker  _et al._ [2018] J. Schuecker, S. Goedeke, and M. Helias, Optimal sequence memory in driven random networks, Physical Review X 8, 041029 (2018).
  * Mastrogiuseppe and Ostojic [2018] F. Mastrogiuseppe and S. Ostojic, Linking connectivity, dynamics, and computations in low-rank recurrent neural networks, Neuron 99, 609 (2018).
  * Landau and Sompolinsky [2018] I. D. Landau and H. Sompolinsky, Coherent chaos in a recurrent neural network with structured connectivity, bioRxiv , 350801 (2018).
  * Landau and Sompolinsky [2021] I. D. Landau and H. Sompolinsky, Macroscopic fluctuations emerge in balanced networks with incomplete recurrent alignment, Physical Review Research 3, 023171 (2021).
  * Toyoizumi and Abbott [2011] T. Toyoizumi and L. Abbott, Beyond the edge of chaos: Amplification and temporal integration by recurrent networks in the chaotic regime, Physical Review E 84, 051908 (2011).
  * Sussillo and Abbott [2009] D. Sussillo and L. F. Abbott, Generating coherent patterns of activity from chaotic neural networks, Neuron 63, 544 (2009).
  * Poole  _et al._ [2016] B. Poole, S. Lahiri, M. Raghu, J. Sohl-Dickstein, and S. Ganguli, Exponential expressivity in deep neural networks through transient chaos, Advances in neural information processing systems 29, 3360 (2016).
  * Keup  _et al._ [2021] C. Keup, T. Kühn, D. Dahmen, and M. Helias, Transient chaotic dimensionality expansion by recurrent networks, Physical Review X 11, 021064 (2021).
  * Amit and Fusi [1992] D. J. Amit and S. Fusi, Constraints on learning in dynamic synapses, Network: Computation in Neural Systems 3, 443 (1992).
  * Amit and Huang [2010] Y. Amit and Y. Huang, Precise capacity analysis in binary networks with multiple coding level inputs, Neural computation 22, 660 (2010).
  * Schuessler  _et al._ [2020] F. Schuessler, A. Dubreuil, F. Mastrogiuseppe, S. Ostojic, and O. Barak, Dynamics of random recurrent networks with correlated low-rank structure, Physical Review Research 2, 013111 (2020).
  * Beiran  _et al._ [2021] M. Beiran, A. Dubreuil, A. Valente, F. Mastrogiuseppe, and S. Ostojic, Shaping dynamics with multiple populations in low-rank recurrent networks, Neural Computation 33, 1572 (2021).
  * Tsodyks [1988] M. Tsodyks, Associative memory in asymmetric diluted network with low level of activity, EPL (Europhysics Letters) 7, 203 (1988).
  * Derrida  _et al._ [1987] B. Derrida, E. Gardner, and A. Zippelius, An exactly solvable asymmetric neural network model, EPL (Europhysics Letters) 4, 167 (1987).
  * Huang  _et al._ [2018] G. Huang, S. Ramachandran, T. S. Lee, and C. R. Olson, Neural correlate of visual familiarity in macaque area v2, Journal of Neuroscience 38, 8967 (2018).
  * Garrett  _et al._ [2020] M. Garrett, S. Manavi, K. Roll, D. R. Ollerenshaw, P. A. Groblewski, N. D. Ponvert, J. T. Kiggins, L. Casal, K. Mace, A. Williford, _et al._ , Experience shapes activity dynamics and stimulus coding of vip inhibitory cells, Elife 9, e50340 (2020).
  * Amit and Brunel [1997] D. J. Amit and N. Brunel, Model of global spontaneous activity and local structured activity during delay periods in the cerebral cortex., Cerebral cortex 7, 237 (1997).
  * Wang [1999] X.-J. Wang, Synaptic basis of cortical persistent activity: the importance of nmda receptors to working memory, Journal of Neuroscience 19, 9587 (1999).
  * Lundqvist  _et al._ [2018] M. Lundqvist, P. Herman, and E. K. Miller, Working memory: Delay activity, yes! persistent activity? maybe not, Journal of Neuroscience 38, 7013 (2018).
  * Aljadeff  _et al._ [2021] J. Aljadeff, M. Gillett, U. Pereira-Obilinovic, and N. Brunel, From synapse to network: models of information storage and retrieval in neural circuits, Current opinion in neurobiology 70, 24 (2021).