A stochastic model of hippocampal synaptic plasticity with geometrical readout of enzyme dynamics

  1. Yuri Elias Rodrigues
  2. Cezar M Tigaret
  3. Hélène Marie
  4. Cian O'Donnell
  5. Romain Veltz  Is a corresponding author
  1. Université Côte d’Azur, France
  2. Institut de Pharmacologie Moléculaire et Cellulaire (IPMC), CNRS, France
  3. Inria Center of University Côte d’Azur (Inria), France
  4. Neuroscience and Mental Health Research Innovation Institute, Division of Psychological Medicine and Clinical Neurosciences,School of Medicine, Cardiff University, United Kingdom
  5. School of Computing, Engineering, and Intelligent Systems, Magee Campus, Ulster University, United Kingdom
  6. School of Computer Science, Electrical and Electronic Engineering, and Engineering Mathematics, University of Bristol, United Kingdom


Discovering the rules of synaptic plasticity is an important step for understanding brain learning. Existing plasticity models are either (1) top-down and interpretable, but not flexible enough to account for experimental data, or (2) bottom-up and biologically realistic, but too intricate to interpret and hard to fit to data. To avoid the shortcomings of these approaches, we present a new plasticity rule based on a geometrical readout mechanism that flexibly maps synaptic enzyme dynamics to predict plasticity outcomes. We apply this readout to a multi-timescale model of hippocampal synaptic plasticity induction that includes electrical dynamics, calcium, CaMKII and calcineurin, and accurate representation of intrinsic noise sources. Using a single set of model parameters, we demonstrate the robustness of this plasticity rule by reproducing nine published ex vivo experiments covering various spike-timing and frequency-dependent plasticity induction protocols, animal ages, and experimental conditions. Our model also predicts that in vivo-like spike timing irregularity strongly shapes plasticity outcome. This geometrical readout modelling approach can be readily applied to other excitatory or inhibitory synapses to discover their synaptic plasticity rules.

Editor's evaluation

Synaptic plasticity is a ubiquitous but also highly complex phenomenon and developing a unifying description has been challenging. This study presents a realistic biophysical model of plasticity induction, with a novel read-out of CaMKII and Calcineurin. It is able to describe a wide range of experimental results and sets a new benchmark for realistic computational models.



To understand how brains learn, we need to identify the rules governing how synapses change their strength in neural circuits. The dominant principle at the basis of current models of synaptic plasticity is the Hebb postulate (Hebb, 1949) which states that neurons with correlated electrical activity strengthen their synaptic connections, while neurons active at different times weaken their connections. In particular, spike-timing-dependent plasticity (STDP) models (Blum and Abbott, 1996; Gerstner et al., 1996; Eurich et al., 1999) were formulated based on experimental observations that precise timing of pre- and post-synaptic spiking determines whether synapses are strengthened or weakened (Debanne et al., 1996; Tsodyks and Markram, 1997; Bi and Poo, 1998; Markram et al., 2011). However, experiments also found that plasticity induction depends on the rate and number of stimuli delivered to the synapse (Dudek and Bear, 1992; Sjöström et al., 2001), and the level of dendritic spine depolarisation (Artola et al., 1990; Magee and Johnston, 1997; Sjöström and Häusser, 2006; Golding et al., 2002; Hardie and Spruston, 2009). The lack of satisfactory plasticity models based solely on neural spiking prompted researchers to consider simple models based on synapse biochemistry (Castellani et al., 2001; Castellani et al., 2005). Following a proposed role for postsynaptic calcium (Ca2+) signalling in synaptic plasticity (Lisman, 1989), previous models assumed that the amplitude of postsynaptic calcium controls long-term alterations in synaptic strength, with moderate levels of calcium causing long-term depression (LTD) and high calcium causing long-term potentiation (LTP) (Shouval et al., 2002; Karmarkar and Buonomano, 2002). However, experimental data suggests that calcium dynamics are also important (Yang et al., 1999; Mizuno et al., 2001; Wang et al., 2005; Nevian and Sakmann, 2006; Tigaret et al., 2016). As a result, subsequent phenomenological models of plasticity incorporated slow variables that integrate the fast synaptic input signals, loosely modelling calcium and its downstream effectors (Abarbanel et al., 2003; Rubin et al., 2005; Rackham et al., 2010; Clopath and Gerstner, 2010; Kumar and Mehta, 2011; Graupner and Brunel, 2012; Honda et al., 2013; Standage et al., 2014; De Pittà and Brunel, 2016). Concurrently, more detailed models tried to explicitly describe the molecular pathways integrating the calcium dynamics and its stochastic nature (Cai et al., 2007; Shouval and Kalantzis, 2005; Miller et al., 2005; Zeng and Holmes, 2010; Yeung et al., 2004). However, even these models do not account for data showing that plasticity is highly sensitive to physiological conditions such as the developmental age of the animal (Dudek and Bear, 1993; Meredith et al., 2003; Cao and Harris, 2012; Cizeron et al., 2020), extracellular calcium and magnesium concentrations (Mulkey and Malenka, 1992; Inglebert et al., 2020) and temperature (Volgushev et al., 2004; Wittenberg and Wang, 2006; Klyachko and Stevens, 2006). This limits the predictive power of this class of plasticity models.

An alternative approach taken by several groups (Bhalla and Iyengar, 1999; Jędrzejewska-Szmek et al., 2017; Blackwell et al., 2019; Chindemi et al., 2022; Zhang et al., 2021) was to model the complex molecular cascade leading to synaptic weight changes. The main benefit is the direct correspondence between the model’s components and biological elements, but at the price of numerous poorly constrained parameters. Additionally, the increased number of nonlinear equations and stochasticity makes fitting to plasticity experiment data difficult (Mäki-Marttunen et al., 2020).

Subtle differences between experimental STDP protocols can produce completely different synaptic plasticity outcomes, indicative of finely tuned synaptic behaviour as detailed above. To tackle this problem, we devised a new plasticity rule based on a bottom-up, data-driven approach by building a biologically-grounded model of plasticity induction at a single rat hippocampal CA3–CA1 synapse. We focused on this synapse type because of the abundant published experimental data that can be used to quantitatively constrain the model parameters. Compared to previous models in the literature, we aimed for an intermediate level of detail: enough biophysical components to capture the key dynamical processes underlying plasticity induction, but not the detailed molecular cascade underlying plasticity expression; much of which is poorly quantified for the various experimental conditions we cover in this study.

Our model is centred on dendritic spine electrical dynamics, calcium signalling and immediate downstream molecules, which we then map to synaptic strength change via a conceptually new dynamical, geometric readout mechanism. It assumes that a compartment-based description of calcium-triggered processes is sufficient to reproduce known properties of LTP and LTD induction. Also, neither spatially-resolved elements (Bartol et al., 2015; Griffith et al., 2016) nor calcium-independent processes are required to predict the observed synaptic change. Crucially, the model also captured intrinsic noise based on the stochastic switching of synaptic receptors and ion channels (Yuste et al., 1999; Ribrault et al., 2011). We report that, with a single set of parameters, the model could account for published data from spike-timing and frequency-dependent plasticity experiments, and variations in physiological parameters influencing plasticity outcomes. We also tested how the model responded to in vivo-like spike timing jitter and spike failures.


A multi-timescale model of synaptic plasticity induction

We built a computational model of plasticity induction at a single CA3-CA1 rat glutamatergic synapse (Figure 1). Our goal was to reproduce results on synaptic plasticity that explored the effects of several experimental parameters: fine timing differences between pre and postsynaptic spiking (Figure 2 and Figure 3); stimulation frequency (Figure 4); animal age (Figure 5); external calcium and magnesium (Figure 6); stochasticity in the firing structure (Figure 7), temperature and experimental conditions variations (Supplemental Information). Where possible, we set parameters to values previously estimated from synaptic physiology and biochemistry experiments, and tuned the remainder within physiologically plausible ranges to reproduce our target plasticity experiments (see Materials and methods).

The model components are schematized in Figure 1 (full details in Materials and methods). For glutamate release, we used a two-pool vesicle depletion and recycling system, which accounts for short-term presynaptic depression and facilitation. When glutamate is released from vesicles, it can bind to the postsynaptic α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid and N-methyl-D-aspartate receptors (AMPArs and NMDArs, respectively), depolarizing the spine head by ∼30 mV (Kwon et al., 2017; Jayant et al., 2017; Beaulieu-Laroche and Harnett, 2018). The dendritic spine membrane depolarization causes the activation of voltage-gated calcium channels (VGCCs) and removes magnesium ([Mg2+]o) block from NMDArs. Backpropagating action potentials (BaP) can also depolarize the spine membrane by up to ∼60 mV (Kwon et al., 2017; Jayant et al., 2017). As an inhibitory component, we modelled a gamma-aminobutyric acid receptor (GABAr) synapse on the dendrite shaft (Destexhe et al., 1998). Calcium ions influx through VGCCs and NMDArs can activate SK potassium channels (Adelman et al., 2012; Griffith et al., 2016), which provide a tightly-coupled local negative feedback limiting spine depolarisation. Upon entering the spine, calcium ions also bind to calmodulin (CaM). Calcium-bound CaM in turn activates two major signalling molecules (Fujii et al., 2013): Ca2+/calmodulin-dependent protein kinase II (CaMKII) and calcineurin (CaN) phosphatase, also known as PP2B (Saraf et al., 2018). We included these two enzymes because of the overwhelming evidence that CaMKII activation is necessary for Schaffer-collateral LTP (Giese et al., 1998; Chang et al., 2017), while CaN activation is necessary for LTD (O’Connor et al., 2005; Otmakhov et al., 2015). Later, we show how we map the joint activity of CaMKII and CaN to LTP and LTD. Ligand-gated ion channels (ionotropic receptors) and voltage-gated ion channels have an inherent random behavior, stochastically switching between open, closed and internal states (Ribrault et al., 2011). If the number of ion channels is large, then the variability of the total population activity becomes negligible relative to the mean (O’Donnell and van Rossum, 2014). However individual hippocampal synapses contain only small numbers of receptors and ion channels, for example they contain ∼10 NMDArs and <15 VGCCs (Takumi et al., 1999; Sabatini and Svoboda, 2000; Nimchinsky et al., 2004), making their total activation highly stochastic. Therefore, we modelled AMPAr, NMDAr, VGCCs and GABAr as stochastic processes. Presynaptic vesicle release events were also stochastic: glutamate release was an all-or-none event, and the amplitude of each glutamate pulse was drawn randomly, modelling heterogeneity in vesicle size (Liu et al., 1999). The inclusion of stochastic processes to account for an intrinsic noise in synaptic activation (Deperrois and Graupner, 2020) contrasts with most previous models in the literature, which either represent all variables as continuous and deterministic or add an external generic noise source (Bhalla, 2004; Antunes and De Schutter, 2012; Bartol et al., 2015).

| The synapse model, its timescales and mechanisms.

(a), Model diagram with the synaptic components including pre and postsynaptic compartments and inhibitory transmission (bottom left). AMPAr, NMDAr: AMPA- and NMDA-type glutamate receptors respectively; GABA(A)r: Type A GABA receptors; VGCC: R-, T- and L-type voltage-gated Ca2+ channels; SK: SK potassium channels. (b), Stochastic dynamics of the different ligand-gated and voltage-gated ion channels in the model. Plots show the total number of open channels as a function of time. The insets show a zoomed time axis highlighting the difference in timescale of the activity among the channels. (c), Dendritic spine membrane potential (left) and calcium concentration (right) as function of time for a single causal (1Pre1Post10) stimulus (EPSP: single excitatory postsynaptic potential, ‘1Pre’; BaP: single back-propagated action potential, ‘1Post’). (d), Left: depletion of vesicle pools (reserve and docked) induced by 30 pairing repetitions delivered at 5 Hz (Sterratt et al., 2011), see Materials and methods. The same depletion rule is applied to both glutamate- and GABA-containing vesicles. Right: BaP efficiency as function of time. BaP efficiency phenomenologically captures the distance-dependent attenuation of BaP (Buchanan and Mellor, 2007; Golding et al., 2001), see Materials and methods. (e), Concentration of active enzyme for CaM, CaN, and CaMKII, as function of time triggered by 30 repetitions of 1Pre1Post10 pairing stimulations delivered at 5 Hz. The vertical grey bar is the duration of the stimuli, 6 s. The multiple traces in the graphs in panels c (right) and e reflect the run-to-run variability due to the inherent stochasticity in the model.

The synapse model showed nonlinear dynamics across multiple timescales. For illustration, we stimulated the synapse with single simultaneous glutamate and GABA vesicle releases (Figure 1b). AMPArs and VGCCs open rapidly but close again within a few milliseconds. The dendritic GABAr closes more slowly, on a timescale of ∼10 ms. NMDArs, the major calcium source, closes on timescales of ∼50 and ∼250 ms for the GluN2A and GluN2B subtypes, respectively.

To show the typical responses of the spine head voltage and Ca2+, we stimulated the synapse with a single presynaptic pulse (EPSP) paired 10 ms later with a single BaP (1Pre1Post10; Figure 1c left). For this pairing, the arrival of a BaP at the spine immediately after an EPSP, leads to a large Ca2+ transient aligned with the BaP due to the NMDArs first being bound by glutamate then unblocked by the BaP depolarisation (Figure 1c right).

Single pre or postsynaptic stimulation pulses did not cause depletion of vesicle reserves or substantial activation of the enzymes. To illustrate these slower-timescale processes, we stimulated the synapse with a prolonged protocol: one presynaptic pulse followed by one postsynaptic pulse 10 ms later, repeated 30 times at 5 Hz (Figure 1d–e). The number of vesicles in both the docked and reserve pools decreased substantially over the course of the stimulation train (Figure 1d left), which in turn caused decreased vesicle release probability. Similarly, by the 30th pulse, the dendritic BaP amplitude had attenuated to ∼85% (∼70% BaP efficiency; Figure 1d right) of its initial amplitude, modelling the effects of slow dendritic sodium channel inactivation (Colbert et al., 1997; Golding et al., 2001). Free CaM concentration rose rapidly in response to calcium transients but also decayed back to baseline on a timescale of ∼500 ms (Figure 1e top). In contrast, the concentration of active CaMKII and CaN accumulated over a timescale of seconds, reaching a sustained peak during the stimulation train, then decayed back to baseline on a timescale of ∼10 and ∼120 s respectively, in line with experimental data (Quintana et al., 2005; Fujii et al., 2013; Chang et al., 2017; Figure 1e).

The effects of the stochastic variables can be seen in Figure 1b–d. The synaptic receptors and ion channels open and close randomly (Figure 1b). Even though spine voltage, calcium, and downstream molecules were modelled as continuous and deterministic, they inherited some randomness from the upstream stochastic variables. As a result, there was substantial trial-to-trial variability in the voltage and calcium responses to identical pre and postsynaptic spike trains (grey traces in Figure 1c). This variability was also passed on to the downstream enzymes CaM, CaMKII and CaN, but was filtered and therefore attenuated by the slow dynamics of CaMKII and CaN. In summary, the model contained stochastic nonlinear variables acting over five different orders of magnitude of timescale, from ∼1 ms to ∼1 min, making it sensitive to both fast and slow components of input signals.

Distinguishing between stimulation protocols using the CaMKII and CaN joint response

It has proven difficult for simple models of synaptic plasticity to capture the underlying rules and explain why some stimulation protocols induce plasticity while others do not. We tested the model’s sensitivity by simulating its response to a set of protocols used by Tigaret et al., 2016 in a recent ex vivo experimental study on adult (P50-55) rat hippocampus with blocked GABAr. We schematized the Tigaret et al., 2016 protocols in Figure 2a. Notably, three leading spike-timing and calcium-dependent plasticity models (Song et al., 2000; Pfister and Gerstner, 2006; Graupner and Brunel, 2012) could not fit well these data (Figure 2b–d). Next, we asked if our new model could distinguish between three pairs of protocols (see Figure 2e–m). For each of these pairs, one of the protocols experimentally induced LTP or LTD, while the other subtly different protocol caused no change (NC) in synapse strength.

| The duration and amplitude of the joint CaN-CaMKII activity differentiates plasticity protocols.

(a), Tigaret et al., 2016 protocols, which inspired this model.(a) is adapted from Figure 2B from Tigaret et al., 2016. (b–d), Standard models for predicting plasticity fail to account for Tigaret et al., 2016 data. Mean weight change for the Tigaret’s data (red), error bars denote ±1 s.d. Plasticity protocols indicated by labels on x-axis. Blue bars show mean plasticity predicted for the same protocols by classic STDP model (Song et al., 2000) (panel b), triplet STDP model (Pfister and Gerstner, 2006) (panel c), or Graupner-Brunel calcium-based STDP (Graupner and Brunel, 2012) model (panel d). (e), Time-course of active enzyme concentration for CaMKII (solid line) and CaN (dashed line) triggered by two protocols consisting of 300 repetitions at 5 Hz of 1Pre2Post10 or 1Pre1Post10 stimulus pairings. Protocols start at time 0 s. Experimental data indicates that 1Pre2Post10 and 1Pre1Post10 produce LTP and no change (NC), respectively. (f), Trajectories of joint enzymatic activity (CaN-CaMKII) as function of time for the protocols in panel e, starting at the initial resting state (filled black circle). The arrows show the direction of the trajectory and filled grey circles indicate the time points at 2, 10, and 60 s after the beginning of the protocol. The region of the CaN-CaMKII plane enclosed in the black square is expanded in panel g. (g), Mean-time (colorbar) spent by the orbits in the CaN-CaMKII plane region expanded from panel f for each protocol (average of 100 samples). For panels g, j and m the heat maps were based on enzyme and 2Post1Pre50 (NC) depicted in the same manner as in panels (e-g). (k-m), CaN-CaMKII activities for the LTD-inducing protocol 2Pre50 (900 repetitions at 3 Hz) and the NC protocol 2Pre10 (300 repetitions at 5 Hz) depicted in the same manner as in panels e-g.

The first pair of protocols differed in intensity. A protocol which caused no plasticity consisted of 1 presynaptic spike followed 10 ms later by one postsynaptic spike repeated at 5 Hz for 1 min (1Pre1Post10, 300 at 5 Hz). The other protocol induced LTP, but differed only in that it included a postsynaptic doublet instead of a single spike (1Pre2Post10, 300 at 5 Hz), implying a slightly stronger initial BaP amplitude. We first attempted to achieve separability by plotting CaMKII or CaN activities independently. As observed in the plots in Figure 2e, it was not possible to set a single concentration threshold on either CaMKII or CaN that would discriminate between the protocols. This result was expected, at least for CaMKII, as recent experimental data demonstrates a fast saturation of CaMKII concentration in dendritic spines regardless of stimulation frequency (Chang et al., 2017).

