# Active oscillatory associative memory

Matthew Du  Department of Chemistry, University of Chicago, Chicago, Illinois
60637, USA  The James Franck Institute, University of Chicago, Chicago,
Illinois 60637, USA  Agnish Kumar Behera  Department of Chemistry, University
of Chicago, Chicago, Illinois 60637, USA  Suriyanarayanan Vaikuntanathan
Department of Chemistry, University of Chicago, Chicago, Illinois 60637, USA
The James Franck Institute, University of Chicago, Chicago, Illinois 60637,
USA

(August 28, 2024)

###### Abstract

Traditionally, physical models of associative memory assume conditions of
equilibrium. Here, we consider a prototypical oscillator model of associative
memory and study how active noise sources that drive the system out of
equilibrium, as well as nonlinearities in the interactions between the
oscillators, affect the associative memory properties of the system. Our
simulations show that pattern retrieval under active noise is more robust to
the number of learned patterns and noise intensity than under passive noise.
To understand this phenomenon, we analytically derive an effective energy
correction due to the temporal correlations of active noise in the limit of
short correlation decay time. We find that active noise deepens the energy
wells corresponding to the patterns by strengthening the oscillator couplings,
where the more nonlinear interactions are preferentially enhanced. Using
replica theory, we demonstrate qualitative agreement between this effective
picture and the retrieval simulations. Our work suggests that nonlinearity in
the oscillator couplings can improve memory under nonequilibrium conditions.

An essential feature of biological organisms is memory, which involves the
learning and retrieval of information. The importance of this cognitive
ability has prompted many studies of memory from a physics perspective [1, 2],
specifically associative, or content-addressable, memory [3, 4, 5, 6].

A well explored model of associative memory is the Hopfield neural network
[7], which resembles a fully connected Ising model with nonuniform and random
interaction strengths. The neurons are modeled as spins and the learned
patterns, which are drawn from a probability distribution, are encoded in the
couplings between spins. Naturally, the memory capabilities of the Hopfield
model can be understood from the perspective of statistical mechanics [8, 9,
10, 11]. Refs. [8, 9, 10] showed that thermal fluctuations in the firing rates
of neurons leads to poorer recall of patterns.

However, it is unclear whether such findings, which characterize the system at
equilibrium, apply to biological systems which are expected to be subjected to
noise sources that do not satisfy the constraints of equilibrium. Indeed, a
proper physical description of living systems requires nonequilibrium aspects
[12, 13]. Notably, the flocking of birds [14, 15] and the dynamics of the
cytoskeleton [16, 17] are modeled as active matter [18, 19, 20], particles
that consume energy from the environment to power their movement. Some active
materials, e.g., epithelial cells [21] and colloids in a bacterial solution
[22], can be described as particles subjected to so-called active noise [23,
24], which has temporal correlations (i.e., is colored noise) and does not
obey detailed balance with respect to the damping. In fact, the membrane
potential of neurons can have temporally correlated fluctuations [25, 26].

Additionally, a feature that the Hopfield model does not exhibit is the
oscillatory activity of neurons [27, 28, 29]. This motivation led to variants
[30, 31, 32, 33, 34, 35, 36] of the Kuramoto model [37] of nonlinearly coupled
oscillators, whose Hopfield-like interactions give rise to associative memory.
These analog systems, due to their simplicity, can be studied in much detail
using analytical methods [31, 32, 36].

Here, we study an oscillatory Hopfield network in the presence of active
noise. Inspired by the Kuramoto-type models of [32, 38], the system here has
not been studied before, to the best of our knowledge, even in the absence of
noise. Our numerical simulations reveal that the recall of information is
enhanced under active noise, relative to under comparable equilibrium
conditions. We then analytically derive, in a perturbative limit, an effective
energy correction associated with the temporal correlations of active noise
[23, 39, 40, 41]. Using replica theory [42], we find that the effective
picture describes qualitatively the improvement of associative memory under
active noise. Our work sheds light on how nonlinearity in the neuron couplings
impacts memory under nonequilibrium conditions.

The paper is structured as follows. We introduce the Kuramoto-type Hopfield
model in Sec. I. Then, we numerically simulate pattern retrieval in Sec. II.
Sec. III describes the theoretical analysis, in which we derive the effective
Hamiltonian (Sec. III.1) and apply replica theory to it (Sec. III.2). We
conclude in Sec. IV.

##  I Model

We consider $N$ periodic oscillators whose interactions endow the system with
associative memory (Fig. 1a). Specifically, we focus on the Hamiltonian

$\displaystyle H$
$\displaystyle=-\frac{1}{2N}\sum_{\mu=1}^{p}\sum_{i,j=1}^{N}\xi_{i}^{\mu}\xi_{j}^{\mu}\cos\theta_{i}\cos\theta_{j}$
$\displaystyle\quad-\frac{\epsilon}{4N}\sum_{i,j=1}^{N}\cos 2\theta_{i}\cos
2\theta_{j},$ (1)

where $\theta_{1},\dots,\theta_{N}$ are the oscillator phases. Analogous in
form to the Hamiltonian of [38], $H$ is a generalization of the Kuramoto model
[37] that features associative memory. Collectively, the phases and their
cosines will be referred to as $\bm{\theta}=(\theta_{1},\dots,\theta_{N}$) and
$\mathbf{c}(\bm{\theta})=(\cos\theta_{1},\dots,\cos\theta_{N})$, respectively.

Figure 1: Oscillatory associative memory. (a) $N$ coupled oscillators,
represented as a graph whose vertices are the oscillator phases $\bm{\theta}$
and the edges are the couplings. Here, $N=5$. (b) A pattern, which is a vector
$\bm{\xi}\in\\{\pm 1\\}^{N}$. Here, $\bm{\xi}=(+1,-1,+1,+1,-1)$. (c) A pattern
configuration, which for pattern $\bm{\xi}$ is any $\bm{\theta}^{\bm{\xi}}$
such that $\cos\theta_{i}^{\bm{\xi}}=\xi_{i}$. Here,
$\bm{\theta}^{\bm{\xi}}=(0,\pi,0,0,\pi)$. (d) Pattern retrieval. The system
starts slightly displaced from a target pattern $\bm{\xi}$
($\bm{\theta}=\bm{\theta}^{\bm{\xi}}+\bm{\delta}$ for nonzero $\bm{\delta}$)
and then approaches the pattern ($\bm{\theta}\approx\bm{\theta}^{\bm{\xi}}$).

The first term of $H$ encodes $p$ patterns (Fig. 1b)
$\bm{\xi}^{1},\dots,\bm{\xi}^{p}\in\\{\pm 1\\}^{N}$ learned by the Hebb rule
[43, 7]. Appearing in the model of Ref. [32], this term is analogous to the
Hopfield Hamiltonian [7], where the variable of spin $i$ is replaced by
$\cos\theta_{i}$ and the couplings are of the form
$\cos\theta_{i}\cos\theta_{j}$. All pattern elements are independent and
identically distributed (i.i.d.) random variables drawn uniformly from $\\{\pm
1\\}$ (i.e., a binomial distribution). When the pattern loading $\alpha=p/N$
is sufficiently small, the system has energy minima near the pattern
configurations (Fig. 1c) $\\{\bm{\theta}^{\bm{\xi}^{\mu}}\\}_{\mu=1}^{p}$ ,
where $\mathbf{c}(\bm{\theta}^{\bm{\xi}^{\mu}})=\bm{\xi}^{\mu}$. Henceforth,
we will sometimes write “pattern” as shorthand for “pattern configuration.”

