Theory of coupled neuronal-synaptic dynamics

In neural circuits, synaptic strengths influence neuronal activity by shaping network dynamics, and neuronal activity influences synaptic strengths through activity-dependent plasticity. Motivated by this fact, we study a recurrent-network model in which neuronal units and synaptic couplings are interacting dynamic variables, with couplings subject to Hebbian modification with decay around quenched random strengths. Rather than assigning a specific role to the plasticity, we use dynamical mean-field theory and other techniques to systematically characterize the neuronal-synaptic dynamics, revealing a rich phase diagram. Adding Hebbian plasticity slows activity in chaotic networks and can induce chaos in otherwise quiescent networks. Anti-Hebbian plasticity quickens activity and produces an oscillatory component. Analysis of the Jacobian shows that Hebbian and anti-Hebbian plasticity push locally unstable modes toward the real and imaginary axes, explaining these behaviors. Both random-matrix and Lyapunov analysis show that strong Hebbian plasticity segregates network timescales into two bands with a slow, synapse-dominated band driving the dynamics, suggesting a flipped view of the network as synapses connected by neurons. For increasing strength, Hebbian plasticity initially raises the complexity of the dynamics, measured by the maximum Lyapunov exponent and attractor dimension, but then decreases these metrics, likely due to the proliferation of stable fixed points. We compute the marginally stable spectra of such fixed points as well as their number, showing exponential growth with network size. In chaotic states with strong Hebbian plasticity, a stable fixed point of neuronal dynamics is destabilized by synaptic dynamics, allowing any neuronal state to be stored as a stable fixed point by halting the plasticity. This phase of freezable chaos offers a new mechanism for working memory.


I. INTRODUCTION
Computations in neural circuits are commonly thought to be implemented through the coordinated dynamics of neurons [1][2][3].Under this view, the role of synaptic connectivity is to sculpt neuronal dynamics to implement computations.In actuality, synapses undergo plasticity over diverse timescales in response to neuronal activity and thus constitute dynamic degrees of freedom in their own right [4].A more accurate picture of computation in neural circuits should involve the coupled dynamics of neurons and synapses.Indeed, it is possible that a network is better described by the states of its synapses than of its neurons [5].Here, we study the consequences of treating neurons and synapses as mutually coupled dynamic variables on equal footing.
Synaptic dynamics are often divided into short-term plasticity, which operates on short timescales and depends on presynaptic activity [6][7][8][9][10][11]; and long-term plasticity, which acts on much longer timescales and depends on both pre-and postsynaptic activity [12][13][14][15].However, short-term forms of Hebbian plasticity exist, suggesting that the timescale distinction is little more than a convention [16][17][18][19][20][21][22][23][24] (see [25] for a review).Hebbian plasticity is more powerful than the presynaptic variety due to its ability to create attractor states of neuronal dynamics, the basis of Hopfield networks [26].We are therefore motivated to introduce ongoing Hebbian plasticity in a recurrent network, without necessarily imposing a separation of timescales between neuronal and synaptic dynamics.This has unexpected, computationally useful consequences, a key example being freezable chaos, a phase in which a stable fixed point of neuronal dynamics is destabilized through synaptic dynamics.By contrast, introducing presynaptic plasticity to this model simply adds an effective constant input to each neuron (Appendix A).Our work thus provides a theoretical impetus for further experimental investigation of ongoing Hebbian plasticity mechanisms.
The view that synapses serve solely to sculpt neuronal dynamics is mirrored in machine learning.In artificial neural networks, weights are trained via gradient descent, then fixed.However, allowing weights to be modulated by the activity of the units has been shown to confer computational advantages [27][28][29][30][31], particularly in tasks requiring short-term memory storage [32].For example, Ba et al. [33] showed that recurrent networks benefit from a combination of "slow weights" trained via backpropagation and "fast weights" that undergo activitydependent updates (the model we study is essentially the continuous-time counterpart of this proposal).In practice, recurrent networks have been superseded by transformers [34].While these models were not neuroscientifically motivated, Schlag et al. [35] showed that linearized transformers [36] are equivalent to fast weight programmers [37], a recurrent network with activity-dependent weight updates [38].
A major impediment to studying neuronal-synaptic dynamics, in both neuroscience and machine learning, is that the analytical methods developed for nonplastic networks often do not translate to plastic networks, particularly when the neuronal and synaptic timescales are not well separated.For example, a common simplification is to study nonplastic networks with linear neuronal dynamics; however, such networks become highly nonlinear when synaptic plasticity is introduced [39].Moreover, synaptic degrees of freedom increase the dimension of phase space from O(N ) to O(N 2 ).In studying and training nonplastic networks, the theory of random networks has played a crucial role [40].In seminal work, Sompolinsky et al. [41] showed that random recurrent networks exhibit a phase transition to highdimensional chaotic activity at a critical coupling variance [42].While this analysis was in firing-rate (nonspiking) networks with fully unstructured couplings, key phenomena, such as the transition to chaos, generalize to spiking networks with more realistic distributions over couplings [43].Here, we extend this approach to plastic networks-that is, develop a theory for such a model in the thermodynamic limit, compute its phase diagram, and characterize its dynamics-thereby providing a foundation for understanding how coupled neuronal-synaptic dynamics could underlie computation.

II. MODEL
We augment the random-network model of Sompolinsky et al. [41] with dynamic couplings.There are N neuronal units with pre-activations x i (t) and activations ϕ i (t) = ϕ(x i (t)), where ϕ(•) is a nonlinearity.Throughout, we use ϕ(•) = tanh(•).Neurons interact through allto-all time-dependent couplings W ij (t) according to the neuronal dynamics Concurrently, W ij (t) displays synaptic dynamics.We first express these couplings as a sum of quenched and fluctuating terms, where J ij ∼ N 0, g 2 /N provides quenched disorder.The fluctuating term A ij (t) follows a local plasticity rule, where k is the sign and strength of the plasticity, which is Hebbian or anti-Hebbian for k > 0 or k < 0, respectively; and p is the synaptic decay timescale in units in which the neuronal decay timescale is unity.We do not require p ≫ 1, though a reasonable constraint from biology is p > 1 since the synaptic timescale is unlikely to be shorter than the neuronal timescale.The couplings W (t) include self-connections (on-diagonals) with the same dynamics as non-self-connections (off-diagonals).However, the effect of these self-connections on each neuron is ∼ 1/ √ N , and thus is negligible as N → ∞.The full set of dynamic variables comprises the N neurons, x i (t), and N 2 fluctuating couplings, A ij (t).We study their collective dynamics as a function of g, k, and p as N → ∞ [Fig.1(a)].
Because A(t) is an average over outer products for each order-one timestep and decays with timescale p, it has rank of order p or smaller.Given that p is order-one (not order-N ), the approximate rank of A(t) is intensive.Note that J ij ∼ 1/ √ N and A ij (t) ∼ 1/N , so plasticity is vanishingly weak at the single-synapse level as N → ∞.Nevertheless, because A(t) is approximately intensiverank, the random and fluctuating couplings both have order-one macroscopic effects.This is because when neuronal activity exhibits alignment to the states encoded in A(t), the J and A(t) terms each contribute an order-one input to neurons, as made clear by the mean-field analysis in Sec.IV.Further intuition for this scaling can be obtained from the spectra of J and A(t); the chosen scaling implies that the eigenvalues of both matrices are order-one, allowing the modes they drive in the network to compete on equal footing.This scaling of low-rank structure relative to randomness at individual synapses is generic in models in which the couplings are the sum of random and low-rank terms [44][45][46] (but see [47] who used a ∼1/ √ N rank-one term given by an outer product of orthogonal vectors).The scaling also appears in spiked matrix or tensor models in statistical physics [48].In experiments that measure changes in synaptic strengths, plasticity on this weak scale could go unnoticed, but nevertheless exert dramatic influence at the network level.