To achieve better separability we set out to test a different approach, which was to combine the activity of the two enzymes, by plotting the joint CaMKII and CaN responses against each other on a 2D plane (Figure 2f). This innovative geometric plot is based on the mathematical concept of orbits from dynamical systems theory (Meiss, 2007). In this plot, the trajectories of two protocols can be seen to overlap for the initial part of the transient and then diverge. To quantify trial-to-trial variability, we also calculated contour maps showing the mean fraction of time the trajectories spent in each part of the plane during the stimulation (Figure 2g). Importantly, both the trajectories and contour maps were substantially non-overlapping between the two protocols, implying that they can be separated based on the joint CaN-CaMKII activity. We found that the 1Pre2Post10 protocol leads to a weaker response in both CaMKII and CaN, corresponding to the lower blue traces in Figure 2f. The decreased response to the doublet protocol was due to the stronger attenuation of dendritic BaP amplitude over the course of the simulation (Golding et al., 2001), leading to reduced calcium influx through NMDArs and VGCCs (data not shown).

Using the second pair of protocols, we explored if this combined enzyme activity analysis could distinguish between subtle differences in protocol sequencing. We stimulated our model with one causal paring protocol (EPSP-BaP) involving a single presynaptic spike followed 50 ms later by a doublet of postsynaptic spikes (1Pre2Post50, 300 at 5 Hz), repeated at 5 Hz for one minute, which caused LTP in Tigaret et al., 2016. The other, anticausal, protocol involved the same total number of pre and postsynaptic spikes, but with the pre-post order reversed (2Post1Pre50, 300 at 5 Hz). Experimentally, the anticausal (BaP-EPSP) protocol did not induce plasticity (Tigaret et al., 2016). Notably, the only difference was the sequencing of whether the pre or postsynaptic neuron fired first, over a short time gap of 50 ms. Although the time courses of CaMKII and CaN activities were difficult to distinguish (Figure 2h), the LTP-inducing protocol caused greater CaN activation, compared to the non LTP-inducing protocol. Indeed, this translated to a horizontal offset in both the trajectory and contour map (Figure 2i–j), demonstrating that this pair of protocols can also be separated in the joint CaN-CaMKII plane.

The third pair of protocols differed in both duration and intensity. We thus tested the combined enzyme activity analysis in this configuration. In line with a previous study (Isaac et al., 2009), Tigaret et al., 2016 found that a train of doublets of presynaptic spikes separated by 50 ms repeated at a low frequency of 3 Hz for 5 min (2Pre50, 900 at 3 Hz) induced LTD, while a slightly more intense but shorter duration protocol of presynaptic spike doublets separated by 10 ms repeated at 5 Hz for 1 min (2Pre10, 300 at 5 Hz) did not cause plasticity. When we simulated both protocols in the model (Figure 2k–m), both caused similar initial responses in CaMKII and CaN. In the shorter protocol, this activation decayed to baseline within 100 s of the end of the stimulation. However the slower and longer-lasting 2Pre50 3 Hz 900 p protocol caused an additional sustained, stochastically fluctuating, plateau of activation of both enzymes (Figure 2k). This resulted in the LTD-inducing protocol having a downward and leftward-shifted CaN-CaMKII trajectory and contour plot, relative to the other protocol (Figure 2l–m). These results again showed that the joint CaN-CaMKII activity can be used to predict plasticity changes.

A geometrical readout mapping joint enzymatic activity to plasticity outcomes

The three above examples demonstrated that plotting the combined CaN-CaMKII activities in a 2D plane (geometrical readout which is abstract e.g. not defined within a physical space) allowed us to distinguish between subtly different protocols with correct assignment of plasticity outcome. We found that the simulated CaN-CaMKII trajectories from the two LTP-inducing protocols (Figure 2e–g and Figure 2h–j) spent a large fraction of time near ∼20 µM CaMKII and 7–10 µM CaN. In contrast, protocols that failed to trigger LTP had either lower (Figure 2h–j, k–m), or higher CaMKII and CaN activation (1Pre1Post10, Figure 2e–g). The LTD-inducing protocol, by comparison, spent a longer period in a region of sustained but lower ∼12 µM CaMKII and ∼2 µM CaN and activation Figure 2k–m. The plots in Figure 2g, j and m, show contour maps of histograms of the joint CaMKII-CaN activity, indicating where in the plane the trajectories spent most time. Figure 2g and j indicate that this measure can be used to predict plasticity, because the NC and LTP protocol histograms are largely non-overlapping. In Figure 2g, the NC protocol response ‘overshoots’ (mostly due to higher CaN concentration) the LTP protocol response, whereas in Figure 2j the NC protocol response ‘undershoots’ (mostly due to lower CaN concentration) the LTP protocol response. In contrast, when we compared the response histograms for the LTD and NC protocols, we found a greater overlap (Figure 2m). This suggested that, in this case, the histogram alone was not sufficient to separate the protocols, and that protocol duration is also important. LTD induction (2Pre50) required a more prolonged activation than NC (2Pre10). We thus took advantage of these joint CaMKII-CaN activity maps to design a minimal readout mechanism connecting combined enzyme activity to LTP, LTD or NC. We reasoned that this readout would need three key properties. First, although the figure suggests that the CaMKII-CaN trajectories corresponding to LTP and LTD could be linearly separable, we will demonstrate later (see Figure 3—figure supplement 3) that the readout requires nonlinear boundaries to activate the plasticity inducing component. Second, since LTD requires more prolonged activity than LTP, the readout should be sensitive to the timescale of the input. Third, a mechanism is required to convert the 2D LTP-LTD inducing signals into a synaptic weight change. After iterating through several designs, we satisfied the first property by designing ‘plasticity regions’: polygons in the CaN-CaMKII plane that would detect when trajectories pass through. We satisfied the second property by using two plasticity inducing components with different time constants which low-pass-filter the plasticity region signals. We satisfied the third property by feeding both the opposing LTP and LTD signals into a stochastic Markov chain which accumulated the total synaptic strength change. Overall, this readout mechanism acts as a parsimonious model of the complex signalling cascade linking CaMKII and CaN activation to expression of synaptic plasticity (He et al., 2015). It can be considered as a two-dimensional extension of previous computational studies that applied analogous 1D threshold functions to dendritic spine calcium concentration (Shouval et al., 2002; Karmarkar and Buonomano, 2002; Graupner and Brunel, 2012; Standage et al., 2014).

We now elaborate on the readout design process (see also Figure 21 of Materials and methods). We first drew non-overlapping polygons of LTP and LTD ‘plasticity regions’ in the CaN-CaMKII plane (Figure 3a). We positioned these regions over the parts of the phase space where the enzyme activities corresponding to the LTP- and LTD-inducing protocols were most different (Materials and methods), as shown by trajectories in Figure 2f, i and l. When a trajectory enters in one of these plasticity regions, it activates LTD or LTP indicator variables (Materials and methods) which encode the joint enzyme activities (trajectories in the phase plots) transitions across the LTP and LTD regions over time (Figure 3b). These indicator variables drove transition rates of a plasticity Markov chain used to predict LTP or LTD (Figure 3c), see Materials and methods. Intuitively, this plasticity Markov chain models the competing processes of insertion/deletion of AMPArs to the synapse, although this is not represented in the model. The LTD transition rates were slower than the LTP transition rates, to reflect studies showing that LTD requires sustained synaptic stimulation (Yang et al., 1999; Mizuno et al., 2001; Wang et al., 2005). The parameters for this plasticity Markov chain (Materials and methods) were fit to the plasticity induction outcomes from different protocols (Appendix 1—table 1). At the beginning of the simulation, the plasticity Markov chain starts with 100 processes (Destexhe et al., 1998) in the NC state, with each variable representing 1% weight change, an abstract measure of synaptic strength that can be either EPSP, EPSC, or field EPSP slope depending on the experiment. Each process can transit stochastically between NC, LTP and LTD states. At the end of the protocol, the plasticity outcome is given by the difference between the number of processes in the LTP and the LTD states (Materials and methods).

Figure 3 with 4 supplements see all
Read-out strategy to accurately model Tigaret et al., 2016 experiment.

(a) Illustration of the joint CaMKII and CaN activities crossing the plasticity regions. Arrows indicate the flow of time, starting at the filled black circle. (b) Region indicator showing when the joint CaN and CaMKII activity crosses the LTD or LTP regions in panel a. For example, the LTP indicator is such that 1LTP(x)=1 if xLTP and 0 otherwise. Leaving the region activates a mechanism with a slow timescale that keeps track of the accumulated time inside the region. Such mechanism drives the transition rates used to predict plasticity (Materials and methods). (c), Plasticity Markov chain with three states: LTD, LTP and NC. There are only two transition rates which are functions of the plasticity region indicator (Materials and methods). The LTP transition is fast, whereas the LTD transition is slow, meaning that LTD change requires longer time inside the LTD region (panel a). The NC state starts with 100 processes. See Figure 23 for more details on the dynamics of the Plasticity Markov Chain. (d) Joint CaMKII and CaN activity for all protocols in Tigaret et al., 2016 (shown in panel f). The stimulus ends when the trajectory becomes smooth. Trajectories correspond to those in Figure 2b, e and h, at 60 s. (e) Region indicator for the protocols in panel f. The upper square bumps are caused by the protocol crossing the LTP region, the lower square bumps when the protocol crosses the LTD region (as in panel d). (f) Synaptic weight (%) as function of time for each protocol. The weight change is defined as the number (out of 100) of states in the LTP state minus the number of states in the LTD state (panel c). The trajectories correspond to the median of the simulations in panel g. (g) Synaptic weight change (%) predicted by the model compared to data (EPSC amplitudes) from Tigaret et al., 2016 (100 samples for each protocol, also for panel h and i). The data (filled grey circles) was provided by Tigaret et al., 2016 (note an 230% outlier as the red asterisk). (h) Predicted mean synaptic weight change (%) as a function of delay (ms) and number of pairing repetitions (pulses) for the protocol 1Pre2Post(delay), where delays are between –100 and 100 ms. LTD is induced by 2Post1Pre50 after at least 500 pulses. The mean weight change along the horizontal dashed line is reported in the STDP curves in panel i. (i) Synaptic weight change (%) as a function of pre-post delay. Each plot corresponds to a different pairing repetition number (color legend). The solid line shows the mean, and the ribbons are the 2nd and 4th quantiles. The filled grey circles are the data means estimated in Tigaret et al., 2016, also shown in panel g.

In Figure 3d, we plot the model’s responses to seven different plasticity protocols used by Tigaret et al., 2016 by overlaying example CaMKII-CaN trajectories for each protocol with the LTP and LTD regions. The corresponding region indicators are plotted as function of time in Figure 3e, and long-term alterations in the synaptic strength are plotted as function of time in Figure 3f. The three protocols that induced LTP in the Tigaret et al., 2016 experiments spent substantial time in the LTP region, and so triggered potentiation. In contrast, the combined response (CamKII, CaN) to 1Pre1Post10 overshoots both regions, crossing them only briefly on its return to baseline, and so resulted in little weight change. The protocol that induced LTD (2Pre50, purple trace) is five times longer than other protocols, spending sufficient time inside the LTD region (Figure 3f). In contrast, two other protocols that spent time in the same LTD region of the CaN-CaMKII plane (2Post1Pre50 and 2Pre10) were too brief to induce LTD. These protocols were also not strong enough to reach the LTP region, so resulted in no net plasticity, again consistent with Tigaret et al., 2016 experiments.

We observed run-to-run variability in the amplitude of the predicted plasticity, due to the inherent stochasticity in the model. To ensure that stochastic components are necessary for adequate model behaviour, we compared stochastic and deterministic versions of the model with and without discrete presynaptic release and found that adding stochastic components indeed modified the model’s behaviour (Figure 3—figure supplement 1). Also, we confirmed that VGCCs are necessary for accurate modelling of Tigaret et al., 2016 data as blocking these channels reproduced the data obtained in VGCC blockers by Tigaret that is no potentiation could be elicited (Figure 3—figure supplement 2). Finally, we stress in Figure 3—figure supplement 3 that the horizontal boundaries (related to CaMKII activity) are indeed necessary.

In Figure 3g, we plot the distributions of the simulation outcomes, along with the experimental data, for the protocols in Tigaret et al., 2016. We find a very good correspondence between the model and experiments. Of note, data fitting of the experiments in Tigaret et al., 2016 (Figure 3g) was more accurate with our model than the fitting obtained with existing leading spike-or calcium-based STDP models (Song et al., 2000; Pfister and Gerstner, 2006; Graupner and Brunel, 2012), see Figure 2b–d.

Experimentally, LTP can be induced by few pulses while LTD usually requires stimulation protocols of longer duration (Yang et al., 1999; Mizuno et al., 2001; Wang et al., 2005). We incorporated this effect into the geometrical readout model by letting LTP have faster transition rates than LTD (Figure 3c). Tigaret et al., 2016 found that 300 repetitions of anticausal post-before-pre pairings did not cause LTD, in contrast to the canonical spike-timing-dependent plasticity curve (Bi and Poo, 1998). We hypothesized that LTD might indeed appear with the anticausal protocol (Appendix 1—table 1) if stimulation duration was increased. To explore this possibility in our model, we systematically Alabi and Tsien, 2012 varied the number of paired repetitions from 100 to 1200, and also co-varied the pre-post delay from –100 to 100 ms. Figure 3h shows a contour plot of the predicted mean synaptic strength change across for the 1Pre2Post(delay) stimulation protocol for different numbers of pairing repetitions. In Figure 3h and a LTD window appears after ∼500 pairing repetitions for some anticausal pairings, in line with our expectation. The magnitude of LTP also increases with pulse number, for causal positive pairings. For either 100 or 300 pairing repetitions, only LTP or NC is induced (Figure 3i). The model also made other plasticity predictions by varying Tigaret et al., 2016 experimental conditions (Figure 3—figure supplement 4). In summary, our geometrical readout mechanism suggests that the direction and magnitude of the change in synaptic strength can be predicted from the joint CaMKII-CaN activity in the LTP and LTD regions.

Frequency-dependent plasticity

The stimulation protocols used by Tigaret et al., 2016 explored how subtle variations in pre and postsynaptic spike timing influenced the direction and magnitude of plasticity (see Appendix 1—table 1 for experimental differences). In contrast, traditional synaptic plasticity protocols exploring the role of presynaptic stimulation frequency did not measure the timing of co-occurring postsynaptic spikes (Dudek and Bear, 1992; Wang and Wagner, 1999; Kealy and Commins, 2010). These studies found that long-duration low-frequency stimulation (LFS) induces LTD, whereas short-duration high-frequency stimulation induces LTP, with a cross-over point of zero change at intermediate stimulation frequencies. In addition to allowing us to explore frequency-dependent plasticity (FDP), this stimulation paradigm also gave us further constraints to define the LTD polygon region in the model since in Tigaret et al., 2016, only one LTD case was available. For FDP, we focused on modelling the experiments from Dudek and Bear, 1992, who stimulated Schaffer collateral projections to pyramidal CA1 neurons with 900 pulses in frequencies ranging from 1 to 50 Hz. In addition to presynaptic stimulation patterns, the experimental conditions differed from Tigaret et al., 2016 in two other aspects: animal age and control of postsynaptic spiking activity (see Appendix 1—table 1 legend). We incorporated both age-dependence and EPSP-evoked-BaPs in our model (Materials and methods). Importantly, the geometrical readout mechanism mapping joint CaMKII-CaN activity to plasticity remained identical for all experiments in this work.

Figure 4a shows the joint CaMKII-CaN activity when we stimulated the model with 900 presynaptic spikes at 1, 3, 5, 10, and 50 Hz (Dudek and Bear, 1992). Higher stimulation frequencies drove stronger responses in both CaN and CaMKII activities (Figure 4a). Figure 4b and c show the corresponding plasticity region indicator for the LTP/LTD region threshold crossings and the synaptic strength change. From this set of five protocols, only the 50 Hz stimulation drove a response strong enough to reach the LTP region of the plane (Figure 4a and d). Although the remaining four protocols drove responses primarily in the LTD region, only the 3 and 5 Hz stimulations resulted in substantial LTD. The 1 and 10 Hz stimulations resulted in negligible LTD, but for two distinct reasons. Although the 10 Hz protocol’s joint CaMKII-CaN activity passed through the LTD region of the plane (Figure 4a and d), it was too brief to activate the slow LTD mechanism built into the readout (Materials and methods). The 1 Hz stimulation, on the other hand, was prolonged, but its response was too weak to reach the LTD region, crossing the threshold only intermittently (Figure 4b, bottom trace). Overall, the model matched well the mean plasticity response found by Dudek and Bear, 1992, see Figure 4e, following a classic BCM-like curve as function of stimulation frequency (Abraham et al., 2001; Bienenstock et al., 1982).

Figure 4 with 2 supplements see all
Frequency dependent plasticity (FDP), Dudek and Bear, 1992 dataset.

(a) Example traces of joint CaMKII-CaN activity for each of Dudek and Bear, 1992 protocol. (b) Region indicator showing when the joint CaMKII-CaN activity crosses the LTD or LTP regions for each protocol in panel a. (c) Synaptic weight change (%) as a function of time for each protocol, analogous to Figure 3c. Trace colours correspond to panel a. The trajectories displayed were chosen to match the medians in panel e. (d) Mean (100 samples) time spent (s) for protocols 1Pre for 900 pairing repetitions at 3, 10, and 50 Hz. (e), Comparison between data from Dudek and Bear, 1992 and our model (1Pre 900 p, 300 sampcomles per frequency, see Appendix 1—table 1). Data are represented as normal distributions with the mean and variance of the change in field EPSP slope taken from Dudek and Bear, 1992. (f), Prediction for the mean weight change (%) when varying the stimulation frequency and pulse number (24x38 × 100 data points, respectively pulse x frequency x samples). The filled grey circles show the Dudek and Bear, 1992 protocol parameters and the corresponding results are shown in panel e. In Figure 4—figure supplement 1, we provide additional graphs of frequency dependent plasticity outcomes, including predictions, when varying experimental parameters in Dudek and Bear, 1992 (external Mg, external Ca, distance from soma, temperature, Poisson spike train during development).

We then used the model to explore the stimulation space in more detail by varying the stimulation frequency from 0.5 to 40 Hz, and varying the number of presynaptic pulses from 50 to 1200. Figure 4f shows a contour map of the mean synaptic strength change (%) in this 2D frequency–pulse number space. Dudek and Bear, 1992 experimental conditions, we found that LTD induction required at least ∼300 pulses, at frequencies between 1 and 3 Hz. In contrast, LTP could be induced using ∼50 pulses at ∼20 Hz or greater. The contour map also showed that increasing the number of pulses (vertical axis in Figure 4f) increases the magnitude of both LTP and LTD. This was accompanied by a widening of the LTD frequency range, whereas the LTP frequency threshold remained around ∼20 Hz, independent of pulse number. This general effect, that increasing pulse number tends to increase the magnitude of plasticity, was also observed in simulation of Tigaret et al., 2016 (see Figure 3h). Ex vivo experiments in Dudek and Bear, 1992 were done at 35°C. However, lower temperatures are more widely usetd for ex vivo experiments because they extend brain slice viability.