The second term of $H$ enhances the stability of the learned patterns, as well
as all other patterns, whose configurations also satisfy
$\mathbf{c}(\bm{\theta})\in\\{\pm 1\\}^{N}$. Characterized by strength
$\epsilon$, this term has couplings of the form $\cos 2\theta_{i}\cos
2\theta_{j}=(2\cos^{2}\theta_{i}-1)(2\cos^{2}\theta_{j}-1)$, making it more
nonlinear than the pattern-encoding term. By varying $\epsilon$ in our
calculations below, we can understand how nonlinearity affects memory in the
presence of noise. We will restrict ourselves to intermediate values of
$\epsilon$, as the learned patterns are further stabilized (compared to
$\epsilon=0$) while the patterns that have not been learned remain relatively
unstable [38].

In pattern retrieval, the system starts near a pattern $\bm{\xi}\in\\{\pm
1\\}^{N}$ and then moves towards it (Fig. 1d). The initial position is
$\bm{\theta}(0)=\bm{\theta}^{\bm{\xi}}+\bm{\delta}$, where
$\mathbf{c}(\bm{\theta}^{\bm{\xi}})=\bm{\xi}$, $\bm{\delta}$ is a vector of
$N$ i.i.d. random variables drawn from $\mathcal{N}(0,\sigma^{2})$, and
$\sigma$ determines the typical size of the initial offset. To represent
neurons isolated from external forces, one would evolve the system via
gradient descent on the energy landscape of Hamiltonian (1). However, real
neurons reside in a biological environment, which can induce nonequilibrium
fluctuations in the firing rates. Thus, we evolve the system according to the
dynamics of active Ornstein-Uhlenbeck particles [23, 24]:

$\frac{d\bm{\theta}}{dt}=-\nabla_{\bm{\theta}}H+\bm{\eta}(t).$ (2)

The noise $\bm{\eta}$ is a stationary Gaussian process that satisfies

$\displaystyle\langle\eta_{i}(t)\rangle$ $\displaystyle=0,$ (3)
$\displaystyle\langle\eta_{i}(t)\eta_{j}(t^{\prime})\rangle$
$\displaystyle=\frac{T}{\tau}\delta_{ij}e^{-|t-t^{\prime}|/\tau},$ (4)

where the temporal correlations exponentially decay with persistence time
$\tau$ and initial value $T/\tau$. For $\tau>0$, the noise is active, and the
system is kept out of equilibrium. In the limit $\tau\rightarrow 0$, the
temporal correlations satisfy

$\lim_{\tau\rightarrow
0}\langle\eta_{i}(t)\eta_{j}(t^{\prime})\rangle=2T\delta_{ij}\delta(t-t^{\prime}),$
(5)

corresponding to thermal noise which drives the system towards equilibrium
with temperature $T$. Hereafter, we refer to this type of noise as passive. To
characterize the proximity of the system to $\bm{\xi}$, the order parameter
that we use is the overlap
$\frac{1}{N}\left[\bm{\xi}\cdot\mathbf{c}(\bm{\theta})\right]$, which ranges
from 0 (farthest from pattern) to 1 (at pattern). After the system has evolved
for a long time, the overlap indicates the accuracy with which the pattern has
been retrieved.

##  II simulations

We simulate pattern retrieval for a system with $N=200$ oscillators. Without
loss of generality, we consider the retrieval of learned pattern $\mu=1$
($\bm{\xi}^{1}$), hereafter referred to as “target pattern.” We choose
$\sigma=0.1$, for which the typical initial overlap with the target pattern is
$\approx 0.995$. The system is evolved from $t=0$ to $t=800$ by numerically
integrating equation of motion (2) using a time step of $\Delta t=0.004$.
Below, the simulation results represent an average over 16 pattern sets and 8
initial configurations. These numbers are high enough such that the error due
to finite sampling can be neglected.

We first focus on the system with just the pattern-encoding term, i.e., with
$\epsilon=0$ (Fig. 2a). In the absence of noise ($T=0$), the learned patterns
can be retrieved at higher pattern loading ($\alpha$) compared to the Hopfield
model, for which retrieval fails at $\alpha\gtrsim 0.138$ [8, 9]. However, the
critical pattern loading quickly decays to zero as passive noise ($T\neq 0$,
$\tau\rightarrow 0$) is added. In fact, the system is less robust to noise
than the Hopfield model, for which patterns can be retrieved up to $T=1$ [8,
9].

Figure 2: Retrieval phase diagram under passive noise, as calculated by
simulation, for stability-enhancing coupling (a) $\epsilon=0$ and (b)
$\epsilon=1$. In the retrieval phase R${}_{\text{L+NL}}$, both learned and
not-learned patterns can be retrieved (i.e., with final overlap > 0.8). In the
retrieval phase R${}_{\text{L}}$, only learned patterns can be retrieved. In
the no retrieval (NR) phase, no patterns can be retrieved.

When the stability-enhancing term [second term of Eq. (1)] is turned on
($\epsilon>0$, Fig. 2b), the retrievability of the learned patterns
deteriorates less rapidly with temperature and pattern loading. In particular,
retrieval is enabled at higher temperatures for which it does not occur when
$\epsilon=0$. These improvements come at the expense of being able to, at
lower temperatures, also retrieve patterns that have not been learned (Fig.
2b, $\text{R}_{\text{L+NL}}$ phase). For this reason, we restrict ourselves to
$\epsilon\leq 1$ throughout this work.

Having characterized the system under passive noise, we now consider the
system under active noise of equal strength ($T$) and with persistence time
$\tau=1$. Compared to passive noise, pattern retrieval is more resistant to
pattern loading and noise strength (Fig. 3, a-b). Notably, retrieval is
successful at noise strengths for which it fails under passive noise. The
critical pattern loading at these noise strengths grows with $\epsilon$ (Fig.
3c). Similarly, for a given pattern loading, active noise raises the critical
noise strength (i.e., the noise strength beyond which retrieval does not
occur), and the relative increase scales with $\epsilon$ (Fig. 3d). These
findings suggest that higher nonlinearity in the oscillator couplings allows
active noise to better stabilize pattern configurations which are unstable
under passive noise.

we

Figure 3: (a, b) Retrieval phase diagram under active noise ($\tau=1$, red),
as calculated by simulation, for stability-enhancing coupling (a) $\epsilon=0$
and (b) $\epsilon=1$. See the caption of Fig. 2 for a description of the
various phases. For comparison, the transition lines under passive noise
(reproduced from Fig. 2) are shown as blue translucent points. (c) Critical
pattern loading (for the $\text{R}_{\text{L}}\rightarrow\text{NR}$ transition)
under active noise ($\alpha_{c}^{(\text{a})}$) as a function of $\epsilon$.
Results are shown for various noise strengths ($T$) at which patterns cannot
be retrieved under passive noise (i.e., $\alpha_{c}^{(\text{p})}=0$, as
indicated by gray dashed line). (d) Relative difference ($\Delta
T_{c}/T_{c}^{(\text{p})}$, where $\Delta
T_{c}=T_{c}^{(\text{a})}-T_{c}^{(\text{p})}$) in the critical noise strength
(for the $\text{R}_{\text{L}}\rightarrow\text{NR}$ transition) between active
($T_{c}^{(\text{a})}$) and passive ($T_{c}^{(\text{p})}$) noise as a function
of $\epsilon$. Results are shown for various pattern loadings ($\alpha$).

