Model

Our model is based on a standard framework of firing-rate neural networks19. It consists of N recurrently connected neurons, with W ij the synaptic connection strength from neuron j to i. Each neuron i transforms its input x i into firing rate via a nonlinearity ϕ(x i ), where the state vector x = (x 1 x 2 ⋯ x N )⊤ evolves as

$$\mathop {{\mathbf{x}}}\limits^. = - {\mathbf{x}} + {\mathbf{W}}\phi \left( {\mathbf{x}} \right) + {\mathbf{b}}\left( t \right),$$ (1)

and b is an external input. Here and below we denote by ϕ(x) the vector obtained by applying ϕ to each coordinate of x.

Connectivity of task-performing networks is often designed to achieve a desired functionality, and assumed to be constant while the network is performing the task18,20. There are models in which connectivity co-evolves with neural dynamics, but changes are usually confined to a training phase, whereas connectivity is kept constant during the test phase21,22. These models are consistent with the expectation of synaptic tenacity in the absence of learning. In our model, to incorporate the recent observations on synaptic fluctuations, the connectivity matrix W continuously co-evolves with neural activity x throughout all task phases, albeit with a slower timescale (Fig. 1, see also ref. 23).

Fig. 1 Co-evolution of neural activity and connectivity. a Illustration of our modeling framework: a recurrently connected neural network, with W ij denoting the connection strength from neuron j to neuron i. The dynamic variables x i evolve by Eq. (1) and connection strengths W ij evolve by Eq. (2). An external signal b i can be added as input to each neuron i. b Two example traces of neural state dynamics (x i (t); top, colors) and of connection strength dynamics (W ij (t); bottom, colors) are shown. Both evolve over time, though on different timescales, and their dynamics are coupled Full size image

In order to study the coexistence of memory with synaptic fluctuations, we let W evolve due to contributions arising from both learning-related and fluctuation-related terms, denoted by Δ L and Δ F respectively:

$$\mathop {{\mathbf{W}}}\limits^. = \eta \left( {{\mathrm{\Delta }}_L + {\mathrm{\Delta }}_F} \right),$$ (2)

with η > 0 the plasticity rate (relative to neural dynamics). The fluctuation term includes stochastic, activity-independent noise in synaptic strength, as well as a homeostatic mechanism to control synaptic and firing-rate stability. These are precisely the processes which endanger the resilience of an acquired memory that is assumed to be stored in synaptic patterns. We first consider how an existing memory is eroded by these processes, and later address the learning part and the interplay between the two.

Homeostatic plasticity erodes real-coded information

We model spontaneous activity-independent synaptic fluctuations by a white noise process ξ ij driving each synapse ij independently. The variance of such a process grows without bound, and thus, without a restraining mechanism these fluctuations would lead to divergence of the synaptic weights W ij . A plausible restraining agent is homeostatic plasticity, modeled here as a network-level mechanism that stabilizes both synaptic weights and neural activities on average over long timescales. To understand this stabilization, we turn to dynamical systems theory that determines stability about a set-point by the spectrum of the appropriate Jacobian matrix (which is the local linear approximation of the dynamics). This spectrum consists of a set of complex eigenvalues—a collection of points in the complex plane, each having a real part and an imaginary part.

In general, the real part of this spectrum defines the system’s stability: a system is only stable if all its eigenvalues have negative real parts. The imaginary part of the spectrum, in contrast, determines the typical timescales of small-amplitude dynamics around this set-point, but not stability itself (Fig. 2a). Therefore, while the real part of the spectrum must be under the control of homeostatic plasticity, its imaginary part is not constrained by the requirement of stability, and is free to store information (Fig. 2b).

