Diverse synaptic plasticity mechanisms orchestrated to form and retrieve memories in spiking neural networks

Synaptic plasticity, the putative basis of learning and memory formation, manifests in various forms and across different timescales. Here we show that the interaction of Hebbian homosynaptic plasticity with rapid non-Hebbian heterosynaptic plasticity is, when complemented with slower homeostatic changes and consolidation, sufficient for assembly formation and memory recall in a spiking recurrent network model of excitatory and inhibitory neurons. In the model, assemblies were formed during repeated sensory stimulation and characterized by strong recurrent excitatory connections. Even days after formation, and despite ongoing network activity and synaptic plasticity, memories could be recalled through selective delay activity following the brief stimulation of a subset of assembly neurons. Blocking any component of plasticity prevented stable functioning as a memory network. Our modelling results suggest that the diversity of plasticity phenomena in the brain is orchestrated towards achieving common functional goals.

(c) Combined input-output relation g of a non-adapting rate neurons from a with STP transfer function from c (black) and excitatory self-feedback w ee = 1.9 (black) and w ee = 2.0 (red). The sigmoidal shape of the effective transfer function gives rise to three intersection points (α 1 , α 2 , α 3 ) with the diagonal (dashed) in the regime < 100 Hz. (d) Lowest possible firing rate at fixed point α 3 as a function of the time constant τ d of recovery from synaptic depression. If synaptic depression vanishes (τ d → 0), the firing rate inside an active assembly is close to 300 Hz 1 . (e) Schematic representation of the rate model. The excitatory network neurons are split into two populations. ν a describes the firing rate of the "assembly" neurons, whereas ν b characterizes the rate of the "background" neurons. Both populations receive external input from external populations ν a ext and ν b ext respectively. Finally, the excitatory populations project to and receive inhibitory input from a shared pool of inhibitory neurons characterized by the firing rate ν inh . All external and excitatory-to-excitatory synaptic connections are modeled as plastic with orchestrated plasticity rules. For simplicity, connections from and to the inhibitory population are taken as static. (f ) Evolution of firing rate (top row) and weight dynamics (bottom row) of the rate model on three different timescales (columns) for constant input ν a ext = ν b ext = 10 Hz. The dashed lines in the top row show the rate dynamics for the same network but with different initial conditions. The insets in the top row show the same network, but without transmitter induced plasticity (δ = 0), which leads to a decay of firing rates to zero. The bottom row shows weight dynamics for one initial state of the network. Solid lines represent the synaptic weights as colored in a. The dashed lines represent the evolution of the respective reference weightsw which were all initialized at zero in this simulation. w ae and w be converge to higher values than w aa , w bb , w ab , or w ba , because the external input pools fire at higher rates (10 Hz) than the stationary rates ν a and ν b of the two excitatory populations. For δ = 0 the firing rates converge to zero.    Fig. 1 f, but for a protocol in which population ν a receives repeated external high frequency stimulation (20 Gaussian rate pulses equally spaced over the interval [15 s,1 h]). The inset in panel P3 illustrates that firing activity is non-zero (≈ 1 Hz) in the absence of stimulation. The weight w aa (solid black) within the assembly and its corresponding reference weightw aa (dashed black) converge to stable hight values of approximately 0.5. (b) Same as above, but during recall dynamics. The assembly population receives brief excitatory stimuli (Gaussian rate profile; σ = 0.5 s) at times 5 s, 300 s and 900 s which are each followed by an equally brief inhibitory stimulus (times: 15 s, 330 s and 930 s). The insets in the top row for a network without heterosynaptic plasticity (β = 0), show that firing rates immediately explode and saturate at the maximum firing rate of 300 Hz, whereas in the normal network retrieval happens at rates of about 100 Hz and is stopped by the inhibitory pulses.   Figure 3. Histograms (from top left to bottom right): 1) Neuronal firing rates determined over the interval indicated below the raster plot. For clarity only neurons that fired at least one spike during the interval are considered for the histogram. 2) Bar plot illustrating the relative numbers of active (at least one spike during one hour) and silent neurons. 3) CV ISI distribution computed over the interval indicated below the raster plot on the left. 4) Excitatory synaptic weights as determined at the end of the simulation (at t = 300 s for b, t = 2h otherwise). (a) Blocked homeostatic metaplasticity of LTD in a network with weak random initial input (cf. Fig. 3). The network fails to develop selectivity. (b) Blocking transmitter-induced plasticity causes many neurons to fall and remain silent. (c) Blocked heterosynaptic plasticity causes rapidly increasing firing rates without learning or emergence of delay activity (simulation stopped after 300 s). Note the bi-modal weight distribution. A weight value of 5 corresponds to the maximally allowed weight. (d) Blocking inhibitory plasticity allows the network to develop delay activity, but at highly elevated firing rates. Compare firing rates (top left histogram) with Fig. 3 Supplementary Fig. 4). (b) Covariance matrix (color coded) of evoked population activity between with previously established stimulus preferences and the novel pattern (grumpy face). The numbers encode the number of units which responded with more than 10 Hz to a single (diagonal) or two patterns (off diagonal). (c) Averaged synaptic weights in matrix representation between the four previously learned patterns and the novel pattern. (d) Schematic representation of different input patterns (left), presented to the network repeatedly after the initial learning session, and how they were typically classified by the network (right).

Supplementary Tables
Supplementary Table 1 -Tabular description of network model.

A Model Summary Populations
Three: excitatory, inhibitory, Poisson input Topology -Connectivity Random sparse connectivity (fixed probability)