##  III Theory

To gain further insight on the improvement of memory due to active noise, we
employ a perturbative expansion with respect to persistence time $\tau$ [23,
39]. This approach yields an effective equilibrium steady state, characterized
by temperature $T$ and a modified Hamiltonian. As we will show below, for the
system studied here, the temporal correlations can be recast as an effective
renormalization of coupling strengths.We then study the associated effective
Hamiltonian using replica theory [42].

###  III.1 Effective equilibrium

Consider the steady-state probability distribution $P(\bm{\theta})$ of system
configurations $\bm{\theta}$ in the limit of small persistence time, $\tau\ll
1$. Following the method of either [23] or [39], we can write $P(\bm{\theta})$
in powers of $\tau$ and keep terms up to $O(\tau)$, yielding a Boltzmann-like
distribution $P(\bm{\theta})\propto\exp(-H_{\text{eff}}/T)$, which is
determined by an effective Hamiltonian with the general form [23]

$H_{\text{eff}}=H+\tau\left(\frac{1}{2}\left|\nabla_{\bm{\theta}}H\right|^{2}-T\nabla_{\bm{\theta}}^{2}H\right).$
(6)

Corresponding to an effective equilibrium steady state, this result can be
alternatively obtained using the unified colored noise approximation [40, 41].
Plugging $H$ [Eq. (1)] in the Laplacian term gives

$\displaystyle H_{\text{eff}}$
$\displaystyle=-\frac{\zeta_{1}}{2N}\sum_{\mu=1}^{p}\sum_{i,j=1}^{N}\xi_{i}^{\mu}\xi_{j}^{\mu}\cos\theta_{i}\cos\theta_{j}$
$\displaystyle\quad-\frac{\zeta_{2}\epsilon}{4N}\sum_{i,j=1}^{N}\cos
2\theta_{i}\cos 2\theta_{j}$
$\displaystyle\quad+\frac{\tau}{2}\left|\nabla_{\bm{\theta}}H\right|^{2}$
$\displaystyle\quad+\tau
T\left\\{\frac{1}{N}\sum_{\mu=1}^{p}\left[\bm{\xi}^{\mu}\cdot\mathbf{s}(\bm{\theta})\right]^{2}+\frac{2\epsilon}{N}\left[\mathbf{1}\cdot\mathbf{s}(2\bm{\theta})\right]^{2}\right\\},$
(7)

where $\zeta_{1}=1+2\tau T$, $\zeta_{2}=1+8\tau T$, and
$\mathbf{s}(\bm{\theta})=(\sin\theta_{1},\dots,\sin\theta_{N})$. While it is
possible to write $\frac{\tau}{2}\left|\nabla_{\bm{\theta}}H\right|^{2}$ more
explicitly in terms of $\bm{\theta}$, we have simply chosen not to do so,
since the current expression is sufficient for physical interpretation (see
below). Previously, we applied this small-$\tau$ expansion to the Hopfield
model under active noise [44]. However, due to the constraint on the norm of
the spin vector, we needed to make additional approximations to arrive at the
effective Hamiltonian.

Note that the $\frac{\tau}{2}\left|\nabla_{\bm{\theta}}H\right|^{2}$ term of
$H_{\text{eff}}$ [third line of Eq. (7)] is the squared norm of the
(conservative) force ($-\nabla_{\bm{\theta}}H$), which is zero at the energy
minima near the patterns. Similarly, the fourth term of $H_{\text{eff}}$
[fourth line of Eq. (7)] is zero at the patterns. Thus, these last two terms
should play a minor role in pattern retrieval, which is determined by the
energy landscape in the neighborhood of the patterns.

We therefore expect pattern retrieval under active noise to be approximately
described by the effective Hamiltonian

$\displaystyle H_{\text{eff}}$
$\displaystyle=-\frac{\zeta_{1}}{2N}\sum_{\mu=1}^{p}\sum_{i,j=1}^{N}\xi_{i}^{\mu}\xi_{j}^{\mu}\cos\theta_{i}\cos\theta_{j}$
$\displaystyle\quad-\frac{\zeta_{2}\epsilon}{4N}\sum_{i,j=1}^{N}\cos
2\theta_{i}\cos 2\theta_{j},$ (8)

which has the same form as the original Hamiltonian [Eq. (1)], except with the
pattern-encoding term (first term) scaled by $\zeta_{1}$ and the stability-
enhancing term scaled by $\zeta_{2}$. Since $\zeta_{1},\zeta_{2}>1$, active
noise stabilizes the learned patterns (compared to passive noise) by deepening
their respective energy wells. In particular, as implied by
$\zeta_{2}/\zeta_{1}>1$, the stability-enhancing term is amplified relative to
the pattern-encoding term. More generally, the pattern-stabilizing
interactions with more nonlinearity (i.e., the stability-enhancing term) are
strengthened over those with less (i.e., the pattern-stabilizing term). This
broader view suggests the importance of nonlinearity for memory under
nonequilibrium settings. The impact of active noise on memory, as described by
$H_{\text{eff}}$ of Eq. (8), are summarized in Fig. 4. The limit of passive
noise is recovered by setting $\tau=0$, which yields the original Hamiltonian
[Eq. (1)].

Figure 4: Schematic illustration of the mechanism, as described by effective
Hamiltonian $H_{\text{eff}}$ [Eq. (8)], by which active noise improves memory
compared to passive noise. The temporal correlations of active noise
effectively strengthen the oscillator interactions, where the more nonlinear
interactions are preferentially enhanced. This strengthening leads to a deeper
energy well at the patterns, and thus their retrieval becomes more robust to
noise intensity and number of learned patterns.

Note that the “active stabilization” described by $H_{\text{eff}}$ of Eq. (8)
arises from the $-\tau T\nabla_{\bm{\theta}}^{2}H$ term of the general form of
$H_{\text{eff}}$ [Eq. (6)] rather than the
$\frac{\tau}{2}\left|\nabla_{\bm{\theta}}H\right|^{2}$ term. The former energy
correction, which depends on the curvature of the potential energy surface,
has been shown to stabilize certain particle configurations in a related
system, whose couplings go as the cosine of relative orientations [39]. The
result here stands in contrast to our previous work [44] on the continuous-
spin Hopfield network of [45] under active noise: there, the interactions
scale with the number $N\gg 1$ of spins in such a way that the the $-\tau
T\nabla_{\bm{\theta}}^{2}H$ term is negligible compared to the
$\frac{\tau}{2}\left|\nabla_{\bm{\theta}}H\right|^{2}$ term.

###  III.2 Replica calculation

