Supplementary Information

Supplementary Figures¶

Supplementary Tables¶

Supplementary Methods¶

Study 4 instructions for upregulation. “During this scan, we are going to ask you to try to imagine as hard as you can that the thermal stimulations are more painful than they are. Try to focus on how unpleasant the pain is, for instance, how strongly you would like to remove your arm from it. Pay attention to the burning, stinging and shooting sensations. You can use your mind to turn up the dial of the pain, much like turning up the volume dial on a stereo. As you feel the pain rise in intensity, imagine it rising faster and faster and going higher and higher. Picture your skin being held up against a glowing hot metal or fire. Think of how disturbing it is to be burned, and visualize your skin sizzling, melting and bubbling as a result of the intense heat.”

Study 4 instructions for downregulation. “During this scan, we are going to ask you to try to imagine as hard as you can that the thermal stimulations are less painful than they are. Focus on the part of the sensation that is pleasantly warm, like a blanket on a cold day. You can use your mind to turn down the dial of your pain sensation, much like turning down the volume dial on a stereo. As you feel the stimulation rise, let it numb your arm, so any pain you feel simply fades away. Imagine your skin is very cool, from being outside, and think of how good the stimulation feels as it warms you up.”

Supplementary Information 1¶

Derivation of the joint steady-state probability of free energy minimizing attractor networks See Spisak & Friston, 2025 for details on the whole framework.

Outline¶

Setup: Continuous–Bernoulli kernel¶

We assume that the single-site conditional distribution is the Continuous–Bernoulli on $[-1,1]$ with canonical parameter

and density (for $x\in[-1,1]$)

with $h(x)=\tfrac12 \mathbf{1}_{[-1,1]}(x)$ and $A(\kappa)=\log(\sinh\kappa)-\log\kappa$ (so the $\kappa\to0$ limit is uniform on $[-1,1]$).

The conditional mean is:

1. Master equation¶

Let’s consider a precise formalization of the asynchronous update procedure: each site $i$ has an independent Poisson clock of rate $\gamma_i$. When it rings, $\sigma_i$ is replaced by a draw $x\sim K(\cdot\mid\kappa_i(\sigma_{-i}))$.

With $\sigma^{(i,x)}$ the configuration equal to $\sigma$ but with coordinate $i$ set to $x$, the transition density is

The master equation for $p(\sigma,t)$ is

2. Probability currents¶

Let’s denote the site-wise flux density as:

At steady state $p^\ast$, by definition we have:

3. Detailed balance condition¶

A special (and maybe the most intuitive) way to satisfy the previous eq. is the case of detailed balance or equilibrium. In this case, every transition is balanced:

It can be shown that:

I.e., the log-density $\Phi$ changes by the local slope at $i$ (see detailed derivation below).

Hence for $j\neq i$:

Mixed partials commute, thus:

Therefore, equilibrium (detailed balance) is possible only when the coupling is symmetric.

4. Non-equilibrium Steady State (NESS)¶

Detailed balance (equilibrium) is only one specific way the steady-state condition can hold:

Let’s consider the general case: an arbitrary (possibly asymmetric) $J$ coupling. $J$ can always be decomposed into:

This leads to a decomposition of the edge currents on each directed edge $(\sigma,\sigma^{(i,x)})$:

so that $F_i = F_i^{\mathrm{sym}} + F_i^{\mathrm{sol}}$ and $F_i^{\mathrm{sym}}(\sigma\to x)=F_i^{\mathrm{sym}}(\sigma^{(i,x)}\to\sigma)$.

Since $F_i(\sigma\to x) = e^{\Phi(\sigma)} w_i(\sigma\to x)$, define rates accordingly:

so $w_i = w_i^{\mathrm{sym}} + w_i^{\mathrm{sol}}$ and $F_i^{\mathrm{sym}} = e^{\Phi}w_i^{\mathrm{sym}}$, $F_i^{\mathrm{sol}} = e^{\Phi}w_i^{\mathrm{sol}}$.

We proceed in two steps: (i) the antisymmetric component is divergence-free and never changes the steady state; (ii) the remaining symmetric component yields exactly the same $p^\ast$ as the equilibrium distribution obtained by setting $J=J^{\mathrm{S}}$. Finally, as a bonus, we show that the induced circulating flow is indeed tangential to the level sets of $\Phi$, in accordance to the Helmholtz decomposition.