III. PHASE DIAGRAM
We summarize the behavior of the model with a phase diagram in (g, k) parameter space for p = 2.5 [Fig.1(b)], noting that constant-p slices of the full (g, k, p) diagram are similar for order-one values of p.We refer to this phase diagram for the rest of this section.
For k = 0, the model reduces to that of Sompolinsky et al. [41] (dashed horizontal line).For g > 1, this nonplastic network is chaotic (i).For g < 1, the trivial neuronal fixed point of the nonplastic network, x(t) = 0 N , is globally stable and the network is quiescent.In analogy with the nonplastic network, the plastic network can produce chaotic activity for g > 1, and the activity is further shaped by synaptic plasticity.Hebbian plasticity, k > 0, slows activity (ii), while anti-Hebbian plasticity, k < 0, quickens activity and generates an oscillatory component (iii).
For g < 1, the trivial neuronal-synaptic fixed point, (x(t), A(t)) = (0 N , 0 N ×N ), is stable.However, in contrast to the nonplastic network, this fixed point is not necessarily globally stable, but coexists with chaotic states for large k (iv).Thus, Hebbian plasticity can induce dynamic activity in an otherwise quiescent network.For g < 1, if k is not large enough to induce dynamic activity, the network is globally quiescent (v).In dynamic states, network activity is shaped by a proliferation of stable fixed points throughout phase space.In particular, if Hebbian plasticity is strong (hatched region), there exist stable fixed points to which finite-size networks settle following transient chaos, accompanied by the rank of A(t) collapsing to unity (vi).Such fixed points are stable with respect to the combined neuronalsynaptic dynamics.Their number, which we compute, is exponential in N .
Two scenarios can lead to transient chaotic states.First, when g < 1 and k is large enough to induce chaotic activity, finite-size networks may collapse to the trivial fixed point.Additionally, when (g, k) is in the hatched region, finite-size networks may collapse to nonzero fixed points.
Finally, we describe freezable chaos.Consider a chaotic state with Hebbian plasticity.Suppose we abruptly halt synaptic dynamics.If k is small, the halted-synapse neuronal dynamics are chaotic, with no trace of the halttime neuronal state in the activity (ii).If k is larger, the halted-synapse neuronal dynamics are chaotic, but neurons fluctuate around the halt-time state, so a memory of this state is retained (vii).If k is sufficiently large, the halted-synapse neuronal dynamics are nonchaotic; neurons flow to a stable fixed point near the halt-time state (viii).In all cases, releasing the synapses returns the network to neuronal-synaptic chaos.As N → ∞, these three cases, which we label nonfreezable, semifreezable, and freezable chaos, respectively, are distinct phases.In freezable chaos, there is, at any instant, a stable fixed point of neuronal dynamics that is destabilized by synaptic dynamics.Thus, a stable fixed point can be created near any neuronal state by halting synaptic plasticity, enabling a new form of short-term memory storage.