Next, we apply replica theory [42] to the effective Hamiltonian
$H_{\text{eff}}$ of Eq. (8). Details of the calculation are presented in
Appendix A. Replica theory allows one to calculate order parameters, such as
the overlap with a pattern, of a system in the thermodynamic limit
($N\rightarrow\infty$) and averaged over randomly generated parameters of the
Hamiltonian (known as “quenched disorder” in the literature [9]), such as the
learned patterns. From the order parameters, one can determine what phase the
system is in, e.g., ferromagnetic (retrieval) and spin glass (no retrieval).
Following the convention of [8, 9, 11], a retrieval phase is defined by
nonzero overlap with the target pattern (Section A.5). In contrast, our
simulated results above (Figs. 2-3) use the criterion that the overlap exceeds
0.8 (caption of Fig. 2). Thus, our analysis below is limited to qualitative
comparisons between the phase diagrams obtained by replica theory and
simulation. Furthermore, our replica calculation is limited to $T>0.1$; for
$T\rightarrow 0$, we anticipate that the saddle-point equations [Eqs.
(27)-(32)] which determine the order parameters can yield unstable solutions
(Sec. A.3).

We first consider the system under passive noise, for which $H_{\text{eff}}$
reduces to $H$ [Eq. (1)], the original Hamiltonian. To our knowledge, $H$ has
not been previously studied using replica theory. Figs. 5a and 5b show the
phase diagram calculated by this method for $\epsilon=0$ and $\epsilon=1$,
respectively. Both plots qualitatively resemble the corresponding phase
diagram for the Hopfield model [8, 9, 11]. At $T\lesssim 0.5$, increasing
pattern loading induces transitions from the ferromagnetic retrieval phase
(i.e., a global free energy minimum has significant overlap with the target
pattern) to the mixed retrieval phase (i.e., the global minimum becomes a
local minimum) to the spin glass phase, which does not support retrieval. The
same phase transitions can be triggered by increasing $T$. At $T\gtrsim 0.5$,
the target pattern cannot be retrieved at all; the system is in the
paramagnetic phase at low pattern loading and the spin glass phase at high
pattern loading.

Figure 5: Retrieval phase diagram under passive noise, as calculated by
theory, for stability-enhancing coupling (a) $\epsilon=0$ and (b)
$\epsilon=1$. In the ferromagnetic (F) retrieval phase, the global minimum of
free energy function $f$ [Eq. (23)] has nonzero overlap with the target
pattern. In the mixed (M) retrieval phase, a local (but not global) minimum of
$f$ has nonzero overlap with the target pattern. No retrieval occurs in the
spin glass (SG) and paramagnetic (P) phases, in which all minima of $f$ have
zero overlap with the target pattern. The gray points indicate where no minima
of $f$ could be found numerically (see Appendix A.4). Note that the plots
start at $T=0.1$, in contrast to the corresponding results from simulations
(Fig. 2).

Consistent with simulations (Fig. 2), replica theory (Fig. 5) reveals that the
stability-enhancing term ($\epsilon>0$) raises the critical temperature and
pattern loading at which the system transitions from a retrieval (mixed) to a
no-retrieval (spin glass) phase. Also increased are the critical values
separating the two retrieval phases.

We proceed to the corresponding phase diagrams for the system under active
noise. Since effective equilibrium theory is correct to $O(\tau)$ (Section
III.1), our application of replica theory to $H_{\text{eff}}$ of Eq. (8)
assumes $\tau=0.2$: this value is small enough for the theory to be a
reasonable approximation but large enough to observe appreciable differences
between active and passive noise. Fig. 6a-b shows phase diagrams obtained by
the theory. Replacing passive with active noise enables the target pattern to
be retrieved at higher noise strengths (for a given pattern loading) and
pattern loadings (for a given noise strength). In addition, the relative
increase in the critical noise strengths scales with $\epsilon$ (Fig. 6c).
Remarkably, these features also appear in the simulated phase diagrams (Fig.
3, a-b, d) for $\tau=1$, which falls outside the range of persistence times
($\tau\ll 1)$ for which effective equilibrium theory is expected to be valid.
This agreement suggests that $H_{\text{eff}}$ of Eq. (8) can qualitatively
describe how active noise improves memory, even for $\tau\geq 1$. As discussed
in Section III.1, the impact of active noise in this effective picture is to
strengthen the interactions that stabilize the patterns, where the
interactions with higher nonlinearity are enhanced over those with less. The
replica calculation further reveals that active noise delays the transition
(with respect to increasing noise strength and pattern loading) from the
ferromagnetic retrieval phase to the mixed retrieval phase.

Figure 6: (a, b) Retrieval phase diagram under active noise ($\tau=0.2$, red),
as calculated by theory, for stability-enhancing coupling (a) $\epsilon=0$ and
(b) $\epsilon=1$. The dotted lines indicate where phase transitions occur
under passive noise (compare to corresponding panels of Fig. 5). See the
caption of Fig. 5 for additional details. (c) Relative difference ($\Delta
T_{c}/T_{c}^{(\text{p})}$, where $\Delta
T_{c}=T_{c}^{(\text{a})}-T_{c}^{(\text{p})}$) in the critical noise strength
(for the $\text{M}\rightarrow\text{SG}$ transition) between active
($T_{c}^{(\text{a})}$) and passive ($T_{c}^{(\text{p})}$) noise, as calculated
by theory and as a function of $\epsilon$. Results are shown for various
pattern loadings ($\alpha$).

Before concluding, let us briefly discuss the quantitative aspects of the
theory. In Appendix A.6, we evaluate how well the theory can reproduce the
final overlap computed by simulations. For the case of passive noise, the
agreement is semi-quantitative. For active noise, we observe similar agreement
between the theory and simulations with the same persistence time
($\tau=0.2$), though only at relatively low $T$. More details can be found in
Appendix A.6.

##  IV Conclusion

We studied how active noise impacts periodic oscillators whose couplings endow
the system with associative memory. One type of interaction encodes patterns,
while the second type, a more nonlinear interaction, enhances the stability of
the patterns. Through simulations, we showed that pattern retrieval under
active noise is more robust to noise strength and number of learned patterns
than under passive noise. This improvement in memory scales with the strength
of the more nonlinear coupling relative to the less nonlinear coupling.

Using effective equilibrium theory, which is perturbative with respect to the
persistence time of active noise, we analytically derived an effective
Hamiltonian [Eq. (8)] that captures how the temporal correlations of active
noise resculpt the potential energy landscape and thereby improve memory.
Specifically, active noise leaves the oscillator interactions unchanged in
form but increased in strength, with the more nonlinear interaction being
preferentially enhanced. As a result of the stronger couplings, the energy
basins at the pattern configurations become deeper. Unlike the Hopfield model
studied in [45, 44], the effect of active noise here arises mainly from the
curvature of the original energy surface rather than the slope.

We then demonstrated qualitative agreement between this effective picture and
the simulated results by applying replica theory to the effective Hamiltonian.
This agreement is remarkable because effective equilibrium theory strictly
holds for persistence time $\tau\ll 1$ while the simulations were run with
$\tau=1$. In addition, the replica analysis reveals that replacing passive
with active noise can increase the critical noise strength at which a target
pattern is converted from a global minimum of the (effective) free energy to a
local minimum.

