## Abstract

Various forms of synaptic plasticity, including spike timing-dependent plasticity, can be accounted for by calcium-dependent models of synaptic plasticity. However, recent results in which synaptic plasticity is induced by multi-spike protocols cannot simply be accounted for by linear superposition of plasticity due to spike pairs or by existing calcium-dependent models. In this paper, we show that multi-spike protocols can be accounted for if, in addition to the dynamics of back-propagating action potentials, stochastic synaptic dynamics are taken into account. We show that a stochastic implementation can account for the data better than a deterministic implementation and is also more robust. Our results demonstrate that differences between experimental results obtained in hippocampus and visual cortex can be accounted for by the different synaptic and dendritic dynamics in these two systems.

## INTRODUCTION

Activity-dependent modification of synaptic connections and strengths is considered a physiological substrate of learning and memory (Bear and Rittenhouse 1999; Morris 2003). Various induction paradigms have been used over the years to induce long-term bidirectional synaptic plasticity (Bliss and Lomo 1973; Dudek and Bear 1992; Feldman et al. 1999). It has recently been shown that the precise timing of pre- and postsynaptic spikes have a significant impact on the sign and magnitude of long-term plasticity (Bi and Poo 1998; Markram et al. 1997). Presynaptic spikes followed closely by postsynaptic spikes generally produce long-term potentiation (LTP), whereas postspikes leading prespikes result in long-term depression (LTD) (Bi and Poo 1998; Feldman 2000; Froemke and Dan 2002; Markram et al. 1997; Sjostrom et al. 2001). Many forms of bidirectional synaptic plasticity are calcium dependent: experimental studies have shown that slow and moderate increase in [Ca] lead to depression, whereas large increase in [Ca] results in potentiation (Cho et al. 2001; Cormier et al. 2001; Cummings et al. 1996; Yang et al. 1999). Various induction protocols of synaptic plasticity, including spike-timing-dependent plasticity, can be qualitatively accounted for by calcium-dependent models of synaptic plasticity (Abarbanel et al. 2003; Shouval and Kalantzis 2005; Karmarkar and Buonomano 2002; Shouval et al. 2002a).

In a natural environment, the pre- and postsynaptic spiking patterns at a typical synapse are much more complex than simple pre-post spike pairs. As a first step toward understanding the underlying mechanisms of plasticity during natural spike trains, experimental studies are using simple extensions of spike pair paradigms, such as spike triplets and quadruplets (Froemke and Dan 2002; Wang et al. 2005). Naïve attempts to use linear superposition of plasticity events from spike pairs cannot account for plasticity induced by multispike paradigms. This is not surprising because the underlying dynamics leading to induction of synaptic plasticity, including calcium transients and the biochemical induction networks, are nonlinear. Additionally, recorded plasticity from multi-spike induction paradigms was significantly different in hippocampal and visual cortical preparations (Froemke and Dan 2002; Wang et al. 2005). On the basis of these experiments, it has been suggested that the timing of the first spike seems to dominate the total response in visual cortical slices, whereas potentiation appears to dominate in hippocampal cultures.

It is reasonable to assume that short-term synaptic dynamics, exhibited in all brain regions, can affect the induction of long-term plasticity. These short-term dynamics, however, can be quite different in different brain regions and even differ between different cell types in the same region. In the visual cortex, neurons usually demonstrate a decrease in the average excitatory postsynaptic potential (EPSP) size over paired presynaptic stimulations. This paired-pulse depression recovers with a time constant of 100–200 ms (Philpot et al. 2001; Tsodyks and Markram 1997; Varela et al. 1997). In hippocampal neurons, such paired-pulse depression (PPD) also exists but with much longer recovery time constant (1–3 s) (Debanne et al. 1996; Dobrunz 2002; Dobrunz and Stevens 1997; Murthy and Stevens 1999; Pyle et al. 2000), and in addition they also exhibit short-term synaptic facilitation, which is characterized by an increase in average response magnitude following multiple presynaptic stimuli (Debanne et al. 1996; Dobrunz and Stevens 1997). Moreover, the stochastic nature of neurotransmitter release and binding and postsynaptic receptor response, results in synaptic transmission which is not only dynamic, but also quite variable.

Spike-timing-dependent synaptic plasticity requires dendritic back-propagating action potentials (BPAPs) (Magee and Johnston 1997); therefore the shape and dynamics of BPAPs can influence plasticity induced by multi-spike protocols. Unlike somatic action potentials, BPAPs exhibit altered magnitudes and shapes during bursts of action potentials. For example, trains of BPAPs elicited in hippocampal neurons exhibit frequency-dependent depression in spike amplitudes (Colbert et al. 1997) due to inactivation of sodium channels (Colbert and Johnston 1998).

Computational modeling, with its ability to precisely monitor and manipulate parameters in complex systems, is well suited to examine the roles that various complex physiological processes might play in synaptic plasticity. For the results to be physiologically relevant, a biophysically based model of plasticity must be used. We use the calcium dependent model of Shouval et al. (2002a) as the basis for the computational study in this paper.

Although calcium-based models can qualitatively account for many properties of synaptic plasticity, they are unable to adequately account for plasticity induced by multi-spike protocols without addition of a questionable physiological assumption (Shah et al. 2006). Previous models used deterministic average behavior to approximate the physiological dynamics; processes like presynaptic transmitter release and postsynaptic binding, however, are more adequately described as stochastic processes. We have previously shown that deterministic approximations of such stochastic processes can be misleading (Shouval and Kalantzis 2005) due to the nonlinear nature of the underlying systems. Thus it is important to examine whether a deterministic model can accurately describe system behavior and to quantify the impacts of a stochastic approach on predicted plasticity.