At this point, having fully described the model, we show the importance of the stochasticity of the different components of the model. We simulated three protocols of Dudek and Bear, 1992 with deterministic equations in Figure 4—figure supplement 2. We show that, for the different protocols, if some of the channels are modelled with deterministic equations, the net effect on synapse weight differs from the expected outcome provided by the original model. The relative contributions of each source of noise differed, depending on the plasticity protocol. We can conclude that all noise sources we introduced in our model are important.

Variations in plasticity induction with developmental age

The rules for induction of LTP and LTD change during development (Dudek and Bear, 1993; Cao and Harris, 2012), so a given plasticity protocol can produce different outcomes when delivered to synapses from young animals versus mature animals. For example, when Dudek and Bear, 1993 tested the effects of low-frequency stimulation (1 Hz) on CA3-CA1 synapses from rats of different ages, they found that the magnitude of LTD decreases steeply with age from P7 until becoming minimal in mature animals >P35 (Figure 5a, circles). Across the same age range, they found that a theta burst stimulation (TBS) protocol induced progressively greater LTP magnitude with developmental age (Figure 5b, circles). Multiple properties of neurons change during development: the NMDAr switches its dominant subunit expression from GluN2B to GluN2A (Sheng et al., 1994; Popescu et al., 2004; Iacobucci and Popescu, 2017), the reversal potential of the receptor (GABAr) switches from depolarising to hyperpolarizing (Rivera et al., 1999; Meredith et al., 2003; Rinetti-Vargas et al., 2017), and the action potential backpropagates more efficiently with age (Buchanan and Mellor, 2007). These mechanisms have been proposed to underlie the developmental changes in synaptic plasticity rules because they are key regulators of synaptic calcium signalling (Meredith et al., 2003; Buchanan and Mellor, 2007). However, their sufficiency and individual contributions to the age-related plasticity changes are unclear and have not been taken into account in any previous model. We incorporated these mechanisms in the model (Materials and methods) by parametrizing each of the three components to vary with the animal’s postnatal age, to test if they could account for the age-dependent plasticity data.

Figure 5 with 1 supplement see all
Age-dependent plasticity, Dudek and Bear, 1993 dataset.

(a), Synaptic weight change for 1Pre, 900 at 1 Dudek and Bear, 1993. The solid line is the mean and the ribbons are the 2nd and 4th quantiles predicted by our model (same for panel b, c and f). (b), Synaptic weight change for theta burst stimulation (TBS - 4Pre at 100 Hz repeated 10 times at 5 Hz given in 6 epochs at 0.1 Hz, see Appendix 1—table 1). (c), Synaptic weight change as a function of frequency for different ages. BCM-like curves showing that, during adulthood, the same LTD protocol becomes less efficient. It also shows that high-frequencies are inefficient at inducing LTP before P15. (d), Synaptic weight change as a function of age. Proposed protocol using presynaptic bursts to recover LTD at ≥ P35 with less pulses, 300 instead of the original 900 from Dudek and Bear, 1993. This effect is more pronounced for young rats. Figure 5—figure supplement 1 shows a 900 pulses comparison. (e), Mean synaptic strength change (%) as a function of frequency and age for 1Pre 900 pulses (32x38 × 100, respectively, for frequency, age and samples). The protocols in Dudek and Bear, 1993 (panel a) are marked with the yellow vertical line. The horizontal lines represent the experimental conditions of panel c. Note the P35 was used for Dudek and Bear, 1992 experiment in Figure 4f. (f), Mean time spent for the 1Pre 1 Hz 900 pulses protocol showing how the trajectories are left-shifted as rat age increases. In Figure 5—figure supplement 1, we provide additional simulations to analyse the synaptic plasticity outcomes, including predictions, of duplets, triplets and quadruplets for FDP, perturbing developmental-mechanisms for Dudek and Bear, 1993, and age-related changes in STDP experiments (Inglebert et al., 2020; Tigaret et al., 2016; Meredith et al., 2003).

We found that elaborating the model with age-dependent changes in NMDAr composition, GABAr reversal potential, and BaP efficiency, while keeping the same plasticity readout parameters, was sufficient to account for the developmental changes in LTD and LTP observed by Dudek and Bear, 1993 (Figure 5a and b). We then explored the model’s response to protocols of various stimulation frequencies, from 0.5 to 40 Hz, across ages from P5 to P80 (Figure 5c and e). Figure 5c shows the synaptic strength change as function of stimulation frequency for ages P15, P25, P35, and P45. The magnitude of LTD decreases with age, while the magnitude of LTP increases with age. Figure 5e shows a contour plot of the same result, covering the age-frequency space.

The 1 Hz presynaptic stimulation protocol in Dudek and Bear, 1993 did not induce LTD in adult animals (Dudek and Bear, 1992). We found that the joint CaN-CaMKII activity trajectories for this stimulation protocol underwent an age-dependent leftward shift beyond the LTD region (Figure 5f). This implies that LTD is not induced in mature animals by this conventional LFS protocol due to insufficient activation of enzymes. In contrast, Tigaret et al., 2016 and Isaac et al., 2009 were able to induce LTD in adult rat tissue by combining LFS with presynaptic spike pairs repeated 900 times at 3 Hz. Given these empirical findings and our modelling results, we observe that LTD induction in adult animals requires that the stimulation protocol: (1) causes CaMKII and CaN activity to stay more in the LTD region than the LTP region and (2) is sufficiently long to activate the LTD readout mechanism. With experimental parameters used by Dudek and Bear, 1993, this may be as short as 300 pulses when multi-spike presynaptic protocols are used since the joint CaMKII-CaN activity can reach the LTD region more quickly than for single spike protocols. We simulated two such potential protocols as predictions: doublet and quadruplet spike groups delivered 300 times at 1 Hz, with 50 ms between each pair of spikes in the group (Figure 5d). The model predicts that both these protocols induce LTD in adults, whereas as shown above, the single pulse protocol did not cause LTD. These simulations suggest that the temporal requirements for inducing LTD may not be as prolonged as previously assumed, since they can be reduced by varying stimulation intensity. See Figure 5—figure supplement 1 for frequency versus age maps for presynaptic bursts.

Dudek and Bear, 1993 also performed theta burst stimulation (Appendix 1—table 1) at different developmental ages, and found that LTP is not easily induced in young rats (Cao and Harris, 2012), as depicted in Figure 5b. The model qualitatively matches this trend, and also predicts that TBS induces maximal LTP around P21, before declining further during development (Figure 5b, green curve). Similarly, the model predicts that high-frequency stimulation induces LTP only for ages >P15, peaks at P35, then gradually declines at older ages (Figure 5e). Note that in Figure 5b, we used six epochs instead of four used by Dudek and Bear, 1993 to increase LTP outcome which is known to washout after one hour for young rats (Cao and Harris, 2012).

In contrast to Dudek and Bear, 1993 findings, other studies have found that LTP can be induced in hippocampus in young animals (<P15) with STDP. For example, Meredith et al., 2003 found that, at room temperature, 1Pre1Post10 induces LTP in young rats, whereas 1Pre2Post10 induces NC. This relationship was inverted for adults, with 1Pre1Post inducing no plasticity and 1Pre2Post10 inducing LTP (as captured by our model in Figure 5—figure supplement 1).

Together, these results suggest that not only do the requirements for LTP/LTD change with age, but also that these age-dependencies are different for different stimulation patterns. Finally, we explore which mechanisms are responsible for plasticity induction changes across development in the FDP protocol (Figure 5—figure supplement 1) by fixing each parameter to young or adult values for the FDP paradigm. Our model analysis suggests that the NMDAr switch (Iacobucci and Popescu, 2017) is a dominant factor affecting LTD induction, but the maturation of BaP (Buchanan and Mellor, 2007) is the dominant factor affecting LTP induction, with GABAr shift having only a weak influence on LTD induction for Dudek and Bear, 1993 FDP.

Plasticity requirements during development do not necessarily follow the profile in Dudek and Bear, 1993 as shown by Meredith et al., 2003 STDP experiment. Our model suggests that multiple developmental profiles are possible when experimental conditions vary within the same stimulation paradigm. This is illustrated in Figure 6—figure supplement 2a–c by varying the age of STDP experiments done in different conditions. We fitted well the data from Wittenberg and Wang, 2006 by adapting the model with appropriate age and temperature.

Effects of extracellular calcium and magnesium concentrations on plasticity outcome

The canonical STDP rule (Bi and Poo, 1998), measured in cultured neurons with high extracellular calcium ([Ca2+]o) and at room temperature, was recently found not to be reproducible at physiological [Ca2+]o in CA1 brain slices (Inglebert et al., 2020). Instead, by varying the [Ca2+]o and [Mg2+]o, Inglebert et al., 2020 found a spectrum of STDP rules with either no plasticity or full-LTD for physiological [Ca2+]o conditions ([Ca2+]o< 1.8 mM) and a bidirectional rule for high [Ca2+]o ([Ca2+]o> 2.5 mM), shown in Figure 6a-c.

Figure 6 with 2 supplements see all
Effects of extracellular calcium and magnesium concentrations on plasticity.

(a), Synaptic weight (%) for a STDP rule with [Ca2+]o=1.3 mM (fixed ratio, Ca/Mg = 1.5). According to the data extracted from Inglebert et al., 2020, the number of pairing repetitions for causal/positive (anti-causal/negative) delays is 100 (150), both delivered at 0.3 Hz. The solid line is the mean, and the ribbons are the 2nd and 4th quantiles predicted by our model (all panels use 100 samples). (b), Same as a, but for [Ca2+]o = 1.8 mM (Ca/Mg ratio = 1.5). (c), Same as a, but for [Ca2+]o = 3 mM (Ca/Mg ratio = 1.5). (d), Mean time spent for causal pairing, 1Pre1Post10, at different Ca/Mg concentration ratios. The contour plots are associated with the panels a, b and c. e, Predicted effects of extracellular Ca/Mg on STDP outcome. Synaptic weight change (%) for causal (1Pre1Post10, 100 at 0.3 Hz) and anticausal (1Post1Pre10, 150 at 0.3 Hz) pairings varying extracellular Ca from 1.0 to 3 mM (Ca/Mg ratio = 1.5). The dashed lines represent the experiments in the panel a, b and c. We used 21x22 × 100 data points, respectively calcium x delay x samples. (f), Predicted effects of varying frequency and extracellular Ca/Mg for an STDP protocol. Contour plot showing the mean synaptic weight (%) for a single causal pairing protocol (1Pre1Post10, 100 samples) varying frequency from 0.1 to 10 Hz and [Ca2+]o from 1.0 to 3 mM (Ca/Mg ratio = 1.5). We used 21 x 18 × 100 data points, respectively calcium x frequency x samples.

We attempted to reproduce Inglebert et al., 2020 findings by varying [Ca2+]o and [Mg2+]o with the following consequences for the model mechanisms (Materials and methods). On the presynaptic side, [Ca2+]o modulates vesicle release probability. On the postsynaptic side, high [Ca2+]o reduces NMDAr conductance (Maki and Popescu, 2014), whereas [Mg2+]o affects the NMDAr Mg2+ block (Jahr and Stevens, 1990). Furthermore, spine calcium influx activates SK channels, which hyperpolarize the membrane and indirectly modulate NMDAr activity (Ngo-Anh et al., 2005; Griffith et al., 2016).

Figure 6a–c compares our model to Inglebert et al., 2020 STDP data at different [Ca2+]o and [Mg2+]o. Note that Inglebert et al., 2020 used 150 pairing repetitions for the anti-causal stimuli and 100 pairing repetitions for the causal stimuli both delivered at 0.3 Hz. At [Ca2+]o=1.3 mM, Figure 6a shows that the STDP rule induced weak LTD for brief causal delays. At [Ca2+]o = 1.8 mM, in Figure 6b, the model predicted a full-LTD window. At [Ca2+]o = 3 mM, in Figure 6c, it predicts a bidirectional rule with a second LTD window for long causal pairings, previously theorized by Rubin et al., 2005.

Figure 6d illustrates the time spent by the joint CaN-CaMKII activity for 1Pre1Post10 using Inglebert et al., 2020 experimental conditions. Each density plot corresponds to a specific Ca/Mg ratio as in Figure 6a–c. The response under low [Ca2+]o spent most time inside the LTD region, but high [Ca2+]o shifts the trajectory to the LTP region. Figure 6—figure supplement 1a presents density plots for the anti-causal protocols.

Inglebert et al., 2020 fixed the Ca/Mg ratio at 1.5, although aCSF formulations in the literature differ (see Appendix 1—table 1). Figure 6—figure supplement 1d shows that varying the Ca/Mg ratio and [Ca2+]o for Inglebert et al., 2020 experiments restrict LTP to Ca/Mg >1.5 and [Ca2+]o>1.8 mM.

Figure 6e shows a map of plasticity as function of pre-post delay and Ca/Mg concentrations and the parameters where LTP is induced for the 1Pre1Post10 protocol. Since plasticity rises steeply at around [Ca2+]o = 2.2 mM (see Figure 6—figure supplement 1b), small fluctuations in [Ca2+]o near this boundary could cause qualitative transitions in plasticity outcomes. For anti-causal pairings, increasing [Ca2+]o increases the magnitude of LTD (Figure 6—figure supplement 1b illustrates this with Inglebert et al., 2020 data). Our model can identify the transitions between LTD and LTP as a function of the ratio between [Ca2+]o and [Mg2+]o, see Figure 6—figure supplement 1.

Inglebert et al., 2020 also found that increasing the pairing frequency to 5 or 10 Hz results in a transition from LTD to LTP for 1Pre1Post10 at [Ca2+]o = 1.8 mM (Figure 6—figure supplement 1c), similar frequency-STDP behaviour has been reported in the cortex (Sjöström et al., 2001). In Figure 6f, we varied both the pairing frequencies and [Ca2+]o and we observe similar transitions to Inglebert et al., 2020. However, the model’s transition for [Ca2+]o = 1.8 mM was centred around 0.5 Hz, which was untested by Inglebert et al., 2020. The model predicts no plasticity at higher frequencies, unlike the data, that shows scattered LTP and LTD (see Figure 6—figure supplement 1c). Another frequency dependent comparison, Figure 3—figure supplement 4c and Figure 6—figure supplement 1h, show that Tigaret et al., 2016 burst-STDP and Inglebert et al., 2020 STDP share a similar transition structure, different from Dudek and Bear, 1992 FDP.

In contrast to Inglebert et al., 2020 results, the model predicts that setting low [Ca2+]o for Tigaret et al., 2016 burst-STDP abolishes LTP, and does not induce strong LTD (Figure 3—figure supplement 4d). For Dudek and Bear, 1992 experiment, Figure 4—figure supplement 1a [Mg2+]o controls a sliding threshold between LTD and LTP but not [Ca2+]o (Figure 4—figure supplement 1b). For another direct stimulation experiment, Figure 6—figure supplement 1f shows that in an Mg-free medium, LTP expression requires fewer pulses (Mizuno et al., 2001).

Despite exploring physiological [Ca2+]o and [Mg2+]o Inglebert (Inglebert et al., 2020) use a non-physiological temperature (30C) which extends T-type VGCC closing times and modifies the CaN-CaMKII baseline (Figure 6—figure supplement 2f). In summary, our model predicts that temperature can change STDP rules in a similar fashion to [Ca2+]o (Figure 6—figure supplement 1a and b). Overall, we confirm that plasticity is highly sensitive to variations in extracellular calcium, magnesium, and temperature (Figure 3—figure supplement 4a, Figure 6—figure supplement 2d–f).

In vivo-like spike variability affects plasticity

In the above sections, we used highly regular and stereotypical stimulation protocols to replicate typical ex vivo plasticity experiments. In contrast, neural spiking in hippocampus in vivo is irregular and variable (Fenton and Muller, 1998; Isaac et al., 2009). Previous studies that asked how natural firing variability affects the rules of plasticity induction used simpler synapse models (Rackham et al., 2010; Graupner et al., 2016; Cui et al., 2018). We explored this question in our synapse model using simulations with three distinct types of additional variability: (1) spike time jitter, (2) failures induced by dropping spikes, (3) independent pre and postsynaptic Poisson spike trains (Graupner et al., 2016).

We introduced spike timing jitter by adding zero-mean Gaussian noise (s.d. σ) to pre and postsynaptic spikes, changing spike pairs inter-stimulus interval (ISI). In Figure 7a, we plot the LTP magnitude as function of jitter magnitude (controlled by σ) for protocols taken from Tigaret et al., 2016. With no jitter, σ=0, these protocols have different LTP magnitudes (corresponding to Figure 3) and become similar once σ increases. The three protocols with a postsynaptic spike doublet gave identical plasticity for σ=50 ms.

Jitter and spike dropping effects on STDP and Poisson spike trains.

(a) Mean weight (%) for the jittered STDP protocols (protocol color legend shown in b). The solid line is the mean, and the ribbons are the 2nd and 4th quantiles predicted by our model using 100 samples (same for panels a, b and g). (b) Mean weight (%) for the same (Tigaret et al., 2016) protocols used in panel a subjected to random spike removal (sparsity %). (c) Effect of jitttering on Mean time (s) spent by joint enzymes trajectories in LTP/LTD regions. Contour plot shows 2Post1Pre50 and 1Pre1Post10 (300 at 5 Hz) without (grey contour plot) and with jittering (coloured contour plot). The circles and squares correspond to the marks in panel a. (d) Effect of sparsity on Mean time (s) spent by joint enzymes trajectories in LTP/LTD regions. Contour plot in grey showing 0% sparsity for 2Post1Pre50 300 at 5 Hz (see Figure 2j). The contour plots show the protocol with spike removal sparsities at 0% (NC), 30% (LTP), and 80% (NC). The triangles correspond to the same marks in panel a. (e) Distribution of the 50 ms jittering applied to the causal protocol 1Pre1Post10, 300 at 5 Hz in which nearly half of the pairs turned into anticausal. The mean frequency is 5 ±13.5 Hz making it to have a similar firing structure and position in the LTP region. The similar occurs for 2Post1Pre50 (panel c). (f) Mean weight change (%) combining both jittering (panel a) and sparsity (panel b) for 2Post1Pre50, 300 at 5 Hz. (g) Mean weight change (%) of pre and postsynaptic Poisson spike train delivered simultaneously for 10 s. The plot shows the plasticity outcome for different presynaptic firing rate (1000/frequency) for a fixed postsynaptic baseline at 10 Hz. The upper raster plot depicts the released vesicles at 40 Hz and the postsynaptic baseline at 10 Hz (including the AP evoked by EPSP). (h) Mean weight change (%) varying the rate of pre and postsynaptic Poisson spike train delivered simultaneously for 10 s. The heat map data along the vertical white dashed line is depicted in panel g.

