## Abstract

Theoretical studies have shown that calcium influx through N-methyl-d-aspartate (NMDA) receptors is a sufficient signal to account for various induction protocols of bidirectional synaptic plasticity, including spike time–dependent plasticity (STDP). The STDP curves obtained by these different models exhibits a form of spike time–dependent long-term depression that occurs when a presynaptic spike precedes the postsynaptic spike (pre-post LTD). We have previously proposed that this novel form of LTD can serve as an experimental test for the validity of these models. These calcium based theoretical models assumed deterministic calcium dynamics that reflect average properties of synaptic calcium currents. In this paper, we show that taking into account the stochastic properties of synaptic transmission significantly alters the form of STDP curves and may significantly reduce the magnitude of pre-post LTD.

## INTRODUCTION

Calcium influx through N-methyl-d-aspartate (NMDA) receptors is necessary for the induction of many forms of bidirectional synaptic plasticity. The hypothesis that the sign and magnitude of synaptic plasticity is determined by calcium concentration in the postsynaptic spines (Lisman 1989) has received some direct and indirect experimental support (Cho et al. 2001; Cormier et al. 2001; Cummings et al. 1996; Yang et al. 1999).

One form of NMDA receptor–dependent synaptic plasticity, spike time–dependent plasticity (STDP), is induced by repetitively firing pre- and postsynaptic neurons with a fixed time lag of Δ*t* between the pre- and postsynaptic spikes. Positive values of Δ*t* (pre-post) result in long-term potentiation (LTP), whereas negative values of Δ*t* (post-pre) result in long-term depression (LTD) (Bi and Poo 1998; Feldman 2000; Froemke and Dan 2002; Markram et al. 1997; Sjostrom et al. 2001; Zhang et al. 1998).

We previously proposed a model of calcium-dependent synaptic plasticity that can account for various induction protocols of synaptic plasticity (Shouval et al. 2002b). This theory is based on three fundamental assumptions. *1*) Calcium influx into spines determines the sign, magnitude, and rate of synaptic plasticity; a moderate elevation in calcium results in LTD, whereas a large elevation in calcium levels results in LTP (Fig. 1*A*). We also assume that the rate of plasticity is a monotonically increasing function of calcium (*Eqs. 1* and *3*). *2*) Calcium influx through NMDA receptors is sufficient to account for the various forms of synaptic plasticity. *3*) The back-propagating action potentials in the postsynaptic neuron have a slow afterdepolarization. Related models were proposed by other groups (Karmarkar and Buonomano 2002; Kitajima and Hara 2000). These calcium-based theories make the testable prediction that, at positive values of Δ*t*, larger than those that produce LTP, LTD will be induced (Fig. 1*B*). This prediction is consistent with results by Nishiyama et al. (2000) in hippocampal slices. Other experiments have not explored this region of Δ*t* extensively; however, from the available data, there is no strong support for a universal form of pre-post LTD (Bi and Poo 1998; Feldman 2000; Froemke and Dan 2002; Markram et al. 1997; Sjostrom et al. 2001).

Both the release of neurotransmitter from presynaptic terminals and the opening of postsynaptic receptors can be characterized as a stochastic process and both can contribute to the variability of calcium transients through NMDA receptors. A presynaptic action potential might or might not cause the release of neurotransmitter, and the amount of neurotransmitter released might vary as well. If there is release, the number of receptors that open also varies. We have calculated the mean and variance of calcium transients as a function of Δ*t*, assuming a simplified stochastic model of synaptic transmission (Yeung et al. 2004). The most significant results of this analysis are that, for a small number of postsynaptic NMDA receptors, the variability is significant, and that the CV of calcium transients increases monotonically as a function of Δ*t* for Δ*t* > 0 (Fig. 1*D*).

Previous theoretical predictions about the existence of a pre-post form of LTD were obtained using deterministic models of synaptic transmission (Karmarkar and Buonomano 2002; Kitajima and Hara 2000; Shouval et al. 2002b). In this paper, we take into account the stochastic nature of synaptic transmission and show that the shape of the STDP curve can be altered by properties of stochastic synaptic transmission. The finding that the shape of plasticity curves can critically depend on the stochastic nature of synaptic transmission is relevant to other nonlinear biophysical theories of synaptic plasticity as well.