In this study, we show that our calcium-based model of synaptic plasticity, with the addition of system-specific physiological parameters, can account for multi-spike plasticity in both visual cortex and hippocampus. We further show that our stochastic implementation of presynaptic transmitter release can more accurately account for available experimental data and is more robust to changes in parameter values than a deterministic implementation. In addition, we show how varying the short-term synaptic dynamics affects the resulting long-term plasticity and identify temporal stimulus regions where the sign of synaptic plasticity (LTP versus LTD) is sensitive to parameter variation. Such differences in short-term synaptic dynamics might play a major role in explaining how apparently different forms of plasticity observed in different systems are related. In contrast to what is ostensibly suggested by present experimental data, our simulations make testable predictions that the overall patterns of plasticity in the two systems (hippocampus and visual cortex) may be qualitatively similar with major quantitative differences arising from variations in short-term synaptic and neuronal dynamics.

## METHODS

### Calcium-dependent plasticity (CaDP) model

The model that serves as the basis for our simulation has been described previously (Shouval et al. 2002a), and here we briefly review the basic model assumptions.

##### ASSUMPTION 1.

The level of [Ca] in the synaptic spine determines the sign and magnitude of synaptic plasticity (1) where *w _{i}* is the synaptic weight of synapse

*i*, and [Ca]

_{i}is the calcium concentration at that synapse, and the Ω function, as depicted in Fig. 1

*A*, has the form (2) where sig(

*x*,β) = exp(β

*x*)/[1 + exp(β

*x*)], α

_{0}= 0.33333, α

_{1}= 0.22, α

_{2}= 039, β

_{1}= 80, and β

_{2}= 40 in both hippocampal and visual cortex simulations. The calcium dependent learning rate function (η) is a monotonically increasing function of [Ca] and has the form (3) where

*p*

_{1}= 1.0,

*p*

_{2}= 0.28,

*p*

_{3}= 3.0,

*p*

_{4}= 1

*e*-5.

##### ASSUMPTION 2.

The primary source of calcium influx is through NMDAR. The kinetics of NMDA current is characterized by a double exponential of the form (4) where *t*_{i} are times of presynaptic APs and *V* is the postsynaptic membrane potential. The voltage dependence of NMDA channels is described by *B*(*V*) (Jahr and Stevens 1990; Shouval et al. 2002a), and parameters used are *I*_{f} = 0.7, *I*_{s} = 1 − *I*_{f}, *t*_{f} = 32 ms, *t*_{s} = 160 ms (Philpot et al. 2001) with approximate calcium reversal potential *V*_{r} = 130 mV and a conductance *G*_{NMDA} of −1.25 mM/[s*mV]. The change of *G*_{NMDA} from the deterministic model (−3.448 in Shah et al. 2006) to the current value is needed to maintain the resulting NMDA current unchanged in the stochastic model because in the stochastic version there is either a vesicle release or no release, thus the *p*_{0} term in the deterministic version (*Eq. 2* in Shah et al. 2006) is incorporated into *G*_{NMDA} in *Eq. 4*. The resting membrane potential is set at −65 mV. Postsynaptic resources are assumed to be completely saturated by a single quantal release of neurotransmitter.

Calcium concentration is determined by a first order ODE (Shouval et al. 2002a) (d[Ca(*t*)]/d*t*) = *I*_{NMDA}(*t*) − [Ca(*+*)]/τ_{Ca}, with time constant τ_{Ca} = 25 ms (Sabatini et al. 2002).

##### ASSUMPTION 3.

The back propagating action potential (BPAP) that contributes to Ca influx has a slow tail. The dendritic BPAPs consist of a fast initial depolarization and a slowly recovering tail and have a double exponential form with *V*_{BPAP} = 80 mV, *I*_{b,f} = 0.7, *I*_{b,s} = 0.3, τ_{b,f} = 2 ms, and τ_{b,s} = 30 ms.

##### ASSUMPTION 4.

Synapses exhibit short-term synaptic plasticity.

###### Synaptic depression.

The probability of presynaptic transmitter release undergoes a depression on a vesicle release and recovers exponentially with a probability of docked vesicle release *p*_{dr0} of 0.3 and depression recovery time constant τ_{r} of 141 ms for visual cortex (Philpot et al. 2001) and p_{dr0} of 0.19 and τ_{r} of 1,000 ms for hippocampus (Dobrunz and Stevens 1997; Murthy and Stevens 1999; Pyott and Rosenmund 2002; Stevens and Tsujimoto 1995) unless otherwise stated. The stochastic implementation of synaptic depression is shown in the pseudo code in the next section.

###### Synaptic facilitation.

In addition to synaptic depression (Shah et al. 2006; Shouval et al. 2002a), we implemented synaptic facilitation to investigate suitability of our model for use in the hippocampus, where short-term facilitation has been observed (Chen et al. 2004; Debanne et al. 1996; Dobrunz and Stevens 1997). We use a simple model of facilitation that allows the probability of release *p*_{r} to increase by a fixed percentage (γ = 0.8) immediately after a presynaptic spike at time *t*_{0}, then recover exponentially to the base value (*p*_{dr0}) with a fixed time constant (τ_{F} = 100 ms) Note that the synapse is facilitated (*p*_{r} increased) even when there is no vesicle release due to a presynaptic spike. This form may lead to release probabilities >1.0 for some choices of *p*_{dr0} and γ, and this potential problem has been take care of during the implementation by setting an upper bound such that *P*_{r} ≤ 1.

Postsynaptic membrane voltage components sum linearly and are dominated by BPAPs. EPSPs are assumed to be relatively small (≈1 mV) with exponential rise and fall times of 5 and 50 ms, respectively (Shah et al. 2006; Shouval et al. 2002a).

##### ASSUMPTION 5.