Fig. 2 Stability and memory-associated connectivity eigenvalues. a Eigenvalues of the Jacobian matrix occupy the complex plane. System stability about a set-point is ensured if all eigenvalues have negative real parts (i.e., reside in the blue half-plane). Eigenvalues with positive real parts (in the red half-plane) correspond to locally unstable directions. b Homeostatic control that prevents noise from accumulating and causing divergence, also affects stored memories represented by the eigenvalues. However, such mechanisms must only control their positive real parts, pushing them to be negative (arrows), while in the left half-plane large-amplitude imaginary eigenvalues can persist. c–e Real (top) and imaginary (bottom) parts of connectivity (W) eigenvalue spectra for systems evolving under the dynamics of Eqs. (1) and (2), with different homeostatic mechanisms in Δ F (N = 128, see “Methods” for simulation details). For each case, a memory item was embedded in W at time t = 2500 (green trajectories), by inducing a real outlier eigenvalue (top panels, generated by the perturbation δW = uu⊤) or an imaginary pair (bottom panels, generated by δW = uv⊤ − vu⊤). c Dissipation of synapses. Both real and imaginary memories decay with the same rate. d Homeostatic rate-control. e Decorrelation homeostasis. In the last two, real-coded memory (top) decays rapidly whereas imaginary-coded (bottom) persists Full size image

The arguments above derive from a general intuition on system stability; they are not a mathematical proof, as they depend on the existence of a set-point and its exact properties. They do, however, provide motivation to test this idea using various homeostatic mechanisms. We perform such tests using the connectivity matrix W as a proxy for the Jacobian. In the case of a linear network, or of linearizing around the origin, the two are equivalent. Our results below indicate that such an approximation is useful also in more general cases.

The embedding of a memory item is often represented in learning theory as a low-rank perturbation to the connectivity matrix, W → W + δW. For example, in the Hopfield model a memory is associated with a particular pattern of activity, a vector u ∈ ℝN, and is embedded in connectivity by adding a perturbation of the form δW = uu⊤. Such a perturbation is symmetric, δW ij = δW ji , and modifies the connectivity spectrum to include a real positive eigenvalue (with an eigenvector in the direction of u). To tap into the resilience of the imaginary part of the spectrum as argued above, one might add an anti-symmetric perturbation (δW ij = −δW ji ) that gives rise to an imaginary conjugate pair of eigenvalues. This defines a different type of memory item, the simplest form being δW = uv⊤ − vu⊤. If the above general arguments on system stability are correct, such memory items should be more resilient to synaptic fluctuations. We test this by comparing the erosion of the two types of memory items under various homeostatic mechanisms. We first embed memories corresponding to either real or imaginary eigenvalues into the connectivity matrix W, and then follow the dynamics of Eqs. (1) and (2) without active learning (Δ L = 0), but with various homeostatic models in Δ F .

Perhaps the simplest implementation of a homeostatic mechanism is by dissipative synaptic dynamics. Together with the noise ξ, this gives a fluctuation term

$$\Delta _F = \xi - \beta {\mathbf{W}},$$ (3)

with β > 0 the rate of dissipation. Figure 2c shows the eigenvalues of the connectivity matrix as a function of time (gray lines), with the eigenvalues corresponding to the memory highlighted in green. It is seen that the memory representation rapidly decays for both real (top) and imaginary (bottom) eigenvalues. This is expected from a dissipative system, where all information decays exponentially with a rate β. Therefore, in the presence of such a mechanism, neither type of memory items can be sustained for longer than the decay time 1/β. However, as will be shown below, this is not the case for more indirect homeostasis mechanisms.

A biologically plausible homeostasis mechanism can be modeled as an activity-dependent rule—where the synaptic matrix is modified to achieve a stable post-synaptic firing-rate24,25,26:

$$\Delta _F = \xi + \left( {\phi _0 - \phi ({\mathbf{x}})} \right)\phi ({\mathbf{x}}^ \top ) \circ {\mathbf{W}},$$ (4)

with ϕ 0 an arbitrary target-rate vector, and \(\circ\) denoting a Hadamard (element-wise) product. Stabilization of firing rates around the set-point ϕ 0 is achieved in this case by scaling the weight of each synapse by a factor dependent on its pre- and post-synaptic neurons. Although not explicitly symmetric (which would lead to a change in real eigenvalues), the stabilization induced by this rule requires control over the real part of the relevant Jacobian. Accordingly, Fig. 2d (top) shows that memories stored as real eigenvalues of W rapidly decay. Imaginary-coded memories, on the other hand, may persist indefinitely without interfering with homeostasis (Fig. 2d, bottom).

Finally, inspired by ref. 27, we consider a homeostasis mechanism that does not have a well-defined firing-rate set-point. Instead, this rule contains an anti-Hebbian term that reduces the connections between correlated neurons, thus pushing towards decorrelated firing rates across the network:

