Memory Maintenance in Synapses with Calcium-Based Plasticity in the Presence of Background Activity

Most models of learning and memory assume that memories are maintained in neuronal circuits by persistent synaptic modifications induced by specific patterns of pre- and postsynaptic activity. For this scenario to be viable, synaptic modifications must survive the ubiquitous ongoing activity present in neural circuits in vivo. In this paper, we investigate the time scales of memory maintenance in a calcium-based synaptic plasticity model that has been shown recently to be able to fit different experimental data-sets from hippocampal and neocortical preparations. We find that in the presence of background activity on the order of 1 Hz parameters that fit pyramidal layer 5 neocortical data lead to a very fast decay of synaptic efficacy, with time scales of minutes. We then identify two ways in which this memory time scale can be extended: (i) the extracellular calcium concentration in the experiments used to fit the model are larger than estimated concentrations in vivo. Lowering extracellular calcium concentration to in vivo levels leads to an increase in memory time scales of several orders of magnitude; (ii) adding a bistability mechanism so that each synapse has two stable states at sufficiently low background activity leads to a further boost in memory time scale, since memory decay is no longer described by an exponential decay from an initial state, but by an escape from a potential well. We argue that both features are expected to be present in synapses in vivo. These results are obtained first in a single synapse connecting two independent Poisson neurons, and then in simulations of a large network of excitatory and inhibitory integrate-and-fire neurons. Our results emphasise the need for studying plasticity at physiological extracellular calcium concentration, and highlight the role of synaptic bi- or multistability in the stability of learned synaptic structures.