The BPAP amplitude is also depressed and recovers exponentially after a postsynaptic spike. This process is characterized by a single exponential with a magnitude of 0.3 (i.e., 30% depression of BPAP magnitude immediately subsequent to each spike) and a time constant of 35 ms in the hippocampal simulations (Fig. 2*D*) and 50% depression and 55 ms recovery constant in visual cortex simulations. The deeper magnitude of BPAP depression in visual cortex is needed to produce a depression region in between the two potentiation regions (cf. ⇓Fig. 4).

### Stochastic implementations

Our stochastic simulations are based closely on the deterministic model described previously. The primary difference is inclusion of a finite pool of presynaptic vesicles. Vesicles from this pool are released and refilled independently, after presynaptic spikes with probability *p*_{dr0} (probability of docked vesicle release). In addition to allowing for short-term depression (a vesicle that has released and not yet recovered cannot be expected to release again), this formulation allows for graded presynaptic release (a single or multiple vesicle releases) as well as complete failures (Fig. 2*C*).

Our basic algorithms for stochastic implementation of PPD and PPF are The function nextRandomNum() returns a random number between 0 and 1.

We did not simulate stochastic binding and opening of the postsynaptic receptors, although our previous work has suggested that a stochastic postsynaptic binding affects the shape of the spiking-timing-dependent plasticity (STDP) curve (Shouval and Kalantzis 2005). Instead, for simplicity, we assumed that a release of a single vesicle saturates all postsynaptic binding sites. Stochastic postsynaptic binding may reduce the size of the second depression in the region of positive pre-post interval in the STDP curve (cf. Fig. 1). However, we do not expect the existence or lack of this depression region to affect the main findings presented in this paper.

We assumed an initial synaptic weight of 1/3 for all simulation results except for those presented in Fig. 3, in which an initial synaptic weight of 0.25 was used (Shah et al. 2006). The choice of initial synaptic weight does not qualitatively change the global synaptic plasticity pattern, it only affects the final results quantitatively.

### Estimates of parameters for stochastic transmitter release

The parameters involved in the stochastic transmitter release are the probability of docked vesicle release (*p*_{dr0}) and the number of vesicles per release site (*N*). In determining the parameters, we first set *p*_{dr0} based on PPD data, then estimated *N* based on *p*_{dr0} and the probability of failure per synaptic contact (*p*_{F}) with the equation (1 –*p*_{dr0})^{N} = *p*_{F}. Because the experimental data are collected under diverse conditions and results vary significantly even in the same experiment (or same preparation), our choice of parameters may be considered more appropriately to be a reasonable starting point, rather than a good estimate of the parameters. For visual cortex simulations, we chose a *p*_{dr0} of 0.3, based on paired-pulse facilitation data from visual cortex neurons (Philpot et al. 2001; Shah et al. 2006), and a failure rate *p*_{F} of 0.51, based on data from L4 to L2/3 connection in the barrel cortex (Feldmeyer et al. 2002). Finally *N* is calculated to be 2 based on values of *p*_{dr0} and *p*_{F}. Increasing the number of vesicles while keeping the same failure rate will produce a pattern that is less random, and has less depression. For hippocampal simulations, we decided to keep the same number of vesicles as used in the visual cortex simulations. Because the effects of depression and facilitation overlap, PPD could not be used to directly set the probability of release. Instead, we chose a *p*_{dr0} value of 0.19, which corresponds to a realistic failure rate of 0.65 (Dobrunz and Stevens 1997) and reflects the fact that in general hippocampal data (Dobrunz 2002; Dobrunz and Stevens 1997; Sudhof 2000) yield lower estimates of *p*_{dr0} than data from the cortex. It should also be noted that in this study we examined plasticity at a single synapse, yet some experimental data were collected from connections with multiple synapses. We have taken this into account when extracting parameters from data obtained in those experiments (e.g., probability of release).

### Simulation and data analysis

All simulations were performed using a Java-based simulation package developed in this laboratory. This package uses a friendly graphic user interface that was originally developed for the object-oriented version of SNNAP [simulator for neural networks and action potentials (Cai et al. 2004)], with improved capacity of data compression and pseudo three-dimensional (3-D) display. The model parameters and simulation conditions are specified in cascaded text-based parameter files, which can easily be changed by using the graphic user interface or a text editor. The results of the simulation are displayed on screen and are updated continuously and can also be saved into files for further analysis. Simulations were performed on computers running Linux and Mac X operating systems. To simulate the plasticity (change in synaptic weight) at a particular triplet condition, we present the same triplet stimuli 100 times, delivered at a rate of 1 Hz. The simulation was usually run for 102 s with a time step of 0.1 ms (cf. Fig. 3) and synaptic weight at the end of the simulation was used to calculate the relative synaptic weight reported in this paper. A synaptic weight >1.0 indicates potentiation and a value smaller than 1.0 signals depression. To generate the pseudo 3-D pattern, we changed *t*_{1} and *t*_{2} at an interval of 3 ms, and run simulation for each condition. To obtain the average response patterns, we run simulations with 12 different random seeds and results are averaged for each data point. As in Shah et al. (2006), to quantify the difference between experimental data and simulation results, we calculated the mean absolute error (MAE) between the experimental and theoretical data points: . The values of theoretical data points used to calculate MAE are averages of 12 simulations.

## RESULTS

### Calcium-based model of synaptic plasticity with and without stochastic transmitter release

The calcium based theory of synaptic plasticity assumes that the concentration of calcium in dendritic spines determines the sign and magnitude of plasticity (Cho et al. 2001; Cormier et al. 2001; Cummings et al. 1996; Yang et al. 1999). The original form of the model (Shouval et al. 2002a) is able to reproduce the basic experimental results for spike pairs: depression when the postsynaptic spike leads the presynaptic spike and potentiation when presynaptic spike leads the postsynaptic spike (Shouval et al. 2002a) (cf. Fig. 1*C*). When simulating more complex spike patterns, like triplets (Fig. 1*D*), the original naïve form of the model is dominated by potentiation and exhibits plasticity that is incorrect in both sign (potentiation in regions of experimental depression or quiescence) and magnitude (overestimating levels in most potentiation regions) (Shah et al. 2006). Such prevalence of LTP might arise from the omission of short-term synaptic and neuronal dynamics. Neurons in neocortex exhibit depression of presynaptic vesicle release (Fig. 2, *A–C*) and are also likely to exhibit depression in the magnitude of BPAPs (Fig. 2*D*). Versions of the model incorporating these two factors, when properly tuned, can account for most of the experimental data (Shah et al. 2006).