$$\Delta _F = \xi + {\mathbf{I}} - \phi _{{\mathrm{post}}}({\mathbf{x}})\phi _{{\mathrm{pre}}}({\mathbf{x}}^ \top ),$$ (5)

where ϕ pre , ϕ post are two sigmoidal functions and I is the identity matrix. This rule leads to persistent and constrained fluctuations in both connectivity and firing rates27. The activity x is dominated by the unstable modes of W, which are then suppressed by the anti-Hebbian term, leading to a new unstable mode in an endless succession.

Once again, we find that the decay of imaginary-coded memories is orders of magnitude slower than that of real-coded ones (Fig. 2e). Note that if the sigmoidal functions are identical, ϕ pre = ϕ post , this rule can only modify the symmetric part of W. In practice, for many non-identical choices of the sigmoid functions, the modification is still mostly symmetric. Nevertheless, the relative decay of imaginary- and real- based memories is similar to the case of the rate-control rule, that does not have any symmetric tendency.

Some evidence for the generality and limits of validity of these results is presented in Supplementary Fig. 1, where sparse networks are considered. The amplitude of embedded imaginary memories decreases smoothly as the network becomes sparser; they remain more resilient than real-coded memories for the entire regime tested.

In light of these results, a natural question arises: can a dynamical learning rule utilize the imaginary subspace to robustly code and store memory representations? Perhaps surprisingly, we find that Spike Time Dependent Plasticity (STDP), a well documented and biophysically plausible learning rule, provides a natural candidate for learning imaginary-coded memories.

STDP stores imaginary-coded information

Symmetric and anti-symmetric matrices give rise to real and imaginary eigenvalues respectively. It is thus reasonable that an anti-symmetric modification to the synaptic weight matrix W would primarily lead to changes in the imaginary part of its spectrum. Local learning rules observed in experiments (e.g., STDP) have a well-defined directionality: consecutive firing of neuron j before i leads to a strengthening of the connection W ij and to the weakening of the reverse connection. The temporal asymmetry of STDP28 leads to an approximately anti-symmetric learning rule when applied to our rate model (see “Methods”); as such, this rule mostly affects the imaginary part of the spectrum. In the case of perfect anti-symmetry, we find the form Δ L = ϕy⊤ − yϕ⊤, which modifies only the anti-symmetric component of W. The vector y is a smoothed version of the activity ϕ(x), and arises in our rate-based formulation when applying an exponential STDP kernel to spike-trains (see “Methods”).

These arguments suggest that a biologically motivated learning rule naturally stores imaginary-coded information, thereby rendering it relatively resilient to the effect of homeostatically controlled synaptic fluctuations. We will next investigate how such a memory can be acquired, retained and retrieved in the presence of synaptic fluctuations. For simplicity, we will use the purely anti-symmetric Δ L .

The acquisition, i.e., the encoding and storage of a new memory trace, is initiated by stimulating the network with an external signal, b(t). A matrix with imaginary eigenvalues is necessarily of (at least) rank 2, corresponding to a two-dimensional space spanning the memory representation. We therefore present the network with a randomly time-varying input evolving on a plane spanned by two arbitrary directions u, v ∈ ℝN (see “Methods”). As the input drives neural activity x onto the (u, v) plane, the activity-dependent learning operator Δ L follows and becomes non-negligible, which in turn causes a change in connectivity.

The learning procedure stores geometric information of the external stimulus, specifically the directions u and v, within the anti-symmetric part of the connectivity matrix. This encoding is manifested as a rank-2 operator uv⊤ − vu⊤ which is embedded into W. To see this, we follow the spectrum of W as a function of time: during stimulus presentation, a complex conjugate eigenvalue pair forms (Fig. 3a), with corresponding eigenvectors overlapping completely with the plane spanned by u, v (Fig. 3b). The strength of the memory representation—corresponding to the magnitude of the imaginary eigenvalue—depends monotonically on stimulation duration and on input amplitude. At later times, additional stimuli may be stored using the same learning protocol (Fig. 3c).

