Emergent Self-Attention from Astrocyte-Gated Associative Memory Dynamics Arnau Vivet1 and Alex Arenas1, 2 1
arXiv:2604.25481v1 [physics.data-an] 28 Apr 2026
Departament d’Enginyeria Informàtica i Matemàtiques, Universitat Rovira i Virgili, 43007 Tarragona, Spain 2 Complexity Science Hub Vienna, Metternichgasse 8, 1030 Vienna, Austria (Dated: April 29, 2026) We introduce a Hopfield-type associative memory in which effective connectivity is multiplicatively modulated by astrocytic gains evolving under an entropy-regularized replicator equation. The coupled neuron–astrocyte dynamics admit a Lyapunov function, ensuring global convergence. At fixed points, astrocytic gains implement a softmax-normalized allocation over pattern similarity scores, yielding a mechanistic realization of self-attention as emergent routing on the gain simplex. In regimes of high memory load and interference, the model significantly improves retrieval accuracy relative to classical Hopfield dynamics and recent neuron–astrocyte baselines. These results establish a dynamical systems framework linking glial modulation, competitive resource allocation, and attention-like computation.
I.
INTRODUCTION
Hopfield networks revolutionized our understanding of how collective neuronal dynamics implement associative memory with biological plausibility [1]. In their original formulation, neurons evolve under Hebbian connectivity that stores a finite set of patterns as energy minima. Although this framework provided a foundational bridge between statistical physics and neural computation, its practical applicability was limited by low storage capacity and the proliferation of spurious attractors under high pattern load or correlation [2]. Interest resurged with generalized energy-based formulations, dense associative memory and modern Hopfield networks, that replace the quadratic Hebbian energy with higher-order or exponential interaction terms, dramatically expanding capacity and retrieval stability [3, 4]. A key insight from this line of work is that the fixed-point readout of modern Hopfield networks is mathematically equivalent to the scaled dot-product attention mechanism at the core of Transformer architectures [5, 6]. This correspondence has reframed associative memory dynamics as a candidate substrate for attention and contextual reasoning in both artificial and biological systems. Neuroscience has concurrently revealed that the brain’s computational resources extend well beyond neurons. Astrocytes, glial cells long considered purely supportive, actively regulate synaptic transmission, integrate signals across local microcircuits, and contribute to memory-related plasticity [7–12]. Through Ca2+ signaling, neurotransmitter uptake, and connexin-mediated gap-junction coupling, astrocytes modulate neuronal excitability and impose structured heterogeneity in effective synaptic strengths [9, 11, 12]. Astrocytes have also been implicated in controlling circuit-level dynamical regimes, including the modulation of gamma-band synchronization associated with memory processing [13– 15]. Critically, astrocytes are now recognized as active participants in memory encoding and consolidation: recent experimental work demonstrates that astrocyte
ensembles selectively encode and stabilize salient experiences, with individual astrocytes reactivated during memory retrieval undergoing noradrenaline-dependent molecular tagging that promotes consolidation of specific traces [16, 17]. These findings extend the classical engram framework [18] by assigning astrocytes a mechanistic role in the selection and stabilization of memory representations, motivating their inclusion in computational models of associative memory. Building on this foundation, recent theoretical work has proposed that astrocytic modulation could implement attention-like computation by providing a flexible, content-dependent reweighting of neuronal ensembles during retrieval [19, 20]. In this view, self-attention corresponds to a gain reallocation across stored patterns driven by their relevance to the current query state. However, a precise dynamical account of how softmaxnormalized gains arise as an emergent, self-organized property of neuron–astrocyte interaction, rather than being imposed by architectural design, has remained absent. Here we close this gap by introducing a Hopfield-type associative memory in which effective connectivity is multiplicatively modulated by astrocytic gains that evolve on the probability simplex under an entropy-regularized replicator equation. The coupled neuron–astrocyte system admits a Lyapunov function, ensuring convergence to stationary points (equilibria) at which gains implement a Gibbs–Boltzmann (softmax) allocation over pattern similarity scores. The stationary readout is thus a convex combination of stored patterns with self-attention weights—emerging purely from the competitive resource dynamics of glial modulation, without any explicit attention mechanism being prescribed. Under high memory load and pattern interference, this mechanism substantially improves retrieval accuracy relative to classical Hopfield dynamics and recent neuron–astrocyte baselines, establishing a dynamical systems framework that links glial gain control, competitive resource allocation, and the emergence of attention-like computation in neural circuits.
2 II.
with effective synaptic matrix
ASTROCYTE-GATED ASSOCIATIVE MEMORY
K
To ground the proposed astrocyte-neuron model, we begin by recalling the continuous-time classical Hopfield network, the minimal framework for associative memory. It is a class of recurrent neural network in which we consider a population of rate-based neurons evolving under a fully connected synaptic matrix. One of the key properties of this framework is that it can be derived from a Lyapunov function that guarantees convergence, enabling a precise characterization of memory retrieval as well as interference and capacity limits [1, 2]. For these reasons, the Hopfield model provides a natural substrate on which biologically motivated extensions can be constructed. More precisely, we consider N rate units x(t) ∈ RN with element-wise activation function ϕ(x) = tanh(σx) ∈ (−1, 1)N and synaptic weights W ∈ RN ×N : X τx ẋi = −xi + W ij ϕ(xj ). (1) j
K
1 X i j ξ ξ . N µ=1 µ µ
(2)
Alternatively in matrix form we have WH = N1 ΞΞ⊺ , where Ξ ≡ (ξ1 , ..., ξK ) ∈ {−1, 1}N ×K is a binary matrix where patterns are stored as columns. Throughout the paper, we use superscripts to label neurons and subscripts to label memory patterns. A.
Astrocyte-gated retrieval as pattern-wise gain control
We extend the classical Hopfield rate dynamics by introducing K astrocytic gain variables p = (p1 , . . . , pK ) that modulate the contribution of each stored pattern to the synaptic matrix. The gains are constrained to the probability simplex, pµ ≥ 0
(∀µ),
K X
pµ = 1,
(3)
µ=1
so modulation is non-negative and globally conserved. The resulting neuronal dynamics read τx ẋi = −xi +
N X j=1
KX K pµ ξµ ξµ⊺ = Ξ diag(p) Ξ⊺ . N µ=1 N
(5)
Thus, pµ acts as a multiplicative gain on the µ-th Hebbian outer product. The classical Hopfield coupling is recovered for uniform gains pµ = 1/K, for which W (p) = WH . We model astrocytic modulation as an adaptive allocation process governed by a softmax-regularized replicator flow (hereafter “s-replicator”): K X τp ṗµ = pµ Fµ − p ν Fν ,
(6)
ν=1
which preserves the simplex constraints. The fitness is defined as Fµ ≡ fµ − T log pµ ,
T > 0,
(7)
where
where τx is the neuronal rate time constant. It sets the intrinsic timescale over which the firing-rate state x(t) relaxes toward the drive W ϕ(x) against the leak term. The synaptic matrix W is constructed using Hebbian storage of K binary patterns ξµ ∈ {−1, 1}N located at the vertices of the hypercube WH =
W (p) ≡
W ij (p) ϕ(xj ),
(4)
fµ ≡
N 2 1 X j ξµ ϕ(xj ) 2N j=1
(8)
is the squared overlap with pattern µ, and the logarithmic term acts as an entropic (temperature-like) regularizer that discourages collapse of the gains onto a single pattern. Throughout, we use the same neuronal nonlin earity as in the classical model, ϕ(xj ) = tanh σxj with σ > 0. We study retrieval dynamics by initializing the neuronal state x(0) as a corrupted version of one of the stored patterns. Unless stated otherwise, the astrocytic gains are initialized uniformly, pµ (0) = 1/K for all µ. During the dynamics, Eq. (6) reallocates gain toward patterns with larger instantaneous overlap fµ (x), thereby reshaping the effective coupling W (p) in Eq. (5). This positive feedback selectively amplifies the synaptic contribution of the most compatible patterns and suppresses competing ones, reducing interference and effectively increasing the basin of attraction of the target memory. The biophysical motivation for Eqs. (4,6) is the following. We view each stored pattern as an engram-like neuronal ensemble, and interpret the gain pµ as an effective proxy for astrocyte-mediated modulation of the synaptic pathways supporting that ensemble. Specifically, pµ summarizes (at a coarse-grained level) how local astrocytic activity can bias presynaptic release probability and thus rescale the effective strength of the synapses recruited by pattern µ. This abstraction is inspired by evidence that synaptic activity can evoke highly localized astrocytic Ca2+ signals in perisynaptic microdomains adjacent to individual synapses [21]. In several experimental settings, astrocytic Ca2+ elevations have been reported to modulate presynaptic release probability and thereby reshape effective synaptic transmission, consistent with
3 an activity-dependent astrocytic gating of synaptic pathways [22]. At the same time, the extent to which Ca2+ -dependent astrocytic signaling modulates excitatory transmission and plasticity under baseline conditions remains debated and depends on preparation and stimulation regime [23]; accordingly, we treat pµ as a phenomenological variable that aggregates multiple microscopic mechanisms rather than as a direct readout of any single molecular process. Finally, we constrain the astrocytic gains to the simPK plex, p ∈ ∆K−1 , i.e., pµ ≥ 0 and µ=1 pµ = 1. In the model, this implements a finite modulatory capacity: increasing gain for one pattern necessarily reduces the gain available to others, thereby enforcing competitive allocation across pathways. This is a deliberate modeling device that operationalizes resource limitation and homeostatic competition at the level of pattern gains, rather than a claim that the brain explicitly maintains a globally conserved scalar “budget” [24]. The neuron–astrocyte coupling is driven by a patternmatch functional fµ (x), which quantifies how strongly the current neuronal state expresses pattern µ (via its overlap with ξµ ). We map this match to a fitness Fµ that controls the gain dynamics: patterns with larger Fµ are preferentially upweighted by Eq. (6). The additional logarithmic term in Fµ acts as an entropic regularizer, penalizing highly concentrated allocations and preventing winner-take-all behavior; consequently, the stationary distribution of gains is softmax-like, with sharpness set by the parameter T .
1.
The system is a global gradient flow
A key property of the classical Hopfield network is the existence of an energy (Lyapunov) function that decreases along the trajectories thus guaranteeing convergence of the dynamics. As shown in the Appendix 1, our system can also be written as a gradient flow whose energy is given by X pµ L(x, p) =KT pµ log − KT log Z π µ µ (9) N X i i i + x ϕ(x ) − L(x ), i
P fµ /T and L(xi ) = where πµ = Z1 efµ /T with Z = µe 1 i The gradient flow condition implies σ log cosh σ x . that this quantity decreases along the trajectories d L(x, p; t) = −∥∇L(x, p; t)∥2 ≤ 0 dt . The last two terms of Eq.9 are similar to the modern Hopfield networks [6, 25], while the first term is the potential for the s-replicator equation also encoding for neuron-astrocyte interaction. Since L is non-increasing and constant only at stationary points, the long-time behavior is fully determined by the coupled equilibria (x∗ , p∗ ), which we analyze next.
2.
Fixed point analysis
The equilibrium (x∗ , p∗ ), satisfies the coupled fixed point condition x∗ = W (p∗ ) ϕ(x∗ ) B.
Analytic results
We now summarize three analytic results that organize the theoretical picture and will be used throughout the remainder of the paper. First, we show that the coupled neuron–astrocyte dynamics admit a global Lyapunov function, so trajectories are dissipative and converge to the set of stationary points rather than exhibiting sustained oscillations or chaos. This structural result justifies focusing on equilibria. Second, we characterize these equilibria explicitly via coupled fixed-point equations for (x∗ , p∗ ), which reveal how astrocytic modulation implements a temperature-controlled, softmax-like allocation over patterns at stationarity. Finally, we connect the framework back to the classical Hopfield model by identifying limiting regimes in which the gain distribution becomes uniform and the effective synaptic matrix reduces to WH . Together, these three results provide (i) a convergence guarantee, (ii) an interpretable description of the attractors, and (iii) a consistency link to the standard associative-memory baseline.
(10)
∗
exp(fµ /T ) fµ (x ) p∗µ = P = softmaxµ T exp(f /T ) ν ν
.
(11)
This can be interpreted as a mixture P of (rank 1) experts [26] commonly expressed as x∗ = µ gµ (x∗ )Eµ (ϕ(x∗ )), where the self-attention like routing gµ (x) ≡ pµ emerges naturally from the s-replicator, and the experts are rank1 matrices Eµ (ϕ(x)) ≡ ξµ ξµ⊺ ϕ(x). In this interpretation, astrocytic modulation provides the routing mechanism that selects the relevant memory patterns in a contextdependent and dynamical manner. These equilibrium relations make clear that the model reduces to the classical Hopfield network whenever the gain distribution is forced to remain (or becomes) uniform, which we discuss next through two limiting regimes.
3.
Recovering the classical Hopfield model
Our model recovers the classical associative memory in two dynamical regimes. For the first case we take
4 the gain timescale to be infinite τp → ∞ such that the astrocytic influence remains constant for all time ṗµ = 0 (see Appendix 4). Since the modulation is initialized uniformly pµ = 1/K, we recover the classical Hopfield mechanism τx ẋ = −x + W (1/K)ϕ(x).
(12)
This should be interpreted as the astrocytic processes being frozen, in which case there may not be a modulatory influence on the neuron dynamics. The second limit in which we recover the classical model is when we take T → ∞, in which case the entropic term dominates and the only stable gain distribution is again uniform pµ = 1/K (Appendix 4). C.
Simulation results
Our focus is on how two control knobs shape retrieval: (i) the relative timescales of neuronal relaxation and astrocytic routing, quantified by τx and τp , and (ii) the selectivity parameter T , which sets the strength of entropic regularization in the gain dynamics. Because our dynamics are Lyapunov (Sec. II B), convergence is guaranteed in principle; numerically, however, we observe a competition between intrinsic convergence time and the finite simulation horizon tf . We therefore report both endpoint observables and empirical convergence times to disentangle genuine failure from slow convergence. We emphasize that these simulations are intended as proof-ofprinciple demonstrations of dynamical regimes controlled by (τx , τp , T ). 1.
Dynamical analysis
In this first analysis, the goal is to understand the long term behavior of the system. To that end, we define two observables, one representative of the neuron component and the other for the astrocytic gains. Our focus is on how two control knobs shape retrieval: (i) the relative timescales of neuronal relaxation and astrocytic routing, quantified by τx and τp , and (ii) the selectivity parameter T , which sets the strength of entropic regularization in the gain dynamics. Because our dynamics are Lyapunov (Sec. II B), convergence is guaranteed in principle; numerically, however, we observe a competition between intrinsic convergence time and the finite simulation horizon tf . We therefore report both end-point observables and empirical convergence times to disentangle genuine failure from slow convergence. We simulate a network of N = 30 neurons storing K = 100 binary (overlapping) patterns {ξµ } via Eq. 5. To probe retrieval, we corrupt a reference memory ξ0 → ξ0η by flipping n randomly chosen entries, defining the noise level η := n/N = 0.2 (six flipped bits). We initialize the neuronal state with the corrupted memory, x(ti ) = ξ0η (see Appendix 5).
We quantify neuronal retrieval at tf using a soft Hamming error, N
Error = ϵ(tf ) :=
1X i ξ − ϕ(xi (tf )) , 2 i=1 0
(13)
which reduces to the standard Hamming distance when ϕ(xi ) ∈ {−1, 1} and satisfies ϵsoft ∈ [0, N ]. To quantify how concentrated the gain vector p(tf ) is over memories, we compute its Shannon entropy H[p] and the associated perplexity ! K X P := exp(H[p]) = exp − pµ log pµ . (14) µ=1
Perplexity satisfies P ∈ [1, K], with P = 1 indicating near winner-take-all routing and P ≈ K indicating nearuniform allocation. Unless stated otherwise, we set τx = 1, τp = 1, and T = 0.01, and vary one parameter at a time. Varying astrocyte timescale τp (keeping τx = 1) Mechanistically, τp controls how rapidly the routing weights pµ (t) track the instantaneous overlaps fµ (x(t)); small τp yields fast reweighting of W (p) toward the correct pattern, whereas large τp delays this reshaping, i.e. τp parameter controls the gain p(t) response time Fig.2. If we let astrocytic modulation be very fast τp → 0, the gains quickly concentrate on the most compatible memories, so the retrieval is fast and we get both low error and low perplexity. As τp increases, we see the gain adaptation become slower and the system needs more time to route the right memory. With a fixed simulation time, it looks like performance gets worse beyond τp ∼ 10, but this behavior is explained simply because the simulation is stopped before convergence. In the extreme τp → ∞ case, since the modulation essentially becomes frozen, the classical Hopfield regime is recovered, which performs poorly in this high memory storage regime (see also Appendix 4). Varying astrocyte timescale τx (keeping τp = 1) τx controls how fast neurons x(t) relax Fig.1. Setting τx → 0 the neuron relaxation is essentially instantaneous, this makes the neuron configuration ”commit” to an attractor too early, before the modulatory signal has time to reshape the landscape. Then the gain modifies the effective connectivity around the wrong attractor leading to high confidence (low perplexity) on the wrong memory. On the other hand, as we increase τx , we see a decrease in the error as well as an increase in perplexity. High perplexity here means multiple patterns share similar overlap, so routing stays diffuse. This seemingly paradoxical region happens because as the neuron dynamics become slower, they do not have enough time to break the symmetry between patterns with similar fitness (see Appendix 3). In fact, if we increase the time limit, symmetry breaks and we get high pattern selectivity again.
5
FIG. 1. End-of-run retrieval error ϵ(tf ) and gain perplexity P (tf ) as functions of the astrocytic timescale τp (with τx = 1), together with the corresponding convergence times. Solid lines show medians across trials; shaded bands denote percentile ranges (5, 95), (10, 90), (20, 80), and (25, 75).
If we keep increasing to τx → ∞, the convergence time increases super-linearly and we see the error stay exactly at ϵ(tf ) = 6 (in our simulations this happens at τx ∼ 102 , simply because the simulation time is bounded) (see Appendix 4).
FIG. 3. Final-time retrieval error ϵ(tf ) and gain perplexity P (tf ) as functions of temperature T (with τx = τp = 1), together with the corresponding convergence times. Solid lines show medians across trials; shaded bands denote percentile ranges (5, 95), (10, 90), (20, 80), and (25, 75).
perplexity are low since the regularization effect is weak and becomes highly selective. When T → ∞ the error increases because the entropy dominates, meaning that p stays close to uniform and once again we recover the Hopfield regime (see Appendix 4). 2.
FIG. 2. Final-time retrieval error ϵ(tf ) and gain perplexity P (tf ) as functions of the neuronal timescale τx (with τp = 1), together with the corresponding convergence times. Solid lines show medians across trials; shaded bands denote percentile ranges (5, 95), (10, 90), (20, 80), and (25, 75).
Varying temperature T (keeping τp = τx = 1): Thus T tunes a selectivity–robustness tradeoff: low T yields sharp routing (small P ) and effective interference suppression, whereas high T enforces near-uniform gains and recovers Hopfield-like behavior, i.e. temperature controls astrocyte selectivity Fig.3. When T → 0, both error and
Comparative retrieval performance
We benchmark retrieval performance against (i) the classical Hopfield network and (ii) the neuron–astrocyte associative-memory model of Kozachkov et al. [20]; see Fig. 4. To enable a dense scan over memory load and corruption, we set N = 20 and vary the number of stored random binary patterns from K = 2 to K = 200. For each (K, n) condition, we generate an independent pattern matrix Ξ, corrupt the target pattern ξ0 by flipping n randomly chosen entries (equivalently, corruption frac(n) tion η = n/N ), and initialize the network as x(ti ) = ξ0 (simulation protocol and model-specific parameters are reported in Appendix 5). All models are evaluated under the same initialization and the same finite simulation horizon tf . We quantify retrieval at time tf by binarizing the final state and computing the Hamming error with respect to the target: N 1X i ϵ(tf ) = ξ0 − sign xi (tf ) , 2 i=1
(15)
so that ϵ ∈ [0, N ] and ϵ = 0 indicates perfect retrieval. For statistical stability, we repeat each condition 50 times with independently sampled Ξ and report the mean error. Across the tested (K, n) grid, our model yields lower mean retrieval error than both baselines, with the largest
6
FIG. 4. Retrieval benchmark across models. Each heat map reports the mean Hamming retrieval error ϵ(tf ) (lighter indicates lower error) as a function of memory load K (x-axis) and corruption level n (y-axis; number of flipped bits in the query pattern). Each entry is averaged over 50 independent random pattern realizations.
gains appearing at high memory loads where interference is strongest (Fig. 4).
III.
DISCUSSION
We introduced an astrocyte-gated extension of a Hopfield-type associative memory in which astrocytic processes allocate a finite modulatory capacity across stored patterns. By dynamically reweighting patternspecific contributions to the effective connectivity, this additional degree of freedom reduces interference during retrieval and enlarges the regime in which accurate recall is achieved at high memory load. We interpret each stored pattern as a memory-specific (engram-like) neuronal ensemble and the gain variable pµ as a coarse-grained proxy for the effective strength of astrocyte-mediated modulation of the synaptic pathway supporting pattern µ. This abstraction is motivated by the tripartite-synapse framework, in which astrocytes integrate local activity and can modulate synaptic efficacy through Ca2+ -dependent signaling and related pathways [7–12]. Importantly, we do not interpret pµ as a direct readout of any single molecular mechanism; rather, it summarizes pathway-level modulation at the scale relevant for associative retrieval. A central modeling assumption is the simplex conP straint µ pµ = 1, which enforces competitive allocation: increasing gain for one pathway necessarily reduces gain available to others. This operationalizes finite modulatory capacity (e.g., limited signaling/metabolic resources distributed across microdomains) and homeostatic competition at the level of pattern gains, without implying that the brain literally implements a globally conserved scalar “budget.” Mechanistically, the gain dynamics implement state-
dependent competition. When the current neuronal state is more compatible with pattern µ, the corresponding match fµ (x) increases the fitness Fµ and thereby upweights pµ , which amplifies the µ-th contribution to W (p) while suppressing competitors. The logarithmic term in Fµ acts as an entropic regularizer: it penalizes highly concentrated allocations and discourages winnertake-all routing. In this view, the temperature parameter T controls selectivity: low T yields sharp, concentrated routing, whereas high T enforces broader, nearuniform modulation and recovers Hopfield-like behavior. The joint neuron–astrocyte system remains Lyapunovconsistent (it admits a global Lyapunov function), ensuring convergence of the coupled dynamics. Modern Hopfield networks and related dense-memory models achieve softmax-like retrieval through neuronal energy descent shaped by higher-order or exponential storage functions [3, 4, 6, 19, 20]. Our model produces a similar normalized reweighting, but via a different mechanism: normalization arises from competitive dynamics on the gain simplex (replicator-type routing) rather than being hard-wired into the neuronal energy landscape. In regimes where gains adapt rapidly, the stationary allocation approaches fµ p∗µ ∝ exp , T so retrieval can be viewed as operating with a dynamically reweighted subset of memories. Compared to the classical Hopfield model, where high load is associated with numerous spurious mixture states, our routing variable tends to concentrate weight on a small set of candidates under ambiguity, consistent with the low-perplexity regimes observed in simulations. In addition, the degree to which attention-like structure is explicit in our framework depends on the choice
7 of pattern score. While the main model uses a squaredoverlap score fµ (x) ∝ ⟨ξµ , ϕ(x)⟩2 (hence invariant under ⟨ξµ , ϕ(x)⟩ 7→ −⟨ξµ , ϕ(x)⟩), one may alternatively define a linear-overlap score fµ (x) = ⟨ξµ , ϕ(x)⟩. In this variant, the neuronal drive becomes a linear combination of stored patterns weighted by the gain vector, i.e. the update takes the form ẋ = −x + Ξp (up to timescale factors), and in the fast-gain limit the stationary allocation becomes p∗ ∝ exp(Ξ⊺ ϕ(x)/T ), yielding the standard attention readout x ← Ξ softmax(Ξ⊺ ϕ(x)/T ) [5, 6]. We do not analyze this variant further here; we include it to emphasize that competitive routing on the gain simplex can mechanistically reproduce attention-like computation under closely related definitions of pattern compatibility. Numerically, this modulatory degree of freedom improves retrieval accuracy across increasing memory load K and increasing corruption level n (equivalently η = n/N ), outperforming both the classical Hopfield model [1] and the neuron–astrocyte associative-memory baseline [20] in the tested regimes (Fig. 4). Intuitively, adaptive gain allocation reshapes the effective connectivity (and hence the attractor structure) during recall by amplifying task-relevant patterns and suppressing competitors, in line with recent evidence that time-dependent modulation can improve robustness of Hopfield-type retrieval [27]. The model also suggests qualitative dependencies that could guide future theoretical and computational work. In particular, changes in astrocytic kinetics or background modulatory tone would be expected to alter the sharpness of pattern selection (captured here by T ), while retrieval speed and stability should depend on the ratio τp /τx setting the competition between routing and neuronal relaxation. More generally, the framework highlights a route by which transient changes in astrocytemediated signaling could bias recall without invoking synaptic plasticity. A limitation of the present formulation is that pµ aggregates multiple biophysical processes and timescales into a single effective variable; identifying which cellular mechanisms and dynamical regimes can realize competitive allocation of this type remains an important direction for future work, and connects naturally to multi-timescale theories of memory stabilization [28]. An important direction is to convert the proposed competitive routing mechanism into a scalable ML architecture, where p acts as a parameterized, contextdependent gating distribution over memory components and T controls gating entropy. This suggests lightweight mixture-of-experts memories with explicit competition and separable routing and state-update timescales (set by τp /τx ). Developing stable, differentiable implementations and benchmarking them on noisy retrieval and continual-learning tasks are natural next steps.
ACKNOWLEDGMENTS
We thank Prof. Luiz Pessoa and Prof. Sergio Gómez for discussions on astrocyte kinetics and associative memory. This work has been supported by spanish Ministerio de Ministerio de Ciencia, Innovación y Universidades PID2024-158120NB-C21. AA also acknowledges the ICREA Academia program of Generalitat de Catalunya. APPENDIX
For the reader’s convenience, we collect here the few definitions that are repeatedly used in the derivations below, so the proofs can be followed without having to refer back to the main text. In particular, the ef⊺ fective synaptic matrix is W (p) = K = N Ξ diag(p) Ξ PK K ⊺ K−1 p ξ ξ . For the gain dynamics, p ∈ ∆ imµ=1 µ µ µ N plies 1⊤ ṗ = 0 (tangent-space constraint), and the projector acts on Tp ∆K−1 . In the gradient-flow formulation we use the Shahshahani (Fisher information) metric on ∆K−1 , G(p) = diag(1/p), and the diagonal metric on neuronal pre-activations, Hx = diag(ϕ′ (x)) with ϕ′ (x) = σ sech2 (σx) for ϕ(x) = tanh(σx). 1.
Gradient-flow structure and Lyapunov function
Here we show that the coupled dynamics in Eqs. (4,6) admit a Lyapunov function and can be written as a (projected) gradient flow. We divide the discussion into the astrocyte domain and the neuron domain, where we derive their respective potentials. Note that the coupled dynamics occur in x, p ∈ RN × ∆K−1 and the metric (needed for the gradient flow result) is given by a blockdiagonal (direct-sum) metric Gx,p = Hx ⊕ Gp . This metric is composed of the astrocyte and neuron domain metrics. The former is given by the Fisher information metric (equivalently the Shahshahani metric), whose components are gµµ′ (p) =
δµµ′ , pµ
so in matrix form G(p) = diag(1/p). The latter is given by Hx = diag(ϕ′ (x)), where in our case ϕ′ (x) = σsech2 (σx). a.
Astrocyte domain
For the astrocyte domain, we first show how the simplex constraints are imposed via the orthogonal projection and then we define the s-replicator potential. With
8 this, we obtain the ODE for the p variable and we discuss its fixed points. Following [29], given that we have a differential 0-form (or scalar potential) Φ : RK → R and a metric G(p), we will define an ODE by computing the gradient of the scalar potential, and then projecting it onto the tangent PK space of the simplex ṗ ∈ Tp ∆K−1 (i.e. 1⊺ ṗ = µ ṗµ = 0). In more detail, to compute the gradient, we get the 1form which lives in the cotangent space dΦ ∈ Tp∗ RK , since the gradient lives in the tangent space, using the sharp map musical isomorphism #G : Tp∗ RK → Tp RK , the gradient becomes z = G−1 (p)dΦ(p) where z ∈ Tp RK . To compute the projection onto the tangent space of the simplex Tp ∆K−1 (since we want the evolution law of the ODE to remain on the simplex), we use the orthogonal projector map ΠG : Tp RK → Tp ∆K−1 , which takes z to the closest vector in Tp ∆K−1 . That is, given a vector z, we solve for ΠG (z) = arg minṗ∈Tp ∆K−1 ||ṗ − z||2G , which returns the vector in Tp ∆K−1 satisfying the minimum distance condition. Written in variational form L(ṗ, λ) =
1 ||ṗ − z||2G + λ1⊺ ṗ, 2
its minimum corresponds to the projected vector. Thus, solving for the minimum condition ∇ṗ L = 0 = G(ṗ − z) + λ1 = 0, we get that: ṗ = z + λG−1 1.
1⊺ z 1⊺ G−1 1
.
0 = (diag(p) − pp⊺ )F, Where F = f − T log p. One (degenerate) way to satisfy this condition is to make (diag(p) − pp⊺ ) = 0 where we get two conditions pµ = p2µ and pµ pν = 0. It can only be satisfied when p is a one-hot vector, thus, the solution lies on the boundary (vertices) of the simplex. Solutions inside the simplex (where pµ > 0 ∀µ) satisfy diag(p)F = pp⊺ F. Expressing it element-wise, since pµ > 0, dividing by pµ 0 = Fµ −
(17)
∀µ,
pν Fν
using that Fµ = fµ − T log pµ and solving for pµ " # K PK X 1 1 log pµ = fµ − pν Fν ⇒ pµ = e T [fµ − ν pν Fν ] . T ν (20) PK Imposing µ pµ = 1 gives 1
1 = e− T
PK ν
pν F ν
K X
1
e T fµ ,
µ
P where taking logarithms, we can identify ν p ν Fν = P 1 fµ T and we can recognize the partition funcT log µ e P 1 tion Z = µ e T fµ . Finally we can rewrite Eq.20 as
P where 1⊺ G−1 1 = i pi = 1. Plugging λ into Eq. (16) and using z = G−1 dΦ: ṗ = [I − G−1 11⊺ ]G−1 dΦ,
K X ν
(16)
The condition of a tangent vector to the simplex is 1⊺ ṗ = 0, we obtain λ=
replicator) as the entropic term T log p prevents winnertake-all dynamics for T > 0. The fixed points of Eq.19 satisfy
pµ =
1 1 fµ eT , Z
where in vector form p∗ = softmax
(21) 1 ∗ Tf
.
and developing the expression, the ODE becomes ṗ = (diag(p) − pp⊺ )dΦ.
Now that we have our gradient flow ODE, we choose 1 our 0-form to be Φ(p) = T DKL (p||π) with π = Z1 e T f , such that Φ(p) = −⟨p, F ⟩ + T log Z, where F = f − T log p and the differential is given by dΦ = T log p + T 1 − f (where neither log Z nor f depend on p). Substituting into Eq.18, we can write the gradient flow as ṗ = −(diag(p) − pp⊺ )(T log p − f ),
b.
(18)
(19)
where T 1 is proportional to 1 and (diag(p) − pp⊺ )1 = 0. We refer to this equation as the s-replicator (for ”soft”
Neuron domain
Here we provide a derivation of the classical Hopfield dynamics starting from the energy. We see that the metric Hx = diag(ϕ′ (x)) emerges as a natural choice. The Hopfield model energy [1] is a scalar function E : RN → R expressed as N
N
X 1X E(x) = − ϕ(xi )W ij ϕ(xj ) + [xi ϕ(xi ) − L(xi )], 2 ij i (22) i
R xi
where L(x ) = ϕ(u)du and ϕ(x) = tanh(k x), thus L(xi ) = k1 log cosh σ xi . It can be shown that this potential induces a global gradient flow similarly to [30]. To show this, we compute the differential of the energy
9 P ∂E i ∗ N and convert it to a dE = i ∂x i dx lying in dE ∈ Tx R gradient flow using the sharp map just like before using the metric Hx = diag(ϕ′ ). Differentiating the first term E1 , we get ∂xk E1 = −
X
=−
X
neuron activation term. From this potential we can derive the dynamics of our system by taking the differential ∗ living in dL(x, p) ∈ Tx,p (RN × ∆K−1 ) such that: dL(x, p) =
k
[∂ x ϕ(xi )W ij ϕ(xj ) + ϕ(xi )W ij ∂xk ϕ(xj )]
N X ∂L i
dxi + ∂xi
K X ∂L µ
∂pµ
dpµ .
(25)
ij
ϕ(xi )W ki ϕ′ (xk ),
i
(23) where we have assumed that W = W ⊺ . For the second term E2 , we get:
K
K∂xk ⟨p, f ⟩ =
∂xk E2 = xk ϕ′ (xk ).
(24)
#
∂xk E = x −
X
W ϕ(x ) ϕ′ (xk ), ik
i
i
or in vector form
K ⊺ Ξ diag(p)Ξ N Putting everything together, the two terms of Eq.25 become: W (p) =
∂x L(x, p) = diag(ϕ′ )[x − W (p)ϕ] ∂p L(x, p) = K(T log p + T 1 − f ).
dE = diag(ϕ′ )[x − W ϕ]. Taking the sharp map #H : Tx∗ RN → Tx RN , we get ẋ = −H −1 dE, where the metric necessarily becomes Hx = diag(ϕ′ (x)). Note that 1/ϕ′ (x) > 0 since in our case ϕ′ = sech2 . Because Hx is diagonal, it also satisfies hx (u, u) = 0 if u = 0 and hx (u, v) = hx (v, u) for all x. This is exactly the same as [30] up to the change of variable z = ϕ(x), where z is their dynamical variable, this way, ż = ϕ′ (ϕ−1 (z))[ϕ−1 (z)−W z] is exactly ϕ′ ẋ = ϕ′ [x−W ϕ].
c.
N
K XX pµ ξµi ξµk ϕ(xi )ϕ′ (xk ), N µ i
P k ′ k where in vector form K µ ξµ pµ ⟨ξµ , ϕ(x)⟩ϕ (x ) = N ′ diag(ϕ )W (p)ϕ(x), where we define the weight matrix as:
Combining both terms, " k
For the ∂p L(x, p) term, it is exactly the same as in Eq.19. For the first derivative, we have the interaction and activation terms. The activation term is exactly the same as in Eq.24, and for the interaction term we have
Joint potential
Finally, we can combine the parts to define the joint potential function L : RN ×∆K → R for the whole system from which we can derive the dynamics as a gradient flow. We define it as L(x, p) = −KT log Z + KT DKL (p||π) + [⟨x, ϕ⟩ − L(x)], which combines the two potentials we have described so far. From Eq. 1 a, we know that T DKL (p||π) = −⟨p, F ⟩+ T log Z and we can express it like
Finally, to get the gradient flow, the metric of the whole space is the direct sum M (x, p) = H(x)⊕G(p) and the orthogonal projection acts only on the astrocyte subspace, thus Π = I ⊕ ΠG = I ⊕ (I − G−1 1⊺ 1G−1 ). The natural gradient descent of L(x, p) can finally be expressed as −1 τx ẋ I 0 H 0 ∂x L(x, z) =− . τp ṗ 0 ΠG ∂p L(x, p) 0 G−1 Thus, the dynamics are as intended τx ẋ = −x + W (p)ϕ(x) τp ṗ = (diag(p) − pp⊺ )(f − T log p),
P i 1 where fµ = 2N . In this form, the interaction i ξµ ϕ P 1 term ⟨p, f ⟩ = 2N µ pµ ⟨ξµ , ϕ⟩2 is in direct analogy to P 1 1 2 µ ⟨ξµ , ϕ⟩ of Eq.22, where now each term 2 ϕW ϕ = 2N comes multiplied by the modulation pµ . The second term is the astrocytic regularization and the last one is the
(28) (29)
where we renormalize the time constant τp /K → τp . This completes the proof that the system is a gradient flow by construction. 2.
Dynamical analysis of the system:
Now that we have defined the evolution of the system, we can explore both its symmetries as well as its longterm behavior in different parameter regimes to better understand its dynamics.
L(x, p) = −K⟨p, f ⟩ + KT ⟨p, log p⟩ + [⟨x, ϕ⟩ − L(x)], i 2
(26) (27)
3. a.
Symmetry remarks
Z2 invariance of the squared-overlap score
In the main model the pattern score is a squared overlap, fµ (x) ∝ ⟨ξµ , ϕ(x)⟩2 . Consequently, the gain dynamics are invariant under the sign flip ⟨ξµ , ϕ(x)⟩ 7→
10 −⟨ξµ , ϕ(x)⟩. Equivalently, defining mµ (x) ≡ ⟨ξµ , ϕ(x)⟩, the gain update depends only on m2µ and cannot distinguish mµ from −mµ . More generally, if two patterns satisfy fµ (x) = fρ (x) at a given state, the gain update has no instantaneous preference between them. b.
Symmetry breaking by neuronal dynamics
Even with the squared-overlap score, degeneracies such as fµ (x) = fρ (x) are typically lifted by the neuronal evolution. To see this, consider the early-time dynamics with uniform gains p = 1/K, for which K
τx ẋi = −xi +
1 X i ξ mµ , N µ=1 µ
Differentiating mµ gives ṁµ = τx ṁµ = −
N X
ξµi ϕ′ (xi )xi +
i=1
mµ ≡
N X
ξµj ϕ(xj ).
j=1 i ′ i i i ξµ ϕ (x )ẋ , hence
P
K N X 1 X mν ξµi ϕ′ (xi )ξνi . N ν=1 i=1
Taking the difference between two overlaps yields N X
K X
N X
1 mν ∆ξ i ϕ′ (xi )ξνi , N ν=1 i=1 i=1 (30) where ∆ṁ = ṁµ − ṁρ and ∆ξ i = ξµi − ξρi . Except for nongeneric trajectories where the right-hand side vanishes identically, ∆ṁ ̸= 0 and the degeneracy is broken over time. Since f˙µ = (mµ /N )ṁµ , this induces a gain reweighting in the astrocytic dynamics. In practice, this mechanism implies that equal-fitness ties are generically transient and are resolved as x(t) evolves away from the initialization. τx ∆ṁ = −
4.
∆ξ i ϕ′ (xi )xi +
Fixed points of the system
Given the system dynamics of Eq.28 fixed point conditions: x∗ = W (p∗ )ϕ(x∗ ) f (x∗ ) ∗ p = softmax , T
(31) (32)
where the second equation comes from Eq.21. Plugging the second equation into the first, it can be rewritten as a single equation on the x domain. In this general regime, there are no analytic solutions, but we can see how the system behaves in different regimes: either hold the temperature fixed and take asymptotic limits of the time constants, or instead fix the time constants and examine the temperature limits. Case τp → 0: In this case we instantly get p∗µ = softmax 2N1 T ⟨ξµ , ϕ(x)⟩2 , so the evolution is written like: τx
dx = −x + W (p) ϕ(x) dt
This corresponds to a biased Hopfield network P in which every memory is weighted by p∗µ like W = µ p∗µ ξµ ξµ⊺ . Since p∗µ is highly heterogeneous when |⟨ξµ , ϕ⟩| is high, we are biasing the Hopfield evolution matrix towards those patterns which are more likely to contain the pattern we are searching for. Case τp → ∞: If we take τp → ∞, the astrocytic modulation is frozen ṗ = 0, that is, p(t) = p̄ is constant. Since the astrocytic influence is initialized by the uniform vector p = 1/K, we recover the classical Hopfield dynamics τx ẋ = −x + W (p̄) ϕ(x),
(33)
where with diag(p̄) = diag(1/K) and thus τx ẋ = −αx + W ϕ(x). Case τx → 0: In this case, the neurons evolve instantaneously to their fixed point x∗ = W (p)ϕ(x∗ ). This way, we have τp ṗ = (diag(p) − pp⊺ )(f − T log p), where fµ = 2N1 T ⟨ξµ , ϕ(x∗ )⟩2 is a transcendental function that depends also on p. Since the initial condition for the modulation equation is the uniform p = 1/K, the neurons initial configuration is already at a fixed point of the classical Hopfield dynamic x∗ = W (1/K)ϕ(x∗ ). For a sufficiently large number of patterns, the initial fixed point will likely be a spurious minima. This means, in the very first instants of time, since x(t0 ) = x∗ , the fitness f (x∗ ) of astrocytic modulation, will select those patterns that are the most similar to the ϕ(x∗ ) configuration, which can be arbitrarily far from x0 and thus even becoming detrimental for retrieval. Case τx → ∞: Here we have ẋ = 0, meaning that x(t) is constant it cannot evolve, so we stay forever in the initial condition even if the astrocyte dynamics converge to a fixed point in the p domain. Case T → ∞: Keeping now the time scales constant, in this regime, the regularization effect on the s-replicator is so strong that the system stays in the p = 1/K solution. Taking the astrocyte component, we can rearrange for the temperature to obtain X X τp dpµ 1 = pµ (fµ − pν fν ) − pµ (log pµ − pν log pν ). T dt T ν ν Taking T → ∞, we see that log pµ = ⟨log p⟩ where we are forced into a maximum entropy distribution. This way we recover the classical Hopfield model. Case T → 0: Similarly, rearranging for temperature, if we take T → 0, the dynamics become τp
X dpµ = pµ (fµ − pν fν ), dt ν
in which case there is no regularization. The fixed points
11 of the system are given by
1 Ξdiag(p∗ )Ξ⊺ ϕ(x∗ ) α 1M (x∗ ) , p∗ = |M (x∗ )|
x∗ =
5.
(34) (35)
where M (x) is the set of pµ elements that satisfy M (x∗ ) = arg maxµ fµ (x∗ ) and 1M (x∗ ) ∈ {0, 1}K . Note that M (x∗ ) will generally have a single element unless in the contrived case discussed before in Eq.(30).
[1] J. J. Hopfield, Proceedings of the National Academy of Sciences 79, 2554 (1982). [2] D. J. Amit, Modeling Brain Function: The World of Attractor Neural Networks (Cambridge University Press, 1989). [3] D. Krotov and J. J. Hopfield, in Advances in Neural Information Processing Systems, Vol. 29 (2016) pp. 1172– 1180. [4] D. Krotov and J. Hopfield, in International Conference on Learning Representations (2021) arXiv:2008.06996. [5] A. Vaswani, N. Shazeer, et al., in Advances in Neural Information Processing Systems (2017). [6] 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), arXiv:2008.02217 [cs.NE]. [7] A. Araque, R. P. Sanzgiri, V. Parpura, and P. G. Haydon, Canadian Journal of Physiology and Pharmacology 77, 699 (1999). [8] G. Perea, M. Navarrete, and A. Araque, Trends in Neurosciences 32, 421 (2009). [9] M. Letellier, Y. K. Park, T. E. Chater, P. H. Chipman, S. G. Gautam, T. Oshima-Takago, and Y. Goda, Proceedings of the National Academy of Sciences of the USA 113, E2685 (2016). [10] M. De Pittà, N. Brunel, and A. Volterra, Neuroscience 323, 43 (2016). [11] C. Giaume, A. Koulakoff, L. Roux, D. Holcman, and N. Rouach, Nature Reviews Neuroscience 11, 87 (2010). [12] A. Verkhratsky and M. Nedergaard, Physiological Reviews 98, 239 (2018). [13] S. S. Purushotham and Y. Buskila, Frontiers in Network Physiology 3, 1205544 (2023). [14] S. Makovkin, E. Kozinov, M. Ivanchenko, and S. Gordleeva, Scientific Reports 12, 6970 (2022). [15] L. Thompson, J. Khuc, M. S. Saccani, N. Zokaei, and M. Cappelletti, Experimental Brain Research 239, 2711 (2021). [16] K.-I. Dewa, K. Kaseda, A. Kuwahara, H. Kubotera, A. Yamasaki, N. Awata, A. Komori, M. A. Holtz, A. Ka-
Simulations details
To integrate Eq. (6), we use an explicit Euler method. Trajectories are computed for 10/dt time steps, with dt = 0.001, using smaller steps for small parameter value regimes (i.e. if α is the parameter then dt = α · 0.05 if α ≤ 0.01). For statistically significant results, we run 50 simulations with different random pattern matrices Ξ for every different parameter configuration. The observables are computed using the values of the last configuration values x(tf ), p(tf ). To plot the percentile bands, we have smoothed out noise using a 1d Gaussian filter for better visualization. The code to reproduce the simulations is available at [31].
sai, H. Skibbe, N. Takata, T. Yokoyama, M. Tsuda, G. Numata, S. Nakamura, E. Takimoto, M. Sakamoto, M. Ito, T. Masuda, and J. Nagai, Nature 10.1038/s41586025-09619-2 (2025), online ahead of print. [17] S. Zbaranska and S. A. Josselyn, Cell Research 35, 241 (2025). [18] S. A. Josselyn and S. Tonegawa, Science 367, eaaw4325 (2020). [19] L. Kozachkov, K. V. Kastanenka, and D. Krotov, Proceedings of the National Academy of Sciences 120, e2219150120 (2023). [20] L. Kozachkov, J.-J. Slotine, and D. Krotov, Proceedings of the National Academy of Sciences 122, e2417788122 (2025). [21] M. A. Di Castro, J. Chuquet, N. Liaudet, K. Bhaukaurally, M. Santello, D. Bouvier, P. Tiret, and A. Volterra, Nature Neuroscience 14, 1276 (2011). [22] G. Perea and A. Araque, Science 317, 1083 (2007). [23] C. Agulhon, T. A. Fiacco, and K. D. McCarthy, Science 327, 1250 (2010). [24] E. Shigetomi, S. Patel, and B. S. Khakh, Trends in Cell Biology 26, 300 (2016). [25] D. Krotov and J. Hopfield, arXiv preprint arXiv:2008.06996 (2020). [26] R. A. Jacobs, M. I. Jordan, S. J. Nowlan, and G. E. Hinton, Neural computation 3, 79 (1991). [27] S. Betteti, G. Baggio, F. Bullo, and S. Zampieri, Science Advances 11, eadu6991 (2025). [28] M. K. Benna and S. Fusi, Nature Neuroscience 19, 1697 (2016). [29] P. Mertikopoulos and W. H. Sandholm, Journal of Economic Theory 177, 315 (2018). [30] A. Halder, K. F. Caluya, B. Travacca, and S. J. Moura, IEEE Transactions on Neural Networks and Learning Systems 31, 4869 (2020). [31] A. Vivet and A. Arenas, Astrocyte-gated associative memory: code repository, GitHub repository (2026), accessed 10 Feb 2026. https://github.com/arnauvivett/Astrocyte-gatedassociative-memory-.