i) Antisymmetric component is divergence-free (does not change $p^\ast$).¶

Starting again from stationarity of the full current,

subtract the same identity written with $F^{\mathrm{sym}}$ in place of $F$. Because $F^{\mathrm{sym}}(\sigma\to x)=F^{\mathrm{sym}}(\sigma^{(i,x)}\to\sigma)$ on every edge, its contribution vanishes pairwise, giving

Thus $F^{\mathrm{sol}}$ has zero divergence and does not change $p^\ast$.

ii) Symmetric component determines $p^\ast$ and matches the equilibrium solution with $J=J^{\mathrm{S}}$.¶

From Section 3, only $J^{\mathrm{S}}$ can appear in a potential gradient (commuting mixed partials). Therefore

which integrates to the same $\Phi=\log p^\ast$ you obtain by enforcing detailed balance with the symmetric coupling $J=J^{\mathrm{S}}$ (i.e., in the absence of $J^{\mathrm{A}}$). Equivalently, the symmetric current $F^{\mathrm{sym}}$ is reversible with respect to $p^\ast$ and alone reproduces the same steady state.

Note: Tangency to level sets of $\Phi$ (no work against the potential).¶

Multiply the preceding identity by a test function $\psi(\sigma)$ and integrate over $\sigma$ (“summation by parts”) to obtain

Taking $\psi=\Phi$ yields

so the solenoidal current is orthogonal to discrete gradients and flows tangent to the isocontours of $\Phi$ (the level sets of $p^\ast$).

3. Link to $J^{\mathrm{A}}$.

Under the linear Continuous–Bernoulli parametrization, $J^{\mathrm{A}}$ cannot contribute to any scalar potential (mixed partials must commute). Therefore, varying $J^{\mathrm{A}}$ at fixed $J^{\mathrm{S}}$ leaves $\Phi$ and $p^\ast$ unchanged; it only affects the circulating, divergence-free component $F^{\mathrm{sol}}$.

Conclusion: explicit steady-state $p^\ast$ for arbitrary $J$¶

By the results above, only the symmetric component $J^{\mathrm{S}} := \tfrac12(J+J^\top)$ can shape the potential $\Phi=\log p^\ast$. Therefore, for arbitrary (possibly asymmetric) $J$, the non-equilibrium steady-state density takes the Boltzmann-like form:

with $J^{\mathrm{S}}_{ii}=0$ and partition function

The antisymmetric component $J^{\mathrm{A}} := \tfrac12(J-J^\top)$ contributes only to divergence-free (solenoidal) probability currents tangent to the level sets of $\Phi$ and thus does not alter $p^\ast$. This expression matches the manuscript’s stationary density, with $J^{\mathrm{S}}$ playing the role of the effective (symmetric) interaction matrix.

Supplementary Information 2¶

Here we provide a concise derivation of how $\sigma$ and $J$ changes as a consequence of free energy minimization resulting in an inference (node update) rule and a local, incremental learning rule, respectively.

For more background and a detailed derivation, see Spisak & Friston, 2025.

Let the instantaneous net input to node $\sigma_i$ be

and let $q(\sigma_i)\propto e^{b_q\sigma_i}$ be the variational marginal with mean $L(b_q)$, where $L$ is the Langevin function (the mean of the Continuous–Bernoulli).

Inference (Hopfield-style)¶

We start by writing up the local variational free energy from the point of view of a single node $\sigma_i$:

Next we express:

with mean $\mathbb{E}_q[\sigma_i]=L(b_q)=A'(b_q)$.

and

As it depends on $\sigma_i$ only through the linear slope $\sum_{j\ne i}J_{ij}\sigma_j$:

Now, we assemble the local free energy (dropping constants independent of $b_q$):

Equivalently, using $u_i=b_i+\sum_{j\ne i}J_{ij}\sigma_j$ and $A'(b_q)=L(b_q)$,

3. Differentiate w.r.t. $b_q$ and use $A'(b_q)=L(b_q)$ to cancel terms:

Setting the derivative to zero gives $b_q^\star=u_i$ and therefore the node-wise update

which is the stochastic Hopfield/Boltzmann-style activation with the Langevin nonlinearity.

Learning (Hebb − anti-Hebb)¶

Make the dependence on $u_i$ explicit by adding and subtracting the CB log-partition $A(u_i)$ (using $\ln p(\sigma_i\mid\sigma_{\setminus i})=u_i\,\sigma_i - A(u_i)+\ln h(\sigma_i)$ and Bayes’ rule):