Fig. 3 Hebbian learning by STDP embeds persistent imaginary-coded memory. a Imaginary part of the spectrum of W before, during and after external stimulation (stim.; marked by the shaded green area) applied between times t = 100 and t = 200. Before learning, the imaginary part of the spectrum is almost constant in time. The learning of a memory item manifests as the growth in imaginary amplitude of one complex conjugate eigenvalue pair (green trajectories). b During stimulus presentation, the learning rule modifies W such that the plane spanned by u, v is invariant. Plotted are the overlaps of this eigenplane of W, corresponding to the largest imaginary eigenvalue pair, with N/2 planes: the u, v (green), and N/2 − 1 orthogonal planes (gray); see “Methods”. c Zoom-out of (a). After the stimulus is removed, the memory representation persists (green). A second stimulus, confined to a second plane, is similarly learned at a later time (blue); both memory items are retained. This figure shows results for the decorrelation rule (N = 128, see “Methods”); qualitatively similar results are obtained for rate-control (see Supplementary Fig. 2) Full size image

The nature of imaginary-coded memories

We have seen that a biologically plausible learning rule can capture the orientation in neuronal state-space of an incoming stimulus, and encode this information as a pair of imaginary eigenvalues in the network connectivity matrix. What is the nature of this memory in terms of network activity? We find that learning creates attractors in state-space, similar in fashion to those in the Hopfield model18. However, rather than fixed points, here the attractors are time-varying stable states—namely, limit cycles. To see this most clearly, we consider a single imaginary-coded memory embedded in the network, and examine neural dynamics while keeping W fixed. Following the Hopfield paradigm, we write W as:

$${\mathbf{W}} = \rho \left( {{\mathbf{uv}}^ \top - {\mathbf{vu}}^ \top } \right),$$ (6)

where the coefficient ρ > 0 represents the strength of the memory representation29.

With one stored memory as in Eq. (6), we find that, from any non-zero initial condition, the dynamics converge to periodic motion concentrated on the ‘memory plane’ spanned by u and v. Figure 4a depicts the projections of neural activity on this plane, for two initial conditions (light gray trajectories), both converging to the limit-cycle attractor (dark closed trajectory). An approximate low-dimensional description of this limit cycle can be obtained for an infinitely steep nonlinearity ϕ (i.e., a step-function). The full dynamics are then well approximated by their projected coordinates on the plane, p u and p v :

$$\begin{array}{l}\dot p_{\mathbf{u}} = - p_{\mathbf{u}} + \rho \,q_{\mathbf{v}}\\ \dot p_{\mathbf{v}} = - p_{\mathbf{v}} - \rho \,q_{\mathbf{u}},\end{array}$$ (7)

where \(q_{\mathbf{v}} \approx {\mathrm{arctan}}\left( {\frac{{p_{\mathbf{v}}}}{{|p_{\mathbf{u}}|}}} \right)\). This low-dimensional system exhibits a stable limit-cycle around the origin (see Supplementary Note 1 and Supplementary Fig. 3). We conclude that imaginary-stored memory items correspond to dynamic attractors, with geometry defined by that of the stimulating input. This behavior stands in contrast to the classic—symmetric—Hopfield model, where memories are represented by fixed-point attractors corresponding to fixed values of neural activity.

Fig. 4 Dynamics of memory retrieval with fixed connectivity. a Overlaps of network activity x with the two directions spanning one embedded memory plane, \(p_{\mathbf{u}}: = \frac{1}{{\sqrt N }}{\mathbf{u}}^ \top {\mathbf{x}}\) and \(p_{\mathbf{v}}: = \frac{1}{{\sqrt N }}{\mathbf{v}}^ \top {\mathbf{x}}\) form a stable limit cycle. Shown are two trajectories, one initiated inside and the other outside the stable orbit, with arrowheads showing the direction of temporal evolution (N = 4096, ρ = 4). b Radial coordinates of overlaps with each of multiple (M = 10, N = 4096) embedded memory planes, \(r_k: = \sqrt {p_{{\mathbf{u}}^{\left( k \right)}}^2 + p_{{\mathbf{v}}^{\left( k \right)}}^2}\). Initiating the network near one plane, u(1), v(1), results in convergence to the associated attractor (green). Overlaps with other planes remain small (black). In this figure we construct connectivity by adding to Eq. (6) a real-coded term (see “Methods”) Full size image

Embedding multiple memory planes \(\{ {\mathbf{u}}^{\left( k \right)},{\mathbf{v}}^{\left( k \right)}\} _{k = 1}^M\) corresponds to setting

