Robustness of STDP to spike timing jitter

In Hebbian plasticity, neural circuits adjust their synaptic weights depending on patterned firing. Spike-timing-dependent plasticity (STDP), a synaptic Hebbian learning rule, relies on the order and timing of the paired activities in pre- and postsynaptic neurons. Classically, in ex vivo experiments, STDP is assessed with deterministic (constant) spike timings and time intervals between successive pairings, thus exhibiting a regularity that differs from biological variability. Hence, STDP emergence from noisy inputs as occurring in in vivo-like firing remains unresolved. Here, we used noisy STDP pairings where the spike timing and/or interval between pairings were jittered. We explored with electrophysiology and mathematical modeling, the impact of jitter on three forms of STDP at corticostriatal synapses: NMDAR-LTP, endocannabinoid-LTD and endocannabinoid-LTP. We found that NMDAR-LTP was highly fragile to jitter, whereas endocannabinoid-plasticity appeared more resistant. When the frequency or number of pairings was increased, NMDAR-LTP became more robust and could be expressed despite strong jittering. Our results identify endocannabinoid-plasticity as a robust form of STDP, whereas the sensitivity to jitter of NMDAR-LTP varies with activity frequency. This provides new insights into the mechanisms at play during the different phases of learning and memory and the emergence of Hebbian plasticity in in vivo-like activity.

where is membrane potential; ! and ! are leak conductance and reversal potential respectively; !"#!$ , !"#$% , !"## and !"#$% are currents through AMPAR, NMDAR, Ltype VSCC (v1.3) and TRPV1, respectively, and AEA stands for anandamide. NMDAR and AMPAR were modeled with two-state kinetic models and 1.0mM Mg 2+ (Destexhe et al., 1995) whereas the model and parameters for the Ca v 1.3 VSCC current was taken from Wolf et al. (2005). The TRPV1 current, including its dependence on AEA, was modeled as: where g TPRV1 is the maximal conductance of TRPV1 and the mathematical expression for the P TRPV1 open was taken from Matta and Ahern (2007).
Biochemical signaling. Free cytosolic calcium is one of the main signaling actors in the model. To model its dynamics, we assumed it can be transferred from/to two main sources: (i) extracellular calcium, via the plasma membrane channels of eq. (SI3) above and (ii) the endoplasmic reticulum (ER), via the IP3-dependent Calcium-Induced Calcium Release (CICR) system. Hence, the concentration of free cytosolic calcium C was computed according to ! ( ) !" !" = !! ! ! − !"#$% + !"#$ + !"#$% + !"## + !"#$% − where the fluxes !"#$ , !"#$% , !"#$ from and to the ER in the CICR system were taken from the model of De Pittà et al. (2009) -see also Cui et al. (2016) and !"#$% , !"## and !"#$% are the calcium fluxes from the plasma membrane channels (eq. SI3), computed as ! = ! ⋅ ! were the ! are constants, ∈ {NMDAR,VSCC,TRPV1}. ! is the basal cytosolic calcium level resulting from equilibration with calcium diffusion out of the cell and ! ! is the corresponding time scale. The CICR model also includes the dynamics of the fraction of inactive IP3 (inositol triphosphate) receptor channels (IP3R), h Likewise, the dynamics of C ER , the calcium concentration in the endoplasmic reticulum (ER) was given by where ρ ER is the ER to cytoplasm volume ratio. In eq.(SI5) and (SI7), ! , x = C or C ER, is a time scaling factor accounting for the presence of endogenous calcium buffers and expressed according to Our model also accounts for the biochemical pathways leading to the production of the where CaMKII* is the amount of activated (phosphorylated) CaMKII (see below); !"#$ is the fraction of active DAG Lipase α and DAGL its total (activated+not activated) concentration. Eq. (SI9-SI11) account for IP3 degradation by IP3 3-kinase (3K) and inositol polyphosphate 5-phosphatase (5P), while DAG and 2AG are degraded by DAG kinase (DAGK) and monoacylglycerol lipase (MAGL), respectively.
2-AG and AEA are retrograde signaling molecules that are produced in the postsynaptic neuron but diffuse to the presynaptic cell, where they activate CB 1 R. We modeled CB 1 R activation by 2-AG and AEA using a simple three-state kinetic model: open (x CB1R ), desensitized (d CB1R ) and inactivated (i CB1R ): where eCB = 2-AG + 0.10 AEA and mass conservation implies CB1R + CB1R + CB1R = 1.
The open fraction x CB1R was then used to compute CB 1 R activation as where ! is a constant that accounts for the modulation of presynaptic plasticity by other pathways and !"#$ quantifies the strength of CB 1 R activation on presynaptic plasticity.
In our model, CB 1 R activation (y CB1R ) controls the presynaptic weight W pre according to the following rule: W pre decreases for intermediate values of y CB1R , i.e. when y CB1R is comprised between two tLTD thresholds ( LTD start < CB1R < LTD stop ) whereas W pre increases when y CB1R is larger than a tLTP threshold ( CB1R > LTP start ). Following Cui et al. (2016), we implemented this rule as Here determines the change of presynaptic plasticity, with a time scale ! !"# set to yield rapid changes of !"# for large CB1R values and very slow changes at very low CB1R : To account for experimental observation that the presynaptic weight ranges from about 50 to 300%, !"# was clipped to 3.0 (hard bound).
Postsynaptic plasticity in the model was based on the activation by calcium of calmodulin and CaMKII and the regulation of this system by PKA, calcineurin and protein phosphatase 1 (PP1). The model of Cui et al (2016) for this subsection is based on the model proposed in Graupner & Brunel (2007). In this model the concentration of the calcium/calmodulin complex with four calcium ions bound (CaM) is considered at equilibrium and computed as: where CaMT is the total calmodulin concentration and K i stands for the equilibrium constant of the binding of the i th calcium ion to calmodulin. Each CaMKII enzyme consists of two subunits, each of which contains 6 phosphorylation sites, thus defining 14 phosphorylation states per subunit. Denoting y i , i=0…13 the concentrations of subunit in phosphorylation state i, the model of Graupner & Burnel (2007) expresses their dynamics as: is the total CaMKII concentration. In eq. (SI21) above the probabilities that CaM binds to phosphorylated ( ) and dephosphorylated ( * ) subunit of CaMKII are computed as = /( + ! ) and * = /( + ! ) . Moreover, the rate of subunit 1 is the concentration of PP1. * is the total concentration of phosphorylated subunits of CaMKII computed accross all possible states of phosphorylation * = where ! is the number of the phosphorylated subunits of CaMKII in state and ! . Finally, the model of Graupner & Brunel (2007) assumes that, in addition to calcineurin (CaN), PKA activity depends on CaM according to a Hill-function activation: The interaction between PP1 and I1 is then described as: where 1 ! and 1 ! are total PP1 and I1 concentrations, respectively.
In agreement with Cui et al. (2016) and Graupner & Brunel (2007), we assumed that postsynaptic plasticity is directly proportional to the calcium-dependent activation of CaMKII and set: is the maximal concentration of activate (phosphorylated).
Finally, the total synaptic weight was taken as the product of the presynaptic and the postsynaptic contributions above: Model implementation and numerics. Our mathematical model comprises 36 ordinary differential equations and close to 150 parameters, among which more than one half is constrained by experimental data. Initial conditions were set to the steady-state of each variable in the absence of stimulation. Numerical solution was obtained with the LSODA solver from the ODEPACK fortran77 library with absolute and relative tolerances both equal to 10 -7 . Initial conditions were set to the steady-state of each variable in the absence of stimulation. Numerical integration proceeded until the synaptic weights reach stable values (typically observed around t ≈ 5min after the stimulation protocol). We used the final values of the pre-and postsynaptic weights to compute the total synaptic weight change due to the stimulations. Importantly, with the exception of the stimulation protocols, we have used the exact same equations and parameter values as in Cui et al. (2016). The list of parameters and their estimated values is given in S2 Table. The current study employs stochastic simulations since the stimulation protocol is stochastic (while the rest of the model is deterministic). Therefore, the model was calibrated using experimental data from deterministic stimulation protocols (Cui et al., 2016) and we tested here whether it can make successful predictions when we applied stochastic stimulation protocols, for which the model was not calibrated.