Overall, our findings highlight that having nonlinear couplings between
neurons can be essential for associative memory under nonequilibrium
conditions and thus biologically relevant contexts. For the Hopfield model in
the absence of noise, highly nonlinear interactions have been shown to
drastically increase the critical number of learned patterns beyond which the
patterns can no longer be retrieved [46, 47]. Our study suggests that such
nonlinearities may confer additional benefits to memory in noisy and
nonequilibrium environments.

###### Acknowledgements.

M. D. and S. V. are supported by DOE BES Grant No. DE-SC0019765. A. K. B. is
supported by a fellowship from the Department of Chemistry at the University
of Chicago. We thank Carlos Floyd, David Martin, and Aditya Nandy for their
feedback on earlier drafts of the manuscript.

##  Appendix A Details of replica calculation

In this section, we apply replica theory [42] to effective Hamiltonian
$H_{\text{eff}}$ of Eq. (8). Our replica calculation largely follows those [8,
11] for the Hopfield Hamiltonian. Note that the first term of $H_{\text{eff}}$
[Eq. (8)] has the same form as the Hamiltonian of the Hopfield model.

###  A.1 Gauge fixing

To simplify the calculation, we first choose a gauge for the patterns. Without
loss of generality, we arbitrarily designate pattern 1 as the target pattern,
$\bm{\xi}\equiv\bm{\xi}^{1}$, as we have done in our numerical simulations
(Sec. II). We also apply an appropriate rotation to all memories such that
$\bm{\xi}=(1,\dots,1)$. This transformation, which is orthogonal, does not
change the final result of the calculation. Since we have fixed the target
pattern, we only need to average over the other $p-1$ patterns.

###  A.2 Replica theory

In replica theory [42], the free energy per oscillator—at (effective)
equilibrium, averaged over the non-target patterns, and in the thermodynamic
limit ($N\rightarrow\infty$)—is expressed in terms of a “replicated” partition
function:

$\lim_{N\rightarrow\infty}\frac{\overline{F}}{N}=-\lim_{N\rightarrow\infty}\frac{1}{\beta
N}\overline{\ln Z}=-\lim_{N\rightarrow\infty}\lim_{n\rightarrow
0}\frac{1}{\beta Nn}\ln\overline{Z^{n}},$ (9)

where $\overline{\cdots}$ denotes the pattern average, $\beta=1/T$ ($k_{B}=1$)
is the inverse temperature, $Z=\int d\bm{\theta}\exp\left[-\beta
H_{\text{eff}}(\bm{\theta})\right]$ is the partition function, and $n$ is the
number of replicas of the system. Due to the periodicity of $H_{\text{eff}}$,
the integral over each $\theta_{i}$ ($i=1,\dots,N$) can be chosen to run over
any (contiguous) interval of length $2\pi$ (e.g., $[0,2\pi]$). Plugging in
$H_{\text{eff}}$ [Eq. (8)] of the chosen gauge, the replicated partition
function can be written as

$\displaystyle\overline{Z^{n}}$ $\displaystyle=\int
d\overleftrightarrow{\bm{\theta}}\underbrace{\exp\left\\{\frac{\beta\zeta_{1}}{2N}\sum_{\gamma}\left[\bm{\xi}\cdot\mathbf{c}(\bm{\theta}^{\gamma})\right]^{2}\right\\}}_{\equiv
E_{1}}$
$\displaystyle\quad\times\underbrace{\overline{\exp\left\\{\frac{\beta\zeta_{1}}{2N}\sum_{\gamma}\sum_{\mu>1}\left[\bm{\xi}^{\mu}\cdot\mathbf{c}(\bm{\theta}^{\gamma})\right]^{2}\right\\}}}_{\equiv
E_{2}}$
$\displaystyle\quad\times\underbrace{\exp\left\\{\frac{\beta\zeta_{2}\epsilon}{4N}\sum_{\gamma}\left[\mathbf{1}\cdot\mathbf{c}(2\bm{\theta}^{\gamma})\right]^{2}\right\\}}_{\equiv
E_{3}},$ (10)

where
$\overleftrightarrow{\bm{\theta}}=(\bm{\theta}^{1},\dots,\bm{\theta}^{n})$
contains the configurations
$\bm{\theta}^{\gamma}=(\theta_{1}^{\gamma},\dots,\theta_{N}^{\gamma})$ of
replica $\gamma=1,\dots,n$. Note that $E_{2}$ includes the pattern average
($\overline{\cdots}$). We next write each $E_{i}$ as an integral, whose
integration variables will become order parameters. The terms $E_{1}$ and
$E_{2}$ are the contributions from the target pattern and non-target patterns,
respectively. In analogy to the replica calculation for the Hopfield model [8,
11], we have

$\displaystyle E_{1}$
$\displaystyle=\left(\frac{N\beta\zeta_{1}}{2\pi}\right)^{n/2}$
$\displaystyle\quad\times\int
d\mathbf{m}\exp\beta\zeta_{1}\left\\{\sum_{\gamma}m^{\gamma}\left[\bm{\xi}\cdot\mathbf{c}(\bm{\theta}^{\gamma})\right]-\frac{N|\mathbf{m}|^{2}}{2}\right\\},$
(11) $\displaystyle E_{2}$ $\displaystyle=\left(\frac{N}{2\pi}\right)^{n^{2}}$
$\displaystyle\quad\times\int
d\mathbf{q}d\hat{\mathbf{q}}\exp\left\\{i\sum_{\gamma\kappa}\hat{q}^{\gamma\kappa}\left[Nq^{\gamma\kappa}-\mathbf{c}(\bm{\theta}^{\gamma})\cdot\mathbf{c}(\bm{\theta}^{\kappa})\right]\right\\}$
$\displaystyle\quad\times\left[\det\left(\mathbf{I}-\beta\zeta_{1}\mathbf{q}\right)\right]^{-p/2},$
(12)

where $\mathbf{I}$ is the $n\times n$ identity matrix. In tensor notation, the
integration variables are $m^{\gamma}$, $q^{\gamma\kappa}$, and
$\hat{q}^{\gamma\kappa}$ for $\gamma,\kappa=1,\dots,n$. Note that Eq. (12) has
been written assuming $N\rightarrow\infty$ [8, 11]. The term $E_{3}$ is the
contribution from the stability-enhancing term [second term of Eq. (1)]. Using
a standard Gaussian integral identity, it is straightforward to show that

$\displaystyle E_{3}$
$\displaystyle=\left(\frac{N\beta\zeta_{2}\epsilon}{4\pi}\right)^{n/2}$
$\displaystyle\times\int
d\mathbf{w}\exp\frac{\beta\zeta_{2}\epsilon}{2}\left\\{\sum_{\gamma}w^{\gamma}\left[\mathbf{1}\cdot\mathbf{c}(2\bm{\theta}^{\gamma})\right]-\frac{N|\mathbf{w}|^{2}}{2}\right\\},$
(13)

where the integration variables are (in tensor notation) $w^{\gamma}$ for
$\gamma=1,\dots,n$. We can then substitute Eqs. (11)-(13) into Eq. (10) and
write the replicated partition function in the form $\overline{Z^{n}}=\int
d\mathbf{m}d\mathbf{w}d\mathbf{q}d\hat{\mathbf{q}}\exp\left(\cdots\right)$.
Finally, we evaluate the free energy per oscillator [Eq. (9)] using the
saddle-point method,