## METHODS

### Unified model of calcium-dependent synaptic plasticity

The calcium-dependent synaptic plasticity equation used here is from Shouval et al. 2002b (1) where *w*_{i} is the synaptic weight of synapse *i*, and [*Ca*]_{i} is the calcium concentration at that synapse, λ is the decay constant, set here to 1, and the Ω function, as depicted in Fig. 1*A*, has the form (2) where sig(*x*,β) = exp(β*x*)/[1 + exp(βx)], and α_{1} = 0.4,α_{2} = 0.65, and β_{1} = β_{2} = 30. The calcium-dependent learning rate function, η, has the form (3) where *p*_{1}=1, *p*_{2}=0.6, *p*_{3}=3, and *p*_{4}=0.00001. These parameters are similar but not identical to those in Shouval et al. 2002b. This learning rate function has a sigmoidal form, which monotonically increases with *x*. The general form of *Eqs.* *1* and *2*, described in our previous paper (Shouval et al. 2002b), has been derived from lower-level biophysical models (Castellani et al. 2001; Shouval et al. 2002a, b). The form of Ω function is based qualitatively on the notion that a moderate rise in calcium produces LTD, whereas a large rise in calcium produces LTP. The actual parameters were set arbitrarily to obtain reasonable plasticity curves, given the calcium transients assumed here. The sigmoidal calcium-dependent learning rate is essential for potentiation during STDP because and equal rate of depression and potentiation would result in oscillations and not robust potentiation. In addition, this monotonically increasing leaning rate is consistent with observations that induction of LTP is faster that induction of LTD (Dudek and Bear 1992; Yang et al. 1999).

The NMDA receptor kinetics, which we have assumed, are characterized by a double exponential curve of the form where *t*_{i} is the time of the presynaptic action potential, *V* is the postsynaptic voltage, and the function *B*(*v*) represents the voltage dependence of the NMDA receptors (Jahr and Stevens 1990; Shouval et al. 2002b). We use here the approximate reversal potential equation instead of the more precise Goldamn-Hodgkin-Katz equation (Johnston and Wu 1995), with *V*_{r} = 130 mV for calcium ions. Here we chose *I*_{f} = 0.75, τ_{f} = 50 ms, and τ_{s} = 150 ms such that the mean time constant is τ_{av} = *I*_{f} × τ_{f} + (1 − *I*_{f}) × τ_{s} = 75 ms, the same time constant as assumed below in the calculations of the CV of calcium transient in the stochastic case. The parameter *G*_{NMDA} is not fixed as in the stochastic case, but is chosen, as explained below, from a distribution to mimic the stochastic calcium dynamics. When we simulate the deterministic case (Fig. 1*B*), we used a fixed value *G*_{NMDA}= 1/325 μM/(ms · V). The voltage dependence of NMDA receptor calcium conductance is as in Shouval et al. 2002b. Calcium concentration is calculated by a first-order differential equation of the form (4) with a calcium time constant of τ*Ca* = 25 ms (Sabatini et al. 2002).

The back-propagating action potential has a double exponential form (Shouval et al. 2002b): *BPAP*(*t*) = *V*_{f}^{bp}exp(−*t*/τ_{f}^{bp}) + *V*_{s}^{bp}exp(−*t*/τ_{s}^{bp}), here we chose *V*_{f}^{bp} = 60 mV, *V*_{s}^{bp} = 25 mV, τ_{f}^{bp} = 2 ms, and τ_{s}^{bp} = 60 ms.

### Stochastic calcium dynamics