S3 Text -Mechanisms of STDP expression in the model
We give here an overview of the mechanisms that give rise to the expression of eCB-tLTP, eCB-tLTD and NMDAR-tLTP in our mathematical model. All the examples given here concern regular (deterministic) STDP protocols with constant IPI and spike timings.
Our model combines the two signaling pathways involved in cortico-striatal STDP: a first signaling pathway leading from NMDAR to calmodulin and CaMKII with a second that couples mGluR and cytosolic calcium to eCB production and the resulting activation of CB 1 R (see Fig. 1f1). In the model, we assume that the total synaptic weight (W Total ) is the product of presynaptic (W Pre ) and postsynaptic (W Post ) weights (see Supporting Information S1).
On the postsynaptic side, W Post is proportional to the amount of CaMKII activated by the NMDAR pathway. CaMKII in our model forms a bistable system that settles at long times either on a "DOWN" state where CaMKII is mostly inactive (no plasticity) or on an "UP" state characterized by high levels of activated CaMKII ( decreases. This decay ER calcium eventually compensates the effect of IP3 accumulation on the calcium-induced calcium release, thus stabilizing the amplitude of the calcium peaks for N pairings >50. Moreover, since the width of the postsynaptic calcium peak in the model is larger with post-pre than pre-post pairings (Fig S2b), the fraction of calcium-activated DAGLα is much larger for post-pre pairings. As a result, the amplitude of eCB peaks and, ultimately, the amplitude of the fraction of activated CB 1 R (y CB1R ), show a more pronounced biphasic profile ( Fig S2d). The biphasic trend is further amplified at the level of CB 1 R activation because of CB 1 R receptor desensitization that amplifies the decay above 20 pairings (Fig S2d). The amplitude of the y CB1R peaks increases sharply for the first 10-20 pairings, and progressively decays afterwards to converge to constant amplitude. As show in the figure, y CB1R reaches large values only for short post-pre pairings (Δt STDP around -15 ms) while even short pre-post pairings (0< Δt STDP <10 ms) do not give rise to such large amplitude peaks.