The results of Shah et al. (2006) were produced with a deterministic model of short-term plasticity. However, many processes underlying neural plasticity, like presynaptic transmitter releases and postsynaptic binding, are stochastic. In Fig. 2, we show schematic representations of deterministic and stochastic PPD. In a PPD experiment, two presynaptic stimuli are presented and the postsynaptic voltage responses are recorded and averaged over many trials. In general, the size (or slope) of an EPSP in response to the second stimulus is smaller than that in response to the first stimulus; this is termed PPD (Fig. 2*A*). The difference in average response does not, however, indicate that the magnitude of successful EPSP events is different for the first and second spikes. As shown in Fig. 2*B*, individual trials will either succeed or fail to produce EPSPs. Successful events produce EPSPs with approximately the same distribution of magnitudes for the first and second presynaptic stimuli, but there are more failures in response to the second stimulus. The smaller average size of the second EPSP (Fig. 2*B*, - - -) deceptively suggests that PPD is a change in response magnitude when it is actually a change in the distribution of success/failure events. This subtle difference has significant implications for system behavior.

In the deterministic approach, parameters are derived from the average behavior. The stochastic approach, on the other hand, simulates synaptic transmission with EPSP dynamics and variability as observed in real single spines. Because the underlying processes are highly nonlinear, there is no reason to expect that deterministic simulations based on the average behavior will produce the same results as the stochastic approach in all regions.

In this paper, while we simulate the variability of the presynaptic stochastic release process, we use a deterministic approximation of the postsynaptic processes, because it is reasonable that the stochasticity of presynaptic release has a stronger effect on multi-spike induction protocols (see discussion).

### Simulations with deterministic versus stochastic synaptic release: visual cortex

To investigate how a stochastic implementation of the presynaptic transmitter release affects multispike plasticity induction, we implemented the stochastic version of the model for the case of pre-post-pre spike triplets (Fig. 3). For each temporal condition (1 combination of *t*_{1}, *t*_{2}), we present the spike triplet event 100 times at a frequency of 1 Hz. A triplet protocol consisting of two presynaptic and one postsynaptic spike (cf. Fig. 1*D*) can be plotted in three regions, depending on the values of *t*_{1} and *t*_{2}: pre-post-pre (*top left quadrant*), pre-pre-post (*top right quadrant*), and post-pre-pre (*bottom left quadrant*). The parameters of the stochastic version were converted directly from the previous deterministic version of the model (Shah et al. 2006), with the minimal changes necessary to the implement stochastic processes. In the deterministic version, the parameter *p*_{0,} which determines the maximal magnitude if PPD, is 0.3 and we can therefore assume the probability of failure is *p*_{F} = 1 − *p*_{0} = 0.7. In the stochastic version, we assumed that there are two vesicles and that each docked vesicle has a probability of release (*p*_{dr0}) of 0.163. This value was chosen so that the resulting probability of failure for a spine (i.e., no vesicle release on a presynaptic spike) of 0.7 is the same as in the deterministic case. We also adjusted the value of the NMDA conductance (*g*_{NMDA}) so as to approximately maintain the same average NMDA current (see methods).

Comparing with experimental data collected from visual cortex neurons by Froemke and Dan (2002) (shown in the bounded squares), the deterministic version of the model produces a plasticity pattern (depression and potentiation) that is broadly consistent with the experimental data (mean absolute error, MAE 1.14; Fig. 4*A* of Shah et al. 2006). However, it fails in the post-pre-pre region (*bottom left quadrant*), producing potentiation where the data show depression. The inability of the model to account for the experimental data in this region may be corrected by including an additional, strong PPD with a very fast recovery (Shah et al. 2006), but there is no experimental evidence supporting such a fast depression and recovery term.

The response pattern produced by the stochastic version of the model (Fig. 3*C*), however, matches the experimental data better. In particular, the region in the post-pre-pre quadrant (*bottom left quadrant*) that showed potentiation in the deterministic version of the model now shows depression consistent with experimental data. The success of the stochastic version of the model in the post-pre-pre region can be attributed to the very nature of the stochastic transmitter release. In the deterministic model, presynaptic release and resulting [Ca] transients are the same for each spike triplet presentation. Accordingly, synaptic weights change monotonically. In the stochastic model, on the other hand, the [Ca] transients have three different characteristic sizes, depending on whether there is a vesicle release on the first, the second presynaptic spike, or both (Fig. 3*A*). In addition, no [Ca] transient is produced when there is no vesicle release on both presynaptic spikes. Accordingly, the changes in the synaptic weight are nonmonotonic, sometimes increasing, and sometimes decreasing (Fig. 3*B*). In simulations of the stochastic model, the [Ca] transients and final synaptic weight also change from trial to trial, depending on the random seed being used even though other parameters are kept the same.

In Fig. 3*C*, *inset*, we show three [Ca] transients taken from a stochastic run of a point in the post-pre-pre region (−6, −12). The largest [Ca] transient corresponds to a triplet presentation that produced transmitter releases from both presynaptic spikes, and the two smaller [Ca] transients correspond to triplet presentations that produced transmitter releases from only one of the two presynaptic spikes. Although the [Ca] transients due to synaptic vesicle releases from both presynaptic events produce increases in synaptic weight (potentiation), the transients due to vesicle releases from only one of the two presynaptic spikes are smaller, producing a decrease in synaptic weight (depression). Due to the relative distributions of the different release events, the effect of depression events outweighs the effect of the potentiation events, producing an overall depression for this particular triplet condition. This example highlights the different mechanisms of plasticity induction in the deterministic and the stochastic versions of the model.