Now

and by the chain rule with $\partial u_i/\partial J_{ij}=\sigma_j$ we obtain

Using a stochastic (sample-based) estimate $\mathbb{E}_q[\sigma_i]\approx\sigma_i$ and descending $F$ gives the local update

This is a predictive coding-based learning rule that strengthens observed correlations and subtracts those already predicted by the current model.

Supplementary Information 3¶

Self-orthogonalization of attractor states

To illustrate how the learning rule in free energy minimizing attractor networks gives rise to efficient, (approximately) orthogonal representations of the external states, suppose the network has already learned a pattern $\mathbf{s}^{(1)}$, whose neural representation is the attractor $\boldsymbol{\sigma}^{(1)}$ and associated weights $\mathbf{J}^{(1)}$. When a new pattern $\mathbf{s}^{(2)}$ is presented that is correlated with $\mathbf{s}^{(1)}$, the network’s prediction for $\sigma_i^{(2)}$ will be $\hat{\sigma}_i = L(b_i + \sum_{k \neq i} J_{ik}^{(1)} \sigma_k)$. Because inference with $\mathbf{J}^{(1)}$ converges to $\boldsymbol{\sigma}^{(1)}$ and $\boldsymbol{\sigma}^{(2)}$ is correlated with $\boldsymbol{\sigma}^{(1)}$, the prediction $\hat{\boldsymbol{\sigma}}$ will capture variance in $\boldsymbol{\sigma}^{(2)}$ that is ‘explained’ by $\boldsymbol{\sigma}^{(1)}$. The learning rule updates the weights based only on the unexplained (residual) component of the variance, the prediction error. In other words, $\hat{\boldsymbol{\sigma}}$ approximates the projection of $\boldsymbol{\sigma}^{(2)}$ onto the subspace already spanned by $\boldsymbol{\sigma}^{(1)}$. Therefore, the weight update primarily strengthens weights corresponding to the component of $\boldsymbol{\sigma}^{(2)}$ that is orthogonal to $\boldsymbol{\sigma}^{(1)}$. Thus, the learning effectively encodes this residual, $\boldsymbol{\sigma}^{(2)}_{\perp}$, ensuring that the new attractor components being formed tend to be orthogonal to those already established.

For more details on self-orthogonalization in these networks, including an empirical demonstration, see Spisak & Friston, 2025.

Supplementary Information 4¶

Reconstruction of the attractor network from the activation timeseries.

We start from Eq. (1), written in matrix notation:
$E_{HN}(\bm{\sigma}) = -\frac{1}{2} \sigma^\top J \sigma + \sigma^\top b \\ = -\frac{1}{2} \sigma^\top J \sigma + \sigma^\top J J^{-1} b \\ = -\frac{1}{2} \sigma^\top J \sigma + \sigma^\top J J^{-1} b - \frac{1}{2} b^\top J^{-1} b + \frac{1}{2} b^\top J^{-1} b$

To complete the square, we added and subtracted $\frac{1}{2} b^\top J^{-1} b$ within the exponent. We recognize that $\sigma^\top J J^{-1} b = \sigma^\top b$. Now add and subtract $\frac{1}{2} b^\top J^{-1} b$

$= -\frac{1}{2} \left[\sigma^\top J \sigma - 2 \sigma^\top J J^{-1} b + b^\top J^{-1} b \right] + \frac{1}{2} b^\top J^{-1} b \\ = -\frac{1}{2} \left[\sigma - J^{-1} b \right]^\top J \left[\sigma - J^{-1} b\right] + \frac{1}{2} b^T J^{-1} b$

Given that $P(\bm{\sigma}) = \exp(-E_{HN}(\bm{\sigma}))$, the expression simplifies to: $P(\sigma) \propto \exp\left(-\frac{1}{2} (\sigma - J^{-1} b)^\top J (\sigma - J^{-1} b) \right)$

Note that the term $\frac{1}{2} b^\top J^{-1} b$ is independent of $\sigma$; we have absorbed it into the normalization constant.
This is exactly the exponent of a multivariate Gaussian distribution with mean $J^{-1} b$ and covariance matrix $J^{-1}$, meaning that the weight matrix of the attractor network can be reconstructed as the inverse covariance matrix of activation timeseries of the lower-level nodes: $\mathbf{J} = -\Lambda = -\Sigma^{-1}$