IV. DYNAMICAL MEAN-FIELD THEORY
The temporal structure of network activity is described in the limit N → ∞ by a dynamical mean-field theory (DMFT) whose main order parameter is the single-unit autocovariance (two-point) function, where we assume statistical stationarity in time.Integrating the synaptic dynamics, Eq. ( 3), and inserting this into the neuronal dynamics, Eq. (1), gives where we have separated terms arising from J (first term on the rhs) and A(t) (second term on the rhs).Taking the limit N → ∞ yields the single-site picture where η(t) is an effective Gaussian field with zero mean and second-order statistics The DMFT is closed by the self-consistency condition In the single-site dynamics of Eq. ( 6), synaptic plasticity introduces a convolutional self-coupling with a kernel that depends self-consistently on C(τ ).For k = 0, the self-coupling vanishes and the DMFT reduces to that of [41], which can be solved analytically.For k ̸ = 0, the nonlinearity of the self-coupling induces a non-Gaussian distribution over x(t)-in particular, a distribution that becomes increasingly bimodal with larger k due to saturation of the tanh function-so we solve the DMFT equations using standard numerical techniques [49][50][51][52][53][54] (Appendix C 1).The DMFT agrees closely with simulations [Fig. 2

(a)].
A. Chaotic states with g > 1 We now examine the solutions C(τ ) of the DMFT.In the regime where the nonplastic network is chaotic, g > 1, Hebbian plasticity slows the activity, broadening C(τ ) (Fig. 2(a), k > 0 solutions).A network with static, symmetric couplings (e.g., the Hopfield network) admits a Lyapunov function that guarantees convergence to fixed points [26].The slow activity generated by Hebbian plasticity results from competition between J , a random, asymmetric matrix that promotes dissipative, chaotic dynamics; and A(t), a symmetric matrix that promotes convergence to a drifting fixed point that trails the neuronal state.
This competition has an interesting dependence on the model parameters.For small k and large p, neurons fluctuate rapidly relative to the synaptic decay timescale and the effect of plasticity is averaged out.As k is increased, synapses slow neurons by attracting the neuronal state toward its history, permitting a stronger synaptic signal to be encoded.This stronger signal causes further neuronal slowing.The interaction of neurons and synapses in this positive-feedback process causes the timescale of neuronal fluctuations to diverge as k increases.For large p and k, finite-size networks can show bistability between a fast state in which plasticity is averaged out and a slow state in which synapses drag neurons (Appendix B).
This slowing behavior can also be understood through the DMFT.We quantify the speed of neuronal fluctuations by defining the dynamic timescale whose dependence on g, k, and p is illustrated in Fig. 2(b).The size of the integral term in the singlesite dynamics [Eq.( 6)] is ∼kτ * /p for τ * ≪ p (and ∼k for τ * ≫ p).Increasing this term slows x(t), increasing τ * .That τ * depends on the size of this term, which itself depends on τ * , produces a positive-feedback loop.Once k is large enough so that the integral term competes with η(t), which occurs when kτ * /p ∼ g, τ * grows rapidly.The inflection of τ * in k is sharpest for large p or g, in which case kτ * /p and g are well separated at small k (Fig. 2(b), p = 20).Under anti-Hebbian plasticity, rather than synapses attracting the neuronal state to its history, this effect is repulsive.This quickens the dynamics and adds an oscillatory component to the activity, tightening C(τ ) and creating oscillations during its decay to zero (Fig. 2(a), k < 0 solutions).In the single-site picture [Eq.( 6)], plasticity modifies the dynamics of x(t) by introducing time-delayed negative feedback, which generically induces oscillations [55,56].While finite-size simulations of the model of [41] show limit cycles, our calculation of C(τ ) for N → ∞ demonstrates that this anti-Hebbian oscillatory component is not merely a finite-size effect.These oscillations are further characterized by the (normalized) power spectrum S(f , where Ĉ(f ) denotes the Fourier transform of C(τ ) [Fig.2(c)].Rather than containing a peak at a nonzero frequency, the power spectra corresponding to the anti-Hebbian autocovariance functions in Fig. 2(a) possess a range of large frequencies (though S(f ) develops a nonzero peak for smaller values of g; not shown).Point estimates of the dominant oscillatory frequency computed from the first three zero-crossings of C(τ ) roughly capture the frequency scale at which S(f ) decays (Fig. 2(c), triangular markers).

B. Chaotic states for g < 1
We next consider the regime g < 1 in which the nonplastic network is globally quiescent.The plastic network has a trivial fixed point, (x(t), A(t)) = (0 N , 0 N ×N ), that is stable when all eigenvalues of J have real part less than unity.For N → ∞, Girko's circular law implies that this requires g < 1 [57].In contrast to the nonplastic network, the trivial neuronal-synaptic fixed point coexists with dynamic states for large k as indicated by DMFT solutions (Fig. 3(a), solid lines).In Sec.V B, we confirm over a restricted parameter regime that these solutions are chaotic (i.e., have positive maximum Lyapunov exponent).
Dynamic states for g < 1 eventually collapse to the trivial fixed point or, if k is large enough, to a stable nonzero fixed point (Sections V B, VI).In this section, we consider values of g < 1 and k leading to chaotic states that are prone to collapsing to the trivial fixed point.Realizing these states in simulations is nontrivial because random initial conditions typically miss the dynamic attractor and decay to zero.A workaround is to deform a chaotic state with g > 1, for which the trivial fixed point is unstable, to a chaotic state for g < 1 by reducing g in the spirit of annealing.Using this method, we verified that simulations agree with the DMFT solutions (Fig. 3(a), dashed lines).
After realizing plasticity-induced chaotic states via this procedure, we observe transient activity with a lifetime that is approximately log-normally distributed (Fig. 3(b), right).The median log-lifetime before collapsing is linear in N over several decades, consistent with the typical lifetime diverging exponentially and becoming infinite in the limit N → ∞ in which the DMFT applies (Fig. 3(b), left).This divergence is faster for larger k.
C. First-order transition to nontrivial DMFT solutions for g < 1 For g < 1, if k is not large enough to induce dynamic activity, the plastic network is globally quiescent.We now analyze the boundary between these phases.For order-one values of p, we find τ * ≫ p for dynamic states for g < 1, which reduces the single-site dynamics to the "slow" form This single-site picture is related to that of a nonplastic network with order-one self-coupling parameter s for which the single-site dynamics read as studied by Stern et al. [53].This network has a continuous transition from quiescence to dynamic activity at g + s = 1.We map solutions of the Stern network onto solutions of the plastic network by enforcing y(s) = s/C(0) = k.For a given g < 1, C(0) becomes nonzero at s = 1 − g, grows with s, and saturates at unity; thus, y(s) descends from infinity at s = 1 − g and grows as s for large s, tracing a U shape [Fig.3(c)].Each k draws a horizontal line intersecting y(s) at selfconsistent solutions of Eq. ( 10).The critical k, defining the boundary between phase-diagram regions iv and v, is the minimum of y(s) (Fig. 3(c), lower dashed line).
As C(0) is finite here, this is a discontinuous, first-order transition.Physically, this is because sustaining dynamic activity for g < 1 requires finite self-coupling.Since the self-coupling depends on C(0), time-dependent solutions with arbitrarily small C(0) are not possible.As g → 0, the dynamic timescale diverges, leaving only fixed points.We show in Sec.VI that fixed points exist for k > 2.02 at g = 0.
Due to the U shape of y(s) for g < 1, y(s) = k has two solutions for k larger than its critical value: one with large s and C(0), the other with small s and C(0) (Fig. 3(c), upper dashed line).As g → 1 − , the small-s solution vanishes (Sec.IV D) while the large-s solution remains finite.For g > 1, y(s) = k has a unique solution whose deformation to g < 1 gives the large-s solution (Fig. 3(c), g = 1.1 curve).

D. Second-order transition at g = 1
Finally, we solve the DMFT for g → 1 + , k < 1.This limit marks a continuous transition from dynamic activity to global quiescence in which C(0) vanishes and τ * diverges.Due to the vanishing variance, the self-coupling in Eq. ( 10) can be linearized, This is equivalent to the single-site dynamics of the nonplastic network with g eff = g/[1 − kC(0)] and time con- to leading order in ϵ.In contrast to the behavior away from this transition, activity becomes faster with increasing k.C(0) diverges as k → 1 − , indicating that solutions with g = 1 for k > 1 have finite variance.Validity of the solution Eq. ( 13) requires k < 1, lest we obtain a nonphysical negative variance.On the other hand, as g → 1 − for k > 1, a positive variance is obtained; this is the small-s solution for g < 1 described in Sec.IV C.

V. HIGH-DIMENSIONAL ANALYSIS
The DMFT describes the temporal structure of network activity through an effective single-site picture.Importantly, the network dynamics result from a complex interaction of high-dimensional neuronal-synaptic modes.We now probe the high-dimensional origin of the dynamics, first through an analytical study of the spectrum of the Jacobian describing the local, linear dynamics, and then through a numerical study of the Lyapunov spectrum describing the global, nonlinear dynamics.Both the Jacobian and Lyapunov spectra show a topological transition at large k to a form with a slow, synapse-dominated band and a fast, neuron-dominated band, with the former driving network network activity.This suggests a flipped view of the network dynamics as being driven by the synaptic couplings, with neurons serving as the connections.

A. Jacobian spectrum
Let us represent Eqs.(1-3) defining the model in the form where a(t) = vec A(t) contains all S = N 2 elements of A(t).We use the notation [∂x/∂y] ij = ∂x i /∂y j for vectors x, y.The Jacobian is a (N + S)-dimensional block matrix, At any instant, the neurons provide N -dimensional input to the S synaptic variables, inducing a low-dimensional structure of the Jacobian.Rather than consider M directly, it is convenient to study the linear dynamics that it generates, The input that δa receives from δx is confined to a Ndimensional subspace spanned by the columns of ∂G/∂x.Perturbations to δa in the (S − N )-dimensional orthogonal complement subspace relax with timescale p. Due to these relaxational modes, M has eigenvalue −1/p with multiplicity S − N .The remaining 2N eigenvalues result from the interaction of δx and the component of δa in the N -dimensional subspace that receives input from δx.Projecting δa into this subspace via δ ã = (∂G/∂x) + δa, where (•) + denotes the pseudoin-verse, the 2N -dimensional dynamics of interest are The associated 2N -dimensional dynamics matrix is In summary, the eigenvalues of M are −1/p with multiplicity S − N together with the 2N eigenvalues of M .
We refer to M as the reduced Jacobian [58].
Each eigenvector of M can be written in terms of its x and a components, v = (v x , v a ).From the lower block row of M , we obtain the relation Letting f a (λ) = ∥v a ∥ 2 /∥v∥ 2 denote the relative weight on the synaptic component, this relation implies which falls off radially with distance from −1/p.Modes become synapse-dominated as λ → −1/p, giving way to a delta function of S − N purely synaptic modes at this point.
The Jacobian analysis thus far holds for any F (•, •), G(•), and S > N .In particular, the simplification to the structure of the Jacobian does not depend on the specific plasticity rule posited in Eq. (3).We now substitute the forms of F (•, •) and G(•) from Eqs. (1-3), yielding M = Mbulk + Mlow-rank , where Here, (x, A) is a point in phase space, ϕ = ϕ(x), ϕ ′ = ϕ ′ (x), and C(0) = ∥ϕ∥ 2 /N .Mbulk generates the bulk of the spectrum of M while Mlow-rank can, in principle, contribute an intensive number of outlier eigenvalues.However, both in chaotic states and at fixed points, we find that the reduced Jacobian does not have outliers [59].We therefore focus on the spectrum of Mbulk .We compute the boundary curve encompassing the compact spectrum of Mbulk for N → ∞ using a theorem of Ahmadian et al.
[60] concerning random matrices expressible as a linear reparameterization of an elementwise independent and identically distributed random matrix (Appendix E).This shows that the limiting spectral density of Mbulk has support at λ ∈ C if where ⟨•⟩ x is an average over the components of x.
We use this result to probe the high-dimensional origin of the dynamics in the plastic network, noting that the spectral density of the Jacobian computed at any point on a dynamic trajectory is time-independent and self-averaging as N → ∞.We use the DMFT to obtain a Monte-Carlo estimate of ⟨•⟩ x .The predicted boundary tightly hugs the Jacobian spectra evaluated from simulations in dynamic states (Fig. 4).These simulation spectra do not contain outliers, justifying our focus on Mbulk .We place a dot at −1/p to indicate the delta function of N 2 − N synaptic relaxational modes present in the spectrum of M .Modes past the stability line, Re(λ) = 0, locally destabilize the network and drive the dynamics.
Setting k = 0 recovers the circularly symmetric spectrum of the nonplastic network, enclosed by |1 + λ| ≤ g ⟨ϕ ′ (x) 2 ⟩ x as predicted by Eq. ( 21) (Fig. 4, k = 0 row).As k is increased, the delta function at −1/p repels eigenvalues leftward, creating a hole (Fig. 4, g = 3, k = 1 spectrum).Further increasing k produces a topological transition to a spectrum with two bands and no holes (Fig. 4, all other spectra with k > 0).Unstable modes are increasingly focused along the real axis, leading to slow activity as seen via the DMFT.The slow, destabilizing band is dominantly synaptic, and the fast, relaxational band is dominantly neuronal.The slowestrelaxing of the fast modes have real part close to −1, the neuronal decay timescale.This two-band topology therefore reflects a dynamic state in which slow network activity is synapse-driven rather than neuron-driven.
As k is decreased from zero into the anti-Hebbian regime, the delta function repels eigenvalues rightward, creating a hole to its right (Fig. 4, k = −1 row).Further decreasing k produces a topological transition to a form with a single band and no holes (Fig. 4, k = −2 row).In contrast to the Hebbian spectra, there are pronounced lobes of dominantly imaginary unstable modes, generating fast, oscillatory activity also seen via the DMFT.The oscillatory frequency computed from zero-crossings of C(τ ), as in Fig. 2(c), corresponds roughly to the locations of the imaginary lobes of the spectra, particularly for smaller values of g (Fig. 4, triangular markers for k < 0 spectra).

B. Lyapunov spectrum
While analytically accessible, the Jacobian spectrum describes only the locally linear dynamics.Rigorously characterizing the nonlinear dynamics requires a calculation of the spectrum of Lyapunov exponents, which are defined as follows.Suppose a ball of radius ϵ is tossed into the flow.As the dynamics unfold, the ball expands and contracts into an ellipsoid.The Lyapunov exponents are the exponential growth or decay rates of the principal axes of this ellipsoid as t → ∞ (simultaneously, ϵ → 0 is taken so that the ellipsoid stays small).A positive Lyapunov exponent implies exponential sensitivity to initial conditions, indicating chaos.The Lyapunov spectrum cannot be derived from the Jacobian spectrum due to both the non-normality of the Jacobian and the timedependence of its eigenvectors.We first study the maximum Lyapunov exponent λ max , which dominates trajectory divergence as t → ∞.We measured λ max in simulations by injecting a small perturbation to the system and measuring the slope (ver-sus t) of the log of the norm of the difference between the perturbed and unperturbed trajectories.For g < 0, we realized dynamic network states using the deformation method described in Sec.IV B. For each setting of (g, k), we simulated 200 random networks.The output of this analysis is displayed as a heatmap in (g, k) parameter space in Fig. 5(a).Parameter values that resulted in convergence to nonzero fixed points (Section VI) within the simulation time for at least 80% of networks are hatched; values that resulted in quiescence of all networks are white.
In regions of parameter space that do not converge to fixed points over the simulation time, λ max is positive, indicating chaos.This includes part of the region g < 1, confirming that plasticity can induce chaos in an otherwise quiescent network.This analysis provides a simulation-based confirmation of the boundary marking the first-order transition to nontrivial DMFT solutions for g < 1 derived in Sec.IV D [Fig.5(a), gray lines].
As k is increased and/or g is decreased, we observe a smaller and smaller Lyapunov exponent that eventually results in simulations reliably collapsing to nonzero fixed points (Fig. 5(a), solid-to-hatched boundary).This crossover occurs in a parameter regime for which phase space is densely filled with stable fixed points (Sec.VI).Using the present finite-N analysis, we are unable to determine whether λ max in the hatched region in Fig. 5(a) is small and positive, or negative.Additionally, solving the DMFT in the hatched region is numerically difficult due to the diverging dynamic timescale (Sec.IV A).As N is increased over a decade, the boundary marking this crossover shifts slightly upward (Fig. 9).
We next study the full Lyapunov spectrum for chaotic states.In general, Lyapunov spectra can be computed numerically by propagating a set of vectors in the tangent space of the flow and periodically orthonormalizing them to prevent their explosion/vanishing and to extract their growth/decay rates, as explained in detail in prior works [49,61,62].The plastic network has N + N 2 variables, so propagating a complete basis is prohibitively computationally expensive for large N .Fortunately, O(N 2 ) exponents concentrate at −1/p, so it suffices to compute the O(N ) largest and smallest exponents.We compute the largest exponents using the aforementioned procedure with an undercomplete set of O(N ) tangent-space vectors.We find the smallest exponents by doing the same for the time-reversed dynamics, noting that the smallest exponents are the largest of the time-reversed system.A complication is that the time-reversed dynamics are unstable.We therefore run time-reversed tangentspace dynamics, tamed by orthonormalization, atop a time-reversed trajectory produced by the forward-time dynamics.We verified that this method accurately computes the largest and smallest Lyapunov exponents of the nonplastic network [41,62].
Histograms of the Lyapunov spectra are shown in Fig. 5(b).For k = 0, the spectrum is the same as that of a nonplastic network with a spike at −1/p [63].We verified that the measurement of λ max obtained using this method matches the perturbation measurement displayed in the heatmap [Fig.5(c)].
For large k, the Lyapunov spectra recapitulate the topological transition to two bands of the Jacobian spectra, further demonstrating a synapse-driven dynamic state [Fig.5(b)].In analogy with the Jacobian spectra, the slow, destabilizing Lyapunov band spans −1/p to λ max , and the fast, relaxational band has an upper limit near −1.
Our calculation of the Lyapunov spectrum allows for further calculation of diffeomorphic-invariant properties of the strange attractor [62].We focus on its dimension, shown in Fig. 5(d), defined by the Kaplan-Yorke conjecture as the number of exponents that must be summed, ranked in descending order, to achieve a cumulative sum of zero [64].We display an intensive version of this quantity that has been divided by N [42,62].Both λ max and the attractor dimension vary nonmonotonically in k [Fig.5(c,d)].This contrasts with the dynamic timescale τ * , which increases monotonically [Fig.2(b)].The eventual decrease of the attractor dimension is reminiscent of the behavior of the participation ratio-based dimension of activity in networks with partially symmetric connectivity, which was recently computed analytically [42].The decline of λ max and the attractor dimension at large k likely reflects the proliferation of stable fixed points throughout phase space, the subject of the next section.

VI. FIXED POINTS
For large k, the dynamics of the plastic network are shaped by a proliferation of stable fixed points, and finitesize networks settle to these fixed points after transient chaos (Fig. 1(b), hatched region and vi).We first probe this settling process by analyzing how A(t) collapses to a rank-one state, measuring the approximate rank of A(t) as the participation ratio of its spectrum, {λ i (t)}, Note that if A(t) encodes P decorrelated neuronal states with equal magnitude (i.e., λ 1 (t), . . ., λ P (t) = const.and λ P +1 (t), . . ., λ N (t) = 0), then PR A(t) = P .Evaluating Eq. ( 22) in the limit N → ∞ gives Thus, PR A(t) is intensive and, as N → ∞, timeindependent.If τ * ≪ p, then T ≈ τ * , so PR A(t) is the ratio of these timescales.On the other hand, if τ * ≫ p, then T is slightly less than p, so PR A(t) is slightly larger than unity.Hebbian plasticity tends to slow chaos.As k is increased, PR A(t) therefore drops closer to unity, with tem- poral fluctuations about the mean-field value in finitesize networks [Fig.6(b)].During chaos, neurons continuously escape the synaptic drag.However, in finite-size networks with sufficiently large k, synapses can "win," namely, the network fluctuates into fixed point of the form (x(t), A(t)) = x, ϕ(x)ϕ(x) T with PR A(t) = 1 [Fig.6(b)].These fixed points are stable with respect to the combined neuronal-synaptic dynamics.
We now compute the number of stable fixed points in phase space.As k → ∞, there is a stable fixed point associated with each of the 2 N binary neuronal states.Thus, any initial condition is pulled to a fixed point with high overlap with the initial neuronal state, consistent with the diverging dynamic timescale τ * in this regime (Sec.IV A).For finite k, we expect exponentially many stable fixed points.This is in contrast to the non-plastic network [41], for which there exist exponentially many fixed points for g > 1, but they are all unstable [65,66].Due to the exponential scaling, we will study the log-number of stable fixed points per neuron, an intensive, self-averaging quantity.Fluctuations around fixed points are governed by the Jacobian spectrum analyzed in Sec.V. Our derivation follows that of Stern et al. [53], differing in its initial steps that enforce stability using this neuronal-synaptic spectrum.
At a fixed point, the DMFT equations reduce to where η ∼ N (0, g 2 C(0)) and C(0) = ϕ 2 (x) η .Eq. ( 24) can be written g(x) = x − kC(0)ϕ(x) = η.The Gaussian distribution over η induces a non-Gaussian distribution over x due to the nonlinearity of g(x).If |η| < η m for some η m , there are three solutions to g(x) = η.The two outer solutions are at points of positive slope of g(x); the central solution has negative slope.The negative-slope solution renders a fixed point unstable for the following reason.The denominator of the averaged quantity in the Jacobian boundary [Eq.( 21)] is Since λ + pλ(1 + λ) varies from zero to infinity as λ varies from zero to infinity, g ′ (x) < 0 causes the denominator to vanish at for some λ > 0, precluding stability.Two solutions remain for |η| < η m .Stern et al. [53] observed that typical fixed points maximize the combinatorial number of solutions subject to stability and, moreover, are marginally stable, meaning that the spectral boundary sits at λ = 0. Setting λ = 0 in Eq. ( 21), marginal stability requires Since p appears in neither Eq. ( 24) nor Eq. ( 26), it drops out of the calculation, indicating that the number of stable fixed points is independent of this timescale [note, however, that the shape of the Jacobian spectrum is pdependent; Fig. 6(e)].Following Sec.IV, we map solutions of the nonplastic network of [53] onto solutions of the plastic network by enforcing s = kC(0).
The structure the solutions is as follows.For each g, there is an onset of fixed points at a critical k [Fig.6(c)].The log-number grows monotonically with k and saturates at log 2 2 N /N = 1.The onset of fixed points is discontinuous for g < 0.76, in which case there are two solutions of the fixed-point mean-field theory for a given (g, k): one with large s and C(0), the other with small s and C(0).The former are exponentially dominant.In Fig. 6(d), we display the log-number of fixed points per neuron as a heatmap in (g, k) parameter space.
We now study the spectra of fluctuations around fixed points.We used the analytical form of the distribution over x from the mean-field analysis to evaluate ⟨•⟩ x in Eq.21, yielding a prediction for the Jacobian spectrum at fixed-points.We also allowed many simulations of chaotic networks to settle to fixed points, then measured their Jacobian spectra.We find that the simulation fixed points are indeed marginal, with spectra described accurately by the theoretical prediction [Fig.6(e)].
The stability of fixed points with this critical value of k is checked by direct evaluation of the eigenvalues of the reduced Jacobian, which are λ = −(1 + p ± (1 + p) 2 − 8p/3)/2p < 0 with multiplicity N − 1 for each sign, comprising the bulk; λ = −1/p − 2/3 < 0 with multiplicity one; and λ = 0 with multiplicity one, indicating marginal stability.Further increasing k gives proper stability.Here, the marginal eigenvalue is an outlier, while the bulk has finite negative real part.This suggests that the mean-field fixed-point calculation, which assumed marginal stability of the bulk, breaks down for small g.Simulations and theory agree down to g = 0.5, k = 1.9, so such a breakdown would have to occur for smaller values of these parameters.

VII. FREEZABLE CHAOS
The previous section characterized stable fixed points of the dynamics, as is typical dynamical-systems studies.Another question is whether a subsystem can have a stable fixed point that is unstable in the context of full system.We now study freezable chaos, a state where a stable fixed point of neuronal dynamics is continuously destabilized through synaptic dynamics, generating chaos.
As described in Sec.III, for chaotic states with Hebbian plasticity, we define nonfreezable, semifreezable, and freezable chaos depending on the neuronal dynamics that ensue after halting synaptic dynamics.In (semi-)freezable chaos, neurons retain a stable memory of the halt-time state as we illustrate in Fig. 7(a).
Picking the halt time to be t = 0, the couplings remain at W (0) = J +A(0).Networks with such "random-pluslow-rank" couplings have been studied in the context of nonplastic networks by Mastrogiuseppe and Ostojic [44], who found chaotic, structured-chaotic, and structuredhomogeneous phases that are qualitatively similar to nonfreezable, semifreezable, and freezable chaos in plastic networks.Crucially, our analysis in plastic networks is complicated by A(0) arising due to synaptic plasticity dependent upon neuronal activity driven by J .More specifically, the neuronal states comprising A(0) are largely confined to a subspace spanned by dominant eigenvectors   2), and freezable (bottom, k = 1.55) chaos for g = 2.25 and p = 2.5.At the first vertical line, A(t) is halted.In (semi-)freezable chaos, a memory of the halt-time neuronal state persists.Stability is demonstrated by perturbing the neurons (second line) and clamping them near zero (third line).In (semi-)freezable chaos, neurons relax back to the stored pattern after these manipulations.At the fourth line, A(t) is released.This process is repeated at the fifth line.(b) Demonstration of the alignment between J and A(t).At the first vertical line, A(t) is halted.At the second line, J , but not the halted A(t), is shuffled, destroying the neuronal fixed point.At the third line, A(t) is released.This process is repeated at the fourth line.g = 2.5, k = 2.55, and p = 2.5 [as in the freezable-chaos plot in (a)].(c) Time-dependent overlap Q(t) from the two-replica DMFT (solid lines) and in simulations (dashed lines) for g = 2, p = 2.5, and various values of k.(d) Halt-time overlap Q(0) from the two-replica DMFT (lines) and in simulations (triangular markers) as a function of k for p = 2.5 and various values of g.Markers are colored according to whether simulations show semifreezable or freezable activity.(e) Same as (d), but using the zero-time correlation coefficient ρ(0).In (c-e), We show only the positive solutions for Q(t) and ρ(t), which are realized upon halting synapses; there is also a symmetric negative solution.
of J [42], inducing alignment between J and A(0).We demonstrate the importance of this alignment by running a simulation in which we halt A(t), storing the halt-time neuronal state as a stable fixed point [Fig.7(b)].We then shuffle J while keeping the halted A(t) fixed.This preserves the statistics of J , but destroys correlations between J and the halted A(t).If these correlations were negligible, the neuronal fixed point would reorganize to a different fixed point.Instead, shuffling destroys the fixed point, causing the network to switch to neuronal chaos.Upon releasing A(t), the network returns to neuronalsynaptic chaos, and A(t) adapts to the permuted version of J .At a later time, we halt synaptic dynamics again.As A(t) has adapted to the shuffled J , a stable neuronal fixed point is created.
We handle this alignment through a replica mean-field analysis involving two networks, A and B, with neuronal states ϕ A (t) = ϕ(x A (t)) and ϕ B (t) = ϕ(x B (t)).Network A has neuronal-synaptic dynamics for all t with quenched random couplings J .Network B has neuronal dynamics for all t with halt-time couplings W (0) constructed from the same J and neuronal states ϕ A (t) as network A, We define order parameters to characterize the phases of interest, starting with the overlap Evaluated at t = 0, this parameter indicates how accurately the halt-time neuronal state is retained as a memory.We also define the autocovariance function of network B, where we assume statistical stationary in time.The autocovariance function of network A is C(τ ), derived previously.Finally, assessing stability of the neuronal fixed point, when one exists, requires the stability matrix , where x B is a fixed point in network B. Stability requires that the spectrum has negative real part.Instability is dominated by the circularly symmetric bulk, so it suffices to compute its radius r, which follows from random matrix theory [60].
Having established these order parameters, we can define the phases of interest quantitatively.Chaos is nonfreezable when Q(t) = 0 is the only solution.Chaos is semifreezable when there is a nontrivial solution for Q(t) associated with a solution for D(τ ) that decays in τ (to a nonzero value).In this case, there may or may not be a fixed-point solution in which there is a distinct nontrivial solution for Q(t) associated with D(τ ) = const.If there is a fixed point, it is unstable, r > 1.Finally, chaos is freezable when only a fixed-point solution for Q(t) and D(τ ) exists.In this case, it is stable, r < 1.Because the fixed point of the halted-synapse system is stable, the lifetime of the memory is infinite.Due to the x → −x symmetry of the network, if Q(t) is a solution, so is −Q(t).Upon halting synapses, the positive solution is realized, barring "flips" that occur in finite-size networks near the onset of semifreezable chaos.
We derive a two-replica DMFT that permits calculation of these order parameters.The high-dimensional equations describing the two replicas are Sending N → ∞ and taking the limit where the time coordinate of network B is much greater than zero, we obtain the single-site picture where η A (t) and η B (t) are Gaussian fields with zero mean and second-order statistics The off-diagonal covariances effectively encode the alignment between J and A(0).These off-diagonals do not depend on the time coordinates of network B in accordance with the limits taken in Eqs.(28) and (29).The system is closed by the self-consistency conditions We solve the DMFT numerically; imposing ∂ t x B (t) = 0 and D(τ ) = const.gives the fixed-point solution (Appendix C 2).We find excellent agreement between theory and simulations [Fig.7(c-e)].
We now examine the solutions of the two-replica DMFT.In the freezable-chaos regime, Q(t) peaks at t > 0, indicating that the fixed point is more aligned with neuronal states that would have unfolded after the halt time than with the halt-time state itself [Fig.7(c)].This reflects a tendency of neurons to continue with "momentum" before becoming trapped in a fixed point [Fig.7(a)].
Increasing k from zero yields a sequence of continuous phase transitions that we analyze by plotting Q(0) and r against k [Fig.7(d)].For small k, Q(0) = 0, indicating nonfreezable chaos.As k is increased, Q(0) develops a nonzero solution associated with a decaying D(τ ), marking the onset of semifreezable chaos.As k is increased further, Q(0) develops an additional nonzero solution associated with D(τ ) = const., marking the onset of an unstable fixed point.Instability is signaled by r > 0, with r computed under the fixed-point solution.Continuing to increase k causes r to decrease and drop below unity, marking the onset freezable chaos, at which point the dynamic and fixed-point solutions converge.The convergence of the dynamic and fixed-point solutions at r = 1 is expected on physical grounds and can also be derived through the DMFT: when r is computed under the dynamic solution, the decay timescale of D(τ ) diverges as ∼1/ √ r − 1 as r → 1 + , giving D(τ ) = const.at r = 1 (Appendix D).The quality of memory retention can be measured by the correlation coefficient which varies between zero and unity [Fig.7(e)].