The stochastic model is more robust than the deterministic model; many different combinations of parameters are able to account for the basic plasticity pattern in the experimental data, and the final results are not as sensitive to the choice of some individual parameters. To show an example, we changed the probability of release in both the deterministic (*p*_{0}) and the stochastic (*p*_{dr0}) cases and examined their effects on the plasticity patterns. Although these different parameters are not equivalent, we can use the resulting probability of failure *p*_{F} to directly compare the models. Decreasing *p*_{0} in the deterministic case reduced both the magnitude of potentiation and the size of the potentiation zone (Fig. 3*D*). Increasing *p*_{0} increased the amplitude of potentiation and broadened the potentiation zone (Fig. 3*E*). In the stochastic case, despite a larger range of p_{F} compared with the deterministic example, changes in *p*_{dr0} have less effect on the magnitude of potentiation and do not change appreciably the area of the potentiation zone (Fig. 3, *F* and *G*). The basic depression-potentiation pattern is preserved, suggesting that the stochastic model is more robust to parameter variations.

### Responses of visual cortex model to spike triplets

To facilitate direct comparisons, parameters in the stochastic model were previously chosen to match parameters from the original deterministic model as closely as possible. This is not, however, the parameter set that best matches experimental data regarding the probability of vesicle release and PPD. Because there is no evidence of facilitation in visual cortex, the magnitude of PPD is equivalent to the probability of vesicle release. To match the experimental data with stochastic synaptic release, we use experimentally observed values of PPD (Markram and Tsodyks 1996; Philpot et al. 2001; Varela et al. 1997) to constrain our value of *p*_{dr0} and then set the number of docked vesicles to obtain a realistic probability of failure (Feldmeyer et al. 2002). The results in Fig. 4 are based on the assumption that there are two vesicles, each with a probability of release (*p*_{dr0}) of 0.3, yielding a probability of failure (no release on a presynaptic spike) of 0.49, (see methods). Although there is a large variability in the experimental data (Feldmeyer et al. 2002), these parameters qualitatively characterize much of the experimental data. This parameter set selected within a reasonable range of biological estimates accounts for experimental data points in all regions (Fig. 4*A*). Averaging 12 stochastic runs (each with a new random seed) produces smooth transitions between regions of potentiation and depression (Fig. 4*B*; MAE 0.3505). These results account for the experimental data better than the deterministic model.

To simulate the post-pre-post induction protocol, we had to include BPAP depression for two proximal post spikes (Shah et al. 2006). We assume that the coefficient of depression (ratio 2nd BPAP magnitude to 1st) is initially 50% and recovers exponentially to 100% with a time constant of 55 ms (see methods). Simulation of the post-pre-post protocol, using our stochastic model, produces a good fit to experimental data (Fig. 4, *C* and *D*, MAE 0.2754). Because postsynaptic BPAPs are calculated deterministically, there is less variation in results between runs with different random seeds than in the pre-post-pre case. There are only two actual spike patterns possible (post-pre-post and post-post) when considering the presence or absence of the pre spike, and the probability distribution for each triplet presentation is the same (assuming stimulation frequency is sufficiently long to allow for full recovery of presynaptic resources in-between stimuli). This explains the minimal difference between weight change calculated by a single run (Fig. 4*C*) and the average over 12 runs (Fig. 4*D*).

We have shown in Fig. 3 that varying the probability of vesicle release in the stochastic model does not significantly change the global plasticity pattern. However, for individual conditions (a *t*_{1}, *t*_{2} pair), changing the probability of vesicle release may considerably impact on resulting plasticity. We varied *p*_{dr0} over a wide range and examined the effects on synaptic plasticity of selected points in three quadrants of the pre-post-pre induction protocol. Depending on the choice of *t*_{1} and *t*_{2} (Fig. 5), altering *p*_{dr0} can produce very little change, change the magnitude of plasticity, or change the sign (from depression to potentiation as *p*_{dr0} is increased). Such dependence of plasticity on *p*_{dr0} (or other parameters) may explain why different neurons exhibit different magnitudes and even different types of plasticity (i.e., depression or potentiation) in response to similar induction stimuli and provides a testable hypothesis of our model.

### Frequency dependence of STDP

Spike timing-dependent plasticity curves are frequency dependent (Markram et al. 1997; Sjostrom et al. 2001). Under some conditions, and in some systems, a minimal frequency might be required for inducing plasticity. Further, as the frequency increases LTP is induced over an increasing range of Δ*t*; eventually, at a high enough frequency, all values of Δ*t* seem to result in LTP. We have previously shown that the calcium-dependent model is frequency dependent but that its STDP curve saturates at ∼10 Hz. This is significantly below the observed experimental saturation point (Shouval et al. 2002a). The addition of stochastic synaptic dynamics and BPAP attenuation can significantly alter the frequency dependence of STDP.

We have simulated spike-pair induction protocols at different frequencies, with the same set of parameters used above (Fig. 6*A*). We normally run our simulation at 1 Hz. At 5 Hz, the STDP curve changed very little. At 10 Hz, the depression region broadened significantly, although the peak levels of depression and potentiation did not change. At 20 Hz, less LTD is observed, and at 40 Hz, we obtain only LTP. At higher frequencies, there is an apparent periodicity in the curves with a period of *T* = 1/frequency. This is due to the symmetry between a stimulation with an interval of Δ*t* and of *T* − Δ*t;* these stimulation paradigms will only be different for the first spike pair. The parameters of our model were set to approximate plasticity induced by pairs and triplet induction protocols in synapses between layer II–III and layer IV in visual cortex slices, and we are not aware of systematic frequency-dependent plots similar to those in Fig. 6*A* for these particular synapses. Nevertheless, these plots are qualitatively similar to the frequency dependent plots of synapses between different excitatory layer V cells in visual cortex (Sjostrom et al. 2001).