In previous papers, calcium dynamics were simulated as described in the section above, assuming deterministic calcium dynamics. However, it is better to model synaptic transmission as a stochastic model since given a presynaptic action potential the presynaptic terminal may fail to release neurotransmitter, and if it does release neurotransmitter, both the number of receptors that open and their duration in the open state can vary. In Fig. 1*C*, we show schematically a histogram of peak calcium transients. The peak at zero is due to release failures from the presynaptic terminal, and the distribution at positive calcium values is due to variability of calcium transients given that there was a release. We have previously calculated analytically both the means and variability of calcium transients invoked by different values of Δ*t*, assuming a simplified model of NMDA receptors (Yeung et al. 2004). These calculations take into account presynaptic variability due to release failures and postsynaptic variability modeled by a simple Markov model of the NMDA receptor. As the number of NMDA receptors increases, the relative postsynaptic contribution to this variability decreases. In Fig. 1*D*, we show the calculated CV of the peak calcium transients, given that there was a presynaptic release, for 10 postsynaptic NMDA receptors. Anatomical experimental evidence of single postsynaptic spines indicates that the number of postsynaptic NMDA receptors is of the order 10–20 (Racca et al. 2000). Here we assume that the time constant of the NMDA receptor is 75 ms. It is important to note that the Markov model we have used for this analysis is significantly simplified with respect to more realistic models that can account for most of the complexity of NMDA receptor dynamics. However, simulations of calcium dynamics using more complex models of NMDA receptors yield qualitatively similar results (see web appendix A^{1} ). Intuitively, this phenomenon arises because the number of NMDA receptors that are still open when the back propagating action potential (BPAP) arrives at the synapse and falls as the time between the presynaptic release and the postsynaptic spike (Δ*t*) increases, resulting in an effectively smaller number of NMDA receptors as Δ*t* increases. This fall in the effective number of NMDA receptors as Δ*t* increases is therefore likely to hold for reasonable models of NMDA receptors. At values of Δ*t* larger than those depicted in Fig. 1*D*, the value of the CV drops abruptly to levels equal to those found at large negative values of Δ*t*. Therefore at large values of positive Δ*t*, larger than those simulated here, the linear fit (Fig. 1*D*, line) fails.

The shape of calcium transients in this paper is determined by the deterministic calcium dynamics outlined above. However, the magnitude of the transients is modulated by a random process by which we choose the maximal chord conductance parameter *G*_{NMDA}. The variance of this random process is determined for each Δ*t* by the results of the theoretical analysis. The stochastic process for choosing *G*_{NMDA} has two phases. First, we emulate the effect of failures of presynaptic release by randomly setting *G*_{NMDA} to zero, with a probability *P*_{failure}= 1 − μ = 0.5, where μ is the probability of release. It must be noted that, in our analysis, we do not differentiate between a release failure or a failure to bind postsynaptic receptors when there was a release. Next, if there is a release, we chose the maximal chord conductance, *G*_{NMDA} randomly from a gamma distribution with a mean of 1/325 μM/(ms · mV), and a variance chosen differently for every Δ*t* and the number of NMDA receptors, *Z*, according to linear fit to the theoretical derivation. The linear fits used in this paper are CV(Δ*t* > 0) = 0.095 + 0.0045Δ*t* and CV(Δ*t* ≤ 0) = 0.095 − 0.00067Δ*t*. We also limit the maximal value of *G*_{NMDA} to be smaller than the mean NMDA receptor conductance times the number of NMDA receptors *Z*; this has a small effect on the results for small *Z*. When this procedure is carried out, the CV of the peak calcium transients matches the linear fit to the calculated CV (Fig. 1*D*, line). Note that we do not change the conductance as a function of *Z*, only the relative variance. This implies that we are scaling the single NMDA receptor conductance parameter inversely to *Z* to approximately preserve the form of the plasticity curves, without changing other parameters.

## RESULTS

