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 ([Mg²⁺]_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): Ca²⁺/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).

Figure 1

Download asset Open asset

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 Ca²⁺, 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 Ca²⁺ 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.

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.

Figure 2

Download asset Open asset

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.

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

Download asset Open asset

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.

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

Download asset Open asset

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.

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

Download asset Open asset

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.

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

Figure 6 with 2 supplements see all

Download asset Open asset

We attempted to reproduce Inglebert et al., 2020 findings by varying [Ca²⁺]_o and [Mg²⁺]_o with the following consequences for the model mechanisms (Materials and methods). On the presynaptic side, [Ca²⁺]_o modulates vesicle release probability. On the postsynaptic side, high [Ca²⁺]_o reduces NMDAr conductance (Maki and Popescu, 2014), whereas [Mg²⁺]_o affects the NMDAr Mg²⁺ 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 [Ca²⁺]_o and [Mg²⁺]_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 [Ca²⁺]_o=1.3 mM, Figure 6a shows that the STDP rule induced weak LTD for brief causal delays. At [Ca²⁺]_o = 1.8 mM, in Figure 6b, the model predicted a full-LTD window. At [Ca²⁺]_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 [Ca²⁺]_o spent most time inside the LTD region, but high [Ca²⁺]_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 [Ca²⁺]_o for Inglebert et al., 2020 experiments restrict LTP to Ca/Mg >1.5 and [Ca²⁺]_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 [Ca²⁺]_o = 2.2 mM (see Figure 6—figure supplement 1b), small fluctuations in [Ca²⁺]_o near this boundary could cause qualitative transitions in plasticity outcomes. For anti-causal pairings, increasing [Ca²⁺]_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 [Ca²⁺]_o and [Mg²⁺]_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 [Ca²⁺]_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 [Ca²⁺]_o and we observe similar transitions to Inglebert et al., 2020. However, the model’s transition for [Ca²⁺]_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 [Ca²⁺]_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 [Mg²⁺]_o controls a sliding threshold between LTD and LTP but not [Ca²⁺]_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 [Ca²⁺]_o and [Mg²⁺]_o Inglebert (Inglebert et al., 2020) use a non-physiological temperature (${30}^{\circ }C$) 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 [Ca²⁺]_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 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. $\sigma$) 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 $\sigma$) for protocols taken from Tigaret et al., 2016. With no jitter, $\sigma =0$, these protocols have different LTP magnitudes (corresponding to Figure 3) and become similar once $\sigma$ increases. The three protocols with a postsynaptic spike doublet gave identical plasticity for $\sigma =50$ ms.

Figure 7

Download asset Open asset

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 $\sigma =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 [Ca²⁺]_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 $\sim 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.

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.

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.

We write ${\displaystyle {\mathbf{1}}_A}$ for the indicator of a set $A$, meaning that ${\displaystyle {\mathbf{1}}_A(x)=1}$ if $x$ belongs to $A$ and zero otherwise.

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 D₀ and R₀ vesicles, respectively. The docked pool loses one vesicle each time a release occurs (Rudolph et al., 2015), with transition $D\to D-1$ (Figure 8). The reserve pool replenishes the docked pool with reversible transition $(R,D)\leftrightarrow (R-1,D+1)$. Finally, the reserve pool is replenished towards initial value $R_{\mathrm{0}}$ over a timescale ${\tau }_D^{r\mathbf{}e\mathbf{}f}$ via the transition $(R,D)\to (R+1,D)$.

Figure 8

Download asset Open asset

In addition to the stochastic dynamics in Table 1, each spike triggers a vesicle release $D\to D-1$ with probability $p_{rel}$:

which is a function of presynaptic calcium ${\displaystyle [C{a}_{pre}]}$ and extracellular calcium concentration, [Ca²⁺]_o, through the threshold $h({[C{a}^{2+}]}_o)$. The function $p_{r\mathbf{}e\mathbf{}l}$ 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 $C{a}_{pre}$ (see Figure 8a) based on a resource-use function (Tsodyks and Markram, 1997):

Upon arrival of the presynaptic spikes, $t\in (t_1,\mathrm{\cdots },t_n)$, we update $C{a}_{pre}$ according to the deterministic jump:

Finally, after ${\displaystyle C{a}_{pre}}$ has been updated, a vesicle may be released with probability $p_{rel}$ (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 [Ca²⁺]_o modifies the release probability dynamics (King et al., 2001), we fixed an initial release probability to 68% for [Ca²⁺]_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 [Ca²⁺]_o from 2.5 mM to 1 mM. Taking these into account, the decreasing sigmoid function in the Figure 8e depicts our [Ca²⁺]_o-dependent release probability model ($p_{rel}$).

Figure 8e shows that our $p_{rel}$ function is in good agreement with a previous analytical model suggesting that $p_{rel}({[C{a}^{2+}]}_o)\propto {({[C{a}^{2+}]}_o)}^{2}m{M}^{-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, [Ca²⁺]_o, [Mg²⁺]_o) to our model of NMDArs, which is the main postsynaptic calcium source. Furthermore, since our model of vesicle dynamics is simple, ${\tau }_{rec}$ in Equation 2 has two roles: to delay the $p_{rel}$ recovery caused by $C{a}_{pre}$ inactivation (enforced by ${\delta }_{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 $p_{r\mathbf{}e\mathbf{}l}$ is purely phenomenological as, in principle, the $\mathrm{[}C\mathbf{}{a}^{\mathrm{2}\mathrm{+}}\mathrm{]}$ jump parameter ${\delta }_{C\mathbf{}a}$ should depend on ${\mathrm{[}C\mathbf{}{a}^{\mathrm{2}\mathrm{+}}\mathrm{]}}_o$, but replacing the model with a more physiological model would not affect the results, since the measured dependence of release probability on ${\mathrm{[}C\mathbf{}{a}^{\mathrm{2}\mathrm{+}}\mathrm{]}}_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.

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.

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

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):

Despite the driving force to the resting concentration, $C{a}_{\mathrm{\infty }}=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 ${\tau }_{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):

Tigaret et al., 2016 experiments used the synthetic calcium-indicator dye Fluo-5f, which is well-modelled by a single Ca²⁺-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, $P_{Ca}$, 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.

where ${\displaystyle {\mathrm{\Phi }}_{Ca}(V_{sp},[{\mathrm{C}\mathrm{a}}^{2+}{]}_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 ${\displaystyle [{\mathrm{C}\mathrm{a}}^{2+}{]}_i}$. Note that for simplicity the calcium external concentration ${\displaystyle [{\mathrm{C}\mathrm{a}}^{2+}{]}_o}$ was kept fixed during the simulation and only altered by experimental conditions given by the aCSF composition.

Figure 12

Download asset Open asset