Input
Poisson input with additional rate coded input patterns Measurements Spike activity, synaptic efficacy (continuously of a subset of synapses and of all synapses at the end of a simulation) for the definition of short term plasticity variables x j and u j see D2 below.
i (t) = 0 except in simulations with adaptation on multiple timescales:

IF neuron, used for all inhibitory neurons Type
The non-adaptive inhibitory neurons were implemented identically except of the omission of the adaptation terms g a i = g b i = 0.

D2 Plasticity Model Name
Short term plasticity following 2 .

Type
Depressing and facilitating short term dynamics Acts on EE,EI,StimE

Plasticity Model Name
Plasticity rule for excitatory synapses Type Triplet STDP rule with pre and post offset terms and slower consolidation dynamics Acts on EE and StimE

Synaptic traces
which is integrated with step size of 1.2s.

Metaplasticity
Where this is mentioned explicitly we allow B i to have a slow time dependence 3,4 :

D4 Plasticity Model Name
Inhibitory Spike Timing Dependent Plasticity (iSTDP) Type Symmetric iSTDP with a constant offset for presynaptic spikes and a global modulation Acts on IE Synaptic traces Global modulation which is constrained to w min = 0 < w ij < 5 = w max . S13 E Input Name Stimulus group Type Rate coded Poisson input Size 4,096 Poisson neurons Firing rates where the ξ µ i are 64×64 grayscale images interpreted as 1D vectors and normalized to [0,1]. Only a single pattern can be active at a time. A pattern stays active during a finite period drawn from an exponential distribution with mean T On . Each pattern activation is followed by a period of inactivity with a duration drawn from another exponential distribution with mean T Off .

F Measurements Type Description
Spike activity for raster and activity plots and weight matrices at the end of the simulation. Supplementary

Supplementary Methods
Population rate model. To better understand the stability properties of our learning rules, we studied a population rate model. In particular, we explicitly modeled one inhibitory population and two excitatory populations which are characterized by their firing rates ν inh , ν a and ν b , respectively. The two excitatory populations correspond to an "assembly" population with firing rate ν a and a "background" population (ν b ). All three populations interact with each other via synaptic couplings w xy , each of which has to be understood as the combined effect of many individual synaptic connections which "see" similar pre-and postsynaptic activity. A schematic representation of the model is given in Supplementary  Fig. 1 e. We modeled the steady state response of a neuronal population to an input h where the maximum firing rate ν max was set to 300 Hz (see inset, Supplementary Fig. 1 a). Each neuronal population followed the external input h on the rapid timescale of population dynamics τ a pop = τ b pop , which we fixed to 0.1 s (τ inh pop = 0.05s for the inhibitory population). The overall temporal evolution of a population with index x, where x stands for a, b, or inh can be described as follows where f = 0.15 characterizes the size of the putative assembly and differences in gain of the different populations were modulated by the gain factors γ exc = 20, γ ext = 10 and γ inh = 2. To capture the effect of short-term plasticity (STP) of excitatory connections in the rate model we introduced the steady-state STP function s based on previous work 2 , defined as s(ν) ≡ x(ν)u(ν)ν with where τ d , τ f and U 0 refer to the parameters of the STP model of Eqs. (9) and (10) in the main text 2 (cf. Supplementary Fig. 1 b). Weights involving the inhibitory population were set statically to w a,inh = w b,inh = w inh,inh = 0.2 and w inh,a = w inh,b = 0.6 like in the main simulation. All other weights were plastic following the equivalent rate model of the orchestrated plasticity rules in which all parameters, if not mentioned otherwise, were identical to those in the plasticity model for the spiking neurons and ξ(ν) is the burst detector, which, in the spiking model, contains the third power of a synaptic trace multiplied by the postsynaptic spike train (see Eq. (13) in the main text). To convert it to a temporal average in the rate model 4 we compute the expectation value of the third moment of this synaptic trace zi (t) 3 assuming Poisson firing. The expression can be obtained conveniently by taking the Laplace transform of the stationary distribution of z − i when interpreting it as a random variable (see the derivation below) which yields ξ(ν y ) = 1 6 2 + 9ν y τ + 6ν 2 y τ 2 ν 2 y . Since neurons in our spiking simulation are refractory, a fact which is not accounted for in the Poisson assumption, the parameter β in Supplementary Eq. (7) was lowered to β = 0.0065 in the rate model. Finally, the consolidation dynamics behindw were the same as in the spiking model (Supplementary Table 2), with a minor parameter change to P = 25. To generate the plots in the Supplementary Figs. 1 and 2, the 15 state variables (three rates, 6 weights and 6 reference weights) were integrated using Python code with the SciPy module. External input to both populations was constant at 10Hz except during external stimulation. To simulate the stimuli shown in Supplementary Fig. 2 a we transiently increased the external input rate to the assembly population ν ext,x . Specifically we added to the spontaneous firing rate of 10Hz a Gaussian temporal profile with a standard deviation of 0.5 s and an amplitude of 80 Hz ( Supplementary Fig. 2 a) and 15Hz for recall and -10Hz (resulting at 0Hz total firing rate) to terminate a recall state ( Supplementary Fig. 2 b).
Derivation of the moments of a synaptic trace for Poisson firing statistics. We can describe the probability flux of a synaptic trace z defined by the ordinary differential equation with time constant τ and spike train S(t) by the following partial differential equation ∂ t p(z, t) = 1 τ ∂ z (zp(z, t)) − λp(z, t) + λp(z − a, t) (10) = 1 τ p(z, t) + z τ ∂ z p(z, t) − λp(z, t) + λp(z − a, t) where the terms involving λ describe the sink and source terms that capture the jumps of size a caused by spikes at rate λ. We now require stationarity ∂ t p = 0, and find