To understand the effects of jittering, we plotted the trajectories of joint CaN-CaMKII activity (Figure 7c). 2Post1Pre50 which ‘undershoots’ the LTP region shifted into the LTP region for jitter σ=50 ms. In contrast, 1Pre1Post10 which ‘overshoots’ (mostly smaller CaN concentration) the LTP region shifted to the opposite direction towards the LTP region.

Why does jitter cause different spike timing protocols to yield similar plasticity magnitudes? Increasing jitter causes a fraction of pairings to invert causality. Therefore, the jittered protocols became a mixture of causal and anticausal pairings (Figure 7c). This situation occurs for all paired protocols. So any protocol with the same number spikes will produce a similar outcome if the jitter is large enough. Note that despite noise the mean frequency was conserved at 5 ±13.5 Hz (see Figure 7e).

Next, we studied the effect of spike removal. In the previous sections, synaptic release probability was ∼60% (for [Ca2+]o = 2.5 mM) or lower, depending on the availability of docked vesicles (Materials and methods). However, baseline presynaptic vesicle release probability is heterogeneous across CA3-CA1 synapses, ranging from 10-90% (Dobrunz et al., 1997; Enoki et al., 2009) and likely lower on average in vivo (Froemke and Dan, 2002; Borst, 2010). BaPs are also heterogeneous with random attenuation profiles (Golding et al., 2001) and spike failures (Short et al., 2017). To test the effects of pre and postsynaptic failures on plasticity induction, we performed simulations where we randomly removed spikes, altering the regular attenuation observed in Tigaret et al., 2016 protocols.

In Figure 7b, we plot the plasticity magnitude as function of sparsity (percentage of removed spikes). The sparsity had different specific effects for each protocol. 1Pre2Post10 and 1Pre2Post50 which originally produced substantial LTP were robust to spike removal until ∼60% sparsity. In contrast, the plasticity magnitude from both 1Pre1Post10 and 2Post1Pre50 showed a non-monotonic dependence on sparsity, first increasing then decreasing, with maximal LTP at ∼40% sparsity.

To understand how sparsity causes this non-monotonic effect on plasticity magnitude, we plotted the histograms of time spent in the CaN-CaMKII plane for 2Post1Pre50 for three levels of sparsity: 0%, 30%, and 80% (Figure 7d). For 0% sparsity, the activation spent most time at the border between the LTP and LTD regions, resulting in no change. Increasing sparsity to 30% caused the activation to shift rightward into the LTP region because there was less attenuation of pre and postsynaptic resources. In contrast, at 80% sparsity, the activation moved into the LTD region because there were not enough events to substantially activate CaMKII and CaN. Since LTD is a slow process and the protocol duration is short (60 s), there was no net plasticity. Therefore for this protocol, high and low sparsity caused no plasticity for distinct reasons, whereas intermediate sparsity enabled LTP by balancing resource depletion with enzyme activation.

Next we tested the interaction of jitter and spike removal. Figure 7f shows a contour map of weight change as a function of jitter and sparsity for the 2Post1Pre50 protocol, which originally induced no plasticity (Figure 3). Increasing spike jitter enlarged the range of sparsity inducing LTP. In summary, these simulations (Figure 7a, b, f and h) show that different STDP protocols have different degrees of sensitivity to noise in the firing structure, suggesting that simple plasticity rules derived from regular ex vivo experiments may not predict plasticity in vivo.

How does random spike timing affect rate-dependent plasticity? We stimulated the model with pre and postsynaptic Poisson spike trains for 10 s, under Dudek and Bear, 1992 experimental conditions. We systematically varied both the pre and postsynaptic rates (Figure 7h). The 10 s stimulation protocols induced only LTP, since LTD requires a prolonged stimulation (Mizuno et al., 2001). LTP magnitude monotonically increased with the presynaptic rate (Figure 7g and h). In contrast, LTP magnitude varied non-monotonically as a function of postsynaptic rate, initially increasing until a peak at 10 Hz, then decreasing with higher stimulation frequencies. This non-monotonic dependence on post-synaptic rate is inconsistent with classic rate-based models of Hebbian plasticity. From this analysis, we can make the prediction that firing variability can alter the rules of plasticity, in the sense that it is possible to add noise to cause LTP for protocols that did not otherwise induce plasticity. For example, we show that protocols inducing LTP can be hindered by jittering, e.g. 1Pre2Post10 in Figure 7a. Also, protocols that are just outside the LTP plasticity region may turn into LTP if jitter is applied, e.g. 2Post1Pre50 and 1Pre1Post10 in Figure 7a. We also investigated how this plasticity dependence on pre- and postsynaptic Poisson firing rates varies with developmental age (Figure 4—figure supplement 1g–i). We found that at P5 no plasticity is induced, at P15 a LTP region appears at around 1 Hz postsynaptic rate, and at P20 plasticity becomes similar to the mature age, with a peak in LTP magnitude at 10 Hz postsynaptic rate.


We built a model of a rat CA3-CA1 hippocampal synapse, including key electrical and biochemical components underlying synaptic plasticity induction (Figure 1). We developed a novel geometric readout of combined CaN-CaMKII dynamics (Figures 24) to predict the outcomes from a range of plasticity experiments with heterogeneous conditions: animal developmental age (Figure 5), aCSF composition (Figure 6), temperature (Supplemental files), and in vivo-like firing variability (Figure 7). This readout provides a simple and intuitive window into the dynamics of the synapse during plasticity. Our model is thus based on the joint activity of these two key postsynaptic enzymes at both fast and slow time scales and considers the stochastic dynamics of their activities dictated by the upstream calcium-dependent components at both the pre- and postsynapse. On this basis alone, our model is akin to biological processes where the outcome is jointly determined by several stochastic signalling components and a combination of multiple enzyme activities, that is, are multi-dimensional. The principal assumption underlying the proposed ‘geometric readout’ mechanism is that all information determining the induction of LTP vs. LTD is contained in the time-dependent spine Ca2+/calmodulin-bound CaN and CaMKII concentrations, and that no extra elements are required. Further, since both CaN and CaMKII concentrations are uniquely determined by the time course of postsynaptic Ca2+ concentration, the model implicitly assumes that the LTP/LTD induction depends solely on spine Ca2+ concentration time course, as in previously compared simplified models (see Introduction).

In addition to providing a new model of CA3-CA1 synapse biophysics, the main contribution of this work is the novel readout mechanism mapping synaptic enzymes to plasticity outcomes. This readout was built based on the concept that the full temporal activity of CaN-CaMKII over the minutes-timescale stimulus duration, and not their instantaneous levels, is responsible for changes in synaptic efficacy (Fujii et al., 2013). We instantiated this concept by analysing the joint CaN-CaMKII activity in the two-dimensional plane and designing polygonal plasticity readout regions (Figure 3a). Here, we used only a two-dimensional readout, but anticipate a straightforward generalisation to higher dimensions. The central discovery is that these trajectories, despite being stochastic, can be separated in the plane as a function of the stimulus (Figure 3). This is the basis of our new synaptic plasticity rule. We generalised previous work with plasticity induction based on single threshold and a slow variable (Badoual et al., 2006; Rubin et al., 2005; Clopath and Gerstner, 2010; Graupner and Brunel, 2012). In contrast, previous models assume that plasticity is explainable in terms of synaptic calcium or enzyme response to single BAP-EPSP pairings (Shouval et al., 2002; Karmarkar and Buonomano, 2002). We expect that future studies using high temporal resolution measurements such as those provided by recent FRET tools available for CaMKII (Chang et al., 2017; Chang et al., 2019) will bring refinements to our model with the possibility to further test our readout predictions.

Let us describe the intuition behind our model more concisely. First, we abstracted away the sophisticated cascade of plasticity expression. Second, the plasticity regions, crossed by the trajectories, are described with a minimal set of parameters. Importantly, their tuning is quite straightforward and done only once, even when the joint activity is stochastic. The tuning of the model is possible thanks to the decoupling of the plasticity process from the spine biophysics which acts as a feedforward input to the plasticity Markov chain and from the distributions of the different trajectories, which are well separated. The separability afforded by the geometrical readout, along with the model flexibility via fitting the plasticity regions, enabled us to reproduce data from nine different experiments using a single fixed set of model parameters. In contrast, we found that classic spike-timing (Song et al., 2000; Pfister and Gerstner, 2006) or calcium-threshold (Graupner and Brunel, 2012) models could not reproduce the range of protocols from Tigaret et al., 2016 (Figure 2b–d). More complicated molecular-cascade models have been shown to account for individual plasticity experiments (Antunes et al., 2016; Jędrzejewska-Szmek et al., 2017; Mäki-Marttunen et al., 2020; Bhalla, 2017), but have not been demonstrated to reproduce the wide range of protocols presented here while considering experimental heterogeneity.

For some protocols, the CaMKII-CaN trajectories overshot the plasticity regions (e.g. Figure 3d). Although abnormally high and prolonged calcium influx to cells can trigger cell death (Zhivotovsky and Orrenius, 2011), the effects of high calcium concentrations at single synapses are poorly understood. Notably, a few studies have reported evidence consistent with an overshoot, where strong synaptic calcium influx does not induce LTP (Yang et al., 1999; Tigaret et al., 2016; Pousinha et al., 2017).

Our model included critical components for plasticity induction at CA3-CA1 synapses: those affecting dendritic spine voltage, calcium signalling, and enzymatic activation. We were able to use our model to make quantitative predictions, because its variables and parameters corresponded to biological components. This property allowed us to incorporate the model components’ dependence on developmental age, external Ca/Mg levels, and temperature to replicate datasets across a range of experimental conditions. The model is relatively fast to simulate, taking ∼1 min of CPU time to run 1 min of biological time. These practical benefits should enable future studies to make experimental predictions on dendritic integration of multiple synaptic inputs (Blackwell et al., 2019; Oliveira et al., 2012; Ebner et al., 2019) and on the effects of synaptic molecular alterations in pathological conditions. In contrast, abstract models based on spike timing (Song et al., 2000; Pfister and Gerstner, 2006; Clopath and Gerstner, 2010) or simplified calcium dynamics (Shouval et al., 2002; Graupner and Brunel, 2012) must rely on ad hoc adjustment of parameters with less biological interpretability.

Intrinsic noise is an essential component of the model. How can the synapse reliably express plasticity but be noisy at the same time (Yuste et al., 1999; Ribrault et al., 2011)? Noise can be reduced either by redundancy or by averaging across time, also called ergodicity (Sterling and Laughlin, 2015). However, redundancy requires manufacturing and maintaining more components, and therefore costs energy. We propose that, instead, plasticity induction is robust due to temporal averaging by slow-timescale signalling and adaptation processes. These slow variables display reduced noise levels by averaging the faster timescale stochastic variables. This may be a reason why CaMKII uses auto-phosphorylation to sustain its activity and slow its decay time (Chang et al., 2017; Chang et al., 2019). In summary, this suggests that the temporal averaging by slow variables, combined with the separability afforded by the multidimensional readout, allows synapses to tolerate noise while remaining energy-efficient.

A uniqueness of our model is that it simultaneously incorporates biological variables such as electrical components at pre and postsynaptic sites, some with adaptive functions such as attenuation, age and temperature, stochastic noise and fast and slow timescales. Some of these variables have been modelled by other groups, for example stochasticity, BaP attenuation or pre-synaptic plasticity (Cai et al., 2007; Shouval and Kalantzis, 2005; Zeng and Holmes, 2010; Miller et al., 2005; Yeung et al., 2004; Shah et al., 2006; Deperrois and Graupner, 2020; Costa et al., 2015), but generally independently from each other. To position the uniqueness of our model in this broader context, we also provide a direct comparison of our model with some of the most recent leading models of excitatory synapse plasticity and the experimental work they reproduce (Appendix 1—table 2 and Appendix 1—table 3).

We identified some limitations of the model. First, we modelled only a single postsynaptic spine attached to a two-compartment neuron (soma and dendrite). Second, the model abstracted the complicated process of synaptic plasticity expression. Indeed, even if this replicated the early phase of LTP/LTD expression in the first 30–60 min after induction, we did not take into account slower protein-synthesis-dependent processes, maintenance processes, and synaptic pruning proceed at later timescales (Bailey et al., 2015). Third, like most biophysical models, ours contained many parameters. Although we set these to physiologically plausible values and then tuned to match the plasticity data, other combinations of parameters may fit the data equally well (Marder and Taylor, 2011; Mäki-Marttunen et al., 2020) due to the ubiquitous phenomenon of redundancy in biochemical and neural systems (Gutenkunst et al., 2007; Marder, 2011). Indeed synapses are quite heterogeneous in receptor and ion channel counts (Takumi et al., 1999; Sabatini and Svoboda, 2000; Racca et al., 2000; Nimchinsky et al., 2004), protein abundances (Shepherd and Harris, 1998; Sugiyama et al., 2005), and spine morphologies (Bartol et al., 2015; Harris and Stevens, 1989), even within the subpopulation of CA1 pyramidal neuron synapses that we modelled here. Fourth, the activation of clustered synapses could influence the plasticity outcome, and the number of synapses activated during plasticity induction can be difficult to control experimentally. Our model concerns plasticity at a single synapse, which is also important during synaptic cluster activation (Ujfalussy and Makara, 2020). We drew from data in Tigaret et al., 2016 where there is little indication of simultaneous clustered synaptic activation. Furthermore, our simulations are in good agreement with plasticity experiments using local field potential recordings (Dudek and Bear, 1993) where the number of activated synapses is uncertain. This indicates that the model proposed here can account for this aspect of synaptic plasticity heterogeneity. Finally, our readout model does not correspond to a specific molecular cascade beyond CaN and CaMKII activations. However, we anticipate that the same mapping could be implemented by simple biochemical reaction networks, with for example, transition rates based on Hill functions for the plasticity boundaries.

Since the model respected the stochasticity of vesicle release (Rizzoli and Betz, 2005; Alabi and Tsien, 2012), NMDAr (Nimchinsky et al., 2004; Popescu et al., 2004; Iacobucci and Popescu, 2017; Sinclair et al., 2016), and VGCC opening (Magee and Johnston, 1995; Sabatini and Svoboda, 2000; Iftinca et al., 2006), the magnitude of plasticity varied from simulation trial to trial. This suggests that the rules of plasticity are inherently stochastic (Bhalla, 2004; Antunes et al., 2016) and that the variability observed in these experiments (Inglebert et al., 2020; Tigaret et al., 2016; Dudek and Bear, 1992; Dudek and Bear, 1993; Mizuno et al., 2001; Meredith et al., 2003; Wittenberg and Wang, 2006) is partly due to stochastic signalling. With our current model, we have been able to reproduce nine experiments with a single set of parameters. By running extensive simulations over the space of protocols beyond those tested experimentally (Figure 3h and i; Figure 4f; Figure 5c, e and f; Figure 6e and f), we made testable predictions for plasticity outcomes that should therefore be of a high level of confidence. For example, Tigaret et al., 2016 did not find LTD when using classic post-before-pre stimulation protocols, but the model predicted that LTD could be induced if the number of pairing repetitions was extended (Figure 3h and i). The model also predicts that the lack of LTD induced by FDP in adults can be recovered using doublets, triplets or quadruplets of spikes in the protocols (Figure 5d). We tested the model’s sensitivity to spike time jitter and spike failure in the stimulation protocols (Figure 7). Our simulations predict that this firing variability can alter the rules of plasticity, in the sense that it is possible to add noise to cause LTP for protocols that did not otherwise induce plasticity.

What do these results imply about the rules of plasticity in vivo? First, we noticed that successful LTP or LTD induction required a balance between two types of slow variables: those that attenuate, such as presynaptic vesicle pools and dendritic BaP, versus those that accumulate, such as slow enzymatic integration (Cai et al., 2007; Mizusaki et al., 2022; Deperrois and Graupner, 2020). This balance is reflected in the inverted-U shaped magnitude of LTP seen as a function of post-synaptic firing rate (Figure 7h). Second, although spike timing on millisecond timescales can in certain circumstances affect the direction and magnitude of plasticity (Figure 3), in order to drive sufficient activity of synaptic enzymes, these patterns would need to be repeated for several seconds. However, if these repetitions are subject to jitter or failures, as observed in hippocampal spike trains in vivo (Fenton and Muller, 1998; Wierzynski et al., 2009), then the millisecond-timescale information will be destroyed as it gets averaged out across repetitions by the slow integration processes of CaMKII and CaN (Figure 7a–d). The net implication is that millisecond-timescale structure of individual spike pairs is unlikely to play an important role in determining hippocampal synaptic plasticity in vivo (Froemke and Dan, 2002; Sadowski et al., 2016; Graupner et al., 2016).

In summary, we presented a new type of biophysical model for plasticity induction at the rat CA3-CA1 glutamatergic synapse. Although the model itself is specific to this synapse type, the study’s insights may generalise to other synapse types, enabling a deeper understanding of the rules of synaptic plasticity and brain learning.

Materials and methods

Experimental datasets

Request a detailed protocol

The datasets at the basis of our model were obtained directly from the authors Tigaret et al., 2016 or extracted from graphs in the references in Appendix 1—table 2 using WebPlotDigitizer v 4.6 software (Rohatgi, A.). The dataset from Tigaret et al., 2016 is freely available upon request.

Code availability

Request a detailed protocol

All simulations were performed in the Julia programming language (version 1.4.2). This choice was dictated by simplicity and speed (Perkel, 2019). The code for the Markov chains is mostly automatically generated from reactions using the Julia package Catalyst.jl (Loman et al., 2022), and could be exported to an SBML representation for porting to other languages. The model is available on GitHub at SynapseElife (copy archived at Veltz, 2023).

Simulating the synapse model is equivalent to sampling a piecewise deterministic Markov process, and this relies on the thoroughly tested Julia package PiecewiseDeterministicMarkovProcesses.jl. These simulations are event-based, and no approximation is made beyond the ones required to integrate the ordinary differential equations by the LSODA method (Livermore Solver for Ordinary Differential Equations). We ran the parallel simulations in the Nef cluster operated by Inria.


Request a detailed protocol

We write 1A for the indicator of a set A, meaning that 1A(x)=1 if x belongs to A and zero otherwise.

Vesicle release and recycling

Request a detailed protocol