$\lim_{N\rightarrow\infty}\frac{\overline{F}}{N}=\min_{\left(\mathbf{m},\mathbf{w},\mathbf{q},\hat{\mathbf{q}}\right)}f\left(\mathbf{m},\mathbf{w},\mathbf{q},\hat{\mathbf{q}}\right),$
(14)

where

$\displaystyle
f\left(\mathbf{m},\mathbf{w},\mathbf{q},\hat{\mathbf{q}}\right)$
$\displaystyle=\lim_{n\rightarrow 0}\frac{1}{\beta
n}\left[\frac{\beta\zeta_{1}|\mathbf{m}|^{2}}{2}+\frac{\beta\zeta_{2}\epsilon|\mathbf{w}|^{2}}{4}\right.$
$\displaystyle\quad-i\sum_{\gamma\kappa}\hat{q}^{\gamma\kappa}q^{\gamma\kappa}+\frac{\alpha}{2}\ln\det\left(\mathbf{I}-\beta\zeta_{1}\mathbf{q}\right)$
$\displaystyle\quad\left.-\frac{1}{N}\ln\int
d\overleftrightarrow{\bm{\theta}}G\left(\mathbf{m},\mathbf{w},\mathbf{q},\hat{\mathbf{q}},\overleftrightarrow{\bm{\theta}}\right)\right]$
(15)

is the function to be minimized and

$\displaystyle
G\left(\mathbf{m},\mathbf{w},\mathbf{q},\hat{\mathbf{q}},\overleftrightarrow{\bm{\theta}}\right)$
$\displaystyle=\exp\left\\{\beta\zeta_{1}\sum_{\gamma}m^{\gamma}\left[\bm{\xi}\cdot\mathbf{c}(\bm{\theta}^{\gamma})\right]\right.$
$\displaystyle\quad+\frac{\beta\zeta_{2}\epsilon}{2}\sum_{\gamma}w^{\gamma}\left[\mathbf{1}\cdot\mathbf{c}(2\bm{\theta}^{\gamma})\right]$
$\displaystyle\quad\left.-i\sum_{\gamma\kappa}\hat{q}^{\gamma\kappa}\left[\mathbf{c}(\bm{\theta}^{\gamma})\cdot\mathbf{c}(\bm{\theta}^{\kappa})\right]\right\\}.$
(16)

As mentioned above, the independent variables of $f$ correspond to order
parameters. The parameter values that minimize $f$ characterize the minima of
the free energy landscape. In particular, the values that globally minimize
$f$ correspond to the (effective) equilibrium values and yield the (effective)
equilibrium free energy through Eq. (14).

The physical meaning of the order parameters can be deduced from the saddle-
point equations. For example, the equations obtained by extremizing $f$ [Eq.
(15)] with respect to $m^{\gamma}$ and $\hat{q}^{\gamma\kappa}$ are

$\displaystyle m^{\gamma}$
$\displaystyle=\frac{1}{N}\left\langle\bm{\xi}\cdot\mathbf{c}(\bm{\theta}^{\gamma})\right\rangle,$
(17) $\displaystyle q^{\gamma\kappa}$
$\displaystyle=\begin{cases}\frac{1}{N}\sum_{i=1}^{N}\left\langle\cos^{2}\theta_{i}^{\gamma}\right\rangle,&\gamma=\kappa,\\\
\frac{1}{N}\sum_{i=1}^{N}\left\langle\cos\theta_{i}^{\gamma}\cos\theta_{i}^{\kappa}\right\rangle,&\gamma\neq\kappa,\end{cases}$
(18)

where we have defined the average $\left\langle
h\left(\overleftrightarrow{\bm{\theta}}\right)\right\rangle=\frac{\int
d\overleftrightarrow{\bm{\theta}}G\left(\overleftrightarrow{\bm{\theta}}\right)h\left(\overleftrightarrow{\bm{\theta}}\right)}{\int
d\overleftrightarrow{\bm{\theta}}G\left(\overleftrightarrow{\bm{\theta}}\right)}$
for any function $h\left(\overleftrightarrow{\bm{\theta}}\right)$. As in the
Hopfield model [8, 11], $m^{\gamma}$ represents the overlap with the target
memory (see Sec. I for definition of overlap), and $q^{\gamma\kappa}$ for
$\gamma\neq\kappa$ is the Edwards-Anderson parameter [48] except with the
variable of spin $i$ replaced by $\cos\theta_{i}$. Thus, minimization of $f$
will provide information about the overlap at free energy minima.

###  A.3 Replica-symmetric solution

To minimize $f$ explicitly, we must impose some constraints on the values of
the order parameters. Let us focus on the replica-symmetric ansatz

$\displaystyle m^{\gamma}$ $\displaystyle=m,$ (19) $\displaystyle w^{\gamma}$
$\displaystyle=w,$ (20) $\displaystyle q^{\gamma\kappa}$
$\displaystyle=Q\delta_{\gamma\kappa}+q(1-\delta_{\gamma\kappa}),$ (21)
$\displaystyle\hat{q}^{\gamma\kappa}$
$\displaystyle=\frac{i\alpha(\beta\zeta_{1})^{2}}{2}\left[R\delta_{\gamma\kappa}+r(1-\delta_{\gamma\kappa})\right].$
(22)

Then $f$ [Eq. (15)] reduces to

$\displaystyle f(m,w,Q,q,R,r)$
$\displaystyle=\frac{\zeta_{1}m^{2}}{2}+\frac{\zeta_{2}\epsilon}{2}\left(\frac{w^{2}}{2}+w\right)$
$\displaystyle\quad+\frac{\alpha}{2}\Biggl{\\{}\beta\zeta_{1}^{2}\left(RQ-
rq\right)$
$\displaystyle\quad+\frac{1}{\beta}\ln\left[1-\beta\zeta_{1}(Q-q)\right]$
$\displaystyle\quad-\frac{\zeta_{1}q}{1-\beta\zeta_{1}(Q-q)}\Biggr{\\}}$
$\displaystyle\quad-\frac{1}{\beta}\int Dz\ln\mathcal{I}_{0}(m,w,R,r,z),$ (23)

where

$\displaystyle\mathcal{I}_{v}(m,w,R,r,z)$ $\displaystyle=\int d\theta
g(m,w,R,r,\theta,z)\cos^{v}\theta,$ (24) $\displaystyle g(m,w,R,r,\theta,z)$
$\displaystyle=\exp\beta\Biggl{\\{}\zeta_{1}\left(m+z\sqrt{\alpha
r}\right)\cos\theta$ $\displaystyle\quad+\left[\zeta_{2}\epsilon
w+\frac{\alpha\beta\zeta_{1}^{2}}{2}(R-r)\right]\cos^{2}\theta\Biggr{\\}}.$
(25)

In evaluating the determinant of Eq. (15), we used the fact that $\mathbf{q}$
is a $n\times n$ matrix of the form $M_{ij}=a\delta_{ij}+b(1-\delta_{ij})$ and
thus has eigenvalues $Q-q$ and $Q+(n-1)q$ with multiplicities $n-1$ and 1,
respectively. In addition, the integral in Eq. (15) was simplified to that in
Eq. (23) using the following facts: $\bm{\xi}=(1,\dots,1)$, all oscillators
are identical, all replicas are identical, and