Introduction
Many experiments have shown that long-lasting changes in synaptic efficacy can be induced by spiking activity of pre-and postsynaptic neurons [1,2]. In hippocampal and neocortical in-vitro preparations, both long-term potentiation and depression can be induced by protocols in which pre-and postsynaptic neurons emit tens to hundreds of spikes in specific temporal patterns [3][4][5][6][7][8][9][10]. In those preparations, plasticity has been shown to depend both on relative timing of pre-and postsynaptic spikes ('spike timing dependent plasticity', or STDP), and the firing rates of pre-and postsynaptic neurons. These induced changes in the connectivity of a neural circuit have then been assumed to maintain or initiate a longterm memory trace of external inputs that triggered these synaptic changes [11]. However, in order for this theory to be valid, the induced synaptic changes need to survive both activity triggered by other inputs, and the ongoing background activity that is pervasive in cortex [12,13]. How changes in synaptic connectivity survive the continuous presentation of other inputs has been the subject of several studies [14,15]. Here, we study the decay of the synaptic memory trace due to background activity using a theoretical approach. Synaptic plasticity has been described using a multitude of different models and approaches [6,7,[16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31]. In early plasticity models, synaptic changes were purely induced by the firing-rates of pre-and postsynaptic neurons [16,17,19]. At the end of the nineties, theorists introduced purely spike-timing based models [21,23]. Finally, more recent models have been striving to capture a wide range of experimental data, and as a result capture both the spike-timing and firing rate dependence of synaptic plasticity [6,[24][25][26][27][28][29][30][31]. These models are natural candidates for studies of the stability of synaptic changes during ongoing activity. In this paper, we choose the model of Graupner and Brunel [28] for the following reasons: (i) the model includes the calcium concentration in the post-synaptic spine, which is known to be a crucial link between pre-and postsynaptic activity and long-term synaptic changes; (ii) the model exhibits bistability of synaptic changes accounting for experimental evidence suggesting that CA3-CA1 synapses in the hippocampus are bistable [32,33]; (iii) the model is simple enough that the dynamics of the synaptic efficacy variable can be computed analytically.
Postsynaptic calcium entry has been identified to be a necessary [34][35][36] and sufficient [37][38][39] signal for the induction of synaptic plasticity (but see ref. [40]). However, most of the in vitro experiments evoking synaptic changes use elevated extracellular calcium concentrations, while in vivo calcium levels are estimated to be around 1.5 mM [41]. The impact of reduced calcium entry due to the lower extracellular calcium concentration in vivo on the time scale of synaptic decay has not been considered heretofore.
In the present paper, we study the persistence of synaptic changes, first in a synapse connecting a pair of independent Poisson neurons, and second in a large network of excitatory and inhibitory leaky integrate-and-fire (LIF) neurons. We show that in the absence of bistability, synaptic changes decay exponentially during ongoing activity and that the time constant exhibits a power-law like behaviour with respect to the present firing rate. We demonstrate that the reduced extracellular calcium concentration in vivo leads to several orders of magnitude longer memory time scales. The introduction of bistability in the synaptic plasticity rule further stabilises synaptic changes at low firing rates and extensively prolongs memory time scales when combined with the in vivo extracellular calcium conditions. Finally, we extend our results to a large recurrent network of LIF neurons, where we demonstrate network firing rate stability under synaptic plasticity, decay of an implanted memory for in vitro parameters and long term memory maintenance for in vivo parameters.

Results
Memories are thought to be stored in the brain thanks to activity-dependent modifications of synaptic connectivity. According to this hypothesis, memories stored by a particular neural circuit are encoded by the state of all the modifiable synapses of the circuit. Synaptic plasticity allows particular patterns of activity to leave a trace in the connectivity matrix, but this trace is then potentially vulnerable to the ongoing activity that follows. An important question is therefore what controls the time scale of the persistence of a particular synaptic state, in the presence of such ongoing activity. To study this question, we initialize the efficacy of a synapse (either an isolated one, or part of a network) at a particular value, and study how this efficacy decays with time in the presence of background activity, using a calcium-based model of synaptic plasticity [28].
In the model, the temporal evolution of the synaptic efficacy variable, r(t)[½0,1, is described by zs ffiffi ffi t p ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi where t is the time constant of synaptic efficacy changes, and c(t) is the calcium concentration. The dynamics of r depends on four terms: 1. The dynamics are governed by a potential U(r) for low calcium concentrations (c(t)vh d ) since all other terms on the right-hand side of Eq. (1) are zero then. In the following we consider two possible scenarios for U(r): (i) a flat potential, U(r)~0 -in this case the synaptic efficacy variable stays constant in time in the absence of calcium transients. This means all possible values of r[½0,1 are stable; (ii) a double well potential, In this case, r evolves towards one of two possible stable fixed points (the minima of U), one at r~0 -the DOWN state -, the other at r~1 -the UP state -, depending on the initial condition.
This corresponds to a bistable synapse. 2. For intermediate calcium concentrations (c(t)wh d ), the synapse is depressed, with a depression rate c d . 3. For large calcium concentrations (c(t)wh p wh d ), the synapse undergoes both potentiation, with a potentiation rate c p , and depression, with the same rate as in the h p wc(t)wh d . Since c p wc d , potentiation dominates over depression in that region. 4. A noise is only active when calcium concentration is above the lowest plasticity threshold h D , and increases in magnitude when the upper plasticity threshold, h P , is also crossed. s defines the amplitude of the noise, and g(t) is a Gaussian white noise process with unit variance. Changes in r are induced by increases in the postsynaptic calcium concentration, c(t) (see Eq. (11), in Methods), evoked by pre-and postsynaptic spikes. The calcium concentration increases by an amount C pre , in response to presynaptic spikes, while it increases by an amount C post in response to postsynaptic spikes. It decays exponentially with a time constant t Ca in between spikes. Calcium transients induced by presynaptic activity are assumed to represent calcium influx through NMDA receptors, while calcium transients induced by postsynaptic spikes are assumed to represent activation of voltage-gated calcium channels [42] (see [28] for more details of the model).
This calcium-based model of synaptic plasticity has been used to successfully fit data from various experimental preparations [28]. Here, we use the data-set that best fits plasticity data obtained in visual cortex slices [6] -hereafter called the 'in vitro' parameter set. In this experiment, the extracellular calcium concentration was set to be 2.5 mM [6], which is significantly higher than the estimated in vivo concentration of about 1.5 mM [41]. Here we assume that a decrease in extracellular calcium concentration leads to a proportional decrease in the calcium influx into the postsynaptic spine. Using this assumption, we can readily predict the effects of decreasing the extracellular calcium concentration on the plasticity rule in the calcium-based model by scaling the amplitudes of the pre-and postsynaptically evoked calcium

Author Summary
Synaptic plasticity is widely believed to be the main mechanism underlying learning and memory. In recent years, several mathematical plasticity rules have been shown to fit satisfactorily a wide range of experimental data in hippocampal and neocortical in vitro preparations.
In particular, a model in which plasticity is driven by the postsynaptic calcium concentration was shown to reproduce successfully how synaptic changes depend on spike timing, specific spike patterns, and firing rate. The advantage of calcium-based rules is the possibility of predicting how changes in extracellular concentrations will affect plasticity. This is particularly significant in the view that in vitro studies are typically done at higher concentrations than the ones measured in vivo. Using such a rule, with parameters fitting in vitro data, we explore how long the memory of a particular synaptic change can be maintained in the presence of background neuronal activity, ubiquitously observed in cortex. We find that the memory time scales increase by several orders of magnitude when calcium concentrations are lowered from typical in vitro experiments to in vivo. Furthermore, we find that synaptic bistability further extends the memory time scale, and estimate that synaptic changes in vivo could be stable on the scale of weeks to months.
transients according to the ratio of calcium concentrations, i.e. 1.5/2.5 = 0.6. This leads to what we call the 'in vivo' parameter set. Values of all parameters for both conditions are indicated in Table 1.
The dynamics of the synaptic efficacy in response to calcium transients under the in vitro and the in vivo conditions are illustrated in Fig. 1. The synaptic efficacy is only modified when the calcium concentration increases above the depression threshold h D (Fig. 1,B-E). For the in vitro case, this happens whenever a postsynaptic spike occurs since C post wh D , but for the in vivo parameter set this happens much more rarely because of the smaller calcium amplitudes (C pre ,C post vh D in the in vivo case; see Table 1). In the latter case, synaptic changes are only induced whenever subsequent spikes occur in short succession such that calcium accumulates and crosses the depression and/or potentiation threshold. Such events are rare at low firing rates (Fig. 1,D-E).

Memory decay for a synapse connecting two independent Poisson neurons
We now proceed to study the time scales of synaptic decay. We start with the case of a synapse connecting two neurons firing according to uncorrelated Poisson processes, and compare the memory time constants in the flat and double-well potential cases. Simulations were performed using an event-based implementation of the synaptic plasticity model, which updates the synaptic efficacy only upon the occurrence of pre-and postsynaptic spikes (see Methods for details).
We initialise the synaptic efficacy to r~1 and investigate the time constant of decay in the presence of an ongoing constant firing rate, initially for the flat potential synapse (Eq. 1, with U(r)~0). Pre-and postsynaptic neurons emit uncorrelated spikes following Poisson statistics, both with a mean rate of 1/s. Under these conditions, a fully potentiated synapse progressively decays and eventually fluctuates around a value of 0.2. On average, this decay is well described by a single exponential function (Fig. 2,A,B). The time constant of this decay is much longer in the case of the in vivo parameter set (Fig. 2,B) than in the in vitro parameter set (Fig. 2,A). The decay time constant is 2.5 minutes for the in vitro case and approximately 2 hours for in vivo in the presence of 1/s pre-and postsynaptic firing (Fig. 2,C).
The dynamics of the synaptic efficacy (Eq. 1) can be described by a truncated Ornstein-Uhlenbeck (OU) process if single calcium transients induce small changes in the synaptic efficacy and if the potential is flat (see [28] for the non-truncated case). Truncation of the process is induced by the bounds at r~0 and r~1. In such a process, the mean synaptic efficacy decays exponentially with a time constant, t eff , which is given by to an asymptotic average efficacy, r r(n) (see Eq. (22) in Methods), where C P (n)~c P a P (n) and C D (n)~c D a D (n) are the net potentiation and depression rates which depend on the rates C p and C d as well as on the average fractions of time spent above the potentiation and depression thresholds, a P (n) and a D (n), respectively. The average fractions of time the calcium traces spend above the potentiation and depression thresholds are given by a P (n)~1{ where P(c) is the probability density function of the calcium variable, that can be computed analytically in the case of independent pre-and postsynaptic Poisson firing [28,43] (see Methods for details). The theory provides an excellent match for the dynamics of the mean synaptic efficacy -compare in Fig. 2,A,B the truncated OU theory (blue and red curves), with the simulation mean (green and light blue curves). Synaptic efficacy decay becomes faster with increasing pre-and postsynaptic firing rates since the calcium trace spends more time above depression and potentiation thresholds (Fig. 2,C). At the same time, the asymptotic value of synaptic efficacy ( r r) increases due to an increase in time spent above the potentiation threshold ( Fig. 2,D). As a result of the smaller in vivo calcium amplitudes, the efficacy decay for the in vivo case is, at all firing rates, much slower than the decay in vitro (Fig. 2,C). The asymptotic efficacy value is lower, at small firing rates (nv1/s), for the in vitro case since isolated postsynaptic spikes always cross the depression threshold (C post wh D ) which results in a large net depression rate C D , compared to in vivo (Fig. 2,D).
To get a deeper understanding of the dependence of the memory time scale on the firing rates of pre-and postsynaptic neurons, we set C pre~Cpost~1 . This simplification allows us to derive a power law relationship between the memory time scale and the firing rate t eff *1=(nt Ca ) k , where k is the number of (pre and/or post-synaptic) spikes required to clear the depression/ potentiation thresholds. To compute the memory time scale, we need to compute the fraction of times spent above the depression and potentiation thresholds, a D and a P . In the case h D v1vh P , one can show that at low rates a P %a D . Consequently it is only necessary to focus our analysis on a D . When h D v1, the time spent above the depression threshold is The in vitro values are obtained in [28] by fitting the model, using a gradient descent method, to cortical plasticity data presented in Figure 8a (frequency dependence of synaptic plasticity for a fixed relative pre-post-spike timing) of [6].
where n is the firing rate of pre-and postsynaptic neurons (n pre~npost~n ), t Ca is the decay constant for the calcium concentration and A~exp ({2nt Ca c)=C(2nt Ca ) (see Eqs. (13) - (17)). This closed form solution allows us to perform an expansion for low firing rates n a D (n)*2nt Ca log Similarly for 1vh D v2vh P we have in the low rate limit, where Li 2 is the dilogarithm, Li 2 (z)~X ? k~1 z k =k 2 . Thus, in both cases we find that the memory time scale depends on the firing rate as where k~Ceil(h D ). We expect this relationship to hold in general. Intuitively, this is due to the fact that we need k spikes arriving simultaneously on a time scale of order t Ca in order for the calcium concentration to cross the depression threshold h D , and that the probability of observing k spikes in a time interval t D is at low rates proportional to (nt Ca ) k . We also expect the result to hold in general for C pre =C post . In this case, we expect that k~Ceil h D max (C pre ,C post ) . The derived power law behaviour for t eff is plotted in Fig. 2,C together with the full analytical solution for t eff . We see that as expected, t eff scales as 1=n for the in vitro parameter set, where a single spike is enough to cross h D , while it scales as 1=n 2 for the in vivo parameter set, where two spikes are needed to cross the depression threshold.
The implication of this theoretical result is that, at low firing rates, there is a direct relationship between the number of spikes required to clear the lower plasticity threshold and the memory time scale. Note that the full synaptic efficacy model with C pre =C post is considered in the following (see Table 1)

Memory decay for a bistable synapse
We now turn to examine the effect of a bistability on memory time scales. The dynamics of the synapse is now described by Eq.
(1), where the potential U(r) is given by Eq. (2). This double well potential leads to a bistable synapse, that can take two possible efficacy states (r~0 and r~1) in the absence of activity. In the presence of background activity, transitions between these two states become possible. We investigate stability and transition  Table 1). Whenever the calcium trace crosses the depression (cyan) or potentiation thresholds (orange), changes in the synaptic efficacy (green) are induced. (D,E) Same as in B,C but with calcium amplitudes corresponding to the in vivo case (C pre~0 :337 and C post~0 :744). The small calcium transients do not cross the depression/potentiation thresholds and no efficacy changes are observed. Note that the flat potential for r is used here and that noise is set to zero for clarity, s~0. doi:10.1371/journal.pcbi.1003834.g001 times for in vitro and in vivo parameter sets as a function of preand postsynaptic firing rates.
The effect of background activity on the dynamics of r can be explained by the fact that it modifies the potential, U(r), leading to an effective potential (see Methods). In Eq. (10), the first term on the r.h.s. corresponds to the 'bare' double well potential U(r) (Eq. (2)); the second term describes the effect of depression on the potential, that tends to strengthen the stability of the lower well (DOWN state) at r*0, at the expense of the other well that tends to disappear when a D (n) increases; finally, the last term describes the effect of potentiation, that shifts the minimum of the only remaining well towards higher values of r when a P (n) increases. Thus, there are two distinct regions of firing rates in the bistable case with respect to the effective potential. For sufficiently low rates, the effective potential still has two minima (see Fig. 3,A, and the effective potentials for 0.1/s and 1/s, indicated by orange and magenta curves in the inset). There is a critical value of the rates at which the high efficacy minimum disappears through a saddlenode bifurcation. Beyond this rate, the synapse is no longer bistable, and synaptic efficacy has one stable state only (Fig. 3,A), equivalent to the asymptotic efficacy value for the flat potential ( Fig. 2,D). Finally, at high firing rates, the 'bare' potential becomes negligible, and the effective potential approaches a quadratic potential with a single stable state whose location depends on the rate (green curve in the inset in Fig. 3,A). The transition from double-well to single well regimes occurs at different firing rates for the in vitro (*0:04/s) and the in vivo (*1:3/s) parameter sets due to the larger calcium amplitudes in the former.
For the in vitro parameter set, adding bistability to the synaptic efficacy has no influence on the decay time constant for firing rates larger than approximately 0.1/s (Fig. 3,B). In contrast, for the in vivo parameter set, bistability considerably prolongs memory decay times with respect to synapses with flat potential at firing rates below ,1.4/s. In the presence of two stable states, the decay of memory occurs only due to synaptic noise fluctuations that push the synaptic efficacy out of the upper well. The influence of the double well potential on the dynamics of the synaptic efficacy traps synapses in the UP state leading to long dwell times before crossing the potential barrier and converging to the low efficacy state (Fig. 3,C). The double-well has a prolongation effect on memory duration up to firing rates of about 3{4/s due to the transition between double-well and single-well regimes. At high firing rates, the potentiation and depression processes dominate and the effects of the double-well becomes negligible for both parameter sets, that is, the decay time constant is indistinguishable between flat and double-well potential synapses (see Fig. 3,B). For low firing rates, we can accurately predict the increase in the decay time constant in the presence of bistability using Kramers escape rate for the mean first passage time across a potential barrier (Fig. 3,B; see Methods Eq. (26)). In this regime, we calculated an effective decay time constant using Kramer's escape theory given by t eff * expf2DU=s 2 eff g, where DU is the height of the effective potential barrier and the noise term, s 2 eff , drives the escape of the efficacy from the upper stable state (see magenta line in Fig. 3,B for the in vivo case). Both terms DU and s 2 eff are dependent on n and are detailed, along with t eff , in Eqs. (23) and (26) (see Methods for more details). In the low rate limit, s 2 eff !1=n k and therefore the memory time scale increases exponentially with the inverse of the rate to a power k, t eff ! exp (a=n k ), where k is again the number of simultaneous spikes needed to cross the depression threshold. This exponential dependence extends the time scale for synaptic decay at 1/s to the order of one month for a bistable synapse with the in vivo parameter set, up from hours for a synapse with flat potential.

Steady-state behaviour of networks of LIF neurons with plastic synapses
We next study the behaviour of the calcium-based synaptic plasticity model in a recurrent network of spiking neurons. We first examine the steady-state of synaptic efficacy and network activity. We again make use of the event-based implementation of the synaptic plasticity rule (described in Methods) allowing us to simulate much longer time scales than are normally attainable by a time stepping simulator.
The recurrent network consists of 8000 excitatory and 2000 inhibitory leaky integrate-and-fire (LIF) neurons. Each neuron receives an external input which consists of a constant (DC) term and a white noise term. External noise is independent from neuron to neuron. Each neuron also receives synaptic inputs from other neurons in the network. The connection probability between any two neurons is 0.05 and uniform in space and across neuron types. Synapses between excitatory neurons are plastic according to the calcium-based plasticity model (Eq. 1), while all synapses involving inhibitory neurons are fixed. Parameters of the network are chosen so that the network settles in a stable asynchronous irregular state [44]. Hence, correlations between neurons are weak. See Methods for more details of the network model.
The fixed point of the network can be determined analytically by solving a set of three self-consistent equations for the excitatory and inhibitory mean rates as well as for the mean excitatory-toexcitatory (E?E) synaptic efficacy (see Methods). Two of these equations give the stationary firing rates of excitatory and inhibitory populations (Eqs. (32) -(33)), as a function of the mean E?E synaptic efficacy [44,45]. The third equation gives the mean E?E synaptic efficacy as a function of the firing rate of the excitatory population, assuming Poisson firing statistics of neurons (Eq. (22)). Starting from the analytically determined initial conditions, the recurrent network converges to a steady-state of constant average firing rates of all neurons in the network, and constant average synaptic efficacy of the plastic connections. Figure 4,A shows how the firing rates observed in the simulations compare with the analytically predicted firing rates. It shows that at sufficiently low rates, the analytical prediction gives a very good estimate of the observed rates; however, for rates above 3Hz the observed rates are significantly lower than the analytical prediction. Likewise, the analytical prediction for the mean E?E synaptic efficacy significantly overestimates the observed efficacies (green dots in 4,B). What is the source of the difference between theory and simulations in predicting the steady network state? When synapses are fixed in the network at the efficacies predicted by the corresponding firing rate, the analytically predicted network firing rates provide a good approximation of the observed activity (blue dots in 4,A). This suggests that the underestimation of firing rates and synaptic efficacy emerges from the mapping of firing rates onto synaptic efficacy, and not due to correlations between spike trains of different neurons. We examine this effect in Fig.4,C, where we show asymptotic mean synaptic efficacy results under three different conditions. First, we simulated a population of disconnected, independent LIF neurons, receiving stochastic independent inputs with the same mean and variance as in the 'real' network simulation (magenta dots). By definition, this simulation led to uncorrelated spike trains. In this simulation, fake synaptic connections between neurons obeyed the plasticity rule, but had no effect on the dynamics of the neurons. Second, we simulated a connected recurrent network with constant synaptic weights. As in the first case, we simulated 'fake' synaptic connections that obeyed the standard plasticity rule, but these fake synapses had no effect on the dynamics (blue dots). Third, we simulated a standard recurrent network in which synaptic weights are plastic according to the plasticity rule (green dots). All three simulations show indistinguishable results, and in all three cases the average (real or fake) synaptic efficacies are consistently smaller compared to the analytical shot noise prediction (4,C, blue line). This suggests that correlations have a negligible effect on mean efficacies and firing rates, and that the differences between simulations and theory are due to differences in spiking statistics between the LIF model and a Poisson process.
To investigate further how the spiking statistics of the LIF model and in particular the interspike-interval (ISI) distribution causes the differences seen in Fig. 4, we varied the ISI distribution of the LIF neuron by changing the reset potential (V reset , see Methods). This change had a strong effect on the average synaptic efficacy (4,D). A reset potential close to threshold (V reset~{ 55 mV, V threshold~{ 50 mV) yields an overrepresentation of short ISI compared to Poisson firing (4,D, inset) and in turn overestimates the average synaptic efficacy (4,D; cyan dots). Conversely, more depolarised reset potentials lead to an under-representation of short ISIs with regard to Poisson firing and consequently to an underestimation of the average synaptic efficacy (4,D; magenta, red and green dots). We use an intermediate value of V reset~{ 60 mV in the following network investigations.
To conclude this section, the calcium-based synaptic plasticity rule does not affect the stability of the asynchronous irregular state in a large recurrent network of LIF neurons. Since LIF neurons in the network exhibit ISI distributions which deviate from those of Poisson neurons, the accuracy of our calculation of the average synaptic efficacy which is based on Poisson firing decreases with increasing firing rates up to a certain point. At high firing rate, calcium remains above the plasticity thresholds most of the time and the fraction of time spent above the thresholds converges to one, irrespective of the underlying neuron model.

Memory decay in a recurrent network of LIF neurons
Finally, we examine the decay of a memory trace in a network for the in vitro and the in vivo parameter set. We initialise all excitatory-to-excitatory synaptic weights at their theoretically predicted asymptotic weights, except for a randomly selected subset of 5% which are set to a weight of 1. With the in vitro parameter set, the potentiated synapses decay relatively quickly to their asymptotic value (Fig. 5,B). The time course of the average decay can be described by a single exponential function and the decay time constant is well approximated by the time constant, t eff , of synaptic decay from the truncated OU process (see Eq. (3); Fig. 5,C). This means that the average dynamics of synaptic decay in the network is equivalent to synapses driven by independent pre-and postsynaptic Poisson neurons firing at the same rate as the excitatory neurons in the network (compare to Fig. 3,B). The addition of the double-well potential does not change the decay time constant for the in vitro parameter set, as for a single synapse driven by independent pre-and postsynaptic Poisson firing (Fig. 5,C orange stars; compare with Fig. 3,B). The lack of short ISIs in LIFs compared to independent Poisson neurons leads to a small increase in observed decay times in the network as compared with the OU theory (see Fig. 5,C).
In contrast, when using the in vivo parameter set with the double-well potential, we observe that the potentiated synapses get locked in the UP state for the duration of the network simulation with an excitatory neuron firing rate of 1/s (Fig. 6,C). None of the synapses in the potentiated subset crosses the unstable fixed point and converges to the DOWN state during a network simulation of 120 min, neither does the reverse transition occur. We expect that the escape from the UP state will be predicted by Kramers escape rate (Eq. (26)) which correctly accounted for escape dynamics of an isolated synapes driven by independent pre-and postsynaptic Poisson processes (Fig. 3B). There, the decay time constant for a firing rate of 1/s is on the order of a month, a time scale that cannot be reached by our network simulation.
Hence, as in case of independent Poisson neurons, the combination of a double-well potential with the in vivo parameter set leads to several orders of magnitude longer memory time constants, compared to the in vitro parameter set and a flat potential.

Discussion
In this paper, we studied the stability of synaptic efficacy, in a plastic synapse subjected to background activity of pre-and postsynaptic neurons. We used a calcium-based plasticity model that has been shown to fit experimental data in hippocampal and neocortical preparations [28]. The model was investigated numerically, using an event-based implementation of the plasticity rule, as well as analytically, using a diffusion approximation. Thanks to this formalism, we derived scaling laws that describe how memory time scale is related to the firing rates of pre-and postsynaptic neurons. At low firing rates, we find that, when synapses are monostable, synaptic efficacies decay to an equilibrium value with a time scale that depends on the firing rates as a power law, t eff *1=n k , where k is the number of simultaneous spikes needed to cross the depression threshold. When synapses are bistable, memory decay is akin to diffusion of a particle out of a potential well, which leads to much stabler memories, with time scales that increase exponentially with the inverse of the firing rates, t eff * exp (a=n), at low rates. We showed that these estimates accurately reproduce the results of simulations, both of a synapse connecting two isolated independent Poisson neurons, and of a large network of LIF neurons.
We have focused here on how changes in the amplitudes of the calcium transients affect memory time scales. A change in other model parameters also affects these time scales. Changing the depression threshold, for example, has a similar pronounced effect, since the exponent in the scaling law between memory time scale and background rate depends on the ratio between this threshold and the amplitudes of the calcium transients (see (9)). On the other hand, changing other parameters of the model (such as the time constants and the potentiation and depression rates) have much milder effects, in the flat potential regime, since the time scale depends algebraically on such parameters, rather than exponentially. In the absence of further data, we have assumed a linear relationship between the external calcium concentration and the calcium influx to explore the in vivo regime. A non-linear relationship between the extracellular calcium concentration and the concentration in postsynaptic microdomains -conjectured to be relevant for synaptic plasticity -should modify our results quantitatively rather than qualitatively.
Previous studies have investigated memory maintenance in networks of neurons connected by synapses endowed with standard spike-timing dependent plasticity rules [46]. Billings and van Rossum (2009) demonstrated that the memory time scale depends dramatically on whether the rule is additive or multiplicative. In a multiplicative STDP rule, in which synaptic change depends on the current value of the weight such that synaptic changes decrease when the weights approach the bounds, distributions of weights are unimodal [46][47][48] and the memory of synaptic changes decay as 1=n 2 , since synaptic changes occur upon coincidence of pre-and postsynaptic spikes in the characteristic time window of the STDP rule. These behaviours are very similar to the behaviour of the calcium-based rule in the flat potential case, in the parameter region in which two spikes are needed to cross the depression threshold. This is due to the fact that the calcium-based rule defined by Eq. (1) is multiplicative. In the calcium-based rule however, the exponent describing the memory decay at low rates can be set to an arbitrary integer number, through an appropriate rescaling of the ratio between the amplitude of the calcium transients and the depression threshold. In additive STDP rules, the picture changes dramatically and the synaptic weight distributions become bimodal, with weights attracted either to the lower or upper bounds through a symmetry breaking mechanism [23,46]. In this situation, the memory time scales are much longer, and decay of synapses is similar to diffusion in a double well potential.
Several studies have shown that synaptic bi-or multi-stability can emerge from a number of mechanisms such as positive feedback loops in extensive protein signaling cascades [49], autophosphorylation of CaMKII [50][51][52][53][54], self-sustained regulation of translation [55], or modulation of receptor trafficking rates [56]. Such mechanims of bistablity are effectively implemented here in the form of the double well potential. Miller et al. (2005) studied the stability of the up state in a model of the bistable calcium/ calmodulin-dependent protein kinase II system with respect to stochastic fluctuations induced by protein turnover [57]. They show that the CaMKII switch composed of a realistic number of CaMKII proteins is stable for years even in the presence of protein turnover, phosphatase as well as free calcium fluctuations. The transitions induced by background activity investigated here impose an upper limit on memory life-time which is typically lower, indicating that in vivo neuronal activity, not protein turnover, will be the limiting factor of memory life-times.
In vivo, memory in synaptic connectivity structures will be affected both by ongoing background activity, but also by changes in network activity induced by external stimuli. How ongoing presentations of external inputs affect memories of past stimuli has been the subject of several studies in recent years (e.g. [15,20,58]), in simpler networks of binary neurons. A detailed exploration of this issue in the model studied here is outside the scope of the present paper, but we anticipate that parameter regions that extend the robustness of synaptic memories in the face of background activity will also tend to protect the network against changes induced by external inputs.
Distributions of synaptic weights have been examined in a number of studies [4,6,[59][60][61][62]. In all of these studies, distributions of synaptic weights appear unimodal and skewed, and peak at a low weight. In some cases, the distribution has been shown to be well fitted by a lognormal distribution [60,62]. This seems at first sight at odds with the distributions of weights shown in the present paper, which are either a truncated Gaussian in the flat potential case, lacking the fatter tail of the lognormal distribution, or bimodal in the double-well case. However, the model in the flat potential case can be made consistent with the data, by choosing synaptic efficacy variables which are an exponential of the r variable, rather than being linearly related to r. In this case, synaptic efficacies themselves become exponentiated Ornstein-Uhlenbeck processes, consistent with [62]. The model with a double-well potential could also be made consistent with a unimodal distribution, provided the synaptic up and down states are highly heterogeneous from synapse to synapse. Finally, we should point out that the distributions we observe are asymptotic distributions under a statistically constant distribution of inputs.
Synapses in vivo are typically subjected to highly non-stationary firing rates of pre and post synaptic neurons. These nonstationarities can also potentially strongly affect distributions of synaptic weights in our model.
A large number of distinct learning rules that capture quantitatively both spike-timing and firing rate effects have been proposed recently [25][26][27][28][29][30][31]. Our rule can be distinguished from most of those rules by the fact that it includes calcium concentration as its primary dynamic variable, which allows us to extrapolate parameters of the rule from in vitro to in vivo conditions, as we have explained here. Scaling laws derived here can be expected to hold also in those other models: at low rates, the time scales of memory decay are expected to be inversely proportional to the rates to a power equal to the number of spikes needed to provoke plasticity. This power should be equal to 2 for standard STDP rules, triplet rules [25], and calcium-based rules in which 2 spikes are needed to cross the depression threshold [24,27]; 1 for spike and voltage based rules [26].
In this work, we have made the hypothesis that synaptic weights are altered during background activity, and that one can treat background activity as being essentially uncorrelated with the synaptic connectivity structure. Memory time scales could in principle be further extended by two factors. A first mechanism would be to gate plasticity by specific neuromodulator(s) that are present only during stimulus presentation. This idea is consistent with a growing body of experimental data showing how plasticity is modulated by dopamine [63], acetylcholine [64,65], noradrenaline [66] (see also [67] and references therein). However, we note that the model we have used here is built from in vitro plasticity data where these neuromodulators were present at very low concentrations, if at all. Hence, we believe that these neuromodulators are likely to enhance plasticity during behaviourally relevant epochs, but that the memory time scales discussed here are likely not to be affected if neuromodulators are not present at high levels during background activity.
A second mechanism that would extend memory time scales would be a scenario in which background activity is in fact strongly correlated with the connectivity structure, and wanders stochastically between network states that are strongly correlated with the states of the network during stimuli presentation. This idea is consistent with a growing experimental literature [68][69][70] showing how spontaneous activity is transiently strongly correlated with sensory responses in visual and auditory cortices, and it is also consistent with the ubiquitous supra-Poissonian variability, potentially due to the doubly-stochastic process of combined rate stochasticity and individual neuronal Poisson spike processes, seen in background activity in cortex [13,71]. Recurrence of activity states resembling the network activity during stimulus presentation could refresh existing memory traces and therefore prolong their lifetimes.
We showed here that the low extracellular calcium concentrations in vivo could have a strong impact on plasticity. A first prediction of calcium-based rules is that plasticity seen in standard protocols should be greatly reduced (and even possibly vanish altogether) at physiological calcium concentrations. While to our knowledge no study has explicitly compared plasticity results at different extracellular calcium concentration, comparisons between different studies using different extracellular concentrations seem to be consistent with this prediction. In hippocampal slices, a standard low frequency STDP protocol produces LTD for all time differences with 2 mM extracellular calcium [9], while it produces the standard STDP curve with 3 mM calcium [10]. A second prediction is that induced synaptic changes should be much more stable in the face of ongoing pre-and postsynaptic activity. These results emphasise the need for experimental studies at physiological calcium concentrations *1:5 mM [72], unlike most published studies that used concentrations in the range 2{3 mM. Our predictions could be easily tested in slice experiments, by providing background activity at a specified rate after the plasticity-inducing protocol. Similar experiments have been performed in the developing Xenopus retino-tectal system in vivo [73], where activity-induced modifications were shown to be erased by subsequent 10 minutes of spontaneous activity. Our model would predict that in cortical slices, at 2.5 mM calcium, induced synaptic changes should disappear on a time scale of minutes, while at 1.5 mM calcium, they should be stable on a time scale of *1 hour.
We provided here an event-based update scheme of plastic synapses which greatly accelerates simulations and should strongly facilitate future studies of the dynamics of recurrent networks with plastic calcium-based synapses. On the theoretical front, it would be interesting to extend the theory to non-Poissonian renewal processes [74] such as for leaky integrate-and-fire neurons used here, which would give a better approximation of average synaptic efficacies, especially at higher firing rates. It would also be of great interest to examine how synaptic connectivity is modulated by non-stationary external inputs, and how such changes in connectivity affect in turn the intrinsic dynamics of the network.
Our investigations show that realistic external calcium concentration and multi-stability of synapses might stabilise memory traces against the potentially deleterious effect of ongoing background activity. These results call for studies of synaptic plasticity induction and maintenance in more realistic conditions and ideally in the intact animal. They provide a glimpse of how plasticity results obtained in vitro might translate to the living organism.

Calcium-based plasticity model
The temporal dynamics of the synaptic efficacy in the calciumbased model are given in Eq. (1) (for details see [28]).
Changes in r are driven by the postsynaptic calcium concentration, c. The calcium dynamics are modelled using instantaneous increases of size C pre and C post in response to pre-and postsynaptic spikes, respectively, followed by an exponential decay dc dt~{ where t Ca is the calcium decay time constant, and C pre , C post the pre-and postsynaptically evoked calcium amplitudes. The sums go over all pre-and postsynaptic spikes occurring at times t i and t j , respectively. The time delay, D, between the presynaptic spike and the occurrence of the corresponding calcium transient accounts for the slow rise time of the NMDAR-mediated calcium influx. We use two parameter sets in this paper. The in vitro parameter set is obtained by fitting the calcium-based plasticity model to plasticity data obtained in cortical slices ( [6]; see [28] for details of the fitting procedure). These experiments were performed with 2.5 mM extracellular calcium concentration. The in vivo calcium amplitudes are obtained by scaling C pre and C post according to the extracellular calcium concentration in vivo, estimated to be 1.5 mM [41] (see Results).

Probability density function of the calcium concentration
We shortly describe here how the probability density function (PDF) of the calcium concentration can be calculated if pre-and postsynaptic neurons fire as independent Poisson processes at rate n~n pre~npost (see [28,43] for more details). In these conditions, the calcium concentration is a shot noise process, whose probability density function is given by the master equation [43], cP'(c)~(2nt Ca {1)P(c){nt Ca P(c{C pre ){nt Ca P(c{C post ):ð12Þ The probability density function P(C) allows us to calculate the fraction of time spent above the depression and potentiation thresholds according to a D (n)~1{ Ð hD 0 P(c)dc and a P (n)~1{ Ð hP 0 P(c)dc.
In the simple case C pre~Cpost~1 , the solution to Eq. (12) is given by where 2 F 1 (a,b,c,z)~X ? n~0 (a) n (b) n (c) n z n n! is the ordinary hypergeometric function [75], c is Euler-Mascheroni constant, c*0:577215665, and C is the gamma function.
Fitting the calcium-based model to cortical plasticity data yields C pre vC post (see Table 1). In this case, the solution of Eq. (12) where B~A=(C pre C post ) nt Ca is a normalisation parameter which assures that Ð P(I)dI~1. The PDF for cw min (2C pre ,C post ) is obtained from a numerical integration of Eq. (12).
Diffusion approximation for the synaptic efficacy with a flat potential Performing a diffusion approximation of the synaptic efficacy r turns Eq. (1) into an Ornstein-Uhlenbeck process (see [28] for details). The temporal evolution of r is then described by for the case of a flat potential (i.e. LU=Lr~0). C P (n)~c P a P and C D (n)~C D a D are the mean potentiation and depression rates, respectively. The bounds at r~0 and r~1 lead to a truncated Ornstein-Uhlenbeck process, whose distribution is a truncated Gaussian, whose mean converges exponentially to r r(n)~C where s 2 r~s 2 (a P (n)za D (n)) is the Gaussian with zero mean and unit variance, and is the complementary cumulative density function of G. The time constant, t eff , of the exponential decay to r r is defined in Eq. (3).

Kramers expected escape time from a double-well potential
In the case of a double-well potential, the diffusion approximation turns Eq. (1) into a Fokker-Planck equation This equation can be rewritten as where the effective potential, U eff , is the sum of the 'bare' potential U and two quadratic terms proportional to the potentiation and depression rates, respectively (see Eq. (10)), and s 2 eff is the amplitude of the effective noise s 2 eff (n)~s 2 (a D (n)za P (n)): When the effective potential is dominated by the double-well term (first term on the rhs of Eq. (10)), the 'escape' rate from the UP state is driven by noise and can be estimated using Kramers theory [76,77]. The height of the potential barrier, DU eff~Ueff (r un ){U eff (r up ), determines the mean dwell time in the UP state, where r up and r un are the local minima and maxima of the effective potential and are obtained solving LU eff Lr~0 . This allows us to calculate the expected escape time from the potential well Numerical methods: Event-based implementation The temporal evolution of individual synaptic weights in the calcium-based model can be calculated in an event-based manner (as opposed to a finite difference method with a fixed time step Dt) in an analytically exact way. This greatly accelerates network simulations since the network variables are updated at the occurrence of spikes only. In the following we describe the practical implementation of the event-based update.
For the event-based update, three variables have to be stored per synapse: the calcium concentration, c, the synaptic efficacy variable, r, and the time of the last spike seen by the synapse, t. The synapse variables c and r must be updated on the occurrence of three events: (1) at the presynaptic spike when the postsynaptic voltage is increased, (2) with delay D after a presynaptic spike when the presynaptically evoked calcium increase occurs (see Eq. (11)), (3) and at the postsynaptic spike when the postsynaptic calcium increase occurs.
Calcium update. The calcium concentration decays exponentially between events and is instantaneously increased by the amplitude C post when a postsynaptic spike occurs. In the case of a presynaptic spike, the calcium increase C pre occurs after the delay D (Eq. (11)). In consequence, we update the calcium concentration as a decay followed by a calcium concentration increase where the amplitude depends on the identity of the spike and the delay D (at time t m for pre-synaptic spikes and t n for post-synaptic spikes). The calcium update for time t iz1 (after the last update at t i ) at the three update events described above reads (1) at the presynaptic spike C pre d(t iz1 {D{t m ) (2) after delay of the presynaptic spike C post d(t iz1 {t n ) (3) at the postsynaptic spike: Synaptic efficacy update. For the propagation of the synaptic efficacy variable, we distinguish between two different regimes, that is, stochastic and deterministic. When the calcium concentration at time t i is lower than both thresholds, c i vh D ,h P , the dynamics of r are described deterministically. When the calcium concentration crosses either or both thresholds, the update is stochastic, for the duration of the suprathreshold period, and the dynamics of the mean and the standard deviation of the synaptic efficacy are described by the corresponding Ornstein-Uhlenbeck process. The different regimes may be updated sequentially in a piecewise manner, accounting for first suprathreshold and then subthreshold behaviour.
Here, we determine the possible potentiation/depression threshold crossings of the calcium trace between events i and iz1 with the inter-event interval Dt~t iz1 {t i . t P is the time the calcium trace spends above the potentiation threshold within that interval, t D is the time the calcium trace spends between the potentiation and the depression threshold, and t 0 the time the calcium trace spends below the depression threshold given by t 0~D t{t P {t D . Here, we consider h P wh D only. The updates for the case h P vh D are equivalent. c end refers to the calcium concentration right before the update at event iz1, that is c end~ci exp ({Dt=t Ca ).
Between events at t i and t iz1 , the following six possible threshold crossings can occur (see Fig. 7 The efficacy, r, is updated in a piece-wise fashion. Stochastic updates are performed when the calcium trace spent time above the potentiation threshold (t P w0, cases: I, II, III), or between the potentiation threshold and the depression threshold (t D w0, cases: II, III, IV, V; and for h P wh D ). A deterministic update is performed if the calcium trace spent time below the depression threshold (t 0 w0, cases: III, V, VI).
In case one or both thresholds are crossed in the interval (t i ,t iz1 ], the temporal evolution of r cannot be solved analytically in the double well potential case, because of the combined presence of the stochastic term and the quartic potential. However, for the parameters used in this paper the double-well potential can be neglected in the suprathreshold regions because c D ,c P & max DLU(r)=LrD (see Table 1). We can therefore approximate the temporal evolution of r by an Ornstein-Uhlenbeck (OU) process, for which the potential of r during stimulation is quadratic with the minimum at r r. We confirmed the validity of this approximation by comparing the event-based simulation based on this approximation, with a simulation which was event-based only in sub-threshold epochs, and based on a forward Euler method in the supra-threshold epochs. The results of both simulations were indistinguishable, which confirms that we can ignore the potential well during the suprathreshold period.
Using this approximations, the updates of the synaptic variable after suprathreshold epochs are given by: where z P and z D are Gaussian random variables of unit variance and zero mean. Note that t D~0 in case I, and therefore only the first update is performed. Equivalently, t P~0 in cases IV as well as V and only the second update rule is performed. When the calcium concentration, c i , is smaller than the potentiation and the depression threshold (t 0~D t{t P {t D w0, cases III, V and VI), the first two terms on the rhs and the noise term of Eq. (1) are zero. That reduces the rhs of Eq. (1) to {dU=dr which can be integrated analytically for the flat-and the double well potential considered here. The update of r is therefore deterministic and depends on the choice of the potential: 1. flat potential 2. double-well potential r(t iz1 )~1 2 x 0~( r(t i zt P zt D ){1=2) 2 =(r(t i zt P zt D )(r(t i zt P zt D ){1)).

The network
We implemented and simulated a recurrent network of 10,000 leaky integrate-and-fire (LIF) neurons, 8,000 of which are excitatory (E) neurons and 2,000 of which are inhibitory (I). Any two neurons have a spatially uniform probability of connection of 0.05. Autapses are specifically disallowed. Synapses between E neurons are plastic and their weight dynamics are described by the calcium-based plasticity model (Eq. 1, [28]). All other synapses have fixed strength w ab (a,b[fE,Ig). A presynaptic spike induces a voltage jump of size w ab in the postsynaptic neuron.
The membrane potential of neuron i of population a evolve according to where is a common external drive to all neurons, comprising a constant input, m aX , and a white noise of amplitude s~5mV. g ai (t) is a Gaussian white noise process with unit variance and zero mean, which is uncorrelated from neuron to neuron. In the absence of synaptic inputs each membrane potential decays exponentially to the leak potential, V leak~{ 70 mV, with a time constant t m~2 0 ms. Spiking occurs when the voltage crosses a threshold, V thr~{ 50 mV, after which it is reset to the reset potential, V reset~{ 60 mV. During all of our simulations, we set the refractory period, during which the membrane potential is fixed at V reset after spiking, to zero. The three sums in the r.h.s. of Eq. (30) go over the two populations b[{E, I}, all presynaptic neurons j, and presynaptic spikes of neuron j in population b, that occur at times t bjk . Each presynaptic spike of neuron j in population b causes a jump of amplitude w aibj in the voltage of neuron i after a delay t L . Here, the delay is chosen to be equal to the integration time step dt~0:01 ms (see below). For all connections involving inhibition (i.e. all (a,b)=(E,E)), the connectivity matrix is set as w aibj~cij w ab where c ij are independent, identically distributed (i.i.d.) Bernoulli variables, c ij~1 with probability 0.05, 0 with probability 0.95, and the fixed synaptic weights are w IE~0 :1mV, w II~{ 0:4mV and w EI~{ 0:4mV. E-E synapses are given by w aibj~cij r ij w EE where c ij are again i.i.d. Bernoulli variables, c ij~1 with probability 0.05, 0 with probability 0.95, r ij obeys Eq. (1), and w EE~0 :2mV. The average value of r is initially, and remains throughout our simulations, much smaller than 0.5, which means that with a 4 : 1 ratio in the E to I populations, for vrwƒ0:5 recurrent inhibition dominates excitation, leading to a stable asynchronous irregular state (see Fig. 8) [44].

Numerical methods: Network simulations
We numerically simulated the recurrent network of LIF neurons using the forward Euler method with a time step of dt~0.01 ms. Synapses were updated using the event-based implementation described above. The simulations were implemented in C and OpenCL and run on general-purpose GPUs. Parallel generation of random numbers, for the Gaussian noise in the LIF equations, was implemented using the Random123 library [78].
In order to initialise the simulations close to their steady-state, with the in-vivo parameter set and the double-well potential, we first calculate the probability distribution function (PDF) for the synaptic weights assuming a 1/s pre-and post-synaptic Poisson firing process. We then use a reverse lookup of the associated cumulative distribution function (CDF) to determine the random initial values for the synaptic efficacies.
Computing analytically mean firing rates and E-E synaptic efficacy In a network of excitatory and inhibitory LIF neurons receiving white noise inputs, the mean firing rates of excitatory and inhibitory neurons are given by [44,45] n E~W (m E ,s E ) ð32Þ n I~W (m I ,s I ) ð33Þ where W is the standard LIF static transfer function [44,45,79,80 and s a is the amplitude of the fluctuations in the inputs to population a[fE,Ig, s E~sext z(C EE r r(n E ) 2 w 2 EE t m,E n E )z(C EI w 2 EI t m,E n I ) ð37Þ s I~sext z(C IE w 2 IE t m,I n E )z(C II w 2 II t m,I n I ) ð38Þ Note that in Eqs. (35,37) r r is given by Eq. (22). Note also that for the parameters studied in this paper the effect of heterogene-ities in numbers of inputs [81,82] have a negligible effect on the mean firing rates of the network.