It should be noted that most of the frequency dependent features in the STDP curves are present in the deterministic model: the higher cutoff frequency at which only potentiation was observed is mainly due to the introduction of short-term synaptic dynamics, a feature already present in our modified deterministic model (Shah et al. 2006). The introduction of stochastic transmitter release mainly affected the final synaptic weight at higher frequencies. The deterministic model, however, results in final synaptic weights at 40 Hz that are much higher than those of the stochastic model (>2.0 for every Δ*t*, not shown).

Because the frequency dependence simulated in Fig. 6A was not tested in the same experimental system in which the multispike protocols were carried out, we decided to examine a different frequency-dependent protocol that was tested in the same visual cortex layer 4 to 2–3 synapses (Froemke et al. 2006). We therefore carried out simulations using trains of bursts containing 5–5 pre-post spike trains (Froemke et al. 2006). Within each burst, the 5 pre- or post spikes are presented at a frequency (referred to as burst frequency) of 10, 20, and 40 Hz, and proximal pre and post spikes are separated from each other by an interval that varies over a range of −20 to 20 ms (experiments used only −6 ms). When the post spike leads the pre spike (negative pre-post spike interval), the synaptic weight increases as the burst frequency was increased, similar to what was observed in Froemke et al. (2006). However, when the pre-post spike interval is positive (pre spike leading the post spike), the model predicted relatively little change in synaptic weight as the frequency was varied (Fig. 6*B*). Clearly these results differ from those observed for pair stimulations in the studies of Markram et al. (1997) and Sojstrom et al. (2001), in which there was no LTP at low frequency. However, because we set our parameters to account for the triplet data of Froemke et al. (2002), where low-frequency LTP was observed, it is no surprise that our results are consistent with the frequency dependence shown by Froemke et al. (2006) in the same system and with the same experimental protocols.

### Simulations of hippocampal data

The plasticity pattern induced by spike pairs in hippocampal cultures (Bi and Poo 1998) is qualitatively similar to the patterns induced in visual cortex (Froemke and Dan 2002; Markram et al. 1997) and somatosensory cortex (Froemke and Dan 2002; Markram et al. 1997). However, plasticity induced in cultured hippocampal neurons by multispike protocols revealed no depression in contrast to the significant depression observed for some parameters in visual cortex. In light of these findings, it was proposed that potentiation dominates synaptic plasticity induction in the hippocampus (Wang et al. 2005). This observation raises the possibility of distinct underlying mechanisms in hippocampus and visual cortex.

In addition to differences in the effects of multispike induction protocols, visual cortex and hippocampus also exhibit different short-term synaptic dynamics. Neurons in the hippocampus, unlike many neocortical neurons, often show paired-pulse facilitation (PPF) (Chen et al. 2004; Debanne et al. 1996; Dobrunz 2002; Dobrunz and Stevens 1997), have lower probability of vesicle release (Dobrunz 2002; Dobrunz and Stevens 1997; Sudhof 2000), and much longer recovery time constants from PPD (1–3 s) (Debanne et al. 1996; Dobrunz 2002; Dobrunz and Stevens 1997; Murthy and Stevens 1999; Pyle et al. 2000) than in visual cortex. Because short-term synaptic dynamics impact significantly induction in visual cortex, we postulate that much of the difference between multispike plasticity in visual cortex and hippocampus could arise from the differences in short-term synaptic dynamics.

Our simple model of synaptic facilitation assumes that the release probability jumps from *p*_{r} = *p*_{dr0} to *p*_{r} = (1 + γ)*p*_{dr0} immediately after a presynaptic spike (Fig. 7*A*) and then recovers with a time constant τ_{F}. By comparing this model to experimental data (Debanne et al. 1996), we extracted plausible parameters for γ and τ_{F} and used them in our simulations. We found that γ = 0.8, and τ_{F} = 100 ms produce a model that is comparable to experimental findings (Debanne et al. 1996). One significant difference is that our facilitation is instantaneous whereas experimental levels require several tens of milliseconds to reach maximal values. Different experimental studies (Chen et al. 2004; Debanne et al. 1996; Pyott and Rosenmund 2002) or even different cells within the same slice can produce significantly different parameters (see methods for discussion of hippocampal parameter selection).

Synaptic facilitation, like synaptic depression, changes the probability of the different types of events and can therefore affect plasticity. We carried out simulations of synaptic plasticity induced by pre-post-pre triplets using different values of γ. Figure 7*B* shows the final synaptic weight as a function of facilitation magnitude for several choices of *t*_{1} and *t*_{2} in a pre-post-pre protocol averaged over 12 stochastic simulations. When the intervals between the two pre spikes are small, increasing the magnitude of facilitation causes plasticity to trend toward depression and, in some cases, even changes the sign from potentiation to depression. When the synapse is facilitated, the probability of post-pre events increases, thus producing more depression. It is not surprising that facilitation has a significant effect for short intervals but no effect for large intervals because of the relatively fast recovery from facilitation. The variation between different stochastic runs using different random seeds is bigger for hippocampal model (Fig. 7*B*) than that for the visual cortex model (Fig. 5) because the average number of available vesicles is smaller due to the long recovery time constant.