$\lim_{n\rightarrow 0}\frac{\ln\int Dz\,h^{n}(z)}{n}=\int Dz\,\ln h(z)$ (26)

for any function $h$ such that the limit can be evaluated using L’Hôpital’s
rule. The extrema of $f$ are solutions of the saddle-point equations

$\displaystyle m$ $\displaystyle=\int Dz\,u_{1}(m,w,R,r,z),$ (27)
$\displaystyle w$ $\displaystyle=-1+2Q,$ (28) $\displaystyle Q$
$\displaystyle=\int Dz\,u_{2}(m,w,R,r,z),$ (29) $\displaystyle q$
$\displaystyle=\int Dz\,u_{1}^{2}(m,w,R,r,z),$ (30) $\displaystyle R$
$\displaystyle=\frac{1}{\beta\zeta_{1}}\frac{1-\beta\zeta_{1}(Q-2q)}{[1-\beta\zeta_{1}(Q-q)]^{2}},$
(31) $\displaystyle r$ $\displaystyle=\frac{q}{[1-\beta\zeta_{1}(Q-q)]^{2}},$
(32)

where

$\displaystyle u_{v}(m,w,R,r,z)$
$\displaystyle=\frac{\mathcal{I}_{v}(m,w,R,r,z)}{\mathcal{I}_{0}(m,w,R,r,z)}.$
(33)

Given that $H_{\text{eff}}$ [Eq. (8)] is similar in form to the Hamiltonian of
the Hopfield model [7], we expect the solution(s) of Eqs. (27)-(32) to
eventually become unstable as $T\rightarrow 0$ for some values of $\alpha$.
Indeed, such instability in the replica-symmetric solution was found for the
Hopfield model [8, 9, 49, 11]. This issue can be addressed by using an ansatz
that breaks replica symmetry [42], which has been done for the Hopfield model
[50, 51].

###  A.4 Numerical details

To numerically solve the saddle-point equations (27)-(32), we optimize with
respect to $m$, $Q$, and $q$ only, while $w$, $Q$, and $q$ are constrained by
Eqs. (28), (31), and (32), respectively. The computation is done in Python
using the optimize.root function (all keyword parameters take their default
value, including method=’hybr’) of the SciPy package. We repeat the
calculation for initial conditions $(m,Q,q)=$ $(0.9,0.9,0.9)$, $(0,0.99,0.9)$,
$(0,0.9,0.9)$, $(0,0.8,0.6)$, $(0,0.6,0)$, $(0,0.4,0.1)$, $(0,0.49999999,0)$.
We accept solutions $(m,Q,q)$ for which the maximum relative error, with
respect to the nonzero elements, is less than $10^{-6}$. For each order
parameter $x=m,Q,q$, the relative error is calculated as
$(x_{\text{solution}}-x_{\text{actual}})/x_{\text{solution}}$, where
$x_{\text{solution}}$ is the value in the solution and $x_{\text{actual}}$ is
the value computed from the right-hand side of the saddle-point equation for
$x$. We note that, for some choices of model parameters ($\alpha$, $\epsilon$,
$T$, $\tau$), no solution could be found using the attempted initial
conditions.

###  A.5 Determination of phases

For a given choice of model parameters ($\alpha$, $\epsilon$, $T$, $\tau$), we
assign a phase according to the sets of order parameters that solve saddle-
point equations (27)-(32) and thus extremize $f$ [Eq. (23)]. Specifically, the
phase is determined by $(m,q)$ of the global minimum of $f$ and $m$ of the
local minima, as done in replica calculations of the Hopfield model [8, 9, 11]
and reviewed in Table 1. The existence of a local minimum with $m\neq 0$
implies that the target pattern can be retrieved, while the absence thereof
implies that retrieval is not possible.

Phase |  Retrieval or no retrieval of target pattern? |  $(m,q)$ of the global minimum of $f$ |  $m$ of the local minima (if any) of $f$
---|---|---|---
ferromagnetic (F) |  retrieval |  $m\neq 0$, $q\neq 0$ |  $m=$ any value
mixed (M) |  retrieval |  $m=0$, $q\neq 0$ |  $m\neq 0$ for $\geq 1$ minimum
spin glass (SG) |  no retrieval |  $m=0$, $q\neq 0$ |  $m=0$ for all minima
paramagnetic (P) |  no retrieval |  $m=0$, $q=0$ |  $m=0$ for all minima
Table 1: Phases of Hopfield-like models as described within replica theory.

###  A.6 Calculation of overlap: theory versus simulation

The highest overlap, with respect to all minima of $f$, is expected to
correspond to the final overlap in a simulation of pattern retrieval. Fig. 7
compares the highest overlap calculated in our theory (lines) with the final
overlap obtained in our simulations (circles).

We first focus on the system under passive noise (Fig. 7, blue). Quantitative
agreement is obtained at low pattern loading, while the theoretical overlap
reasonably approximates the decay of the simulated overlap as the pattern
loading increases. Note, however, that the theory overestimates the overlap
during the decay and incorrectly predicts that the overlap approaches zero. We
attribute these discrepancies to the fact that the simulations are done for a
finite number $N$ of oscillators while the theory assumes
$N\rightarrow\infty$.

Figure 7: Overlap with target pattern versus pattern loading $\alpha$, as
calculated by simulation (circles) and theory (lines). Results are shown for
passive (purple) and active (blue, $\tau=0.2$) noise; noise strength (a-c)
$T=0.2$, (d-f) $T=0.4$, and (g-i) $T=0.5$; and stability-enhancing coupling
(a, d, g) $\epsilon=0$, (b, e, h) $\epsilon=0.5$, and (c, f, i) $\epsilon=1$.

Having validated the theory for the case of passive noise, we now turn to the
system under active noise with $\tau=0.2$ (Fig. 7, red). Like the case of
passive noise, we see good agreement between simulation and theory at
relatively low pattern loading and noise strength. In contrast, at
sufficiently high noise strength, the theory predicts significant overlap
under active noise while the simulations show little overlap (Figs. 7g-7h). We
believe this discrepancy is related to the following: the effective
Hamiltonian [Eq. (8)] used in the replica calculation does not include the
$\frac{\tau}{2}\left|\nabla_{\bm{\theta}}H\right|^{2}$ term of the full
effective Hamiltonian [Eq. (6)] derived by effective equilibrium theory
(Section III.1). At higher noise strengths, the system is kicked away from the
patterns, and so $-\nabla_{\bm{\theta}}H$, the force pushing the system
towards the patterns, becomes significant.