VIII. DISCUSSION
We characterized the dynamics of N neurons coupled to N 2 dynamic synapses.Strong Hebbian plasticity causes the timescales of the system, measured through the Jacobian or Lyapunov spectra, to segregate into a slow, synapse-dominated band and and a fast, neurondominated band.The synapse-dominated band drives the dynamics.It is possible that this two-band structure could be detected through in-vivo recordings of neuronal activity.Takens' embedding theorem implies that it is possible, in principle, to extract the spectrum of neuronal-synaptic timescales from neuronal activity alone [67].This segregation of timescales could also be examined during task execution.If the dynamics are synapse-driven, neurons may revert to their trialaverage trajectories upon optogenetic or electrophysiological perturbation.Prior studies have attributed such robustness to neuronal mechanisms such as excitatoryinhibitory balance [68], but our study invites reevaluation of such data with an emphasis on synaptic dynamics.Indeed, Hebbian plasticity can enhance the robustness of an attractor manifold against distractors [69].
Increasing the strength of Hebbian plasticity initially enriches network dynamics, indicated by an increased maximum Lyapunov exponent and attractor dimension.Beyond a certain plasticity strength, these metrics decrease, likely due to the increased presence of stable fixed points throughout phase space.This implies that there may be an optimal level of plasticity for task performance-one that is robust enough to enrich the dynamics compared to a nonplastic network but not so overpowering that it simplifies the dynamics through the overabundance of fixed points.This could be investigated by training plastic networks to solve tasks, e.g., using FORCE or backpropagation [40].
Our analyses point to a region of parameter space with large k and/or small g where it is possible that the behaviors of the fixed-point density and maximum Lyapunov exponent are more interesting than what we have explored.In particular, as (g, k) → (0, 2.02), fixed points are marginally stable by virtue of outliers, not the bulk (Sec.VI).Additionally, our Lyapunov analysis leaves open the possibility that the maximum exponent is negative for sufficiently large k even as N → ∞ (Sections V B, F; Fig. 9).Given that phase space is densely filled with stable fixed points in this parameter regime, it is possible that these features signal a novel glassy phase of the model.Such phases in dynamic networks remain poorly understood and are an important direction for future research [70,71].
Our study addresses chaotic networks that either generate activity autonomously or produce complex responses to inputs [72].A potentially desirable alternative property is stability, defined by the ability of a network to generate input-driven trajectories that are robust against perturbations; however, such networks cannot generate rich activity autonomously.Kozachkov et al. [73] analyzed a plastic network with the same governing equations as our model, but without quenched disorder, establishing conditions for stable dynamics [74].A key finding was that anti-Hebbian plasticity can promote stable dynamics, pointing to a possible function for ongoing plasticity in input-driven computations.
We considered plastic synapses with strengths weaker than random synapses by a factor of 1/ √ N .In experiments measuring synaptic strength changes, such plasticity may easily be overlooked despite its orderone network-level impact.Our model assumes all-toall connectivity; if neurons receive K < N inputs, the structure-to-randomness scaling is 1/ √ K, making detection of plasticity more feasible if K is not too large.
Humans and animals can remember a stimulus over a delay period, implying a form of rapid information storage in neural circuits, i.e., working memory (WM).Freezable chaos provides a new WM mechanism that we now compare to prior models.Most WM models rely on either cell-intrinsic or network-level mechanisms that support self-sustained activity.These "persistent activity" models are supported by some experimental studies, but undermined by others showing "activity-silent" WM [75].These latter studies suggest that information can be rapidly stored in synapses, requiring fast synaptic plasticity.Synaptic WM models typically use short-term facilitation (STF) due to the convention that fast plasticity is presynaptic [76,77].Because such plasticity cannot create attractor states of neuronal dynamics, STF models require existing symmetric structure in the synapses, potentially formed through prior Hebbian plasticity.A prototypical example is the model of Mongillo et al. [9] in which clusters of excitatory neurons with broad inhibition prime a network to function in a metastable regime.Due to STF, an activity pattern can be selectively sustained by providing transient external input to one of the clusters.A key requirement of this class of proposals is that the possible neuronal states to be stored are known in advance.
The inability of STF models to store novel patterns suggests the existence of fast Hebbian plasticity.This is at odds with conventional wisdom, but supported experimentally [16][17][18][19][20][21][22][23][24].For examples of Hebbian WM models, see [69,[78][79][80][81][82][83][84][85].Due to its Hebbian nature and ability to store novel patterns, freezable chaos aligns more with these proposals than with STF models.A crucial feature distinguishing freezable chaos from both STF and Hebbian WM models is that plasticity is deactivated, rather than activated, to store a pattern.Whereas other models require an external input carrying the pattern to be stored, this feature allows our model to store the neuronal state while it is engaged in strongly recurrent dynamics (in our random-network model, chaos).
Hinton and collaborators considered the possibility of a network performing a computation, saving its state in synapses, using neurons to perform a subroutine, and resuming computation from the saved state.This was termed "true recursion" [29,33,86].Freezable chaos provides a minimal example of this: the neuronal state can be saved by halting plasticity, allowing neurons to engage in arbitrary dynamics before returning to the saved state.In our model, halting plasticity leaves a globally stable fixed point, so neuronal dynamics during the subroutine must be driven by external inputs.An interesting question is whether halting plasticity can leave the network with a fixed point that coexists with a dynamic regime that can be used for recurrent computation.This could be implemented in an ad-hoc manner by turning on feedback loops upon halting plasticity.
This feature of freezable chaos suggests a method of detecting it in-vivo, namely, by "interupting" a task requiring strong recurrent dynamics, such as evidence integration, for variable periods of time [87].Finding that neurons involved in the computation show continuity in their activities at the beginning and end of the interruption period would be suggestive of freezable dynamics.This conclusion would be further supported if task performance degrades when synaptic plasticity is disabled by genetic or pharmacological manipulations [22].The activity expressed by the neurons during the interruption period would depend on whether and how they are recruited in this interval.
In other Hebbian WM models, an external neuronal or neuromodulatory signal is typically required to erase information stored in the synapses and return the network to a dynamic state.Freezable chaos avoids this requirement by leveraging fixed points that are stable with respect to neuronal, but not neuronal-synaptic dynamics.The presence or absence of a resetting signal could enable experimental disambiguation of our proposal.
Our model offers a mechanism for WM, but lacks a mechanism for long-term memory.This could be addressed by introducing slow Hebbian dynamics into J .For example, prior DMFT studies have taken J to be a static result of associative plasticity, resulting in various combinations of chaos and long-term memory retrieval [44,[88][89][90].Models like these could be extended by incorporating fast Hebbian synaptic dynamics atop this static structure.It would be particularly interesting if the short-and long-term dynamics could be made to interact, e.g., to implement memory consolidation.For example, during frozen chaos, if long-term plasticity was activated while the shorter-term plasticity of A(t) was disabled, the frozen state could be consolidated into J .
Humans have a WM capacity of ∼ 7 items, but freezable chaos can store just one item because synaptic plasticity is halted once a pattern is stored.One way of overcoming this limitation would be to use multiple A(t) matrices that can be independently modulated, corresponding either to different biophysical plasticity mechanisms or disjoint sets of synapses.Organisms might benefit from these different sets of synaptic variables possessing a hierarchy of timescales.
Large language models display a remarkable capacity for in-context learning: producing output that incorporates information or algorithms contained in the input [91].This is surprising because, in both neuroscience and machine learning, such learning is generally thought to require weight updates.One possibility is that these models have such capabilities because they emulate weight dynamics through attention layers [35].Nevertheless, explicit weight dynamics could benefit machine-learning models, for example by enabling in-context learning with fewer model parameters [33] or on longer sequences [36].An impediment to work in this direction is that weight dynamics are computationally costly.An important direction of work therefore pertains to ameliorating this burden, e.g., by formulating new forms of weight dynamics that are both expressive and have low computational complexity, or by exploring low-power neuromorphic architectures with dynamic weights.
Manuel Beiran, and Patricia Cooney for comments on a manuscript.The authors were supported by the Gatsby Charitable Foundation and NSF NeuroNex award DBI-1707398.dτ e −τ /p C(τ ).
(A7) 0 2000 4000 6000 8000 10000 x 800 (t) x 6 (t) x 5 (t) x 4 (t) x 3 (t) x 2 (t) x For large p and k, the plastic network has two approximate dynamic solutions: a fast solution in which plasticity is averaged out, with dynamic timescale τ * = O(1); and a slow solution in which synapses drag neurons, with τ * ≫ p.In finite-size networks, both behaviors can be realized depending on initial conditions.Moreover, for appropriately tuned values of the parameters, the system can switch between these behaviors in a bistable manner (Fig. 8).A QEP can be solved by linearization, which entails forming a 2N -dimensional linear eigenvalue problem whose eigenvalues coincide with those of the QEP [92].Applying the aforementioned identity to the determinant defining the characteristic polynomial of the reduced Jacobian M shows that M is one such linearization of the QEP.Thus, we recover the originally derived spectrum.There are infinitely many linearizations given by matrices related to M by similarity transformation, so in this sense, the reduced Jacobian is not unique.
[59] In chaotic states, this is a result of the network being locally destabilized by a continuous band of eigenvalues in analogy with the nonplastic network [41].At fixed points, this is a result of the marginal stability of typical fixed points, which have a continuous band of Jacobian eigenvalues brushing against the stability line, Re(λ) = 0.

FIG. 1 .
FIG. 1.(a) Dynamics of a pair of neurons (top panel) and of the synapses through which they are reciprocally coupled (bottom panel).Synapses fluctuate about quenched random strengths (dashed lines) in response to pre-and postsynaptic activity according to a Hebbian rule.(b) Left: phase diagram of the plastic network for p = 2.5.Right: example neuronal traces xi(t) from simulations of each phase-diagram region, with parameters given by the location of the associated square marker.

FIG. 2 .
FIG. 2. Chaotic states with g > 1.(a) C(τ ) from the DMFT (solid lines) and in simulations (dashed lines) for g = 2, p = 2.5, and various values of k (indicated by the lower color bar).(b) Dynamic timescale τ * [Eq.(9)] as a function of k for various values of g.Dotted line indicates p. (c) Power spectra (normalized such that S(0) = 1) for the autocovariance functions shown in (a).For anti-Hebbian (k < 0) power spectra, triangular markers indicate an oscillatory frequency computed from the zero-crossings of C(τ ).

FIG. 3 .
FIG. 3. Chaotic states for g < 1.(a) C(τ ) from the DMFT (lines) and in simulations for g = 0.9, p = 2.5, and various values of k.(b) Left: median log-lifetime of transient activity, before collapsing to the trivial fixed point, as a function of N for g = 0.9 and values of k from (a).Right: histograms of the log-lifetime of transient chaos, corresponding to stars in the left plot.Simulations were terminated at time Tsim = 10 7 .(c) Curves y(s) for solutions of the model of Stern et al. [53] for various values of g.Dashed horizontal lines correspond to different values of k, intersecting y(s) at self-consistent solutions of Eq. (10).

FIG. 4 .
FIG. 4. Spectra of the Jacobian for p = 2.5 and various values of g (running horizontally) and k (running vertically).Lines: predicted boundary curves from random matrix theory and DMFT.Dots: spectra measured in simulations of chaotic plastic networks.Modes are colored by fa(λ), the weight on the synaptic part of the corresponding eigenvector of the reduced Jacobian [Eq.(19)].The red dot at λ = −1/p indicates a delta function of N 2 − N synaptic modes.For anti-Hebbian (k < 0) spectra, triangular markers indicate an oscillatory frequency computed from the zero-crossings of C(τ ).

FIG. 5 .
FIG. 5. (a) Maximum Lyapunov exponent λmax, computed by a perturbation method, throughout (g, k) parameter space with p = 2.5, N = 4000.White: quiescence.Hatched: convergence to nonzero fixed points.(b) Histograms of Lyapunov spectra, computed using tangent-vector propagation with N = 900, for p = 2.5 and various values of g (running horizontally) and k (running vertically).Black outline for k = 0 histograms: spectra of nonplastic network.The red triangle marks −1/p, where there are O(N 2 ) exponents in the full spectrum.(c) λmax as a function of k for various values of g.Solid lines: estimate from tangent-vector propagation.Dashed lines: estimate from perturbation method.(d) Attractor dimension divided by N as a function of k for various value of g.

FIG. 6 .
FIG. 6. Stable fixed points.(a) DMFT participation ratio [Eq.(23)] of A(t) as a function of k for g = 3 and various values of p.(b) Top: empirical participation ratio [Eq.(22)] of A(t) with g = 2, k = 2.25, and p = 2.5.The network settles to a stable fixed point.Dashed line: DMFT value.Bottom: neuronal traces during the same settling event.(c) Log-number of stable fixed points per neuron as a function of k for various values of g.(d) Log-number of stable fixed points per neuron throughout (g, k) parameter space.(e) Jacobian spectra at fixed points.Lines: predicted boundary curve from mean-field analysis.Dots: spectra after settling to a fixed point in many simulations.Top: g = 2, k = 2.25, p = 2.5 [as in (b)].Bottom: g = 0.5, k = 1.9, p = 2.5.

2 FIG. 7 .
FIG. 7. Freezable chaos.(a)Neuronal traces in nonfreezable (top, k = 0.6), semifreezable (middle, k = 1.2), and freezable (bottom, k = 1.55) chaos for g = 2.25 and p = 2.5.At the first vertical line, A(t) is halted.In (semi-)freezable chaos, a memory of the halt-time neuronal state persists.Stability is demonstrated by perturbing the neurons (second line) and clamping them near zero (third line).In (semi-)freezable chaos, neurons relax back to the stored pattern after these manipulations.At the fourth line, A(t) is released.This process is repeated at the fifth line.(b) Demonstration of the alignment between J and A(t).At the first vertical line, A(t) is halted.At the second line, J , but not the halted A(t), is shuffled, destroying the neuronal fixed point.At the third line, A(t) is released.This process is repeated at the fourth line.g = 2.5, k = 2.55, and p = 2.5 [as in the freezable-chaos plot in (a)].(c) Time-dependent overlap Q(t) from the two-replica DMFT (solid lines) and in simulations (dashed lines) for g = 2, p = 2.5, and various values of k.(d) Halt-time overlap Q(0) from the two-replica DMFT (lines) and in simulations (triangular markers) as a function of k for p = 2.5 and various values of g.Markers are colored according to whether simulations show semifreezable or freezable activity.(e) Same as (d), but using the zero-time correlation coefficient ρ(0).In (c-e), We show only the positive solutions for Q(t) and ρ(t), which are realized upon halting synapses; there is also a symmetric negative solution.