The detailed exocytosis model we implemented was motivated by taking into account the following minimal requirements: synaptic failures (which impact STDP protocols), vesicle depletion (for frequency/pulse number repetition dependent protocols), external calcium (motivated by Inglebert et al., 2020). This is best modeled by counting the released vesicles, hence our choice of a stochastic model. In biological synapses, an action potential arriving at the presynaptic terminal can trigger the release of a neurotransmitter–filled vesicle, which activates postsynaptic receptors. We derived a vesicle–release Markov chain model based on a deterministic approach described in Sterratt et al., 2011.

Vesicles can be in two states, either belonging to the docked pool (with cardinal D) with fast emptying, or to the reserve pool (with cardinal R) which replenishes D (Rizzoli and Betz, 2005). Initially the docked and reserve pools have D0 and R0 vesicles, respectively. The docked pool loses one vesicle each time a release occurs (Rudolph et al., 2015), with transition DD-1 (Figure 8). The reserve pool replenishes the docked pool with reversible transition (R,D)(R-1,D+1). Finally, the reserve pool is replenished towards initial value R0 over a timescale τDref via the transition (R,D)(R+1,D).

Presynaptic release model validation.

(a) Presynaptic calcium in response to the protocol 1Pre, 300 pulses at 5 Hz displaying adaptation. (b) Release probability for the same protocol as panel A but subjected to the docked vesicles availability. (c) Number of vesicles in the docked and reserve pools under depletion caused by the stimulation from panel a. (d) Plot of the mean (300 samples) release probability (%) for different frequencies for the protocol 1Pre 300 pulses at [Ca2+]o = 2.5 mM. (e) Release probability (%) for a single presynaptic spike as a function of [Ca2+]o. Note that King et al., 2001 model was multiplied by the experimentally measured release probability at [Ca2+]o = 2 mM since their model has this calcium concentration as the baseline.

In addition to the stochastic dynamics in Table 1, each spike triggers a vesicle release DD-1 with probability prel:

(1) 910prel([Ca]pre,[Ca2+]o,D)=([Ca]pre)s([Ca]pre)s+h([Ca2+]o)s1D>0     h([Ca2+]o)=0.654+1.3491+e4([Ca2+]o1.708 mM) 

which is a function of presynaptic calcium [Capre] and extracellular calcium concentration, [Ca2+]o, through the threshold h([Ca2+]o). The function prel was fitted in Figure 8e against data from Tigaret et al., 2016; Hardingham et al., 2006. To decide whether a vesicle is released for a presynaptic spike, we use a phenomenological model of Capre (see Figure 8a) based on a resource-use function (Tsodyks and Markram, 1997):

(2) {dCapredt=CapreτpreCapre(0)=0dCajumpdt=1CajumpτrecδCaCajumpCapreCajump(0)=1.
Table 1
Stochastic transitions used in the presynaptic vesicle pool dynamics.

Note that the rates depend on the number of vesicles in each pool (Pyle et al., 2000).

TransitionRateInitial Condition

Upon arrival of the presynaptic spikes, t(t1,,tn), we update Capre according to the deterministic jump:


Finally, after Capre has been updated, a vesicle may be released with probability prel (Figure 8b).

Parameters for the vesicle release model are given in Table 2. We constrained these parameters against the data reported by Hardingham et al., 2006 and Tigaret et al., 2016. Because [Ca2+]o modifies the release probability dynamics (King et al., 2001), we fixed an initial release probability to 68% for [Ca2+]o = 2.5 mM as reported by Tigaret et al., 2016 (initial value in Figure 8b and d). Additionally, Hardingham et al., 2006 reports a 38% reduction in the initial release probability when changing [Ca2+]o from 2.5 mM to 1 mM. Taking these into account, the decreasing sigmoid function in the Figure 8e depicts our [Ca2+]o-dependent release probability model (prel).

Table 2
Parameter values for the presynaptic model.

Our model does not implement a larger pool of ∼ 180 vesicles (CA3-CA1 hippocampus) sometimes named ‘reserve pool’ (Südhof, 2000) or ‘resting pool’ (Alabi and Tsien, 2012). Furthermore, what we term the ‘reserve pool’ in our model is sometimes called the ‘recycling pool’ in other studies, for example Rizzoli and Betz, 2005; Alabi and Tsien, 2012.

Vesicle release model (stochastic part)
Initial number of vesicles at DD0=255–20 (Rizzoli and Betz, 2005; Alabi and Tsien, 2012)
Initial number of vesicles at RR0=3017–20 vesicles (Alabi and Tsien, 2012)
Time constant R→ D (D recycling)τD=5s1 s (Rizzoli and Betz, 2005)
time constant D → R (R mixing)τR=45 s20 s (when depleted) to 5 min (hypertonic shock) (Rizzoli and Betz, 2005; Pyle et al., 2000)
Time constant 1→ R (R recycling)τRref=40s20–30 s (Rizzoli and Betz, 2005)
Release probability half-activation curvehsee Equation 1
Release probability sigmoid slopes=2fixed for all [Ca2+]o
Vesicle release model (deterministic part)
Capre attenuation recoveryτpre=0.02 s50-500 ms with dye (Maravall et al., 2000) therefore < 50 to 500 ms without dye
Deterministic jump attenuation recoveryτrec=20s∼20 s (Rizzoli and Betz, 2005)
Deterministic jump attenuation fractionδca=0.0004Forsythe et al., 1998

Figure 8e shows that our prel function is in good agreement with a previous analytical model suggesting that prel([Ca2+]o)([Ca2+]o)2mM-2 (King et al., 2001). Our model also qualitatively reproduces the decrease of calcium dye fluorescence levels after 20 s of theta trains from Tigaret et al., 2016. We interpret their fluorescence measurements as an effect of short-term depression (see Figure 8b).

Despite these agreements between our model and previous work, it is a simplified presynaptic model that does not encompass the highly heterogeneous nature of vesicle release. Vesicle release dynamics are known to be sensitive to various experimental conditions such as temperature (Fernández-Alfonso and Ryan, 2004), the age for some brain regions (Rudolph et al., 2015) or magnesium concentration (Hardingham et al., 2006). Although we did not consider this for the presynaptic model, note that we do incorporate such experimental parameters (age, temperature, [Ca2+]o, [Mg2+]o) to our model of NMDArs, which is the main postsynaptic calcium source. Furthermore, since our model of vesicle dynamics is simple, τrec in Equation 2 has two roles: to delay the prel recovery caused by Capre inactivation (enforced by δCa in Equation 2) and to prevent vesicle release after HFS-induced depression (King et al., 2001; Rizzoli and Betz, 2005). Our presynaptic model for prel is purely phenomenological as, in principle, the [Ca2+] jump parameter δCa should depend on [Ca2+]o, but replacing the model with a more physiological model would not affect the results, since the measured dependence of release probability on [Ca2+]o is already satisfied by this phenomenological model. Also, multi-vesicular release (MVR) at SC-CA1 synapses was shown to be prominent after manipulations that increase release probability, for example during the facilitation seen with paired-pulse stimulations (Christie and Jahr, 2006; Oertner et al., 2002). Yet, we chose not to incorporate this mechanism in the pre-synaptic model because we do not hold enough information on how MVR participates to plasticity outcomes of the different protocols we used in this study.

Model compartments

Request a detailed protocol

Our model has three compartments. 1) a spherical dendritic spine linked by the neck to 2) a cylindrical dendrite connected to 3) a spherical soma. The membrane potential of these compartments satisfy the equations below (parameters in Table 3). Since the dendrite is a single compartment, the precise spine location is undefined. For more detailed morphological simulations to predict plasticity, see Ebner et al., 2019, Chindemi et al., 2022 and Jędrzejewska-Szmek et al., 2017. The distance from the soma to the spine functionally mimics the BaP attenuation as shown in Golding et al., 2001, and it is set to 200 µm for all simulations, except in Figure 3—figure supplement 4e and Figure 4—figure supplement 1c. In these panels, we modified this distance parameter as described in the graph y-axis obtaining results similar to Ebner et al., 2019. The different currents in the soma, dendrite and spine are described in Equations 3–5.

Table 3
Parameters for the neuron electrical properties.

* The membrane leak conductance in the spine is small since the spine membrane resistance is so high that is considered infinite (>106MΩ) (Koch and Zador, 1993). The current thus mostly leaks axially through the neck cytoplasm. The dendrite leak conductance is also small in order to control the distance-dependent attenuation by the axial resistance term gBaPadapt in Equation 4 and Equation 5. The table provides the parameters associated with these equations which were adjusted (see comparison with Reference value) to fit the dynamics seen in Golding et al., 2001 and Buchanan and Mellor, 2007 experiments as in Figure 9a and b.

Passive cable
Leak reversal potentialEleak=-70mV69 mV (Spigelman et al., 1996)
Membrane leak conductance per area (for spine and passive dendrite)gleak=410-6nS/μm2* see table legend (Koch and Zador, 1993)
Membrane leak conductance per area (only soma)gsoma=5.31103 nS/μm23104 to 1.3103nS/μm2 (Fernandez and White, 2010)
Membrane capacitance per areaCm=610-3pF/μm2110-2pF/μm2
(Hines and Carnevale, 1997)
Axial resistivity of cytoplasmRa=110-2GΩμm210-3GΩμm
(Golding et al., 2001)
Dendrite diameterDdend=2μmsame as Yi et al., 2017
Dendrite lengthLdend=1400μmapical dendrites, 1200–1600 μm
López Mendoza et al., 2018
Dendrite surface areaAdend=8.79103μm2πDdendLdend
Dendrite volumeVoldend=4.4103μm3π(Ddend/2)2Ldend
Dendritic membrane capacitanceCdend=52.77pFCmAdend
Dendrite leak conductancegLdend=3.5110-2nSgleakAdend
Dendrite axial conductancegdiff=50nSRaAdend
Soma diameterDsoma=30μm21 μm (Stuart et al., 2016) page 3
Soma area (sphere)Asoma=2.82103μm2(4π/3)(Dsoma/2)3;
(Zhuravleva et al., 1997)
Soma membrane capacitanceCsoma=16.96pFCmAsoma
Soma leaking conductancegLsoma=15nSgsomaAsoma
(Fernandez and White, 2010)
Dendritic spine
Spine head volumeVolsp=0.03μm3Bartol et al., 2015
Spine head surfaceAsp=4.6610-1μm24π(3Volsp/4π)2/3
Spine membrane capacitanceCsp=2.810-3pFCmAsp
Spine head leak conductancegLsp=1.8610-6nSgleakAsp
Dendritic spine neck
Spine neck diameterDneck=0.1μm0.05–0.6 μm Harris et al., 1992
Neck lengthLneck=0.2μm0.7±0.6μm (Adrian et al., 2017)
Neck cross-sectional areaCSneck=7.8510-3μm2π(Dneck/2)2
Neck resistancegneck=3.92nS255.1MΩCSneck/(LneckRa)
50 to 550 MΩ (275±27 MΩ) (Popovic et al., 2015)

Membrane potential and currents

Request a detailed protocol

The membrane potential of these compartments satisfy the equations below (parameters in Table 3). The different currents are described in the following sections.

(3) CspdVspdt=gneck(VdendVsp)+gLsp(ErevVsp)+IT+IL+IR+INMDA+IAMPA+ISK
(4) CdenddVdenddt=gBaPadapt(VsomaVdend)+gneck(VspVdend)+gLdend(ErevVdend)+IGABA
(5) CsomadVsomadt=gBaPadapt(VdendVsoma)+gLsoma(ErevVsoma)+λage(IBaP+INa)+IK

Action-potential backpropagation (BaP)

Postsynaptic currents

Request a detailed protocol

The postsynaptic currents are generated in the soma, backpropagated to the dendritic spine via a passive dendrite which filters them. The soma generates BaPs using a version of the Na+ and K+ channel models developed by Migliore et al., 1999. The related parameters are described in Table 4 (the voltage unit is mV). We used the following the units: αx(V/mV)[ms-1] and βx(V/mV)[ms-1].

Table 4
The Na +and K+conductances intentionally do not match the reference because models with passive dendrite need higher current input to initiate action potentials (Levine and Woody, 1978).

Therefore, we set it to achieve the desired amplitude on the dendrite and the dendritic spine according to the predictions of Golding et al., 2001 and Kwon et al., 2017.

Soma parameters for Na +and K+ channel
Sodium conductanceγNa=8102nSgeneric value, see legend commentary
Potassium conductanceγK=40nSgeneric value, see legend commentary
Reversal potential sodiumErevNa=50mVMigliore et al., 1999
Reversal potential potassiumErevK=-90mVMigliore et al., 1999
BaP attenuation parameters
Attenuation step factor (age)δagesee Equation 7 and Figure 9b and c bottom Buchanan and Mellor, 2007; Golding et al., 2001
Attenuation step factorδdecay=1.72710-5adjusted to fit Buchanan and Mellor, 2007; Golding et al., 2001
Auxiliary attenuation step factorδaux=2.30410-5adjusted to fit Buchanan and Mellor, 2007; Golding et al., 2001
Recovery time for the attenuation factorτrec=2sadjusted to fit Buchanan and Mellor, 2007; Golding et al., 2001
Recovery time for the age attenuation factorτrecage=0.5sadjusted to fit Buchanan and Mellor, 2007; Golding et al., 2001
AP evoked by EPSP
Decay time for VevokeτV=0.04sHines and Carnevale, 1997
Delay AP evoked by EPSPδdelay-AP=0.015sFricker and Miles, 2000
Sodium channelαm(Vsoma)=0.4Vsoma+301eVsoma+307.2βm(Vsoma)=0.124Vsoma+30eVsoma+307.21minf(Vsoma)=αm(Vsoma)αm(Vsoma)+βm(Vsoma)mτ(Vsoma)=1αm(Vsoma)+βm(Vsoma)αh(Vsoma)=0.01Vsoma+45eVsoma+451.51βh(Vsoma)=0.03Vsoma+451eVsoma+451.5dhdt=αh(Vsoma)(1h)βh(Vsoma)hdmdt=minfmmτINa=γNam3h(ErevNaVsoma).Potassium channelαn(Vsoma)=e0.11(Vsoma13)βn(Vsoma)=e0.08(Vsoma13)ninf(Vsoma)=11+αn(Vsoma)nτ(Vsoma)=max(50βn(Vsoma)1+αn(Vsoma);2)dndt=ninfnnτIK=γKn(ErevKVsoma)

To trigger a BaP, an external current IBaP is injected in the soma at times t{t1,,tn} (postsynaptic input times) for a chosen duration δinj with amplitude Iamp(nA). Taking H as the Heaviside function, this is expressed as:


The currents underlying the BaP in the soma are filtered in a distance-dependent manner by the dendrite before reaching the dendritic spine. Biologically, BaP adaptation is caused by the inactivation of sodium channels and variations in sodium and potassium channel expression along the dendrite (Jung et al., 1997; Golding et al., 2001). We used a phenomenological model, implementing distant-dependent BaP amplitude attenuation by modifying the axial resistance gBaPadapt (see Equation 4 and Equation 5) between the dendrite and the soma as follows (Figure 9c top):

(6) gBaPadapt=λgdiffϕdist(dsoma),ϕdist(dsoma)=0.1+1.41+e0.02(dsoma230.3μm)

where δdecay is the distance of the spine from the soma and where the factor is dynamically regulated based on a resource-use equation from Tsodyks and Markram, 1997 with a dampening factor λaux changing the size of the attenuation step δdecay:


The BaP attenuation model is based on Golding et al., 2001 data for strongly attenuating neurons. Therefore, the second type of attenuation (weakly attenuating) in neurons is not considered (dichotomy in Figure 9a). Figure 9a compares Golding’s data to our model and illustrates the effect of BaP attenuation in the upper panels of Figure 9a, b.

AP Evoked by EPSP.

(a) Model and data comparison for the distance-dependent BaP amplitude attenuation measured in the dendrite and varying the distance from the soma. The stimulation in panel a is set to reproduce the same stimulation as Golding et al., 2001. Golding described two classes of neurons: those that are strongly attenuated and those that are weakly attenuated (dichotomy mark represented by the dashed line). However, in this work we consider only strongly attenuated neurons. (b) Attenuation of somatic action potential from Buchanan and Mellor, 2007 and model in response to five postsynaptic spikes delivered at 100 Hz. The value showed for the model is the spine voltage with distance from the soma set to zero (scale 25 ms, 20 mV). (c) Top panel shows the λsoma used in Equation 6 to modify the axial conductance between the soma and dendrite. Bottom panel shows the age-dependent changes in the step of the resource-use equation (Equation 7) that accelerates the BaP attenuation and decreases the sodium currents in Equation 5. (d) Probability of evoking an AP multiplied by the successfully evoked AP, pAP(Vevoked)1(evoked), for the protocol 1Pre, 300 at 5 Hz (2.5 mM Ca). (e) Two-pool dynamics with the same stimulation from panel D showing the vesicle release, the reserve and docked pools, and the evoked AP. (f) Probability of evoking an AP for the protocol 1Pre 300 pulses at different frequencies (3 and 5 Hz have the same probability).

Table 4 shows the BaP attenuation parameters. The plasticity outcomes as function of the dendritic spine distance from the soma are shown in Figure 3—figure supplement 4e and Figure 4—figure supplement 1c.

Age-dependent BaP adaptation

Request a detailed protocol

Neuronal bursting properties are altered during development through development through the maturation and expression of potassium and sodium channels (Gymnopoulos et al., 2014), which change the interaction of hyperpolarizing and depolarizing currents (see Figure 9b; Grewe et al., 2010; Jung et al., 1997). We fitted the data of the age dependent somatic attenuation profiles from Buchanan and Mellor, 2007 (Figure 9b) with our model by including an age-dependent BaP amplitude attenuation factor. We define the attenuation factor λage (Figure 9c bottom), as follows.

(7) dλauxdt=1IageτrecageδageλageIBaP(t),δrecage=1.3911041+e0.135(age16.482 days).

In Equation 5, the age effects are introduced by multiplying the sodium INa and the external IBaP currents by the attenuation factor λage.

AP evoked by EPSP

Request a detailed protocol

Biologically, a presynaptic stimulation triggers a BaP if sufficient depolarization is caused by the EPSPs reaching the soma (Stuart et al., 2016). To model this effect for some LFP recordings protocols, we included an option to choose whether an EPSP can evoke an AP using an event generator resembling the previous presynaptic release probability model prel, as in Equation 1. Like prel, the BaPs evoked by EPSPs are calculated offline, before the postsynaptic simulation. We use a variable Vevoke which is incremented by 1 at each presynaptic time t(t1,,tn) and has exponential decay:

(8) {dVevokedt=VevokeτvVevoke(0)=0VevokeVevoke+1.

Since the BaPs evoked by EPSPs are triggered by the afferent synapses and are limited by their respective docked pools (D), we use the previous prel to define the probability of an AP to occur. We test the ratio of successful releases from 25 synapses to decide if a BaP is evoked by an EPSP, setting a test threshold of 80%. Therefore, we express the probability of evoking an AP, pAP(Vevoke), with the following test:

251(rand<prel(Vevoked,[Ca2+]o,D))25>80% .

In real neurons, the EPSP summation dynamics in the soma and dendrites depend on the complex neuron morphology (Etherington et al., 2010; Ebner et al., 2019) which was not implemented by our model. Instead our ‘AP evoked by EPSP test’ is a simplified way to produce BaPs, similar to an integrate-and-fire model (Sterratt et al., 2011).

Previous work (Mayr and Partzsch, 2010) suggests that low-frequency stimulation, as used in Dudek and Bear, 1992 ([Ca2+]o = 2.5 mM), can evoke BaPs with a ∼5% probability. Our model accounts for this, but also allows the probability of eliciting an AP to increase with stimulation frequency (Etherington et al., 2010). This is captured by Vevoke as shown in Figure 9f. Figure 9d, e show how a 5 Hz stimulation can evoke APs. The delay between the EPSP and the evoked AP is set to δdelay-AP=15ms, similar to the EPSP-spike latency reported for CA1 pyramidal neurons (Fricker and Miles, 2000).


Markov chain

Request a detailed protocol

The AMPArs are modeled as a Markov chain (Figure 10) described by Robert and Howe, 2003 and Coombs et al., 2017 and adapted to temperature changes according to Postlethwaite et al., 2007. Here, we introduce the additional parameters ρfAMPA,ρbAMPA to cover AMPAr temperature-sensitive kinetics (Postlethwaite et al., 2007). The corresponding parameters are given in Table 5.

AMPAr Markov chain with three sub-conductance states and two desensitisation levels.

It includes parameters ρfAMPA, ρbAMPA (binding and unbinding of glutamate) which depend on temperature. Open states are O2, O3, and ; closed states are C0, C1, C2, C3, and C4; desensitisation states are D0, D1, D2, D3, and D4; deep desensitisation states are D22, D23, and D24.

Table 5
Parameter values for the AMPAr Markov chain and glutamate release affecting NMDAr, AMPAr.

Properties of GABA release are the same as those for glutamate.

Glutamate parameters
Duration of glutamate in the cleftgluwidth=0.001sSpruston et al., 1995
Concentration of glutamate in the cleftgluamp=103μMSpruston et al., 1995
Glutamate variability (gamma distribution Γ)glucv=Γ(1/0.52,0.52)Liu et al., 1999
Glutamate signalGluglucvgluamp for AMPAr, NMDAr and copied to GABA neurotransmitter
AMPAr parameters
Number of AMPArsNAMPA=120Bartol et al., 2015
Reversal potentialErevAMPA=0mVBartol et al., 2015
Subconductance O2γA2=0.0155nS0.0163nS (Coombs et al., 2017)
Subconductance O3γA3=0.026nS0.0287nS (Coombs et al., 2017)
Subconductance O4γA4=0.0365nS0.0378nS (Coombs et al., 2017)
glu bindingk1=16μM-1s-1Robert and Howe, 2003
glu unbinding 1k-1=7400s-1Robert and Howe, 2003
glu unbinding 2k-2=0.41s-1Robert and Howe, 2003
Closingα=2600s-1Robert and Howe, 2003
Openingβ=9600s-1Robert and Howe, 2003
Desensitisation 1δ1=1500s-1Robert and Howe, 2003
Desensitisation 2δ2=170s-1Robert and Howe, 2003
Desensitisation 3δ0=0.003s-1Robert and Howe, 2003
Re-desensitisation 1γ1=9.1s-1Robert and Howe, 2003
Re-desensitisation 2γ2=42s-1Robert and Howe, 2003
Re-desensitisation 3γ0=0.83s-1Robert and Howe, 2003

The AMPAr conductance is given as the sum of the occupancies of the three subconductance states O2, O3 and O4 of the Markov chain in Figure 10. The AMPAr current is then:


The adaptation of the Markov chain from Robert and Howe, 2003 is made by changing the forward ρfAMPA and backward ρbAMPA rates in a temperature-dependent manner matching the decay time reported by Postlethwaite et al., 2007:


The effects of temperature change on AMPAr dynamics are presented in Figure 11, which also shows that the desensitisation is not altered by temperature changes (Figure 11b and c). The recovery time from desensitisation is the same as at room temperature (Robert and Howe, 2003). Desensitisation measurements are required to account for a temperature-dependent change in the rates of the ‘vertical’ transitions in Figure 10, see Postlethwaite et al., 2007. This could be relevant for presynaptic bursts as indicated in Figure 11b–c.

Effect of temperature on AMPArs.

(a) Probability of AMPAr opening (O2+O3+O4NAMPA) and the decay time at different temperatures, in response to 1 mM glutamate applied for 1 ms (standard pulse). Postlethwaite et al., 2007 data (our model) suggests that AMPAr decay time at 35°C is 0.5ms (0.6ms) and at 25°C is 0.65ms (0.95ms). This shows a closer match towards more physiological temperatures. (b) Desensitisation profile of AMPAr at 35°C showing how many AMPAr are open in response to repeated glutamate saturating pulses (5 mM Glu for 20 ms) separated by an interval (x-axis). (c) Same as in panel b but for 25°C.

Postsynaptic Ca2+ influx

Request a detailed protocol

The effects of experimental conditions on the calcium dynamics are due to receptors, ion channels and enzymes. A leaky term models the calcium resting concentration in the Equation 9. The calcium fluxes from NMDAr and VGCCs (T, R, L types) are given in Equation 10. The diffusion term through the spine neck is expressed in Equation 11. Finally, the buffer, the optional dye and the enzymatic reactions are given in Equation 12 (parameter values given at the Table 6):

(9) dCadt=CaCaτCa+
(10) CaNMDA+IT+IR+IL2FAsp+
(11) max(Ca,Ca/3)CaτCaDiff
(12) dBuffCadtdDyedt+enzymes.
Table 6
Postsynaptic calcium dynamics parameters.

Note that the buffer concentration, calcium diffusion coefficient, calcium diffusion time constant and calcium permeability were considered free parameters to adjust the calcium dynamics.

Buffer and dye
Association buffer constantkonBuff=247μM-1s-1Bartol et al., 2015
Dissociation buffer constantkoffBuff=524s-1Bartol et al., 2015
Generic buffer concentrationBuffcon=62μM76.7μM (Bartol et al., 2015)
Calcium dynamics
Calcium baseline concentrationCa=0.05μM0.037±0.005 to 0.054±0.005 μM
(Maravall et al., 2000)
Calcium decay timeτCa=10-2s0.05–0.5 s with dye (Maravall et al., 2000) therefore < 0.05–0.5 s undyed (unbuffered)
Calcium diffusion coefficientDCa=333.8μm2s-1200 to 400 μm2s1 (Bartol et al., 2015; Holcman et al., 2005)
Calcium diffusion time constantτCaDiff=Volsp2DCa2Dneck+Lneck22DCa=5104 s8 ms for a 
Volsp=0.7μm3 (Holcman et al., 2005)
GHK equation
TemperatureT=35Cconverted to Kelvin in the Equation 14 given the protocol
Faraday constantF=96.485Cmol-1Hille, 1978
Gas constantR=8.314JK-1mol-1Hille, 1978
Calcium permeabilityPCa=45 μm s1adjusted to produce 3μM Calcium in response to a Glu release supplementary file from Chang et al., 2017
Calcium ion valencezCa=2Hille, 1978

Despite the driving force to the resting concentration, Ca=50nM, the tonic opening of T-type channels causes calcium to fluctuate making its mean baseline value dependent on temperature, extracellular calcium, and voltage. The effects of this tonic opening in various experimental conditions are shown in Figure 6—figure supplement 2f. To avoid modelling dendritic calcium sources, we use a dampening term as one-third of the calcium level since calcium imaging comparing dendrite and spine fluorescence have shown this trend (Segal and Korkotian, 2014). Equation 11 implements the diffusion of calcium from the spine to the dendrite through the neck. The time constant for the diffusion coefficient τCaDiff, is estimated as described in Holcman et al., 2005. The calcium buffer and the optional dye are described as a two-state reaction system (Sabatini et al., 2002):

(13) dBuffCadt=konBuff(BuffconBuffCa)CakoffBuffBuffCadDyedt=konFluo5(Fluo5fconDye)CakoffFluo5Dye.

Tigaret et al., 2016 experiments used the synthetic calcium-indicator dye Fluo-5f, which is well-modelled by a single Ca2+-dye binding reaction (Maravall et al., 2000; Bartol et al., 2015). Although we include a detailed model of Calmodulin, which is a major endogenous calcium buffer, the other types are poorly quantified experimentally. Instead, we used a parsimonious generic buffer model that represents an aggregate of these largely unknown endogenous buffers. Future iterations of the model could include more detailed versions of these endogenous buffers, for example calbindin (Bartol et al., 2015).

Unlike other calcium-based plasticity models (Graupner and Brunel, 2012) using the dye fluorescence decay as an approximation to calcium decay, our model is based on receptor and ion channel kinetics. Additionally, our model can simulate the dye kinetics as a buffer using (Equation 13) when appropriate. Figure 12 highlights differences between calcium and dye dynamics which is affected by the laser-induced temperature increase (Wells et al., 2007; Deng et al., 2014). We estimated the calcium reversal potential for the calcium fluxes using the Goldman–Hodgkin–Katz (GHK) flux equation described in Hille, 1978. The calcium ion permeability, PCa, was used as a free parameter adjusting a single EPSP to produce a calcium amplitude of ∼3 µM (Chang et al., 2017). This free scaling is needed to compensate for the fact that that GHK equation is derived for a model that assumes ionic currents pass through the membrane as a distributed and continuous flux, rather that the ion channels we modelled as having discrete conductance levels. Although this adaptation implies that we are using the magnitude of the GHK flux in a phenomenological way, it nevertheless captures the nonlinear dependence of relative calcium flux on extracellular calcium concentration.

(14) ϕ(Vsp,T)=zCaVspF/(R(T+273.15K))ΦCa(Vsp,[Ca2+]i)=PCazCaFϕ(Vsp,T)[Ca2+]i[Ca2+]oeϕ1eϕ,

where ΦCa(Vsp,[Ca2+]i) is used to determine the calcium influx through NMDAr and VGCC in the Equation 15, Equation 16, Equation 17 and Equation 18 using the spine membrane voltage and calcium internal concentration [Ca2+]i. Note that for simplicity the calcium external concentration [Ca2+]o was kept fixed during the simulation and only altered by experimental conditions given by the aCSF composition.

Differences between dye measurements and simulated calcium.

(a), Pre and postsynaptic stimuli as used in Tigaret et al., 2016. (b), Calcium imaging curves (fluorescence ΔF/A) elicited using the respective stimulation protocols above with Fluo5 200 μM (extracted from Tigaret et al., 2016). Scale 100 ms, 0.05 ΔF/A. (c), Dye simulation using the model. The dye is implemented by increasing temperature to mimic laser effect on channel kinetics and decreases the interaction between NMDAr and voltage elicited by BaP. Temperature effects over NMDAr are shown in Korinek et al., 2010. Also, the effects of temperature on calcium-sensitive probes shown in Oliveira et al., 2012 (baseline only, likely related to T-type channels). Other examples of laser heating of neuronal tissue are given in Deng et al., 2014. Such a dye curve fitting was obtained by increasing temperature by 10°C to mimic laser-induced heating (Wells et al., 2007; Deng et al., 2014). We achieved a better fit by decreasing the amplitude of the BaP that reaches the dendrite. Additionally, for fitting purposes, we assumed that a temperature increase lead to a decrease in BaP amplitude. Scale 0.6 µM dye, 100 ms. (d), Calcium simulation without dye. Scale 0.85 µM Ca2+, 100 ms.

NMDAr - GluN2A and GluN2B

Markov chain

Request a detailed protocol

In hippocampus, NMDArs are principally heteromers composed of the obligatory subunit GluN1 and either the GluN2A or GluN2B subunits. These N2 subunits determine the kinetics of these receptors, with the GluN1/GLUN2B heteromers displaying slow kinetics (∼250 ms) and the GluN1/GluN2A heteromers displaying faster kinetics (∼50 ms). We modeled both NMDA subtypes. The NMDAr containing GluN2A is modeled with the following Markov chain (Popescu et al., 2004):


where we have introduced the additional parameters ρfNMDA,ρbNMDA to account for temperature dependence (below).

The NMDAr containing GluN2B is modeled with a Markov chain based on the above GluN2A scheme. We decreased the rates by ∼75% in order to match the GluN2B decay at 25°C as published in Iacobucci and Popescu, 2018.


The different rates are given in Table 7.

Table 7
NMDAr parameters.

The existing model of NMDAr (GluN2A) was adapted to obtain the NMDAr (GluN2B) model. The decay time of NMDAr (GluN2B) was fitted to match decay time in Iacobucci and Popescu, 2018 and the temperature dependence uses the EPSP decay time from Korinek et al., 2010.

NMDAr (GluN2A)
Glutamate bindingka=34μM-1s-1Popescu et al., 2004
Glutamate bindingkb=17μM-1s-1Popescu et al., 2004
Forward ratekc=127s-1Popescu et al., 2004
Forward ratekd=580s-1Popescu et al., 2004
Opening rateke=2508s-1Popescu et al., 2004
Opening ratekf=3449s-1Popescu et al., 2004
Closing ratek-f=662s-1Popescu et al., 2004
Closing ratek-e=2167s-1Popescu et al., 2004
Backward ratek-d=2610s-1Popescu et al., 2004
Backward ratek-c=161s-1Popescu et al., 2004
Glutamate unbindingk-b=120s-1Popescu et al., 2004
Glutamate unbindingk-a=60s-1Popescu et al., 2004
NMDAr (GluN2B)
Glutamate bindingsb=0.25kbadapted from GluN2A (Popescu et al., 2004; Iacobucci and Popescu, 2018)
Glutamate bindingsc=0.25kcadapted from GluN2A (Popescu et al., 2004; Iacobucci and Popescu, 2018)
Forward ratesc=0.25kcadapted from GluN2A (Popescu et al., 2004; Iacobucci and Popescu, 2018)
Forward ratesd=0.25kdadapted from GluN2A (Popescu et al., 2004; Iacobucci and Popescu, 2018)
Opening ratese=0.25keadapted from GluN2A (Popescu et al., 2004; Iacobucci and Popescu, 2018)
Opening ratesf=0.25kfadapted from GluN2A (Popescu et al., 2004; Iacobucci and Popescu, 2018)
Closing rates-f=0.23k-fadapted from GluN2A (Popescu et al., 2004; Iacobucci and Popescu, 2018)
Closing rates-e=0.23k-eadapted from GluN2A (Popescu et al., 2004; Iacobucci and Popescu, 2018)
Backward rates-d=0.23k-dadapted from GluN2A (Popescu et al., 2004; Iacobucci and Popescu, 2018)
Backward rates-c=0.23k-cadapted from GluN2A (Popescu et al., 2004; Iacobucci and Popescu, 2018)
Glutamate unbindings-b=0.23k-badapted from GluN2A (Popescu et al., 2004; Iacobucci and Popescu, 2018)
Glutamate unbindings-a=0.23k-aadapted from GluN2A (Popescu et al., 2004; Iacobucci and Popescu, 2018)
Other parameters
Total number of NMDArNNMDA=155–30 (Spruston et al., 1995; Bartol et al., 2015; Nimchinsky et al., 2004)
Distribution of GluN2A and GluN2Bdefined by  ragefitted from Sinclair et al., 2016, see Figure 13b and e
NMDAr conductance depending on calciumγNMDAfitted from Maki and Popescu, 2014, see Figure 13c
NMDAr reversal potentialErev =NMDA0mVDestexhe et al., 1994
Fraction of calcium carried by NMDArfCa=0.1Griffith et al., 2016

NMDAr and age switch

Request a detailed protocol

The age-dependent expression ratio of the subtypes GluN2A and GluN2B (rage) was obtained from experimental data of mouse hippocampus (Sinclair et al., 2016). We added noise to this ratio causing ∼1 NMDAr subunit to flip towards GluN2A or GluN2B (see Figure 13e). The population of 15 total NMDArs is divided in the two subtypes according to the ratio plotted in Figure 13b, as a function of age. The ratio to define the number NMDAr subtypes as function of age reads:

rage=0.507+0.9641+e0.099(age25.102 days)+N(0,0.05)NGluN2B=round(NNMDAragerage+1)NGluN2A=round(NNMDArage+1).
NMDAr changes caused by age, temperature and extracellular and magnesium concentrations in the aCSF.

(a) Decay time of the NMDAr-mediated EPSP recorded from neocortical layer II/III pyramidal neurons (grey) (Korinek et al., 2010) compared to the decay time from the GluN2B channel estimated by our model (yellow) and data from Iacobussi’s single receptor recording (purple) (Iacobucci and Popescu, 2018). (b), Comparison of our model of the GluN2B:GluN2A ratio and the GluN2B:GluN2A ratio from the mouse CA1 excitatory neurons. (c), Comparison of our model of NMDAr conductance change as a function of extracellular calcium, against data (Maki and Popescu, 2014). (d), Forward and backwards temperature factors implemented to approximate NMDAr subtypes decay times at room temperature (Iacobucci and Popescu, 2018) and temperature changes observed in Korinek et al., 2010. (e), NMDAr subtype number in our model as a function of animal age. We added noise to have a smoother transition between different ages. (f), Calcium concentration changes for causal and anticausal protocols in response to different aCSF calcium and magnesium compositions with fixed Ca/Mg ratio (1.5). Scale bars 50 ms and 5 μM.

The round term in the two previous equations ensures that we have an integer value for the NMDAr subtypes, making the stair shaped curve seen in Figure 13e.

NMDAr and temperature

Request a detailed protocol

We adjusted the GluN2A and GluN2B forward and backward rates to follow the temperature effects on NMDAr-mediated EPSPs (Korinek et al., 2010), see Figure 13a and d. Because GluN2B dominates the NMDAr-mediated EPSP, we fit the GluN2B decay time to data on the NMDAr-mediated EPSP,w as function of temperature as reported by Korinek et al., 2010 using logistic functions ρfNMDA and ρbNMDA. The decay time comparison is shown in Figure 13a. Then, we applied the same temperature factor ρfNMDA and ρbNMDA for GluN2A. The decay times of GluN2A and GluN2B are similar to those reported by Iacobucci and Popescu, 2018. The forward and backward factors are described as follows:


NMDAr current and Ca2+-dependent conductance

Request a detailed protocol

NMDAr conductance is modulated by external calcium and is modelled according to the next equations using NMDAr subconductances AO1 and AO2 (GluN2A), and BO1 and BO2 (GluN2B).

γNMDA=33.949pS+58.3881+e4([Ca2+]o2.701 mM)pSB(Vsp,[Mg]o)=11+[Mg]o3.57mMe0.062Vsp/mVNMDA=(BO1+BO2+AO1+AO2)B(Vsp,[Mg]o)γNMDAINMDA=(ErevNMDAVsp)NMDA

We modified the conductance γNMDA as a function of extracellular calcium from that reported by Maki and Popescu, 2014. The reported NMDAr conductance at [Ca2+]o = 1.8 mM is 53±5pS. Here, we used the higher conductance 91.3pS for NMDAr (for both subtypes) at [Ca2+]o = 1.8 mM to compensate for the small number of NMDArs reported by Nimchinsky et al., 2004. Hence, we adjusted the Maki and Popescu, 2014 data to take into account this constraint: this caused a rightward-shift in the NMDA-conductance curve (Figure 13c). The calcium influx CaNMDA is modulated by the GHK factor, Equation 14, as a function of the internal and external calcium concentrations and the spine voltage:


The combined effect of extracellular Magnesium (Jahr and Stevens, 1990) and Calcium concentration are displayed in Figure 13f.

GABA(A) receptor

Since the precise delay of GABA release relative to glutamate is not known, we assumed GABA and glutamate release are synchronized for simplicity (see Table 5). We used the GABA(A) receptor Markov chain (Figure 14) presented in Busch and Sakmann, 1990; Destexhe et al., 1998 and we estimated temperature adaptations using the measurements reported by Otis and Mody, 1992. Table 8 presents the GABAr model parameters.

GABAr Markov chain model.

Closed states (C0, C1 and C2) open in response to GABAr and can go either close again or open (O1 and O2).

Table 8
GABAr parameters.

The GABAr number and conductance were modified to fit GABAr currents as in Figure 15b and e.

GABA(A) receptor
Number of GABArNGABA=3430 Edwards et al., 1990
Chloride reversal potentialsee age-dependent equationfitted from Rinetti-Vargas et al., 2017
GABAr conductanceγGABA=0.036nS0.027 nS (Macdonald et al., 1989)
Bindingrb1=20μM-1s-1Busch and Sakmann, 1990
Unbindingru1=4.6103s-1Busch and Sakmann, 1990
Bindingrb2=10μM-1s-1Busch and Sakmann, 1990
Unbindingru2=9.2103s-1Busch and Sakmann, 1990
Opening raterro1=3.3103s-1Busch and Sakmann, 1990
Opening raterro2=10.6103s-1Busch and Sakmann, 1990
Closing raterc2=400s-1temperature changes to fit Otis and Mody, 1992; Busch and Sakmann, 1990
Closing raterc2=9.8103s-1temperature changes to fit Otis and Mody, 1992; Busch and Sakmann, 1990

GABA(A)r and temperature

Request a detailed protocol

Because the amplitude of GABA(A) current is altered by the GABAr reversal potential shift during development (Rinetti-Vargas et al., 2017), we applied temperature changes only to the closing rates using a logistic function for ρbGABA, estimated by fitting to the measurements from Otis and Mody, 1992 (data comparison in the Figure 15b and e).

GABA(A)r current, kinetics and chloride reversal potential.

(a) States of GABA(A)r Markov chain at 25°C in response to a presynaptic stimulation. Opened = O1+O2, closed = C0+C1+C2. (b) Model and data comparison (Otis and Mody, 1992) for GABA(A)r current at 25°C. Even though data were recorded from P70 at 25°C and P15 at 35°C, we normalize the amplitude to invert the polarity and compare the decay time. This is done since the noise around P15 can either make GABAr excitatory or inhibitory as shown by Ecl data in panel c. (c) Chloride reversal potential (ErevCl) fitted to Rinetti-Vargas et al., 2017 data. Note that we used both profiles from axon and dendrite age-depended ErevCl changes since exclusive dendrite data is scarce. (d) States of simulated GABA(A)r Markov chain at 35°C in response to a presynaptic stimulation. (e) Model and data comparison (Otis and Mody, 1992) for GABA(A)r current at 25°C (same normalization as in panel b). (f) Change in the polarization of GABA(A)r currents given the age driven by the ErevCl.

GABA(A)r current and age switch

Request a detailed protocol

The GABA(A)r-driven current changes during development (Meredith et al., 2003) passing from depolarizing (excitatory) to hyperpolarizing (inhibitory) (Chamma et al., 2012). The reversal potential of chloride ions permeating GABA(A)r shifts from above the membrane resting potential (inward driving force - excitatory) to below the membrane resting potential (outward driving force - inhibitory; Rinetti-Vargas et al., 2017). This effect mediated is associated with the KCC2 pump (K Cl co-transporter) which becomes efficient in extruding chloride ions during maturation (Rinetti-Vargas et al., 2017). To account for the GABA(A)r age-dependent shift, we fit a function for the chloride reversal potential (ErevCl) to the data published by Rinetti-Vargas et al., 2017 (Figure 15c):

ErevCl=92.649+243.5151+e0.091(age0.691 days)IGABA=(O1+O2)(ErevClVdend)γGABA.

VGCC - T, R, and L type

Markov chain

Request a detailed protocol

A stochastic VGCC model was devised using the channel gating measurements from rat CA1 (2–8 weeks) pyramidal neurons by Magee and Johnston, 1995 at room temperature. Our model has three different VGCC subtypes described by the Markov chains in Figure 16: the T-type (low-voltage), the R-type (medium-to-high-voltage) and the L-type (high-voltage).

From left to right, R-, L-, and T-type VGCCs Markov chain adapted from Magee and Johnston, 1995.

The R- (left scheme) and T- type (right scheme) have a single open state (red colour), respectively, Or and OT. The L-type VGCC (middle) has two open states, OL1 and OL2.

The VGCC Markov chain models were derived from voltage activation and inactivation profiles reported in Magee and Johnston, 1995. The T- (Figure 17a and d) and R-type (Figure 17b and e) models are composed of independent activation (m) and inactivation (h) gating variables, while the L-type (Figure 17c) model has one closed state but two open states, to capture the two timescales of channel closing kinetics reported by Magee and Johnston, 1995. The VGCC model equations are given below. We used the following the units: αx(V/mV)[ms-1] and βx(V/mV)[ms-1].

VGCC rates and temperature factors.

(a), Activation (αm(Vsp)) and deactivation rates (βm(Vsp)) for the T-type m-gate. (b), Activation (αm(Vsp)) and deactivation rates (βm) for the R-type m-gate. (c), Activation (αm(Vsp)) and both deactivation rates (β2L(Vsp) and β21(Vsp)) for the L-type VGCC. (d), Activation (αh(Vsp)) and deactivation rates (βh(Vsp)) for the T-type h-gate. (e), Activation (αh(Vsp)) and deactivation rates (βh(Vsp)) for the R-type h-gate. (f), Temperature factor applied to all the rates, forward change (ρfVGCC) for the α rates and backward change (ρbVGCC) for the β rates.

R-type h-gate ratesτhR=100hinfR(Vsp)=11+eVsp+399.2αhR(Vsp)=hinfRτhRβhR(Vsp)=1hinfRτhRL-type ratesαL(Vsp)=0.831+e13.7Vsp6.1β1L(Vsp)=0.531+eVsp11.56.4β2L(Vsp)=1.861+eVsp18.86.17
T-type h-gate ratesτhT=50hinfT(Vsp)=11+eVsp+706.5αhT(Vsp)=hinfTτhTβhT(Vsp)=1hinfTτhT
R-type m-gate ratesβmR=40minfR=11+e3108αmRr=βmRminfR1minfRτmR=1αmR+βmRminfR=11+e3Vsp8αmR(Vsp)=minfRτmRβmR(Vsp)=1minfRτmRT-type m-gate ratesβmT=1minfT=11+e32+207αmTr=βmTminfT1minfTτmT=1αmT+βmTminfT=11+e32Vsp7αmT(Vsp)=minfTτmTβmT(Vsp)=1minfTτmT

VGCC and temperature

Request a detailed protocol

We used the same temperature factor for every VGCC subtype, respectively ρfVGCC and ρbVGCC (see Figure 17f), as follows:


The VGCC subtypes have different sensitivities to temperature, with temperature factors for decay times ranging from 2 (Iftinca et al., 2006) to 50-fold (Peloquin et al., 2008). T-type isoforms further complicate the issue: different isoforms can have temperature factors that either speed up or slow down the kinetics. For instance, when passing from room to physiological temperatures, the Ca3.3 isoform has a closing rate ∼50% faster (Iftinca et al., 2006), but the Ca3.1 isoform becomes ∼15% slower. For simplicity in our model, the same temperature factor was adopted for all VGCC subtypes.

VGCC currents

Request a detailed protocol

The VGCC currents are estimated using the GHK (Equation 14), as follows:

(16) IT=γTΦCaOT
(17) IR=γRΦCaOR
(18) IL=γLΦCa(OL1+OL2)

Table 9 presents the parameters to model the VGCC channels. VGCC rates and temperature factors are shown in Figure 17.

Table 9
VGCC parameters.

The number of VGCC was set to 3 to reproduce the calcium dynamics measured with a dye as in Figure 12 (Tigaret et al., 2016).

VGCC T-type conductanceγCaT=12103nSsame as Magee and Johnston, 1995
VGCC R-type conductanceγCaR=17103nSsame as Magee and Johnston, 1995
VGCC L-type conductanceγCaL=27103nSsame as Magee and Johnston, 1995
number of VGCCs3 for each subtype1–20 Higley and Sabatini, 2012

SK channel

Request a detailed protocol

The small potassium (SK) channel produces hyperpolarizing currents which are enhanced in the presence of intracellular calcium elevations. We included SK channels to incorporate a key negative feedback loop between spine calcium and voltage due to the tight coupling that exists between SK channels and NMDArs (Adelman et al., 2012; Griffith et al., 2016). Although there are a few publications on the single channel recording of SK channels (Hirschberg et al., 1998; Hirschberg et al., 1999) and at least one stochastic model of SK channel (Stanley et al., 2011), we chose to model SK channels deterministically. In tests, we found that this assumption had only a negligible impact on the outcomes of plasticity protocols (data not shown). Although SK channels can additionally be regulated by metabotropic glutamate receptors and muscarinic receptors (Tigaret et al., 2016), we did not include these regulatory steps in the model. The SK channel current was based on the description from Griffith et al., 2016 as follows:


There is little information on how temperature effects SK channel function, but van Herck et al., 2018 suggests a left-ward shift in the SK half-activation when changing from 37°C (hSK=0.38±0.02μM) to 25°C (hSK=0.23±0.01μM); that is a 65% decrease. Thus, to mimic temperature dependence of SK, we decrease the decay time of the SK hyperpolarizing current by a factor of two when passing from physiological to room temperature.


Table 10 presents the parameters to model the SK channel.

Table 10
SK channel parameters.
SK channel
Number of SK channelsNSK=15Lin et al., 2008
SK conductanceγSK=104nSMaylie et al., 2004
SK reversal potentialErevSK=90mVGriffith et al., 2016
SK half-activationhSK=0.333 μMGriffith et al., 2016
SK half-activation slopeσ=64 Griffith et al., 2016
SK time constantτSK=0.0063sGriffith et al., 2016

Enzymes - CaM, CaN, and CaMKII

To model enzyme dynamics, we adapted a monomeric CaM-CaMKII Markov chain from Chang et al., 2019 which was built on the model by Pepke et al., 2010. Our adaptation incorporates a simplified CaN reaction which only binds to fully saturated CaM, i.e. CaM bound to four calcium ions on its N and C terminals (see Markov chain in Figure 18). A consequence of the Pepke coarse-grained model is that calcium binds and unbinds simultaneously from the CaM terminals (N,C). We assumed a lack of dephosphorylation reaction between CaMKII and CaN since Otmakhov et al., 2015 experimentally suggested that no known phosphatase affects CaMKII decay time which is probably caused only by CaM untrapping (Otmakhov et al., 2015). This was previously theorized in the Michalski’s model Michalski, 2013, and it is reflected in Chang data (Chang et al., 2019; Chang et al., 2017). The structure of the corresponding Markov chain is shown in Figure 18.

Coarse-grained model of CaM, CaMKII, and CaN adapted from Chang et al., 2019 and Pepke et al., 2010.

Figure 18 is adapted from Figure 5 from Pepke et al., 2010. Reaction from the CaM-Ca reactions (first layer) are attributed to 2Ca release and binding from different CaM saturation states CaM2C (2Ca bound to terminal C), CaM2N (2Ca bound to terminal N), CaM0 (no calcium bound), CaM4(Ca bound to both C and N terminal). Note that CaN is allowed to bind only to fully saturated CaM. Activated CaN is represented by the state CaNCaM4. Reactions between the first (CaM-Ca reactions) and the second layer (KCaM-Ca reactions) represent the binding of free/monomeric CaMKII (mKCaM) (Pepke et al., 2010) to different saturation levels of CaM. Reactions within the layer KCaM-Ca represent the binding of calcium to Calmodulin bound to CaMKII (KCaM0, KCaM2C, KCaM2N, KCaM4). Transition of layer KCaM-Ca reactions to layer KCaM-phosphorylation represents CaMKII bound to CaM that became phosphorylated (PCaM states) (Pepke et al., 2010; Chang et al., 2017; Chang et al., 2019). PCaM can become self-phosphorylated (Autonomous layer with P and P2) and release CaM. Once the KCaM deactivates from autonomous states, it returns to free monomeric CaMKII (mKCaM). The CaMKII activity in this work represent the states (KCaM +PCaM + P + P2). See Chang et al., 2019 for further explanation on this system. CaNCaM4 represents the CaN activity. For graphical reasons, we could not show the complete list of reactions as given in Table 12.

Chang et al., 2019 data provides a high-temporal resolution fluorescence measurements for CaMKII in dendritic spines of rat CA1 pyramidal neurons and advances the description of CaMKII self-phosphorylation (at room temperature). We modified Chang’s model of CaMKII unbinding rates k2,k3,k4,k5 to fit CaMKII dynamics at room/physiological temperature as shown by Chang et al., 2017 supplemental files. Previous modelling of CaMKII Chang et al., 2019; Pepke et al., 2010 used a stereotyped waveform with no adaptation to model calcium. Our contribution to CaMKII modelling was to use calcium dynamics sensitive to the experimental conditions to reproduce CaMKII data, therefore, allowing us to capture physiological temperature measurements from Chang et al., 2017. Note that the CaMKII dynamic has two time scales and we capture only the fastest timescale which ends after stimulation ceases (at 60 s). The slowest dynamic occurs at the end of the stimulus, close to the maximum (Figure 19a). This may be caused by the transient volume increase in the dendritic spine as measured by Chang et al., 2017. Table 11 shows the concentration of the enzymes and Table 12 shows the parameters to model enzymes reactions in shown in Figure 18.

CaMKII temperature changes in the model caused by 1Pre, 30 at 0.49 Hz with glutamate uncaging (no failures allowed), 1 mM Ca, 2 mM Mg, P4-7 organotypic slices from mouse hippocampus.

(a) CaMKII fluorescent probe lifetime change measured by Chang et al., 2017 for 25°C (blue) and 35°C (red). The decay time (τ) was estimated by fitting the decay after the stimulation (30 pulses at 0.49 Hz) using a single exponential decay, y=aetb.(b) Simulation of the CaMKII concentration change (with respect to the baseline) at 25°C in response to same protocol applied in the panel a. The simulations on the panels b, c, e, f show the mean of 20 samples. (c) Same as in panel b but for 35°C. (d) Estimated temperature change factor for the dissociation rates k2, k3, and k5 in the Markov chain in Figure 18. (e) Change in the concentration of the CaMKII states (25°C) which are summed to compose CaMKII change in the panel b. (f) Same as in panel e for 35°C with reference to the panel c.

Table 11
Concentration of each enzyme.
Enzyme concentrations
Free CaM concentration (spine)CaMcon=30 µMKakiuchi et al., 1982
Free KCaM concentration (spine)mKCaMcon=70 µMFeng et al., 2011; Lee et al., 2009
Free CaN spine concentration (spine)mCaNcon=20 µM>10 μM (estimation from Kuno et al., 1992)
Table 12
Parameters for the coarse-grained model published in Pepke et al., 2010 and adapted by Chang et al., 2019 and this work.

Pepke et al., 2010 rate adaptation for the coarse-grained model adapt(a,b,c,d,Ca)=abc+dCa. Refer to Figure 18 for definition of variables.

Coarse-grained model, CaM-Ca reactions
CaM0+2 Ca⇒ CaM2C CaM2N+2 Ca⇒ CaM4kf2C=adapt(kon1C,kon2C,koff1C,kon2C,Ca)Pepke et al., 2010
CaM0+2 Ca⇒ CaM2N CaM2C+2 Ca⇒ CaM4kf2N=adapt(kon1N,kon2N,koff1N,kon2N,Ca)Pepke et al., 2010
CaM2C⇒ CaM0+2 Ca CaM4⇒ CaM2N+2 Cakb2C=adapt(koff1C,koff2C,koff1C,kon2C,Ca)Pepke et al., 2010
CaM2N⇒ CaM0+2 Ca CaM4⇒ CaM2C+2 Cakb2N=adapt(koff1N,koff2N,koff1N,kon2N,Ca)Pepke et al., 2010
kon1C=5 μM1s11.2 to 9.6 μM-1s-1 (Pepke et al., 2010)
kon2C=10 μM1s15 to 35 μM-1s-1 (Pepke et al., 2010)
kon1N=100 μM1s125 to 260 μM-1s-1 (Pepke et al., 2010)
kon2N=200 μM1s150 to 300 μM-1s-1 (Pepke et al., 2010)
koff1C=50 s110 to 70 s-1 (Pepke et al., 2010)
koff2C=10 s18.5 to 10 s-1 (Pepke et al., 2010)
koff1N=2000 s11 . 103 to 4 . 103 s-1 (Pepke et al., 2010)

koff2N=500 s1
0.5 . 103 to > 1.103 s-1 (Pepke et al., 2010)
Coarse-grained model, KCaM-Ca reactions
KCaM0+2 Ca⇒ KCaM2C KCaM2N+2 Ca⇒ KCaM4kfK2C=adapt(konK1C,konK2C,koffK1C,konK2C,Ca)Pepke et al., 2010
KCaM0+2 Ca⇒ KCaM2N KCaM2C+2 Ca⇒ KCaM4kfK2N=adapt(konK1N,konK2N,koffK1N,konK2N,Ca)Pepke et al., 2010
KCaM2C⇒ KCaM0+2 Ca KCaM4⇒ KCaM2N+2 CakbK2C=adapt(koffK1C,koffK2C,koffK1C,konK2C,Ca)Pepke et al., 2010
KCaM2N⇒ KCaM0+2 Ca KCaM4⇒ KCaM2C+2 CakbK2N=adapt(koffK1N,koffK2N,koffK1N,konK2N,Ca)Pepke et al., 2010
konK1C=44μM-1s-1Pepke et al., 2010
konK2C=44μM-1s-1Pepke et al., 2010
konK1N=76μM-1s-1Pepke et al., 2010
konK2N=76μM-1s-1Pepke et al., 2010
koffK1C=33 s1Pepke et al., 2010
koffK2C=0.8 s10.49 to 4.9 s-1 (Pepke et al., 2010)
koffK1N=300s-1Pepke et al., 2010
koffK2N=20 s16 to 60 s-1 (Pepke et al., 2010)
Coarse-grained model, CaM-mKCaM reactions
CaM0+mKCaM⇒ mKCaM0kfCaM0=0.0038μM-1s-1Pepke et al., 2010
CaM2C+mKCaM⇒ mKCaM2CkfCaM2C=0.92μM-1s-1Pepke et al., 2010
CaM2N+mKCaM⇒ mKCaM2NkfCaM2N=0.12μM-1s-1Pepke et al., 2010
⇒ mKCaM4
kfCaM4=30 μM1s114 to 60 μM-1s-1 (Pepke et al., 2010)
mKCaM0⇒ CaM0+mKCaMkbCaM0=5.5s-1Pepke et al., 2010
mKCaM2C⇒ CaM2C+mKCaMkbCaM2C=6.8s-1Pepke et al., 2010
mKCaM2N⇒ CaM2N+mKCaMkbCaM2N=1.7s-1Pepke et al., 2010
⇒ CaM0+mKCaM
kbCaM4=1.5 s11.1 to 2.3 s-1 (Pepke et al., 2010)
Coarse-grained model, self-phosphorylation reactions
KCaM0⇒ PCaM0 KCaM2N⇒ PCaM2N KCaM2C⇒ PCaM2C KCaM4⇒ PCaM4k1=12.6s-1Chang et al., 2019
Fraction of activated CaMKIIF=CaMKII/mKCaMconsee Equation 19 (Chang et al., 2019)
PCaM0⇒ P+CaM0 PCaM2N⇒ P+CaM2N PCaM2C⇒ P+CaM2C PCaM4⇒ P+CaM4k2=0.33s-10.33s-1 ; adapted from Chang et al., 2019
P⇒mKCaMk3=40.17s-10.17s-1 adapted from Chang et al., 2019
P⇒P2k4=40.041s-10.041s-1 adapted from Chang et al., 2019
P2⇒Pk5=80.017s-10.017s-1 adapted from Chang et al., 2019
Calcineurin model, CaM-CaM4 reactions
CaM4+mCaN⇒mCaNCaM4kfCaN=10.75μM-1s-146μM-1s-1 (Quintana et al., 2005)
mCaNCaM4⇒CaM4+mCaNkbCaN=0.02s-1fit (Fujii et al., 2013) see Figure 20

We provide an example of equation describing the binding reaction associated to the state CaM0. Note that these equations are automatically generated by the code which implements Table 12.


The CaN concentration was chosen as the total concentration used in a previous model (Stefan et al., 2008) (1.6 μM) scaled by a factor of 12 due to a higher CaN concentration in dendritic spines (Goto et al., 1986; Baumgärtel and Mansuy, 2012) and taking into account the discrepancy between different CaN concentration studies (Kuno et al., 1992; Goto et al., 1986). Kuno et al., 1992 proposes 9.6 μg/mg (7.0+2.6 μg/mg for Aα and Aβ isoforms) for the catalytic subunit A of CaN (CaNA) in the hippocampus, while Goto et al., 1986 proposes 1.45 μg/mg (presumably for both isoforms). There is therefore a lack of consensus on CaN concentration in neurons, which seems to range between 1 and 10 μg/mg. However, models of CaN in spines (Stefan et al., 2008) use low values of CaN concentration (e.g. 1.6 μM), without adjusting for the fact that these values were estimated from measurements from the entire neuropil. There is little information on CaN concentration in spines, but Kuno et al., 1992 note that the concentration of CaN is 50% to 84% higher in synaptosomes than in neuronal nuclei. With this information in mind, we set CaN spine concentration 20 μM in our model. CaN was entirely activated through CaM for the following reason: CaNA is activated by calcium-CaM in a highly cooperative manner (Hill coefficient 2.8–3), whereas the activation of CaN by calcium (via CaNB) is at most 10% of that achieved with CaM (Stemmer and Klee, 1994). In other words, CaNA affinity for CaM is 16 nM to 26 pM (Creamer, 2020), while CaNB affinity for calcium ranges from 15 μM to 24 nM (Kakalis et al., 1995). CaN decay time was modeled using experimental spine CaN activity dynamics measured in Fujii et al., 2013.

The lack of reactions between CaN and CaMKII

Request a detailed protocol

The protein phosphatases responsible for CaMKII dephosphorylation have not been established unequivocally (Lisman, 1989). Our model of CaMKII is based directly on a quantitative model fit to FRET imaging data (Chang et al., 2017; Chang et al., 2019), which implicitly account for the effects of any ‘hidden’ phosphatases, absorbing their contribution into the decay rates of the CaMKII activity. As pointed out by Otmakhov et al., 2015, FRET sensor imaging of CaMKII activity unfortunately does not capture the identity of the phosphatases involved in the dephosphorylation of CaMKII. More specifically, Otmakhov et al., 2015 observed no significant changes in the decay constant of their CaMKII FRET sensor when selectively inhibiting PP1 and PP2A. Given that these two phosphatases are widely used in models to determine plasticity, we believe that our model is more aligned with data of CaMKII activity in vivo.

Yet, our decision to include CaN in the model was determined by the evidence supporting CaN as the strongest candidate for calcium-sensitive protein phosphatase in the brain (Baumgärtel and Mansuy, 2012). Furthermore, the central role of CaN in synaptic plasticity has been demonstrated both pharmacologically and with genetic manipulation (Onuma et al., 1998; Malleret et al., 2001).

Temperature effects on enzymatic activity

Request a detailed protocol

We included temperature factors in the coarse-grained model using Chang’s data (Chang et al., 2019), as shown in Figure 19. For CaMKII, we fit the modified dissociation rates of the phosphorylation states k2, k3, and k5 to match the data on relative amplitude and decay time using the following logistic function:


For CaN, we fit the Fujii et al., 2013 data at 25°C as seen in Figure 20a. However, since CaN-CaM dissociation rates at physiological temperatures were not reported, we set the temperature factor to CaN that fits the outcomes of the protocols that we aimed to reproduce. A reference value from the CaN-AKAP79 complex (Li et al., 2012) showed a Q10=4.46=(2.19s-1/9.78s-1), which is nearly the temperature factor used in our model for CaN. Therefore, both the association and dissociation rates are modified using the following logistic functions:

CaN temperature changes in our model caused by 1Pre, 100 at 20 Hz with glutamate uncaging (no failures allowed), 2 mM Ca, Mg-free, 11–13 days in vitro.

(a), Simulated CaN change (blue solid line) in response to the same stimuli of the CaN measurement from Fujii et al., 2013 RY-CaN fluorescent probe (green solid line). The decay time (τ) estimated from data (y=ae-tb) is 94.83 s (dashed purple line) and 82.66 s for our model (solid purple line). (b), Simulated CaN change for physiological temperature with decay time of 54.44 s. (c), Temperature change, ρfCaN and ρbCaN, applied to CaN association and dissociation rates.

Positioning of the plasticity regions

Request a detailed protocol

Tigaret et al., 2016 LTP protocols were used to set the LTP region and as a first approximation of the LTD region. See Figure 21 TopDudek and Bear, 1992, Dudek and Bear, 1993 and Inglebert et al., 2020 were used to further define the LTD region. See Figure 21 Middle and Bottom. We highlight further a few points. For simplicity, we positioned the right border of the LTD region at the left border of the LTP region. The bottom part of the geometrical readout, under 4 μM of CaMKII, does not code for any dynamics. Note that some protocols may also enter and leave the plasticity regions multiple times, for example, TBS in Figure 5b, protocols in between LTD/LTP region in Figure 3d and age related LTD in Figure 5f. Because of this, we created an integrate and leak variable instead of using the time spent for predicting plasticity (see next section). This way, only after a certain time spent in a region would the synaptic weights start to change, as in Figure 23h. The coordinates of the plasticity regions are given in the last two rows of Table 13.

Positioning the plasticity regions.

The figure shows how Tigaret et al., 2016, Dudek and Bear, 1992, Dudek and Bear, 1993 and Inglebert et al., 2020 contributes to define the plasticity regions. In summary, Tigaret et al., 2016 data was used to define the LTP region, and Dudek and Bear, 1992, Dudek and Bear, 1993, Inglebert et al., 2020 data were used to define the LTD region.

Table 13
Parameters of the plasticity readout.

The variables in this table were fitted as described in the section Positioning of the plasticity regions.

Leaking variable (a.u.)
Rise constant inside the LTD regionaD=100a.u.s-1
Rise constant inside the LTP regionaP=200a.u.s-1
Decay constant outside the LTD regionbD=210-2a.u.s-1
Decay constant outside the LTP regionbP=0.1a.u.s-1
Plasticity Markov chain
LTD rate time constanttD=18s
LTP rate time constanttP=13s
Half occupation LTPKP=1.3104a.u.
Half occupation LTDKD=8104a.u.
Plasticity regions (vertices determining the polygons)
LTP region (CaN, CaMKII)[6.35,1.4], [10,1.4], [6.35,29.5], [10,29.5]
LTD region (CaN, CaMKII)[6.35,1.4], [6.35,23.25], [6.35,29.5], [1.85,11.32]
[1.85,23.25], [3.76,1.4], [5.65,29.5]

Geometrical readout

Request a detailed protocol

We describe here the geometrical readout mechanism which allows for plasticity outcome assignment. First, we define the following variables which are representative of ‘active CaMKII’ and ‘active CaN’:


The first two equations of (20) represent the total activation of K and P components of CaMKII associated with CaM as defined in the code that accompanies Chang et al., 2019. The last equation is from Chang et al., 2019.

Calcium entry in the spine initiates a cascade of events that ultimately leads to long term plasticity changes. Specific concentrations of CaMKII and CaN trigger activation functions actD and actP when they belong to one of the two polygonal regions (P and D), termed plasticity regions in the main text:


The variables actD and actP act as low pass filters of CaMKII and CaN activities with some memory of previous passages in the respective plasticity regions. To specify the LTP/LTD rates, termed Drate and Prate, we use the activation functions, actD and actP, as follows:


The Markov plasticity chain (see Figure 22) starts with initial conditions NC=100, LTD=0 and LTP=0.

Plasticity Markov Chain.

Table 13 provides the parameters that define the boundaries of the plasticity regions (see Figure 21 and Figure 23b).

Plasticity readout for the protocol 1Pre2Post10, 300 at 5 Hz, from Tigaret et al., 2016.

(a) CaMKII and CaN activity in response to protocol 1Pre2Post10. (b) Enzymatic joint activity in the 2D plane showing LTP and LTD’s plasticity regions. The black point marks the beginning of the stimulation, and the white point shows the end of the stimulation after 60 s. (c) Region indicator illustrating how the joint activity crosses the LTP and the LTD regions. (d) The leaky activation functions are used as input to the LTP and LTD rates, respectively. The activation function has a constant rise when the joint-activity is inside the region, and exponential decay when it is out. (e) The LTD rate in response to the leaky activation function, actD, in panel d. Note that this rate profile occurs after the stimulation is finished (60 s). The joint-activity is returning to the resting concentration in panel A. (f) The LTP rate in response to the leaky activation function, actP, in panel D. (g) Outcome of the plasticity Markov chain in response to the LTD and LTP rates. The EPSP change (%) is estimated by the difference between the number of processes in the states LTP and LTD, LTP-LTD. (h) Normalized LTP and LTD rates (multiplied to their respective time constant, tD, tP) sigmoids. The dashed line represents the half-activation curve for the LTP and LTD rates. Note in panel d that the leaky activation function reaches the half-activation Kp=1.3e4.

Figure 23 shows how the readout works to predict plasticity for a single orbit. Figure 23a shows the enzyme’s activity alone which is combined to form an orbit as shown in Figure 23b. The region indicator of the respective orbit is shown in Figure 23c. Simultaneously, Figure 23d depicts the leaky activation actP and actD, which will define the rate of plasticity induction in Figure 23e and f. The rates in the plasticity Markov chain will not reset to 0 if the orbit leaves the readout. The plasticity Markov chain is shown in Figure 23g with the prediction outcome represented as a weight change (%). Figure 23h shows the rate, Prate and Drate, activation profile. The LTP activation rate is steep, meaning that orbits do not need to spend a long time inside the readout to promote LTP induction, while the LTD region requires five-fold longer activation times.

Appendix 1

Appendix 1—table 1
Synaptic plasticity protocol parameters.

To fit the data from publications displaying a parameter interval (e.g. 70–100), we used a value within the provided limits. Otherwise, we depict in parentheses the value used to fit to the data. Further information is available in the github code and Appendix 1—table 3. Some of these experiments did not control AP generation following EPSP stimulation: Mizuno et al., 2001, Dudek and Bear, 1992 Dudek and Bear, 1993. We modeled this effect, described below. In addition, Tigaret et al., 2016 used GABA(A)r blockers, which we modelled by setting the GABAr conductance to zero. Also, Mizuno et al., 2001 LTD protocol used a partial NMDA blocker, which we modelled by reducing NMDA conductance by 97%.

ExperimentPaperRepetitionsFreq (Hz)Age (days)Temp. (C)[Ca2+]o(mM)[Mg2+]o(mM)
STDPTigaret et al., 2016300556352.51.3
STDPInglebert et al., 2020100, positive delays0.321(30.45)1.3—3Ca/1.5
STDPInglebert et al., 2020150, negative delays0.314301.3—3Ca/1.5
STDPMeredith et al., 2003200.29—4524—2822
STDPWittenberg and Wang, 200670—100514—21(22.5–23)21
pre-burstTigaret et al., 2016300 and 9003 and 556352.51.3
FDPDudek and Bear, 19929001—5035352.51.5
FDPDudek and Bear, 199390017—35352.51.5
TBSDudek and Bear, 19933—4 (5) epochs4Pre at 100 Hz, 10 x at 5 Hz6, 14 and 17352.51.5
LFSMizuno et al., 20011—600112—28(26.5–31)2.40
Appendix 1—table 2
Comparison of recent computational models for plasticity highlighting the experimental conditions implemented and the experiments in the hippocampus and cortex they reproduce.

See Appendix 1—table 3 for additional details on experimental conditions of experimental works.

ModelGraupner and Brunel, 2012Ebner et al., 2019Jędrzejewska-Szmek et al., 2017Inglebert et al., 2020Chindemi et al., 2022This paper
Model frameworkExtension of Shouval et al., 2002Extension of Clopath et al., 2010 and modified from Hay et al., 2011Modified from Evans et al., 2013Extension of Graupner and Brunel, 2012Extension of Graupner and Brunel, 2012
TemperatureAbsentAbsentTemperature corrected ion channels (but not receptors)No temperature control needed (experiments covered are at 30° C)Only in the GHK equationTemperature is selectable on the dendritic spine level for ion channels, receptor and the calcium cascade
DevelopmentAbsentAbsentAbsentAbsentAbsentAge is selectable and implemented by GABAr and NMDAr switch and BaP maturation
aCSFAbsentAbsentAbsentPhenomenological changes in pre and post amplitudes to mimic extracellular calcium effectsIn vivo or in vitro changes for release probability, calcium reversal potential on NMDAr-induced calcium influxExternal Ca and Mg are selectable and affect release probability, reversal potential, NMDAr and VGCCs calcium current driving force
Plasticity experiments (quant. comparisons only)
Sjöström et al., 2001XXX
Wittenberg and Wang, 2006XX
Wang et al., 2005X
Sjöström and Häusser, 2006XXX
Nevian and Sakmann, 2006X
Letzkus et al., 2006X
Weber et al., 2016X
Fino et al., 2010X
Pawlak and Kerr, 2008X
Shen et al., 2008X
Inglebert et al., 2020XX
Markram et al., 1997X
Rodríguez-Moreno and Paulsen, 2008X
Egger et al., 1999X
Tigaret et al., 2016X
Dudek and Bear, 1992X
Dudek and Bear, 1993X
Mizuno et al., 2001X
Meredith et al., 2003X
O’Connor et al., 2005 (not included due to space)X
Bittner et al., 2017 (not included due to space)X
Appendix 1—table 3
Comparison of the experimental conditions for the differentdatasets in Appendix 1—table 2 covering experiments from neocortex, hippocampus and striatum.
Experimental workAge (days)[Ca2+]o (Mm)[Mg2+]o (Mm)Temperature (°C)
Sjöström et al., 200112–212.5132–34
Wittenberg and Wang, 200614–212124–30 or 30–34
Wang et al., 2005embryonic day 17–1832room
Sjöström and Häusser, 200614–212132–35
Nevian and Sakmann, 200613–152132–35
Letzkus et al., 200621–422134–35
Weber et al., 201649–771.251.3 or 0.132–35
Fino et al., 201015–212134
Pawlak and Kerr, 200819–222.5231–33
Shen et al., 200819–2621room
Inglebert et al., 202014–201.3–3.0Ca/1.530
Markram et al., 199714–162132–34
Rodríguez-Moreno and Paulsen, 20089–1422room
Egger et al., 199912–142134–36
Tigaret et al., 201650–552.51.335
Dudek and Bear, 1992352.51.535
Dudek and Bear, 19937–352.51.535
Mizuno et al., 200112–282.4Mg-Free (most experiments)30
Meredith et al., 20039–452224–28
O’Connor et al., 200514–212127.5–32
Bittner et al., 201742–632135

Data availability

The datasets at the basis of our model were obtained directly from the authors Tigaret et al., 2016 or extracted from graphs in the references in Appendix 1—table 1 using WebPlotDigitizer v 4.6 software (Rohatgi,A.)The dataset from Tigaret et al., 2016 is freely available upon request. Concerning the dataset Tigaret et al., 2016, Jack Mellor the person to contact. One can contact Cezar M. Tigaret as well. The model is available on GitHub at https://github.com/rveltz/SynapseElife (copy archived at Veltz, 2023).


  1. Book
    1. Destexhe A
    2. Mainen ZF
    3. Sejnowski TJ
    Kinetic models of synaptic transmission
    In: Koch C, Segev I, editors. Methods in Neuronal Modeling. MIT Press. pp. 1–25.