To observe the complete triplet induced plasticity map, we set γ = 0.8, a realistic value according to experimental data, and varied *t*_{1} and *t*_{2}. Simulation results are displayed in Fig. 8*A*. These results are in good agreement with the limited experimental data points from hippocampal neurons (Wang et al. 2005) for the pre-post-pre (MAE 0.0550). To simulate post-pre-post data, we used the same release model and assumed that the attenuation of the postsynaptic spikes is 30% (less than in visual cortex) with a recovery time constant of 35 ms. These assumptions are consistent with experimental data (Colbert et al. 1997; C. M. Colbert, personal communication). The post-pre-post simulations (Fig. 8*B*) are also in good agreement with the available experimental data (MAE 0.1132). The results for the specific data points obtained experimentally by Wang et al. (2005) are quite different from results for visual cortex using the same values of *t*_{1} and *t*_{2}. However, we observe that the global plasticity patterns are not that different from those of visual cortex for the pre-post-pre protocol (compare Fig. 8 and Fig. 4).

The small number of experimental data points from hippocampal neurons is not sufficient to fully constrain the model, and there are many choices of parameters that would match the experimental data reasonably well. Because we do not have very good estimates of the different physiological parameters, we investigated how changes to the parameters of short-term plasticity affect the resulting multispike plasticity. One such parameter is synaptic facilitation. Because facilitation was not observed in visual cortex neurons but is evident in hippocampal data, we thought it must play some role in the formation in the plasticity, especially in the pre-post-pre protocol. However, turning off facilitation did not change the pattern significantly (supplemental Fig. 1^{1} ). This is not surprising considering the relatively short time constant of recovery from facilitation (100 ms in our model) and the robustness of the stochastic model. Decreasing the time constant for vesicle recovery produced bigger (although still small) changes in responses (supplemental Fig. 1), especially when the probability of presynaptic vesicle release was increase from 0.19 to 0.3. However, their effect on the global plasticity pattern is still small. In fact, all the variations we tried yielded a good match to the limited number of data points (MAEs ranged from 0.04 to 0.2).

For the post-pre-post protocol, the general patterns in these two brain regions are also similar (compare Figs. 8*B* and Fig. 4). The main difference between the patterns of plasticity in hippocampus and visual cortex simulations is the depression region that separates two regions of potentiation in visual cortex (Fig. 4*C*), whereas in the hippocampus, there is one continuous potentiation region (Fig. 8*B*). This difference between hippocampal and visual cortex results arises from different assumptions about the depression of the second BPAP: deeper depression in BPAP results in better separation between the two potentiation regions. Changing BPAP depression in either the visual cortex or hippocampal parameter set has confirmed this observation. Obviously, in visual cortex, the plasticity data (Froemke and Dan 2002) support a separation between the two potentiation regions, but in hippocampus, the limited data points do not preclude a continuous potentiation pattern. Thus differences in simulation results of visual cortex and hippocampus may suggest that the BPAP in the visual cortex is depressed more than in hippocampus; a hypothesis that can be tested experimentally. In our simulation, the BPAP depression parameters for hippocampal neurons were estimated from experimental data (C. M. Colbert, unpublished data), but there are no data available in visual cortex regarding the magnitude of BPAP depression.

## DISCUSSION

Plasticity induced in different brain regions by similar induction protocols can lead to significantly different results. In this paper, we demonstrated how differences in the synaptic transmission statistics and the dynamics of BPAPs affect the induction of bidirectional synaptic plasticity and how this might account for the differences in the plasticity induced in visual cortex versus hippocampus. We also demonstrated that deterministic and stochastic simulations produce significantly different results and that the stochastic simulations better approximate experimental data. Interestingly, stochastic simulations were also less sensitive to precise parameter values. Our results also demonstrated that, although at specific inter-spike intervals, plasticity induced in visual cortex and hippocampus seems significantly different, the overall qualitative patterns of plasticity induced by triplets are similar. It is, however, problematic drawing conclusions about differences in underlying mechanisms on the basis of a small number of experimental data points.

### Stochastic synaptic transmission

One of the primary results of this study is that predictions from the stochastic version of the plasticity model match the experimental data more accurately than the deterministic version of the model. The stochastic version could match data points from visual cortex slices in all quadrants (Fig. 3*C*), whereas the deterministic version could do so only with strong depression and very fast recovery of presynaptic transmitter release (Shah et al. 2006), an assumption that is not supported by available physiological data. The better results obtained from the stochastic version are encouraging because the stochastic simulations are a much better approximation of what really occurs during synaptic transmission.

Most models of neuronal properties and of synaptic plasticity are deterministic. Deterministic models are usually justified when one models whole-neuron properties because of the large number of channels and other relevant molecules. However, signals that control synaptic plasticity in a single synapse are highly variable due to the small number of vesicles available for release, the small number of postsynaptic receptors, and the small volume of the spine. The small volume of a synaptic spine implies that some of the key molecules have a very small copy number. Previously, we have shown the effect of fluctuations due to the small number of postsynaptic NMDA receptors (Shouval and Kalantzis 2005). Here we demonstrated the distinct effects of stochasticity arising from presynaptic variables and showed that taking it into account improves our model's predictions. We chose here a number of docked vesicles *n* = 2, which is consistent with experimental results. An increase in the number of docked vesicles would result in less variability. We expect that different synapses in the brain will have a different number of docked vesicles, different synaptic dynamics, and different depression rates. We did not include here the stochastic opening of postsynaptic NMDA receptors. Including stochastic postsynaptic synaptic transmission would alter the triplet plasticity map (Fig. 3.), but significant changes would occur primarily outside the region for which there exist many triplet data points. Such additional complexity would not impact the central claim of this paper about how the stochasticity of presynaptic vesicle release affects the resulting triplet plasticity curves.