Our simulations of STDP with stochastic synaptic transmission are based on previous theoretical analysis (Yeung et al. 2004). The key assumptions derived from this analysis aredepicted in Fig. 1*D*. Most significantly, for positive values of Δ*t*, the CV significantly and monotonically increases as Δ*t* increases. This is quite surprising since two different values of Δ*t*, one positive and one negative, that produce the same mean peak transient result in significantly different variability. An intuitive explanation for this result is that as Δ*t* increases (for Δ*t* > 0), the number of NMDA receptors that are still bound to the neurotransmitter when the BPAP arrives decreases; thus the effective number of NMDA receptors is smaller, and the relative level of fluctuations is larger. This intuitive explanation suggests that the qualitative form of the dependence of CV on Δ*t* does not depend on the exact details of the Markov model used for the NMDA receptors, and indeed, simulations using more complex Markov models of NMDA receptors produce qualitatively similar results (See web appendix A^{1}). The CV of the calcium transients is inversely proportional to the square root of the number of NMDA receptors at each spine (CV ∝ 1/), and as the number of NMDA receptors increases, the relative postsynaptic contribution to the variance of calcium transients decreases. We have used these analytical results about the dependence of CV on Δ*t* and the number of postsynaptic NMDA receptors at each spine (*Z*) as a basis for a stochastic implementation of the unified plasticity model.

Simulations of STDP with the unified model (Shouval et al. 2002b), using deterministic calcium dynamics, result in STDP curves that exhibit a pre-post form of LTD (Fig. 1*B*). In Fig. 2, we show the effect of stochastic synaptic transmission on the shape of the STDP curve. For a small number of NMDA receptors (*Z* = 10, Fig. 2*A*) the pre-post form of LTD is nearly eliminated, whereas the post-pre form of LTD is not. This happens because the CV at values of Δ*t*, which produce pre-post LTD in the deterministic model (dashed black line), are significantly larger than for values of Δ*t* that produce post-pre LTD. For example: CV(60)/CV(−10), given that there is a release, is 3.6. Since the omega function (Fig. 1*A*) is nonmonotonic in the LTD range, a large variability of calcium transients implies that some of the transients can actually produce LTP. Due to the calcium dependent learning rate η, these instances of LTP have a more significant effect than the instances of LTD. Experimental evidence (Racca et al. 2000) indicates that the number of NMDA receptors in postsynaptic spines is 10–20. To compare the STDP curve for *Z* = 10 to experimental results, we have fit these data points to exponential curves; we obtain time constants τ = 14 ms for Δ*t* > 0 and τ = 57 ms for Δ*t* < 0. These results are comparable with experimental results (Feldman 2000; Froemke and Dan 2002), although we have not made an attempt to use biophysical parameters appropriate for each system.

As the number of postsynaptic NMDA receptors is increased, the variability, arising from postsynaptic sources, is decreased. Therefore as *Z* increases (Fig. 2, *B–D*), the shape of the plasticity curve approaches that of the deterministic curve. However, even at a very large *Z* (Fig. 2*D*, *Z* = 3,200), a case in which the postsynaptic variability is negligible, the shape of the STDP curve is not identical to the deterministic STDP curve. This happens because of the presynaptic release failures, which have the primary effect of reducing the actual number of presynaptic events, thus reducing the total magnitude of synaptic plasticity.

## DISCUSSION

Several biophysical models of calcium-dependent synaptic plasticity make the prediction that there must be an LTD window at positive values of Δ*t* (Karmarkar and Buonomano 2002; Kitajima and Hara 2000; Shouval et al. 2002b). This seems like a robust prediction, independent of the precise parameters of the model, since it stems from continuity; as Δ*t* increases, the magnitude of calcium influx through NMDA receptors monotonically and continuously decreases, thus it seems that it must pass through the LTD region before plasticity is completely eliminated at large values of Δ*t*. However, these calcium-based models assume deterministic calcium dynamics that represent average values of calcium transients.

We found that, due to stochastic properties of NMDA receptors, the pre-post form of LTD predicted by these calcium-dependent models is not a robust prediction of such theories. Instead we propose that, in systems that have a large number of NMDA receptors, and consequently a smaller relative variability of calcium transients, the pre-post form of LTD should be observed. However, if the number of postsynaptic NMDA receptors is small and the variability large, the pre-post form of LTD should have a significantly reduced magnitude.

Our prediction that the CV of peak calcium transients increases as a function of Δ*t* (for Δ*t* > 0) could be tested in experiments in which calcium transients in spines are imaged (Koester and Sakmann 1998; Nevian and Sakmann 2004; Nimchinsky et al. 2004; Sabatini et al. 2002). Since the interesting quantity to be experimentally examined is the variability of calcium transients, results would be meaningful if the variability of the method of measurement is small, the number of trials for each synapse large, and estimated variability should be over a single synapse rather than over several synapses, because that will obscure the single spine variability by adding between spine variability.