## References

  * Keim  _et al._ [2019] N. C. Keim, J. D. Paulsen, Z. Zeravcic, S. Sastry, and S. R. Nagel, Reviews of Modern Physics 91, 035002 (2019).
  * Sokolov  _et al._ [2021] A. S. Sokolov, H. Abbas, Y. Abbas, and C. Choi, Journal of Semiconductors 42, 013101 (2021).
  * Amit [1989] D. J. Amit,  _Modeling brain function: The world of attractor neural networks_ (Cambridge university press, 1989).
  * Nakano [1972] K. Nakano, IEEE Transactions on Systems, Man, and Cybernetics SMC-2, 380 (1972).
  * Kohonen [1972] T. Kohonen, IEEE Transactions on Computers C-21, 353 (1972).
  * Amari [1972] S.-I. Amari, IEEE Transactions on Computers C-21, 1197 (1972).
  * Hopfield [1982] J. J. Hopfield, Proceedings of the National Academy of Sciences 79, 2554 (1982).
  * Amit  _et al._ [1987] D. J. Amit, H. Gutfreund, and H. Sompolinsky, Annals of Physics 173, 30 (1987).
  * Amit  _et al._ [1985] D. J. Amit, H. Gutfreund, and H. Sompolinsky, Physical Review Letters 55, 1530 (1985).
  * Feigelman and Ioffe [1986] M. V. Feigelman and L. B. Ioffe, Europhysics Letters 1, 197 (1986).
  * Coolen [2001] A. C. C. Coolen, in  _Handbook of Biological Physics_, Neuro-Informatics and Neural Modelling, Vol. 4, edited by F. Moss and S. Gielen (North-Holland, 2001) pp. 553–618.
  * Gnesotto  _et al._ [2018] F. S. Gnesotto, F. Mura, J. Gladrow, and C. P. Broedersz, Reports on Progress in Physics 81, 066601 (2018).
  * Fang  _et al._ [2019] X. Fang, K. Kruse, T. Lu, and J. Wang, Reviews of Modern Physics 91, 045004 (2019).
  * Vicsek  _et al._ [1995] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Physical Review Letters 75, 1226 (1995).
  * Toner  _et al._ [2005] J. Toner, Y. Tu, and S. Ramaswamy, Annals of Physics Special Issue, 318, 170 (2005).
  * Jülicher  _et al._ [2007] F. Jülicher, K. Kruse, J. Prost, and J. F. Joanny, Physics Reports 449, 3 (2007).
  * Banerjee  _et al._ [2020] S. Banerjee, M. L. Gardel, and U. S. Schwarz, Annual Review of Condensed Matter Physics 11, 421 (2020).
  * Ramaswamy [2010] S. Ramaswamy, Annual Review of Condensed Matter Physics 1, 323 (2010).
  * Marchetti  _et al._ [2013] M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Reviews of Modern Physics 85, 1143 (2013).
  * Bechinger  _et al._ [2016] C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe, and G. Volpe, Reviews of Modern Physics 88, 045006 (2016).
  * Sepúlveda  _et al._ [2013] N. Sepúlveda, L. Petitjean, O. Cochet, E. Grasland-Mongrain, P. Silberzan, and V. Hakim, PLOS Computational Biology 9, e1002944 (2013).
  * Maggi  _et al._ [2014] C. Maggi, M. Paoluzzi, N. Pellicciotta, A. Lepore, L. Angelani, and R. Di Leonardo, Physical Review Letters 113, 238303 (2014).
  * Fodor  _et al._ [2016] E. Fodor, C. Nardini, M. E. Cates, J. Tailleur, P. Visco, and F. van Wijland, Physical Review Letters 117, 038103 (2016).
  * Martin  _et al._ [2021] D. Martin, J. O’Byrne, M. E. Cates, E. Fodor, C. Nardini, J. Tailleur, and F. van Wijland, Physical Review E 103, 032607 (2021).
  * Chacron  _et al._ [2000] M. J. Chacron, A. Longtin, M. St-Hilaire, and L. Maler, Physical Review Letters 85, 1576 (2000).
  * Fourcaud and Brunel [2002] N. Fourcaud and N. Brunel, Neural Computation 14, 2057 (2002).
  * Selverston and Moulins [1985] A. I. Selverston and M. Moulins, Annual Review of Physiology 47, 29 (1985).
  * Ward [2003] L. M. Ward, Trends in Cognitive Sciences 7, 553 (2003).
  * Wang [2010] X.-J. Wang, Physiological reviews 90, 1195 (2010).
  * Schuster and Wagner [1990] H. G. Schuster and P. Wagner, Biological Cybernetics 64, 77 (1990).
  * Arenas and Vicente [1994] A. Arenas and C. J. P. Vicente, Europhysics Letters 26, 79 (1994).
  * Fukai and Shiino [1994] T. Fukai and M. Shiino, Europhysics Letters 26, 647 (1994).
  * Fukai [1994] T. Fukai, Biological Cybernetics 70, 427 (1994).
  * Park and Choi [1995] K. Park and M. Y. Choi, Physical Review E 52, 2907 (1995).
  * Aonishi [1998] T. Aonishi, Physical Review E 58, 4865 (1998).
  * Aoyagi and Kitano [1998] T. Aoyagi and K. Kitano, Neural computation 10, 1527 (1998).
  * Kuramoto [1975] Y. Kuramoto, in  _International Symposium on Mathematical Problems in Theoretical Physics: January 23-29, 1975, Kyoto University, Kyoto/Japan_ (Springer, 1975) pp. 420–422.
  * Nishikawa  _et al._ [2004] T. Nishikawa, Y.-C. Lai, and F. C. Hoppensteadt, Physical Review Letters 92, 108101 (2004).
  * Maitra and Voituriez [2020] A. Maitra and R. Voituriez, Physical Review Letters 124, 048003 (2020).
  * Jung and Hänggi [1987] P. Jung and P. Hänggi, Physical Review A 35, 4464 (1987).
  * Maggi  _et al._ [2015] C. Maggi, U. M. B. Marconi, N. Gnan, and R. Di Leonardo, Scientific Reports 5, 10742 (2015).
  * Mézard  _et al._ [1987] M. Mézard, G. Parisi, and M. A. Virasoro,  _Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications_ , Vol. 9 (World Scientific Publishing Company, 1987).
  * Hebb [2005] D. O. Hebb,  _The organization of behavior: A neuropsychological theory_ (Psychology press, 2005).
  * Behera  _et al._ [2022] A. K. Behera, M. Rao, S. Sastry, and S. Vaikuntanathan, arXiv:2203.03024 (2022).
  * Bollé  _et al._ [2003] D. Bollé, T. M. Nieuwenhuizen, I. P. Castillo, and T. Verbeiren, Journal of Physics A: Mathematical and General 36, 10269 (2003).
  * Krotov and Hopfield [2016] D. Krotov and J. J. Hopfield, in  _Advances in Neural Information Processing Systems_, Vol. 29 (Curran Associates, Inc., 2016).
  * Ramsauer  _et al._ [2021] H. Ramsauer, B. Schäfl, J. Lehner, P. Seidl, M. Widrich, T. Adler, L. Gruber, M. Holzleitner, M. Pavlović, G. K. Sandve, V. Greiff, D. Kreil, M. Kopp, G. Klambauer, J. Brandstetter, and S. Hochreiter, Hopfield Networks is All You Need (2021).
  * Edwards and Anderson [1975] S. F. Edwards and P. W. Anderson, Journal of Physics F: Metal Physics 5, 965 (1975).
  * Naef and Canning [1992] J.-P. Naef and A. Canning, Journal de Physique I 2, 247 (1992).
  * Crisanti  _et al._ [1986] A. Crisanti, D. J. Amit, and H. Gutfreund, Europhysics Letters 2, 337 (1986).
  * Steffan and Kühn [1994] H. Steffan and R. Kühn, Zeitschrift für Physik B Condensed Matter 95, 249 (1994).