Our stochastic model was also more robust than the deterministic version. The deterministic and stochastic versions of the model responded differently to changes in parameters, for example, the probability of transmitter release. The deterministic version was sensitive to changes in release probability, resulting in drastic changes of the boundary of depression-potentiation regions (Fig. 3, *D* and *E*, and Fig. 4*A* of Shah et al. 2006), whereas in the stochastic version, the basic plasticity pattern is preserved in response to similar changes in release probability (Fig. 3, *C, F*, and *G*). Because estimates of release probability vary widely across different experimental studies, and vary with different experimental parameters such as temperature, the reduced sensitivity of the stochastic model might have a functional significance. The relative stability of the plasticity pattern predicted by the stochastic model suggests that neurons may have the ability to ensure their functional stability under diverse conditions.

### Can results in visual cortex and hippocampus be accounted for by the same model?

Plasticity induced experimentally by multispike induction protocols seems to produce significantly different results in visual cortex and in hippocampus. In visual cortex, plasticity induced by spike triplets seems to be dominated by the first spike pair (Froemke and Dan 2002). In hippocampus, neither the first nor second spike pair dominates, but plasticity induced by triplets exhibits no LTD, suggesting that LTP dominates (Wang et al. 2005). These results might indicate that the underlying plasticity mechanisms are significantly different in these different systems. However, in this paper, we showed that these different experimental results can all be consistent with a common plasticity mechanism. We demonstrated that the different experimental results might arise from different assumptions about synaptic and neuronal dynamics. In excitatory connections between layers 4 and 2–3 neurons in visual cortex, PPD seems to dominate (Varela et al. 1997), whereas in hippocampus there is also significant paired-pulse facilitation (Debanne et al. 1996). In addition, the recovery time constant is much longer in the hippocampus (order of seconds) than in the visual cortex (hundreds of milliseconds). These differences in synaptic dynamics between visual cortex and hippocampus can account for the different results in pre-post-pre induction protocols, as shown when comparing Figs. 4 and 8. Moreover, in Fig. 7, we showed specifically how increasing the magnitude of facilitation can convert LTP into LTD.

To account for synaptic plasticity induced by post-pre-post triplets in hippocampus and visual cortex, we postulated attenuation of BPAP magnitude in the second post spike. Although there is limited information about BPAP depression in hippocampus (Colbert et al. 1997), we could not locate similar data for visual cortex slices. In our simulations, we used 30% BPAP depression for the hippocampal model based on recording from hippocampal neurons (Colbert et al. 1997 and unpublished data) but had to use stronger BPAP depression (50%) to match the visual cortex plasticity data (Fig. 4). Although 50% attenuation might seem extreme, it is slightly less than the 60% attenuation assumed by Froemeke et al. (2006) in their revised suppression model. Because these parameters are not based on direct measurements, but were set to account for the plasticity data, our different assumptions about BPAP depression could be viewed as a prediction of this model. Thus we predict that in visual cortex the second BPAP in a pair is more significantly depressed than in hippocampus.

### Other models

Froemke and co-workers (Froemke and Dan 2002; Froemke et al. 2006) presented phenomenological spike-suppression models to account for plasticity observed in multispike induction protocols. This model assumes addition of effects from individual spike pairs but assumes that the contribution of subsequent spikes to plasticity are suppressed immediately after the first spike and then recover with different time constants for pre- versus postspikes. These different spike suppression mechanisms for pre- and postspikes could be identified with PPD and BPAP attenuation. However, the recovery parameters in their work were not based on measured parameters of PPD and BPAP attenuation but were chosen to optimize agreement between the model and the experimental results. Two additional models that can account for STDP in visual cortex (Abarbanel et al. 2003; Senn 2002) have also been used to simulate plasticity elicited by spike triplets in visual cortex. More recently, a phenomenological learning rule, which is based on pair and triplet interactions only, has been developed and can account well for plasticity induced by various induction protocols at different frequencies (Pfister and Gerstner 2006). These models were able to accomplish this without additional assumption regarding short-term synaptic and dendritic dynamics and therefore have no direct equivalents of PPD and BPAP attenuation. Another calcium-dependent model that was developed to account for spike timing-dependent plasticity in hippocampus incorporates several molecular components, which compete to induce either LTP or LTD (Rubin et al. 2005). This model incorporates a veto mechanism to eliminate the pre-post form of LTD present in most calcium-dependent models, and its assumptions are sufficient to account for plasticity induced by multispike protocols in hippocampal culture (Wang et al. 2005).

Our stochastic model is able to qualitatively reproduce results of multi-spike induction protocols in many systems and under various conditions. However, even if we can account for experimental results in two different types of neurons using the same fundamental plasticity model, this does not imply that these two different systems must employ a similar underlying mechanism. Indeed, experiments that directly address the underlying mechanisms of synaptic plasticity have shown that there are significant differences between induction of LTD in the hippocampus and in cortex (Bender et al. 2006; Sjostrom et al. 2003) and that there are even different underlying mechanisms in different layers within the same cortical region (Crozier et al. 2005). Moreover, in hippocampus, two forms of LTD have been established: an NMDAR-dependent and an mGluR-dependent form (Oliet et al. 1997). Nevertheless, this does not invalidate the central methodological points of this paper, as follows: synaptic and neuronal dynamics have a significant impact on the outcome of different induction protocols; stochastic nature of synaptic transmission has a significant impact on plasticity, independent of the precise mechanism, and therefore the statistics of the different types of events must be taken into account; stochastic models are less sensitive to changes in system parameters, such as probability of release, and are therefore more likely to produce robust predictions; and we cannot claim simply by observing different plasticity curves in different types of neurons that the underlying plasticity mechanisms are fundamentally different because such differences can arise from differences in other physiological parameters such as synaptic and neuronal dynamics.

## GRANTS

This work was funded by National Science Foundation Colabortive Reserve in Computational Neuroscience Grant 0515285 and by National Institute of Neurological Disorders and Stroke Grant 2 P01-NS-038310.

We thank Y. Dan and R. Froemke for generously sharing experimental data. We also thank B. Blais for many useful discussions.

## Footnotes

↵1 The online version of this article contains supplemental data.

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 © 2007 by the American Physiological Society