In the extreme cases (*Z* < 10) where the pre-post form of LTD is nearly eliminated, most pairs of spikes still fall within the LTD region of the Ω function. However, the effect of those spikes that fall in the LTP region is larger because we have assumed that the learning rate increases with intracellular calcium levels. This assumption arises from lower level biophysical models (Castellani et al. 2001; Shouval et al. 2002a) and is necessary for avoiding oscillations during the induction of spike time–dependent LTP. It is also consistent with several experimental observations, first, rate-based induction of LTD is typically induced by 900 (Dudek and Bear 1992) spikes over 15 min, whereas LTP is much faster and requires less spikes; furthermore, when synaptic plasticity is obtained by direct control of intracellular calcium levels (Yang et al. 1999), LTP is obtained by a brief large increase in calcium levels, whereas LTD is induced by a moderate and sustained increase in calcium levels.

In the original deterministic model, we assumed calcium influx only through NMDA receptors. This was done because influx through NMDA receptors was sufficient to qualitatively account for the different induction mechanisms (Shouval et al. 2002b), not because this is the only calcium source. However, imaging of calcium transients when the postsynaptic cell is depolarized (Sabatini et al. 2002) shows that NMDA receptors are the major source for calcium in spines. Furthermore, when the voltage dependence of NMDA receptors is eliminated using a magnesium-free solution, the difference in calcium transients between pre-post and post-pre conditions was eliminated (Nevian and Sakmann 2004). Nevertheless, inclusion of other calcium sources such as voltage-gated channels or calcium-induced calcium release from internal stores will have a quantitative effect on results of the model.

The parameters assumed in this paper are generic parameters and are not based on measurements in a specific experimental system. Therefore the consequences of the model should not be taken as quantitative predictions. For example, the LTP and LTD time windows in Fig. 2*A* depend on specific assumptions about the NMDA receptor and the BPAP time constants, and therefore should not be taken as quantitative predictions of what these time constants should be. This theory, if shown to be qualitatively correct, could be taken to the next quantitative level by using measured parameters of a specific system and making quantitative predictions for that specific system.

We have based our estimate of the number of postsynaptic NMDA receptors of anatomical measurements (Racca et al. 2000). Recent physiological measurements (Nimchinsky et al. 2004) indicate that the median number of NMDA receptors open during synaptic transmission is ∼3, given that there is release. In our analysis, we have not distinguished between probability of release and probability of receptors binding; however, for 10 postsynaptic NMDA receptors (*Z* = 10) and probability of release of 0.5, our model assumes a mean number of open receptors of 5, above the physiological estimate of 3. The physiological estimate of three bound receptors would imply even larger fluctuations than we have assumed, which would result in an even larger difference between the deterministic and stochastic plasticity curves (see Web appendix B^{1}).

Although the consequences of the unified plasticity model are qualitatively consistent with experimental results, recent experiments indicate that the basic assumptions of this theory might not be. It has been shown, in different systems, that spike time–dependent LTD requires additional receptors such as cannabinoid receptors in neocortex (Sjostrom et al. 2003) and metabotropic glutamate receptors in hippocampus (Nishiyama et al. 2000). These findings might imply that the basic assumptions of the unified plasticity model need to be revised.

The central finding of this paper is that the stochastic properties of synaptic transmission not only add variability to resulting plasticity curves, but can also alter their shape. The significance of this result is independent of whether the exact assumptions of the unified plasticity model are correct, and it potentially applies to all nonlinear biophysical models of synaptic plasticity.

## GRANTS

This work was supported in part by the Institute for Brain and Neural Systems at Brown University.

## Acknowledgments

We thank D. Feldman for reading an early version of the paper and the members of the Institute for Brain and Neural Systems at Brown for useful discussions.

## Footnotes

↵1 The Supplementary Material for this article (Appendix A and Appendix B is available online at http://jn.physiology.org/cgi/content/full/00504.2004/DC1.

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2005 by the American Physiological Society