$${\mathbf{W}} = {\mathbf{UDU}}^ \top$$ (8)

where the columns of U are the memory patterns (interleaved u(k) and v(k)), and D is a 2M × 2M block-diagonal matrix, with the k-th block reading \(\left( {\begin{array}{*{20}{c}} 0 & {\rho _k} \\ { - \rho _k} & 0 \end{array}} \right)\). Now, a locally stable limit-cycle lies on each embedded plane, and the network functions as an associative memory: initiating the dynamics within the basin of attraction of one plane—providing the network with partial information of the memory to be retrieved—leads to the recovery of the full memory item (Fig. 4b). Similar to the Hopfield model, the memory capacity is found to be proportional to system size30. Numerical simulations presented in Supplementary Fig. 4 show that in fact the proportionality constant is slightly higher compared with that of the symmetric Hopfield model (when normalized by a factor of two, since each memory resides on a plane; see Supplementary Note 2).

Life cycle of a memory trace

We next consider the entire life-cycle of a memory in the presence of synaptic fluctuations and homeostasis, starting from learning, through retention and to retrieval. During a learning event, implemented by presenting a stimulus in the two-dimensional memory plane, a memory representation is formed by the Hebbian learning rule. Figure 5a (left) shows the overlaps of neural activity onto the two planes, r 1 (green) and r 2 (blue), together with the stimulus which drives learning (shades). These projections are elevated during stimulation, which—via the Hebbian learning rule—modifies the synaptic matrix to store each of the planes in connectivity. The projection r 3 onto a third plane, which was not learned, is negligible as shown in the bottom line (orange).

Fig. 5 Stable memory with unstable synapses. a Left: Learning phase. Two distinct stimuli are sequentially presented to the network (green and blue shaded areas). Right: Retrieval phase. The embedded memories are read out of the network dynamics via a brief stimulation bearing partial information of the original stimulus. Overlaps of the network state x with the two learned planes are shown in corresponding color. During retrieval, activity follows the oscillatory trajectory of the dynamic memory item (blue zoom). A novel cue, on the other hand, does not elicit a significant response in neural activity (bottom trajectory, orange). b Relative contributions of activity-dependent (left bar) and spontaneous (right bar) synaptic fluctuations are estimated to be of similar magnitude during all phases of the memory life-cycle (see “Methods”). Error bars denote one standard deviation from the mean; overlaid data points are marked by green circles. c Effect of retrieval on existing memories can differ, with the green memory degrading, and the blue one strengthening. d Weights of four out of N2 synapses (gray shades) as a function of time. The spontaneous and homeostatic contributions to plasticity drive perpetual fluctuations of synaptic weights, which occur during and after learning. This figure shows results for the rate-control homeostasis rule (N = 128, see “Methods”) Full size image

After learning, the two memory items are stored as pairs of imaginary eigenvalues, remaining stable over time, until they are retrieved at times t 1 and t 2 , respectively. At retrieval, activity is transiently attracted to the respective memory planes, as indicated by the spikes in the overlaps (Fig. 5a, right). During retrieval, activity follows the stored dynamic trajectory, exhibiting its typical oscillations (Fig. 5a, right, blue zoom). At the same time, the projection onto an arbitrary plane shows no temporal structure (orange zoom). Finally, stimulating the network with a novel cue does not elicit a significant response in neural activity in any of the projections (t 3 in Fig. 5a).

The effect of retrieval on the connectivity, namely on the stored memory itself, is somewhat unpredictable and depends on the exact state of the network and on the memory properties. As an example, in Fig. 5c it is seen that the green memory is damaged by retrieval, namely the magnitude of the corresponding imaginary eigenvalue is decreased. This may be caused by the homeostatic mechanism that constrains activity, in particular the component projected onto the memory plane by the retrieval event. In contrast, the blue memory is slightly strengthened by retrieval, as seen by the increased magnitude of the eigenvalue pair. In other cases the memory remains unaffected.

Throughout this entire cycle, synapses fluctuate under the effect of activity-independent noise and homeostasis. Figure 5d shows a few example synapses tracked across time, during both phases. We may disentangle the two effects, spontaneous and activity-dependent, and estimate their relative contribution to synaptic fluctuations; their magnitude is found to be similar (Fig. 5b), in agreement with the experimentally observed phenomenon14.