Abstract

In stochastic models of synaptic plasticity based on a random walk, the control of fluctuations is imperative. We have argued that synapses could act as low-pass filters, filtering plasticity induction steps before expressing a step change in synaptic strength. Earlier work showed, in simulation, that such a synaptic filter tames fluctuations very well, leading to patterns of synaptic connectivity that are stable for long periods of time. Here, we approach this problem analytically. We explicitly calculate the lifetime of meta-stable states of synaptic connectivity using a Fokker-Planck formalism in order to understand the dependence of this lifetime on both the plasticity step size and the filtering mechanism. We find that our analytical results agree very well with simulation results, despite having to make two approximations. Our analysis reveals, however, a deeper significance to the filtering mechanism and the plasticity step size. We show that a filter scales the step size into a smaller, effective step size. This scaling suggests that the step size may itself play the role of a temperature parameter, so that a filter cools the dynamics, thereby reducing the influence of fluctuations. Using the master equation, we explicitly demonstrate a bifurcation at a critical step size, confirming this interpretation. At this critical point, spontaneous symmetry breaking occurs in the class of stochastic models of synaptic plasticity that we consider.

1.  Introduction

Accumulating evidence suggests that synaptic plasticity may involve discrete, all-or-none changes at single synapses (Yasuda, Sabatini, & Svoboda, 2003; Bagal, Kao, Tang, & Thompson, 2005; Sobczyk & Svoboda, 2007). Other evidence indicates that single synapses may, furthermore, occupy only a very limited number of discrete states of synaptic strength (Petersen, Malenka, Nicoll, & Hopfield, 1998; Montgomery & Madison, 2002, 2004; O'Connor, Wittenberg, & Wang, 2005a, 2005b; but see Enoki, Hu, Hamilton, & Fine, 2009), although it is possible that experimental protocols may saturate synapses at upper or lower strength limits (Elliott, 2010a). Discrete state synapses, whether there are many or few states, are theoretically appealing and have received increasing attention (Willshaw, Duneman, & Longuet-Higgins, 1969; Abarbanel, Talathi, Gibb, & Rabinovich, 2005; Senn & Fusi, 2005; Baldassi, Braunstein, Brunel, & Zecchina, 2007; Romani, Amit, & Amit, 2008; but see Satel, Trappenberg, & Fine, 2009, for a critique). Despite this possibility of discrete, step-like changes in synaptic strength, biophysical models of spike-timing-dependent plasticity (STDP) (Caporale & Dan, 2008) mostly assume graded changes in synaptic strength at single synapses (Castellani, Quinlan, Cooper, & Shouval, 2001; Castellani, Bazzani, & Cooper, 2009; Senn, Markram, & Tsodyks, 2001; Senn, 2002; Shouval, Bear, & Cooper, 2002; Rubin, Gerkin, Bi, & Chow, 2005; but see Graupner & Brunel, 2007). We proposed instead a model of STDP based on fixed-amplitude changes in synaptic strength at single synapses in which the probability, but not amplitude, of change depends on spike timing (Appleby & Elliott, 2005). Such a model can be shown to constitute a non-Markovian random walk in synaptic strength (Elliott, 2010b).

Controlling fluctuations in synaptic strength in such models represents a formidable challenge. The corrosive influence of fluctuations in binary-state synapse models is also a serious problem. Proposed solutions to this problem include slow, stochastic updating of synaptic strength (Senn & Fusi, 2005; Fusi & Senn, 2006) and the introduction of cascades of internal states in synapses that essentially act as probabilistic shift registers, storing the history of plasticity induction events (Fusi, Drew, & Abbott, 2005; Leibold & Kempter, 2008). We have proposed instead that single synapses could act as low-pass filters, filtering plasticity induction processes before expressing synaptic plasticity (Elliott, 2008; Elliott & Lagogiannis, 2009). Because we are particularly interested in the emergence of segregated states of neuronal connectivity through activity-dependent competition during neuronal development (Purves & Lichtman, 1985; for reviews of models, see, e.g., Swindale, 1996; van Ooyen, 2001), (Purves & Lichtman, 1985) (Swindale, 1996) (van Ooyen, 2001) we have used our model to investigate, in simulation, the stability against fluctuations of these segregated states over long periods of time (Elliott & Lagogiannis, 2009), and we have computed the mean time to express synaptic plasticity in a class of three explicit models of the filtering process (Elliott, 2011). In this letter, our principal aim is to address stability against fluctuations not in simulation but rather analytically. We strive to maintain as much generality as possible by considering generic plasticity induction and plasticity expression processes, restricting to our tristate switch model of STDP only when we require an explicit, concrete model of plasticity induction. We use a Fokker-Planck equation formalism to perform these calculations.

Our results for the mean transition time between different segregated states of synaptic connectivity, driven by fluctuations, reveal, however, some extremely suggestive properties, which we also explore. We find that the presence of a synaptic filter scales the plasticity step size, generating a smaller, effective plasticity step size. This observation hints at the possibility that the plasticity step size may play the role of a temperature or temperature-like parameter in stochastic, random-walk models of synaptic plasticity. To show this, we move from a Fokker-Planck equation formalism to a master equation formalism and explicitly demonstrate the existence of a bifurcation at a critical value of the plasticity step size. The presence of a bifurcation indicates spontaneous symmetry breaking in these models at the critical step size.

The outline of the remainder of this letter is as follows. In section 2, we summarize our previous results for filtering models of plasticity expression. In section 3, we formulate our approach to synaptic dynamics in terms of the Fokker-Planck and master equations, making clear any approximations and simplifications that are necessary or convenient. We consider the Fokker-Planck formalism in section 4, first examining the fixed-point structure induced by the deterministic Liouville dynamics, and then explicitly calculating the lifetimes of meta-stable states of segregated patterns of synaptic connectivity in the presence of fluctuations induced by the diffusion term in the Fokker-Planck equation. In order to address the possibility of a bifurcation in the dynamics, we consider the master equation approach in section 5, demonstrating the presence of a bifurcation both numerically and analytically, and calculating explicitly the critical plasticity step size. Finally, in section 6, we discuss our results, examine the relationship between our work and that of others, and consider our approach to the control of fluctuations more generally.

2.  Summary of Previous Results

We consider a model of synaptic plasticity in which changes in synaptic strength at single synapses occur in fixed-amplitude, discrete steps of sizes T+ and T for potentiation and depression, respectively. We distinguish between the induction of synaptic plasticity and the expression of synaptic plasticity, and suppose that expression is governed by a process that filters induction steps in order to suppress or control fluctuations in synaptic strength (Elliott, 2008, 2011; Elliott & Lagogiannis, 2009). Our major assumption, for the purposes of analytical tractability, is that potentiating and depressing induction steps at single synapses are renewal processes, which limits the history dependence of induction processes to the last induction step and transforms a typically intractable non-Markovian process into a tractable, semi-Markovian process with, in general, nonexponentially distributed waiting times between successive induction steps. Let the rates of potentiating and depressing induction steps be r±in(t; λπ, λp), where t denotes the time elapsed since the last induction step occurred and λπ and λp are the pre- and postsynaptic firing rates, respectively, at the synapse of interest. We will consider below a model for r±in based on our stochastic model of STDP (Appleby & Elliott, 2005), but where possible, we leave the forms of r±in unspecified for the purposes of generality. We define the derivatives of these rates as ρ±in(t; λπ, λp) = dr±in(t; λπ, λp)/dt. For notational ease, we will usually not denote explicitly the dependence of r±in(t) and ρ±in(t) on λπ and λp unless it is necessary to do so for clarity. Since r±in(t) = ∫t0dt′ρ±in(t′), we see that , where the notation denotes the Laplace transform of the function f(t) with transformed variable p. Hence, are the asymptotic rates of plasticity induction when the limit p → 0 exists. The probability density functions (PDFs) for induction steps, excluding the possibility that an induction step of the opposite type has already occurred, are given by Elliott (2010b):
2.1
These distributions g±in(t) are in general defective, meaning that their integrals are not unity, ∫0dtg±in(t) < 1, so that an induction event of a given type (either potentiating or depressing) is not inevitable. We refer to such distributions as excluding PDFs (ePDFs). Their sum g+in(t) + gin(t) is, however, the PDF for a single induction step of any type, and we assume this not to be a defective distribution, so that induction is inevitable. Given g±in(t), we can determine ρ±in(t) through the inverse relation:
2.2
We may expand as a Taylor series in p, writing
2.3
Because the Laplace transforms are, up to a sign convention, just the moment generating functions of the two plasticity induction step processes, are the probabilities of these induction steps, with by assumption, and are the conditional mean induction step times, with being the mean time for any induction step, regardless of its type, to occur. From equation 2.2, we see that
2.4
so that if we calculate the two induction step probabilities and the (unconditional) mean induction step time, then we can immediately deduce the asymptotic induction rates.
For the simplest version of our stochastic model of STDP, using independent Poisson pre- and postsynaptic spiking processes, we found that (Elliott, 2010b)
2.5
2.6
where λ± are the rates of the two exponentially distributed stochastic decay processes postulated in our tristate switch model of STDP to return the plasticity meta-states POT (with rate λ+) and DEP (with rate λ) to the OFF meta-state in the absence of spiking (Appleby & Elliott, 2005). Since the waiting times between independent Poisson spikes are exponentially distributed, the memorylessness of the exponential process guarantees the renewal of spiking processes and, therefore, of the induction processes in our model.

We suppose that these plasticity induction steps drive changes in the state of a synaptic filter, denoted by X(t), with a potentiating induction step incrementing the filter by ξ+ and a depressing induction step decrementing it by ξ (Elliott & Lagogiannis, 2009). For simplicity, we assume that the magnitudes of change in X(t) in both cases are equal, so we may scale these steps to be of unit size without loss of generality, so that ξ± = 1. Starting at X(t0) = 0, some t0, X(t) for t>t0 therefore performs a random walk driven by induction steps, and we take potentiation expression to occur when X(t) first exceeds an upper threshold +Θ+, and depression expression when X(t) first falls below a lower threshold −Θ, with Θ±>0. For convenience, we set these two thresholds equal: Θ± = Θ. Following plasticity expression, X(t) is reset to 0, and we ignore the possibility of a subsequent refractory period in which X(t) does not respond to induction steps. The filter X(t) may possess additional dynamics involving an X(t)-state-dependent decay to 0 at rate η|X|, decaying either deterministically and continuously in a manner similar to an Ornstein-Uhlenbeck process, or jumping stochastically in steps of unit size so that X(t) is always an integer. Before, we called the filter model without decay the linear summation model; the version with deterministic, continuous decay the low-pass filter model; and the variant with stochastic, discrete decay the molecular filter model (Elliott, 2011). Figure 1 summarizes schematically our models for the expression of synaptic plasticity based on the integration of synaptic plasticity induction signals. For these three filter models, we calculated the PDFs for potentiating and depressing expression steps, from which we can obtain the rates of plasticity expression steps, denoted by r±ex(t; λπ, λp), and their derivatives, ρ±ex(t; λπ, λp) = dr±ex(t; λπ, λp)/dt. Of course, an expressed plasticity step will change the postsynaptic firing rate λp if postsynaptic firing is synaptically coupled to presynaptic firing, as is usually the case. Typical experimental protocols, however, often hold neurons in voltage or current clamp, uncoupling spike generation from synaptic transmission, so that λp could remain unchanged even following the expression of synaptic plasticity. We take care below to distinguish between these two scenarios.

Figure 1:

A schematic illustration of the integration and filtering of synaptic plasticity induction signals leading to the expression of synaptic plasticity. The synaptic filter variable X(t), denoted here along the x-axis, undergoes increments or decrements of ξ± in response to potentiating or depressing induction signals, respectively. The filter starts at some initial time t0 at value X(t0) = 0. The PDF for a potentiating induction signal excluding the possibility that a depressing induction signal has already occurred (an ePDF) is g+in(t), where t is the time elapsed since the last induction step (of any type); the ePDF for the induction of depression is gin(t). In two of the models of plasticity expression that we consider, in addition to these induction-driven steps ξ± in X(t), there is also a decay process, represented by the dashed lines, of X(t) to zero, the overall rate of decay governed by a rate parameter η and proportional to the magnitude of the current value of X(t). This decay can be either a continuous, deterministic process or a discrete, stochastic process. When X(t) reaches the upper threshold +Θ+, a single potentiating step of size T+ is expressed. Conversely, when X(t) reaches a lower threshold −Θ, a single depressing step of size T is instead expressed. In either case, the filter X(t) is reset to zero. The ePDFs for single expression steps, g±ex(t), involving the transition of X(t) from its initial value 0 to either of its two final values ±Θ± in time t, may be calculated in terms of the induction ePDFs g±in and the possible decay process.

Figure 1:

A schematic illustration of the integration and filtering of synaptic plasticity induction signals leading to the expression of synaptic plasticity. The synaptic filter variable X(t), denoted here along the x-axis, undergoes increments or decrements of ξ± in response to potentiating or depressing induction signals, respectively. The filter starts at some initial time t0 at value X(t0) = 0. The PDF for a potentiating induction signal excluding the possibility that a depressing induction signal has already occurred (an ePDF) is g+in(t), where t is the time elapsed since the last induction step (of any type); the ePDF for the induction of depression is gin(t). In two of the models of plasticity expression that we consider, in addition to these induction-driven steps ξ± in X(t), there is also a decay process, represented by the dashed lines, of X(t) to zero, the overall rate of decay governed by a rate parameter η and proportional to the magnitude of the current value of X(t). This decay can be either a continuous, deterministic process or a discrete, stochastic process. When X(t) reaches the upper threshold +Θ+, a single potentiating step of size T+ is expressed. Conversely, when X(t) reaches a lower threshold −Θ, a single depressing step of size T is instead expressed. In either case, the filter X(t) is reset to zero. The ePDFs for single expression steps, g±ex(t), involving the transition of X(t) from its initial value 0 to either of its two final values ±Θ± in time t, may be calculated in terms of the induction ePDFs g±in and the possible decay process.

When the induction processes correspond to our tristate switch model of STDP, equations 2.5 and 2.6, it is important to notice that the expression process itself also exhibits the same dependence on spike timing. Consider, for example, a scenario in which a presynaptic-then-postsynaptic spike pair is repeatedly presented, as in typical experimental protocols, with either shorter or longer interspike intervals. In the case of the pairs with shorter interspike intervals, the probability of plasticity induction is higher, thus driving the synaptic filter more readily to the threshold for plasticity expression. However, when the interspike intervals are longer, the probability of induction is lower, making expression less likely. Of course, in the presence of an expression filter, multiple stimulus pairs are required before plasticity is expressed, consistent with experimental results (Petersen et al., 1998; Bagal et al., 2005; O'Connor et al., 2005a, 2005b; Wittenberg & Wang, 2006), so the overall magnitude of plasticity is reduced compared to that in which an expression filter is absent. Despite this, a model of synaptic plasticity expression based on an STDP induction process still exhibits expressed plasticity depending on spike time differences. Although conventional STDP models typically fail to distinguish between plasticity induction and plasticity expression, and thus we could be tempted to regard a filterless model of STDP as in some sense “standard,” it is critical to recognize that experimental data do not, in fact, mandate such an easy position. Nevertheless, we will indicate results obtained in the absence of an expression filter, corresponding to the simple parameter selection, Θ = 1, although we must be extremely cautious in interpreting these results as in any sense “standard” STDP. What is “standard” is defined by the experimental system, not by models based on simplifying (and often entirely unexamined) assumptions.

Given the induction ePDFs g±in(t), we calculated the ePDFs for single expression steps in the three filter models described above and deduced the probabilities of expression steps and the mean time for the expression of a plasticity step of any type, assuming that λπ and λp are held constant during the integration period (Elliott, 2011). Denoting these ePDFs by g±ex(t), in the linear summation model (without decay), we found that
2.7
2.8
where
2.9
The probabilities of expression steps, denoted by , and the mean time for the expression of a plasticity step of any type, denoted by , were found to be
2.10
2.11
for general induction processes defined by ρ±in(t). These lowest-order moments of the expression process are therefore sensitive only to the asymptotic rates of plasticity induction steps. The asymptotic rates of the expression of plasticity steps are given by an argument identical to that leading to equation 2.4, so that . So,
2.12
Notice that we may write
2.13
2.14
for a single function ε(u, v), that is, the expression process treats the two induction processes symmetrically. Such a symmetry must be present even with arbitrary Θ± ≠ Θ provided that we make the simultaneous replacements and Θ+ ↔ Θ. In general, the expression ePDFs g±ex(t) must also exhibit this symmetry, so that g+ex(t) = Γ(ρ+in(t), ρin(t); Θ+, Θ) and gex(t) = Γ(ρin(t), ρ+in(t); Θ, Θ+) for a single function Γ(x, y; α, β), for any induction processes ρ±in(t) and thresholds Θ±.
Results for the expression ePDFs in the low-pass filter model are considerably more complicated than those in the linear summation model and may be found elsewhere (Elliott, 2011). Such results can be derived only in the diffusion limit. Here we give the results only for expression probabilities and the mean expression time in this limit. In this filter model, in order to regularize the mean expression time as a function of λπ and λp, we previously set the filter decay rate η to be (Elliott, 2011),
2.15
where η′ is a dimensionless constant. The decay rate η then scales with the summed asymptotic induction rates, these latter rates themselves depending on pre- and postsynaptic firing rates (Elliott, 2011). We give results for this choice because in this case, they take on a particularly compelling and natural form. Defining
2.16
and
2.17
we see that with the choice in equation 2.15, z± and z0 depend on only the relative drift between the asymptotic induction rates:
2.18
The expression probabilities are then
2.19
and the mean expression time is
2.20
where erfi(z) is the imaginary error function, and
2.21
where 2F2 is a generalized hypergeometric function. The asymptotic expression rates may therefore be written as
2.22
2.23
Thus, the scale of the asymptotic expression rates is set by the overall asymptotic induction rate, , but otherwise the expression rates depend on only the relative drift Δρ between the asymptotic induction rates via z± and z0. These results hold for general induction processes ρ±in(t), and these lowest-order moments of the expression process once again depend on only the asymptotic induction rates. The asymptotic expression rates can, as with the linear summation model, be written in terms of one function, and , where an additional argument is required to represent the dependence on the dimensionless decay parameter η′.

We remark briefly that these results for the low-pass filter model do not reduce to those for a filterless model when we set Θ = 1. This is because the above results are derived by employing the diffusion limit. To reduce the low-pass filter model, in the diffusion limit, to a filterless model, we must take the limit of the low pass filter model to the linear summation model by taking η′ → 0, and set Θ = 1. When we compare results below in the presence of an expression filter to those in the absence of an expression filter by writing Θ = 1, for the low-pass filter model, it should be implicitly understood that the η′ → 0 limit is also taken. Of course, were it possible to solve the low-pass filter model exactly, this additional η′ → 0 limit would not be necessary.

Below, it will be necessary to consider the behaviors of in the low-pass filter model in the limit of large η′ in order to understand some of our results. We find that
2.24
and
2.25
By writing for , we have somewhat abused standard asymptotic notation, but we mean that these forms are asymptotically suppressed even compared to the other forms. We see that these other forms are also suppressed in the large η′ limit, as we should expect. Although plasticity expression is inevitable in this filter model (), increasing η′ makes reaching threshold harder and harder, and therefore plasticity expression takes longer and longer. For large enough η′, the suppression switches to an exponential dependence, which means that both expression rates collapse to very close to zero, for any sign of the relative induction drift Δρ. This transition to exponential suppression as η′ increases occurs over a quite narrow range of η′, and we will see its impact on our results below.
For the molecular filter model, the expression ePDFs g±ex(t) can in general be determined exactly, without the need to take the diffusion limit, but only up to a matrix inversion. If are two (2Θ − 1)-dimensional vectors of expression ePDFs with components g±Y(t) corresponding to expression with the filter X starting from state Y ∈ {−(Θ − 1), …, − 1, 0, + 1, …, + (Θ − 1)}, then are the solutions of the matrix equations
2.26
where is the identity matrix, and the matrix and the vectors may be found elsewhere (Elliott, 2011). The mean expression probabilities and the mean expression time are then obtained by the usual moment expansion of the zeroth components, . Of course, as usual, the asymptotic expression rates are given by . Unlike the linear summation and low-pass filter models, in which and are sensitive only to the asymptotic induction rates , in the molecular filter model, these lowest-order moments depend on the full dynamics of ρ±in(t) for finite t and are thus sensitive to any non-Markovian (or semi-Markovian) characteristics of the induction processes. For Markov induction processes, or for non-Markovian induction processes that attain their asymptotic behaviors swiftly, so that we can to a good approximation replace the full, non-Markovian induction processes by Markov induction processes with rates , we can avoid the matrix inversion problem above and instead use standard methods to calculate splitting (here expression) probabilities and mean first passage (here expression) times. Writing and hY = π+Yh+Y + πYhY, we may determine π±Y and hY from the solutions of the recurrence relations,
2.27
2.28
for Y ⩾ 0, subject to boundary conditions π+−Θ = 0, π+ = 1; π−Θ = 1, π = 0; and h±Θ = 0. Although messy, explicit solutions for π±Y and hY, and therefore for and , may be given. We again regularize mean expression times as functions of pre- and postsynaptic firing rates by setting the decay rate η according to equation 2.15, as in the low-pass filter model. With this choice, we tabulate the first few asymptotic expression rates for Θ ⩽ 3 in order to give a flavor for how the terms develop. In general, we may write
2.29
where and are polynomials in Δρ and η′. The lowest few of these are given by
and
While the computation of these polynomials is routine, they rapidly become unwieldy for larger Θ, and their evaluation is best left to a symbolic algebra package. We again observe that equation 2.29 exhibits the required symmetry under the replacement , under which +Δρ ↔ −Δρ. Although the molecular filter model in the presence of general rather than specifically Markov induction processes ρ±in(t) is considerably more complicated, it too must treat the induction processes symmetrically.

3.  Synaptic Dynamics: General Approach

We now turn to an explicit formulation of the dynamical evolution of synaptic strengths driven by induction processes g±in(t) that are filtered to generate expression processes g±ex(t) leading to discrete steps in synaptic strength of sizes T±. Just as the induction ePDFs g±in(t) lead to a (possibly decaying) random walk of the filter variable X with step sizes ξ± = 1, the expression ePDFs g±ex(t) lead to a directed, filtered random walk in synaptic strength with basic step sizes T±. Our aim is to formulate the dynamics so that we can obtain the evolution of the probability distribution for synaptic strengths and, in particular, extract the stationary distribution, if it exists. A maximum of the stationary probability distribution corresponds to a state of synaptic strength that is an attractor of the synaptic dynamics. If there are multiple maxima, then these strength states will only be meta-stable because fluctuations in synaptic strength, even though suppressed by the filtering of induction steps, will drive transitions between them. We can then calculate the lifetimes of such meta-stable states explicitly as functions of the parameters defining the induction and expression processes. In this way, we can determine parameters for which these lifetimes are so large that for all practical purposes, they may be regarded as infinite and the strength states effectively completely stable against fluctuations. Of course, we do not expect to be able to perform all these calculations exactly, so approximations will be required. Our intention is to understand the qualitative features of simulations of synaptic plasticity in the presence of filtered induction processes such as those presented in Elliott and Lagogiannis (2009). Although we do not demand exact, quantitative agreement, we at least expect our approximations to be good enough to obtain order-of-magnitude agreement. We shall find that this is indeed the case.

We showed before that the derivatives of the rates ρ±(t), ignoring for the moment the distinction between plasticity induction and plasticity expression, constitute memory kernels in the non-Markovian random walk defining synaptic strength dynamics (Elliott, 2010b) or synaptic filter dynamics (Elliott, 2011). In those earlier calculations, we held the pre- and postsynaptic firing rates λπ and λp fixed in ρ±(t; λπ, λp). If the dynamics are in synaptic strength, then this corresponds to the assumption that the postsynaptic neuron is held in voltage or current clamp and that its firing rate is therefore uncoupled from synaptic transmission (Elliott, 2010b). Changes in synaptic strength under such a scenario report the underlying mechanisms of synaptic plasticity, but there is no feedback of strength changes into postsynaptic firing rate changes. On the other hand, if the dynamics are in the state of the synaptic filter, then the postsynaptic firing rate will indeed remain constant (subject to processes such as spike rate adaptation, for example) until the filter reaches threshold, provided that the presynaptic firing rate remains constant. In this second scenario, it is not necessary to regard the postsynaptic neuron as being held in voltage or current clamp (Elliott, 2011). Here, however, we wish to study the full, natural dynamics of synaptic strength with postsynaptic firing coupled to synaptic transmission in order to determine the probability distribution of synaptic strength and thus the meta-stable strength states to which these dynamics are attracted under normal stimulus paradigms. The postsynaptic firing rate λp must therefore not be held constant, but must change as synaptic strengths change. The expression kernels ρ±ex(t; λπ, λp) will therefore depend on synaptic strength through λp. This transforms a homogeneous random walk into an inhomogeneous one, even in the absence of boundaries or constraints on synaptic strength.

Consider for the moment a single synapse of strength S. Let P(S; t) denote the probability distribution of S at time t. We consider only the simplest model of the postsynaptic firing rate, so that the postsynaptic firing rate is a synaptic-strength-weighted sum of presynaptic firing rates. For a single presynaptic neuron synapsing on a single postsynaptic neuron, we therefore write λp = Sλπ. The master equation governing the evolution of P(S; t) can then be written as
3.1
The strength S should not be allowed to become negative. It is therefore convenient to erect a reflecting boundary at S = 0 as argued elsewhere (Elliott, 2010b). This is most easily achieved by setting the strength step sizes T± equal, so we write T± = T. The boundary form of the master equation then becomes
3.2
For general but sufficiently small T±, we can take the diffusion limit of equation 3.1 and obtain the non-Markovian Fokker-Planck equation,
3.3
For small T± = T, we can also obtain the continuum version of equation 3.2.

With a model of induction processes ρ±in based on our tristate switch model of STDP, it is not necessary to erect a reflecting boundary corresponding to an upper strength limit because the model is intrinsically self-stabilizing in synaptic strength (Appleby & Elliott, 2006; Elliott, 2008). This stands in contrast to many conventional models of STDP in which soft or hard upper bounds in synaptic strength are imposed to prevent synaptic strengths from diverging (see, e.g., Song, Miller, & Abbott, 2000; Rubin, Lee, & Sompolinsky, 2001) or in which the STDP update rule is strength dependent, motivated either experimentally (see, e.g., van Rossum, Bi, & Turrigiano, 2000) or theoretically in order to introduce fixed-point dynamics (see, e.g., Gütig, Aharonov, Rotter, & Sompolinsky, 2003; Burkitt, Meffin, & Grayden, 2004). The interactions between bounds on synaptic strength and strength-dependent STDP update rules in conventional models can dramatically influence the final distribution of synaptic strengths (Rubin et al., 2001). In our tristate switch model of STDP, however, stable nonzero fixed points emerge not via any of these methods, but rather because of the stabilizing influence of multispike interactions (Appleby & Elliott, 2006; Elliott, 2008). Of course, it is legitimate to ask whether the value of synaptic strength to which synapses tend in our model is consistent with experimental data. However, this is a question of an entirely different character from that concerning the implementation of bounds on synaptic strength in conventional models of STDP that would otherwise be unstable. The fact is that at least theoretically, upper bounds in synaptic strength are not required by our tristate switch model of STDP, and we consider this to be a virtue, not a vice. Hence, we do not examine the impact of upper-strength limits here, instead preferring to require that any models that we consider are intrinsically self-stabilizing.

Although the strength dynamics of the connection between a single afferent cell and a single target cell are illustrative of more general dynamics, we are actually interested in the interactions between multiple afferent cells synapsing on a target cell. In particular, we are interested in the regimes in which stable, segregated patterns of afferent connectivity emerge on the target cell, driven by activity-dependent competition between multiple afferents (Purves & Lichtman, 1985). Consider, therefore, m afferents with synaptic strengths Si and activities , i = 1, …, m. We write the vectors of strengths and presynaptic activities as and , respectively. The postsynaptic firing rate is then , where “·” denotes the dot product. Let be an m-dimensional vector in which all components except the ith component are zero and the ith component is unity, so that eij = δij, the Kronecker delta. The set forms the basis set of an m-dimensional Cartesian space. Then it is extremely tempting to replace equation 3.1 for m = 1 afferent by the generalized form
3.4
which is a master equation based on a possible one-step potentiating or depressing expression step for each of the m afferents. We will explain the meaning of the modified equality symbol “” shortly. In writing this equation, we have tacitly assumed that changes in strength over different synapses are independent. We will scrutinize and justify this assumption in section 6. Equation 3.1 for one afferent is an equation that ultimately owes its origin to the assumption of renewal. Without renewal, we could not tame the intractable complexity of the full, non-Markovian processes. In considering the dynamics of a single synapse, renewal induction processes limit the history dependence of the induction processes to the last induction step. Renewal of induction processes also ensures renewal of expression processes. However, other than commonality of the postsynaptic firing rate, the induction processes in different synapses are independent, and except for exponentially distributed and therefore memoryless processes, independent processes cannot be renewal processes. Thus, in general, an induction step in one synapse will not renew the induction processes in all the other synapses. A dependence on the induction history before the last induction step will thus inevitably leak back into the dynamics over all synapses. Equation 3.4 is therefore generally incorrect, except when the expression processes g±ex are Markov processes and therefore memoryless and renewing. The modified equality sign “” means equality only for Markov expression processes. In this case, we have that , where the rates of the Markov processes are no longer explicitly time dependent; δ(τ) is the Dirac delta function. Of course, the Markov rates continue to be implicitly time dependent through the time dependence of the strengths , but lose their explicit support away from τ = 0 when g±ex are Markov processes. The integral in equation 3.4 therefore collapses because the memory kernels develop amnesia. Equation 3.4 then transforms into
3.5
The Fokker-Planck equation in equation 3.3 becomes
3.6
for Markov expression processes.

In order to proceed, we must make the approximation that the expression processes are Markov processes, so that we can employ equation (3.5), or its equivalent, equation (3.6), in the continuum limit. Of course, even with Markov induction processes, any nontrivial filtering process will lead to non-Markovian expression processes. The simplest and most natural approximation to make is to replace the non-Markovian expression processes having asymptotic rates by Markov expression processes having rates . The severity of this approximation will depend on the timescales on which the non-Markovian expression processes with rates r±ex(t) attain their asymptotic rates relative to the mean expression times . So as not to distract from the overall flow of our discussion, in the appendix we show that the timescale governing the approach to asymptotic behavior is typically smaller than the mean expression time, so that the approximation is not a bad one.

So far we have been careful to distinguish between time intervals separating the occurrence of plasticity expression steps, during which the postsynaptic firing rate remains constant provided that the presynaptic firing rates remain constant, and longer time intervals over which changes in synaptic strength must be allowed to feed back directly into changes in the postsynaptic firing rate. We must now consider, however, the fact that presynaptic firing rates will change on timescales typically much shorter than the characteristic mean plasticity expression time. In the linear summation model, can be on the order of 10 s or higher (see the appendix), and in the low-pass and molecular filter models, the introduction of the decay parameter η′ increases compared to that in the linear summation model for the same threshold, Θ (Elliott, 2011). The rate of plasticity expression is thus around 0.1 Hz or lower, while the rate of plasticity induction is, in our tristate switch model of STDP, around 20 Hz (see the appendix). A typical stimulus regime, however, might consist of, for example, a few saccades per second, so that activity patterns might be expected to change a few times per second—perhaps at least once per second—or at a stimulus rate of 1 Hz. Thus, induction processes, at least in our tristate switch model of STDP, occur on timescales an order of magnitude faster, while expression processes occur on timescales an order of magnitude slower. For induction processes, the stimulus may therefore be regarded as constant, but for expression processes, the stimulus will change several times between expression steps, causing the rate of induction to change several times between expression steps too. In principle, we could establish a master equation formalism governing the probability distribution of stimulus patterns, but the resulting analysis of synaptic dynamics, involving multiple stochastic processes on multiple timescales, would become hopelessly complicated. We therefore make what is in effect an adiabatic approximation. We have seen that the asymptotic rates of plasticity expression can be written in the form . It is these rates that we will use in equations 3.5 and 3.6 when approximating the non-Markovian expression processes by Markov processes. We now make the further, adiabatic approximation,
3.7
where the angle brackets 〈〉 denote an average of the induction rates over the ensemble of stimulus patterns defined by the presynaptic activity vector . In making this second approximation, we of course lose a source of fluctuations by neglecting the fluctuations in the dynamics due to changing activity patterns. We make two remarks. First, we expect the fluctuations due to the induction and expression processes to dominate those due to the stimulus process, especially when the plasticity step sizes T± are sizable. Indeed, it was precisely the problems with fluctuations in synaptic strength arising from our tristate switch model of STDP that led us to consider plasticity expression as a filtering process that suppresses fluctuations in synaptic strength (Elliott, 2008; Elliott & Lagogiannis, 2009). Second, for the purposes of analytical tractability, we typically consider that the presynaptic activity vector fluctuates around its stationary mean in a manner that is controlled by a parameter that we take to be small, so that we need to expand up to second-order moments only (Appleby & Elliott, 2006; Elliott, 2008). We can then make the fluctuations due to activity (rates) as small as we please while retaining the qualitative features of the synaptic dynamics of interest to us. These remarks of course depend on the coding scenario. For nonstationary, non-Poisson presynaptic firing processes, rate fluctuations may be as important as the intrinsic fluctuations due to the induction and expression of synaptic plasticity. However, by restricting to stationary Poisson processes, we expect the intrinsic fluctuations to dominate the stimulus-driven fluctuations.
As we are interested in activity-dependent, competitive dynamics during normal rather than abnormal development (Purves & Lichtman, 1985), we set
3.8
so that all afferents have a common mean activity μ, and the αi are parameters governing the fluctuations around the common mean, so that 〈αi〉 = 0, ∀i. For analytical purposes, we will consider |αi| ≪ 1, but for numerical results, this assumption is not required. We also assume a common variance ς2 in activity for all afferents, so that the afferents are treated in an unbiased manner. Defining q = ς/μ, we then may write
3.9
Finally, defining the correlation coefficient, denoted by Cij, between two afferents i and j in the usual way, we have
3.10
We will for convenience treat all afferents symmetrically in terms of their second-order moments, so we finally assume that Cij = C for ij, so that there is a common pairwise correlation coefficient between all pairs of afferents. The two assumptions of unbiased, symmetrical second-order moments are not vital to our subsequent analysis but considerably simplify it. The assumption of common means is vital to considering normal development, in which we would not expect any single afferent to be privileged.
We may finally consider averaging a two-argument function over such an activity ensemble. Because
3.11
where the vector has components xi, and ∇i = ∂/∂xi, we have
3.12
where we have retained Cij for generality, and means that the function and its derivatives are evaluated at xi = μ. Since μ2q2 = ς2 is by assumption small, at least for analytical purposes, the averaging involves a Taylor expansion that can be truncated at second order. Suppose that the ensemble in fact involves only the 2m activity patterns with αi ∈ {−q, + q}, ∀i. Then, taking into account the covariance, the average can be explicitly written as
3.13
as above in equation 3.12, although of course the higher-order terms in general will differ. Thus, provided that the fluctuations about the mean are small, |αi| ≪ 1, explicitly averaging over the ensemble αi ∈ {−q, + q}, ∀i is entirely equivalent to a formal, analytical ensemble average in which the first- and second-order statistical properties of the ensemble are defined by 〈αi〉 = 0, 〈α2i〉 = q2, and 〈αiαj〉 = Cijq2, but for which its elements are not explicitly specified. An identical conclusion holds for a general multiargument function, not just a two-argument function. This observation regarding ensemble averaging allows us to consider simulations involving only a small set of possible activity patterns with the assurance that their results are not restricted to this set but also apply to generic, statistically defined ensembles, provided that the activity fluctuations are adequately small.

4.  Synaptic Dynamics: Fokker-Planck Equation Approach

We now consider the dynamics of synaptic strength by employing the Fokker-Planck equation in equation 3.6. We have seen that the asymptotic rates of expression processes in the linear summation and low-pass filter models are sensitive only to the asymptotic rates of induction processes. This is not true in the molecular filter model, but for this model, we will work for simplicity in the approximation in which the non-Markovian induction processes are replaced by Markov processes of rates . This allows all three filter models to be treated simultaneously in a common framework. Analogous to , we define . Writing
4.1
4.2
where these quantities do not depend explicitly on because of the ensemble averaging, the Fokker-Planck equation in equation 3.6 becomes
4.3
If at some initial time t0, then a standard calculation shows that and determine the subsequent evolution of the mean and variance in synaptic strength, and define the so-called jump moments of the process. For completely deterministic dynamics in which we discard the second-order terms in the Fokker-Planck equation, we have the Liouville equation,
4.4
where is a vector with components , and ∇i = ∂/∂Si. The synaptic strength then evolves according to the deterministic equation
4.5

4.1.  Deterministic Dynamics.

In previous work, we have extensively analyzed the deterministic dynamics governed by equation 4.5 in our tristate switch model of STDP, but only for the case in which each induced plasticity step is immediately expressed without filtering, so that ε(u, v, η′) ≡ u, corresponding to Θ = 1 (Appleby & Elliott, 2006; Elliott, 2008). In such a case, these dynamics induce the existence of fixed points in corresponding to final, segregated states of afferent connectivity in which all but one afferent has zero strength. Using the basis vectors , i = 1, …, m, with components eij = δij defined earlier, these fixed points form the set , where is the common, overall strength of the surviving afferent. Strictly, these are quasi-fixed points, because in previous work, we truncated Si, any i, at zero (but did not hold it there) if the dynamics try to drive it negative. However, if T± = T, then this procedure corresponds precisely to the erection of reflecting boundaries along hyperplanes in which any component of is zero. With T± = T, the vectors are therefore genuine fixed points of the dynamics with reflecting boundaries. The stability (and, indeed, existence) of these fixed points depends on the parameters defining our tristate switch model of STDP. By stability, here for deterministic dynamics, we do not mean stability against fluctuations, but stability only in the deterministic dynamical systems sense.

We have characterized the parameter regimes in our model that lead to stable, segregated fixed points in the nonnegative hyperquadrant (Appleby & Elliott, 2006; Elliott, 2008). As typical examples of these competitive dynamics in the unfiltered model (Θ = 1), in Figure 2 we show the flows in the S1S2 plane for m = 2 afferents for a choice of parameters guaranteeing the presence of competition. This choice corresponds to the standard set used elsewhere (Elliott, 2008, 2011; Elliott & Lagogiannis, 2009), for which we take μ = 50 Hz, λ = 50 Hz, T/T+ = 0.9, and γ ≡ (T+λ)/(Tλ+) = 0.6 or λ+ = 50/(0.6 × 0.9) ≈ 93 Hz. We have taken three different values for the correlation coefficient, C = 0 and . We have also taken the nonstandard value , and rather than employing a statistically-defined ensemble using equation 3.12, we have explicitly averaged over the four activity patterns defined by αi ∈ {−q, + q}, i = 1 and i = 2, as in equation 3.13, but without the expansion to second order in q as q is not small. We have taken a large value of q in order to illustrate the fact that in general, the separatrix partitioning the nonnegative quadrant (hyperquadrant for m>2) into two distinct regions is not a straight line (hyperplane for m>2). This is especially clear in Figure 2A with . For the larger values of C, however, the separatrix is close to a straight line. For our “standard” choice of q, which is , the separatrix is always close to a hyperplane for general m, even for strongly negative C.

Figure 2:

The flows in equation 4.5 for two afferents of strengths S1 and S2 for the deterministic dynamics of the tristate switch model of STDP in the absence of an expression filter. Flows are shown for three choices of the correlation coefficient between the two afferents’ activities, corresponding to (A) , (B) C = 0, and (C) . Most other parameters correspond to our hitherto “standard” choices. In particular, we set λ = 50 Hz, T/T+ = 0.9, and γ ≡ (T+λ)/(Tλ+) = 0.6, which fixes λ+ ≈ 93 Hz. Activity patterns correspond to mean firing rates μ = 50 Hz, but the standard deviation is set by the nonstandard value , compared to the standard value . Averaging is therefore performed as an explicit average over four activity patterns rather than by statistical averaging via a Taylor expansion in q.

Figure 2:

The flows in equation 4.5 for two afferents of strengths S1 and S2 for the deterministic dynamics of the tristate switch model of STDP in the absence of an expression filter. Flows are shown for three choices of the correlation coefficient between the two afferents’ activities, corresponding to (A) , (B) C = 0, and (C) . Most other parameters correspond to our hitherto “standard” choices. In particular, we set λ = 50 Hz, T/T+ = 0.9, and γ ≡ (T+λ)/(Tλ+) = 0.6, which fixes λ+ ≈ 93 Hz. Activity patterns correspond to mean firing rates μ = 50 Hz, but the standard deviation is set by the nonstandard value , compared to the standard value . Averaging is therefore performed as an explicit average over four activity patterns rather than by statistical averaging via a Taylor expansion in q.

In the presence of a nontrivial filter expression function ε(u, v; η′) ≠ u (Θ ≠ 1), the deterministic dynamics are somewhat, although not substantially, modified, so we must examine the fixed-point structure and stability arguments for a general filter. We expand to second order in q. First, we have that
4.6
where by ∂/∂x and ∂/∂y, we mean derivatives with respect to the first and second arguments of ϱ±in(x, y), and the derivatives are evaluated at x = μ and . The vector is an m-dimensional vector all of whose components are unity. We then have
4.7
where here, by ∂/∂u and ∂/∂v, we mean derivatives with respect to the first and second arguments of ε(u, v; η′), and the derivatives are evaluated at and . Consider the leading order, O(q0) terms in . We have
4.8
independent of i. The fixed points are solutions of . We have seen that the functions ε(u, v; η′) in the various filter models can be very complicated. In the linear summation model, we can, in fact, solve exactly, but for other filter models, in general we require the solutions to transcendental equations to obtain the fixed points even at O(q0). Notice, however, that if we set the plasticity step sizes equal, T± = T, then at O(q0) when . This is the equation that governs the location of the fixed points in the absence of an expression filter, that is, when ε(u, v; η′) = u. Thus, for T± = T, an expression filter (with Θ>1) does not modify the location of the fixed points, at least at O(q0), induced by the induction processes ϱ±in (the Θ = 1 case). Henceforth, we therefore set T± = T. This is also very conveniently the condition that enables the erection of reflecting boundaries in the master equation formulation to prevent synaptic strengths from becoming negative, and ensures that the truncation of synaptic strengths at zero in simulations corresponds precisely to the imposition of reflecting boundaries.

At O(q0), we see from equation 4.8 that only the quantity appears in . An entire hyperplane therefore corresponds to the fixed-point manifold. This is because at O(q0), all afferents have the same firing rates and are indistinguishable. In equation 4.7, however, which controls the O(q2) corrections to the locations of the fixed points, we see terms controlled by , which vanishes when C = 1. Of course, when C = 1, the afferents are again indistinguishable since their firing rates are perfectly correlated. For C < 1, the terms will deform the hyperplane. The separatrix is therefore governed in general by both and , but for small q, the deformation of the separatrix from a hyperplane will be small. We also see in equation 4.7 the appearance of pure Si terms, again multiplied by 1 − C. When C < 1, these Si terms break the otherwise perfect symmetry between the afferents, and the fixed-point manifold collapses to a set of m points corresponding to m segregated states for appropriate choices of parameters. These observations were made previously for the tristate switch model of STDP (Appleby & Elliott, 2006), and therefore with Θ = 1, but here they are extended to dynamics in the presence of a general expression process ε(u, v; η′), with Θ>1.

It is now convenient, although not necessary, to restrict to an analysis of competition between m = 2 afferents. Competitive dynamics between m>2 afferents are not fundamentally different from those between m = 2 afferents, but for small q, we must consider the hyperplane while for m = 2, this is just a line in the S1S2 plane. More important, we will see in section 5 that the dominant mode by which fluctuations drive a change in meta-stable states involves only two afferents: the afferent that is currently strong and just one of the m − 1 currently weak afferents. Assuming that well-segregated states exist and are sufficiently stable to be identified as such, the development of fluctuations in the weak afferents will occur with low probabilities and, roughly speaking, be approximately independent of each other. Thus, we expect fluctuations involving multiple weak afferents to arise with much lower probability than those involving just a single weak afferent. Of course, if the firing rates between a pair of weak afferents are strongly correlated, then this argument will break down, because the development of fluctuations in these two afferents would not be expected to be approximately independent. However, in this case, we in effect have a single afferent. If all afferents’ firing rates are strongly correlated, then we would expect fluctuations to dominate, and indeed, the lifetimes of segregated states are so short that such well-defined states cannot in practice be identified or observed.1 In summary, provided that well-segregated states exist and are adequately stable, fluctuations involving multiple weak afferents are expected to be suppressed compared to those involving just one weak afferent. We will verify this explicitly in section 5 when we construct stationary solutions of the master equation for m = 3 afferents. We may conclude, then, that in order to examine stability against fluctuations in section 4.2, we need only consider m = 2 afferents. Our strategy is to change synaptic variables from to and the m − 1 mutually orthogonal variables , i = 1, …, m − 1, where , ∀i and , for ij. For small q, the competitive dynamics reside almost exclusively in the variables Si. For m = 2, there is only one such orthogonal variable, SS1S2. Thus, we restrict to a consideration of the dynamics in S+ = S1 + S2 and S = S1S2.

At O(q0), the fixed-point manifold is given by the solution of the equation ϱ+in(μ, μS+) = ϱin(μ, μS+) in only the variable S+, independent of S. Let S+0 be this solution. We then need to find the O(q2) correction to S+0 to locate the segregated fixed points to O(q2). These fixed points are at for some , and because of the nature of the reflecting boundaries, must be determined by solving, for example, . Writing , we can expand to O(q2) using equation 4.7. After some algebra, we find that in a compact notation,
4.9
where the derivatives are evaluated at x = μ and y = μS+0. The locations of the fixed points therefore depend on the induction processes ϱ±in(x, y) (the case Θ = 1) but not on the expression process ε(u, v; η′) (when Θ>1), at least to O(q2).
To determine the stability of these fixed points, we examine the dynamics of S. We find that to O(q2), S evolves according to the equation
4.10
where the derivatives of ϱ±in(x, y) are as usual evaluated at x = μ and y = μS+0, and the derivatives of ε(u, v; η′) are evaluated at u = ϱ+in(μ, μS+0) and v = ϱin(μ, μS+0). Because ϱ+in(μ, μS+0) = ϱin(μ, μS+0) from the definition of S+0, both derivatives of ε(u, v; η′) are evaluated at a common value. To O(q2), we see that in the evolution of S, the dependence on the induction processes ϱ±in(x, y) and the dependence on the expression process ε(u, v; η′) factorize. We can therefore write dS/dt in the form
4.11
where Gin(S+0) and Gex(S+0) are two induction and expression functions, respectively. Consider for the moment a filterless process, so that ε(u, v; η′) = u and (∂/∂u − ∂/∂v)ε = 1. This corresponds to setting Θ = 1. For competition to be present and for the fixed points to be stable, we require that the line S = 0, or S1 = S2, is unstable. Hence, for C < 1, we require
4.12
so that the induction processes must drive competition. Notice that equation 4.10 is linear in S, and if S = 0 is unstable, then S could in principle grow unboundedly. However, the reflecting boundaries will prevent |S| from growing greater than , at which point will have reached one of the fixed points . With a general expression function ε(u, v; η′) (thus now with Θ>1) and competition driven by the induction processes, we must have
4.13
Indeed, for the expression processes ε(u, v; η′) in the three filter models considered above, this requirement is always satisfied. In general, the expression term (∂/∂u − ∂/∂v)ε slows the approach toward the fixed points, but does not fundamentally alter the deterministic dynamics governed by the induction processes.
The advantage of the generality of our approach above is that equations 4.9 and 4.10 hold for general induction processes ϱ±in(x, y) and a general expression process ε(x, y; η′). In particular, these results are not specific to our tristate switch model of STDP. For this model, however, using equations 2.5 and 2.6, we find that , determining the location of the segregated fixed points, is given at O(q2) by
4.14
At O(q0), we require λ + μ>λ+ for the fixed hyperplane to intersect the positive hyperquadrant. For μ sufficiently small and λ+, therefore, the dynamics can transform into being entirely depressing, as discussed elsewhere (Elliott, 2008). The function Gin(S+0), with S+0 = (λ + μ − λ+)/μ, is given by
4.15
One choice of parameters to ensure that Gin(S+0)>0 is λ+, corresponding to γ ≡ (T+λ)/(Tλ+) = λ+ < 1 as elsewhere (Appleby & Elliott, 2005, 2006; Elliott, 2008), (Appleby & Elliott, 2005) (Appleby & Elliott, 2006) (Elliott, 2008) and also 2(λ + μ)>λ+. If λ + μ>λ+ for S+0>0, then this latter condition is automatically satisfied.
We also give explicit results for the expression function Gex(S+0) for the three filter models discussed above. For the linear summation model, we have the particularly simple result,
4.16
For the low-pass filter model, we have
4.17
For the molecular filter model, results involve messy polynomials in η′. We give the first few:
4.18
4.19
4.20
The general forms are unwieldy, but for Θ large enough we expect that continuum limit low pass filter model to be a good approximation to the discrete molecular filter model, as discussed elsewhere (Elliott, 2011). As Θ increases in all three filter models, Gex(S+0) decreases, as we should expect. However, as η′ increases in both the low pass and molecular filter models, Gex(S+0) undergoes a transition from being roughly constant for small η′ to being rapidly suppressed. This behaviour coincides with the collapse of the expression rates ϱ±ex in equations 2.24 and 2.25 for the low pass filter model, and is evident in the trend in equation 2.29 for the molecular filter model. For η′ large enough in these two models, plasticity essentially stops, and dS/dt is therefore negligibly small.

4.2.  Stability Against Fluctuations.

Having characterized the deterministic dynamics in the presence of an expression process ε(x, y; η′), we can now turn to an examination of the impact of turning on fluctuations in synaptic strength by considering the full Fokker-Planck equation in equation 4.3 rather than the Liouville equation in equation 4.4. Turning on fluctuations will have the effect of transforming stable fixed points of the deterministic dynamics into maxima of the probability distribution. Fluctuations can then drive transitions between meta-stable segregated states of synaptic strength. In this subsection, we therefore calculate the mean first passage time for a transition between two segregated states. As discussed earlier, this is the dominant mode by which fluctuations drive a change of state, so it is sufficient to restrict to m = 2 afferents. We can then use the derived formula for the mean transition time to determine the dependence of this time on the step size parameter T. Previously, we have found that in an unfiltered model (corresponding to Θ = 1), T must be unrealistically small to guarantee an adequately long, developmentally appropriate transition time (Appleby & Elliott, 2006; Elliott, 2008). Introducing a filter allows more realistic, larger plasticity steps T to be employed, but we were content previously with results only from simulation (Elliott & Lagogiannis, 2009). Using the Fokker-Planck equation in equation 4.3 allows a more thorough and revealing treatment.

When states are already segregated, fluctuations drive a change to a different state of segregation along the separatrix. Trajectories that veer off the separatrix are suppressed compared to those that remain on it. For m>2 afferents, the separatrix is an (m − 1)-dimensional manifold, but in this case, the fluctuation trajectories do not cover the entire manifold. Rather, they are restricted to the one-dimensional paths on the separatrix that connect all possible pairs of segregated fixed points. For m = 2 afferents, these paths and the separatrix precisely coincide. For small q, the separatrix is very close to a line, or a hyperplane, and we may characterize it as S+ = S+0 at O(q0) or S+ = S+0 + q2ΔS at O(q2), with ΔS given by equation 4.9. Along a fluctuation trajectory, S+ therefore remains (approximately) constant, and the fluctuations drive a change only in S, for m = 2. Even for m>2, the dynamics are essentially one-dimensional, with S corresponding to the difference in strength between any pair of afferents that swap control of the target cell.

In the variables S+ and S, the Fokker-Planck equation in equation 4.3 transforms into
4.21
where Mi, Vi and P on the right-hand side are now functions of S± (and t), and the Jacobian factor of is absorbed into P(S+, S; t). However, because small q fluctuations drive a change in S and not in S+, we need only consider the one-dimensional Fokker-Planck equation in S, obtained by dropping any terms with derivatives in S+:
4.22
The first-order derivative term governs the deterministic dynamics in the absence of fluctuations, and we already know its form, from equation 4.11. If S = 0 is an unstable line, then equation 4.11 describes motion in an inverted quadratic potential,
4.23
so that M1M2 = −dU(S)/dS in equation 4.22. Because of the reflecting boundaries, we may erect infinite potential barriers to confine S in the interval , so that U(S) = +∞ for .
The variable S must climb over this inverted quadratic potential in order to change sign and thus move from one segregated state to another. These fluctuations are controlled by the diffusion term on the right-hand side of equation 4.22. Since the potential is already O(q2), in the calculation of the transition time it suffices to evaluate the V1 + V2 term to O(q0). It is given, using equation 4.7, simply by
4.24
using the fact that ϱ+in(μ, μS+0) = ϱin(μ, μS+0). Thus, V1 + V2 is, to O(q0), constant on the O(q0) fixed hyperplane , and independent of S. We therefore write V1 + V2 = σ2. Because ε(x, y; η′) is evaluated here in the driftless induction rate case x = y, the overall scale of ε(x, y, η′) must be set by the common induction rate x = y. In the driftless limit, ε(x, x; η′) should therefore factorize into an overall scale set by x and a function whose form is entirely determined by the details of the expression process. Hence, we write σ2 = 4T2ϱ+in(μ, μS+02ex. In the tristate switch model of STDP, we have
4.25
and ϱ±in(μ, μS+0)>0 precisely when S+0>0, or λ + μ>λ+. For σ2ex, in the linear summation model, we have
4.26
and for the low-pass filter model,
4.27
Results are, as usual, messy for the molecular filter model. For small Θ, we have
4.28
4.29
4.30
Notice that while an expression filter reduces the overall size of the potential in equation 4.23 because of the presence of the additional term Gex(S+0) < 1 compared to a filterless (Θ = 1) process, an expression filter suppresses even more strongly the diffusion term in equation 4.22, controlled by σ2ex. Roughly speaking, while Gex(S+0) ∝ 1/Θ, we have that σ2ex ∝ 1/Θ2. The reduction in the size of the potential means that we approach the fixed points of the deterministic dynamics more slowly, but the even greater reduction in the diffusion term by the presence of the synaptic filter means that these fixed points possess increased stability against fluctuations.
Writing , where Λ = T(1 − C2q2Gin(S+0)Gex(S+0), we may determine the mean first passage time for a transition between S = ±S+0 and S = ∓S+0. We have
4.31
where erf is the error function and ΔU is the size of the potential barrier opposing the transition, given by . Notice that this result differs somewhat from the standard formula for the mean first passage time between meta-stable states, which involves the Arrhenius factor exp(2ΔU2). This is because the standard formula is derived by making a parabolic approximation of the potential around its extrema, while the above evaluation of the integrals is exact since, to O(q2), U(S) is strictly quadratic in q, and the reflecting boundaries at |S| = S+0 allow τtrans to be evaluated in terms of special functions without approximation. Equation 4.31 is derived for m = 2 afferents, so that there is only one fluctuation mode that drives a change in state. For m>2 afferents, for any segregated state, there are m − 1 possible, dominating modes that drive changes in state to the other m − 1 segregated states. Thus, in general, the right-hand side of equation 4.31 should be divided by m − 1 to generalize it to m>2 afferents, if all afferents have a common, pairwise correlation coefficient C.
In simulations of our tristate switch model of STDP, we have not explicitly measured the mean lifetimes of segregated states of connectivity, but rather have quantified their stability using a segregation index (Appleby & Elliott, 2006; Elliott & Lagogiannis, 2009). For m = 2 afferents, this is just SI ≡ (S1S2)/(S1 + S2). Although the mean lifetimes can be estimated well when they are small, large mean lifetimes will be very hard to estimate in simulation because of the need to run many very long-lasting simulations to obtain decent averages. Using SI, we can instead average SI over a simulation for a fixed period of time after allowing an initial period during which segregation occurs. If segregation is stable against fluctuation-induced transitions during this later probe period, then its average will be close to +1 or −1, but if segregation is unstable, then its average will be close to 0. We form a mean segregation index (MSI) by averaging SI over the probe period during a simulation and then by averaging the moduli of these averages over many simulations. If the MSI is close to 1, then we can at least say that τtrans exceeds the duration of this probe period. In fact, we can obtain a better lower bound on τtrans. Suppose that the probe period is of duration τprobe and consider a sequence of alternating intervals of length τtrans in which we assign to even-numbered intervals a value SI = +1 and to odd-numbered a value SI = −1. Considering that our probe interval may fall anywhere in this sequence, we easily find that MSI for τtrans ⩾ τprobe. If τtrans ≫ τprobe, then the MSI will be very close to 1, and it will be difficult in simulations to measure accurately its deviation from 1 for the same reason that it is difficult to measure τtrans when it is very large. Thus, for measured MSIs close to 1, we may conclude that
4.32
but not strict equality. For MSI =0.95, we may infer that τtrans is at least an order of magnitude longer than τprobe.

Figure 3 illustrates the dependence of the MSI on the plasticity step size T in our tristate switch model of STDP in the presence of the linear summation expression filter, for different thresholds Θ. As usual, the Θ = 1 case corresponds to a filterless process. Although T± = T (contrary to our hitherto standard choice of T/T+ = 0.9), all other parameter choices are standard. Increasing the threshold Θ enables larger plasticity steps to be employed without compromising the stability of segregated states. In these simulations, we have used a probe period τprobe = 107 s. We compare these simulation results for the MSI with the mean transition time in equation 4.31 by reading off, for a given value of Θ, the value of T at which the MSI =0.95. Although τprobe = 107 s, the mean lifetimes for MSI =0.95 must, by the above argument, be at least 10 times longer.

Figure 3:

The mean segregation index (MSI) for two afferents plotted against the plasticity step size parameter, T, in the tristate switch model of STDP in the presence of the linear summation expression filter, for different thresholds, Θ. From left to right, the thresholds correspond to Θ = 1, 2, 3, 5, 10, 15, 20, 25, and 30. Means are obtained by averaging over 80 spike-based simulations. The probe period in which the segregation index is measured is of duration τprobe = 107 s, and an initial period of the same duration is allowed for segregation to occur. Simulations are run with our new “standard” set of parameters, in which we take T/T+ = 1 instead of T/T+ = 0.9. This new choice means that simulations implement correctly reflecting boundaries when an afferent has zero synaptic strength. The parameters are: λ = 50 Hz; T/T+ = 1; γ ≡ (T+λ)/(Tλ+) = 0.6; μ = 50 Hz; C = 0; . Because q is small, the use of only four activity patterns μ(1 ± q) for two afferents in simulation is equivalent to employing a statistically defined ensemble average.

Figure 3:

The mean segregation index (MSI) for two afferents plotted against the plasticity step size parameter, T, in the tristate switch model of STDP in the presence of the linear summation expression filter, for different thresholds, Θ. From left to right, the thresholds correspond to Θ = 1, 2, 3, 5, 10, 15, 20, 25, and 30. Means are obtained by averaging over 80 spike-based simulations. The probe period in which the segregation index is measured is of duration τprobe = 107 s, and an initial period of the same duration is allowed for segregation to occur. Simulations are run with our new “standard” set of parameters, in which we take T/T+ = 1 instead of T/T+ = 0.9. This new choice means that simulations implement correctly reflecting boundaries when an afferent has zero synaptic strength. The parameters are: λ = 50 Hz; T/T+ = 1; γ ≡ (T+λ)/(Tλ+) = 0.6; μ = 50 Hz; C = 0; . Because q is small, the use of only four activity patterns μ(1 ± q) for two afferents in simulation is equivalent to employing a statistically defined ensemble average.

Thus, s for an MSI =0.95 in Figure 3. We therefore solve equation 4.31 to obtain the step size T that gives τtrans = 108 s. These results from equation 4.31 and those obtained by reading off T from Figure 3 are shown in Figure 4. The trend in the analytical results follows the trend in the simulation results perfectly. Thus, in spite of the two approximations made above, our analysis has captured the trends exhibited by the full, unapproximated system. This perfect qualitative agreement is both surprising, given the approximations, and satisfying. We may thus use the analysis developed above to examine the meta-stability of patterns of synaptic connectivity rather than rely on time-consuming simulations as in previous work. Further, we have remarkably good quantitative agreement, given both the approximations and the ambiguities in defining and measuring τtrans. Indeed, also shown in Figure 4 are analytical results in which we simply divide the derived analytical result for T by a constant numerical factor of 1.2. The agreement with the simulation results is almost exact. There is a discrepancy, then, of an overall, constant factor of just 20% between analytical results and simulation results, which is well within the order of magnitude target that we set ourselves above.

Figure 4:

The step size T that gives a minimum lifetime for segregated states of 108 s, as a function of Θ. Results are based on the tristate switch model of STDP and the linear summation expression filter. The solid line denotes the value of T read off from Figure 3 that gives an MSI = 0.95, leading to a lifetime in excess of 108 s from equation 4.32. The dashed line represents the analytical result for T obtained from the mean transition time between segregated states in equation 4.31. Because of the approximations required to derive this result, we would not expect exact agreement. However, dividing the analytical result by a constant factor of 1.2, shown by the dotted line, gives almost perfect agreement with the results from simulation.

Figure 4:

The step size T that gives a minimum lifetime for segregated states of 108 s, as a function of Θ. Results are based on the tristate switch model of STDP and the linear summation expression filter. The solid line denotes the value of T read off from Figure 3 that gives an MSI = 0.95, leading to a lifetime in excess of 108 s from equation 4.32. The dashed line represents the analytical result for T obtained from the mean transition time between segregated states in equation 4.31. Because of the approximations required to derive this result, we would not expect exact agreement. However, dividing the analytical result by a constant factor of 1.2, shown by the dotted line, gives almost perfect agreement with the results from simulation.

Figure 4 shows the value of the step size T that produces an analytical mean transition time of 108 s in equation 4.31 for the linear summation model, as a function of the threshold Θ. In comparison, Figure 5A shows equivalent analytical results for the low-pass filter model for different choices of the decay parameter η′. We see the same initial trend in Figure 5A as in Figure 4 for values of Θ below around 30 for η′ = 0.01, 20 for η′ = 0.02, 15 for η′ = 0.04, and 10 for η′ = 0.08. For larger values of Θ, however, the behaviors differ, with the required step size T diverging. Of course, any large values for T derived from formulas relying on the Fokker-Planck equation, and therefore on the diffusion limit that requires T to be small, must be interpreted with caution. However, the onset of the divergence in T owes its origin to the onset of the collapse of the expression rates in the low-pass filter model as η′ grows. Examining equation 4.31, we see an overall factor of Λ−1, which is inversely proportional to the product TGex(S+0). From equation 4.17, we see that Gex(S+0) collapses for large η′, and this collapse itself reflects the collapse in the expression rates in equations 2.24 and 2.25. These rates, for the driftless induction case, are shown in Figure 5B, where we clearly see a narrow tipping region in Θ for fixed η′ (and therefore vice versa) in which the expression rates fall rapidly to very close to zero. As Gex(S+0) collapses, T must therefore grow in order to prevent τtrans in equation 4.31 from exceeding any given mean lifetime, such as the value 108 s used here. “Large” η′ is evidently not very large at all, and from equations 2.24, 2.25 and 4.17, it must also be Θ-dependent. Roughly speaking, this tipping point occurs when η′Θ2 ≈ 10. Above this value, expression rates rapidly fall to zero, and plasticity therefore essentially stops. The very notion of the mean transition time between segregated states becomes meaningless because the dynamics will never even reach the segregated states in this parameter regime: all synapses are frozen in their initial states as the timescales for any dynamics become irrelevantly large. The low-pass filter model must therefore operate in a regime in which the parameters η′ and Θ satisfy . Similar considerations apply to the molecular filter model: if the decay rate η′ becomes too large for a given threshold Θ in either the low-pass or molecular filter models, then reaching threshold, although theoretically inevitable, takes longer than the lifetime of the universe, let alone that of an ephemeral animal.

Figure 5:

Analytical results illustrating the lifetimes of segregated states in the low-pass filter model. (A) The step-size parameter T that gives a lifetime of 108 s for segregated states, as a function of Θ, for different choices of the decay parameter η′, as shown. The initial trends for smaller values of Θ resemble that in Figure 4 for the linear summation model. However, the required step sizes switch to diverging behavior as Θ increases. (B) Such divergence occurs because the plasticity expression rates, shown here for driftless induction processes for the same choices of η′ used in A, collapse for larger values of Θ. The expression rates as a function of Θ have been scaled by these rates at Θ = 1, so that the curves do not depend on the overall induction rate and may be plotted on the same scale.

Figure 5:

Analytical results illustrating the lifetimes of segregated states in the low-pass filter model. (A) The step-size parameter T that gives a lifetime of 108 s for segregated states, as a function of Θ, for different choices of the decay parameter η′, as shown. The initial trends for smaller values of Θ resemble that in Figure 4 for the linear summation model. However, the required step sizes switch to diverging behavior as Θ increases. (B) Such divergence occurs because the plasticity expression rates, shown here for driftless induction processes for the same choices of η′ used in A, collapse for larger values of Θ. The expression rates as a function of Θ have been scaled by these rates at Θ = 1, so that the curves do not depend on the overall induction rate and may be plotted on the same scale.

It is now instructive, after examining the analytical results for T and comparing them to those obtained in simulation, to examine the dependence of the factor 2ΔU2 in equation 4.31 on the step-size parameter T and the expression process, ignoring all aspects of the induction processes. We have
4.33
For the linear summation model, this becomes
4.34
while for the low-pass filter model, we have
4.35
We notice, first, a dependence on the inverse plasticity step size, T. This result is highly suggestive, indicating that the plasticity step size T is perhaps playing the role of a temperature parameter. It is in fact clear that T controls the size of the fluctuations on the mean dynamics, and a physical temperature indeed sets the overall scale of fluctuations. If T were playing this role in our dynamics, then we might expect a bifurcation in the mean dynamics as T is reduced below some critical level, TC. Figure 3 already hints at such a critical value. However, from equation 4.1, with T± = T, T appears merely as an overall factor, so a bifurcation at a critical value TC in the mean dynamics of the Fokker-Planck equation for our problem is not possible. Technically the derivation of the Fokker-Planck equation requires taking the diffusion limit, so it has already assumed that T is small. It is therefore likely that the approximation made in reducing the master equation to the Fokker-Planck equation destroys even the very possibility of a bifurcation in T. To resolve this question of a bifurcation in T, we will consider the full synaptic dynamics induced by the master equation in section 5. There we will find a bifurcation in the dynamics as T is reduced, confirming the interpretation of the plasticity step size as playing precisely the role of a temperature parameter.
Our second observation is that in the linear summation model, the threshold Θ scales the plasticity step size T in 2ΔU2, so that we have an effective step size Teff = T/Θ. In equation 4.31, for τtrans, there is also an overall factor of Λ−1, but this too depends only on Teff. The introduction of an expression filter therefore “cools” the dynamics by reducing the effective plasticity step size in the linear summation model. This scaling is expected to break down for larger step sizes, however, because its derivation depends on the diffusion limit and therefore the assumption that T is small. In particular, the presence of the bifurcation in the master equation formalism in section 5 must lead to only approximate scaling for larger T. In the low-pass filter model, similar observations hold, but the precise form of the scaling of T by the filter is complicated by the additional, η′-dependent factor in equation 4.35. This additional factor is monotonic increasing as a function of η′. The parameter η′ itself therefore scales the threshold Θ in a rather nonlinear fashion, but by no more than a factor of 2, in the large η′ limit. This scaling of Θ by η′ does not, however, mean that the low-pass filter model with a given threshold Θ and η′ → ∞ is equivalent to the linear summation model with a threshold 2Θ, since we know that when η′ → ∞, the expression rates go to zero and plasticity switches off. Finally, the molecular filter model is somewhat similar to the low-pass filter model in that η′ induces a scaling of Θ, but matters are more complicated because the scaling of Θ by η′ is nonmonotonic. For Θ = 3, for example, we have
4.36
with a maximum at . This suggests an optimal value of η′ for each value of Θ to achieve the strongest control of fluctuations in the molecular filter model, but ηmax typically exceeds the range of η′ over which expression rates in the molecular filter model collapse.

This effective scaling of T by the expression mechanism is, in fact, already clear in Figure 3. If instead of plotting the MSI against T, we plot it against Teff, then we would expect all the MSI curves to be superimposed, provided that T is sufficiently small. This is illustrated explicitly in Figure 6A for the linear summation model in which Teff = T/Θ. The scaling in the low-pass filter model is more complicated, giving a more elaborate form for Teff, but this is also illustrated in Figure 6B. Here we examine the additional scaling involving the parameter η′ for a fixed value of Θ. Again, we see that all of the MSI curves are superimposed when we plot them against Teff. In both cases, we see that the scaling is more or less exact. This is because the considered ranges for T are well below the bifurcation point TC that we will see in section 5.

Figure 6:

Scaling behavior in the linear summation (A) and low-pass filter (B) models. (A) The MSI in Figure 3 for the linear summation model is plotted not as a function of the step size T but as a function of the effective step size Teff = T/Θ. All the MSI curves in Figure 3 for different Θ then coincide, showing explicitly the scaling of the step size induced by the linear summation expression filter. (B) In the low-pass filter model, the scaling of T by the filter mechanism is rather more complicated, as seen in equation 4.35. Nevertheless, the MSIs obtained from simulations of the low-pass filter model with Θ = 10 and different values of η′ (η′ = 0.00, 0.01, 0.02, 0.03, 0.04) all coincide when plotted against the effective step size.

Figure 6:

Scaling behavior in the linear summation (A) and low-pass filter (B) models. (A) The MSI in Figure 3 for the linear summation model is plotted not as a function of the step size T but as a function of the effective step size Teff = T/Θ. All the MSI curves in Figure 3 for different Θ then coincide, showing explicitly the scaling of the step size induced by the linear summation expression filter. (B) In the low-pass filter model, the scaling of T by the filter mechanism is rather more complicated, as seen in equation 4.35. Nevertheless, the MSIs obtained from simulations of the low-pass filter model with Θ = 10 and different values of η′ (η′ = 0.00, 0.01, 0.02, 0.03, 0.04) all coincide when plotted against the effective step size.

We stress that the above observations hold regardless of the induction processes. In particular, they do not depend on any feature of our preferred, tristate switch model of STDP. In general, we expect that if there is a fundamental unit (here, T) that governs changes in synaptic strength, then this unit acts as a temperature, or temperature-like parameter, even if changes in strength can occur in multiples of this unit rather than in single units as considered here. This unit then acts as an irreducible source of fluctuations in synaptic strength. Any mechanism that filters induction processes before expression reduces the effective size of this fundamental unit, reducing fluctuations, and therefore effectively “cooling” the dynamics of synaptic plasticity.

5.  Synaptic Dynamics: Master Equation Approach

We have seen evidence that the plasticity step size parameter T may be playing the role of a temperature parameter, but that there is no bifurcation in the mean dynamics as a function of T. Since the derivation of the Fokker-Planck equation requires the assumption of small T, this likely precludes the existence of a bifurcation. We must return to the discrete dynamics governed by the master equation in equation 3.5. We are therefore interested in the solution of the master equation, and specifically whether the structure of this solution exhibits a qualitative change at some critical value TC as T is reduced. In particular, we require the existence of maxima of away from , the zero vector, in order for segregated states to exist. If there is a bifurcation, then we would expect a unique maximum at for T>TC, but the sudden appearance of nonzero maxima as T passes downward through TC.

It would be convenient to address this question by examining the stationary solution of the master equation. Segregated states appear over time from an initially unsegregated state, so it is reasonable to require that segregated states should correspond to maxima of the stationary solution of the master equation. Unfortunately, the definition of the postsynaptic firing rate ensures that the state is an absorbing state of the dynamics, because when λp = 0, synaptic plasticity stops, at least under STDP. With , the stationary solution of the master equation can only be . However, the critical question is the timescale on which approach to stationarity occurs, and whether on shorter timescales, a quasi-stationary solution exists that exhibits the required segregated states, but that subsequently evolves, on a much longer timescale, to the stationary solution. Consider, for example, a very small step size T. We know from simulation that segregated states exist, are attainable, and are stable against fluctuations for very long periods of time. The fluctuations required to drive such states to are therefore enormous and must develop on timescales that are irrelevant not only biologically but also cosmologically. Such a stationary solution is for all practical purposes irrelevant, but extracting the relevant, quasi-stationary solution requires separating out dynamics that develop on vastly different timescales.

A much simpler approach to this technical awkwardness is merely to prevent from being an absorbing state at all. We could achieve this by writing , where λspont>0 is some low, spontaneous postsynaptic firing rate, or even more simply by setting λp = λspont only when . We opt for the latter, since it is the minimal change, and we set λspont = 0.1 Hz. In section 5.2, we will see that the exact value of λspont is unimportant, provided it is nonzero. The above simulations are in fact run with λspont = 0.1 Hz but are of course identical to those running with λspont = 0 Hz because the timescale for the approach to stationarity in this latter case is so great that never gets anywhere near , let alone reaches it. We may regard setting λp = λspont>0 for as either a method that introduces spontaneous activity and that incidentally changes the structure of the asymptotic solution to that required for our analysis, or a method that enables us to examine the structure of the quasi-stationary solution while ignoring the biologically irrelevant stationary solution. In either case, the method does not fundamentally alter and, more important, does not introduce the bifurcation in T that we now explore. We first consider numerical results for m>1 afferents and then examine the bifurcation analytically for m = 1 afferent.

5.1.  Numerical Results for m>1 Afferents.

For one-dimensional dynamics, constructing the stationary solution of the master equation in the presence of a reflecting boundary is trivial. Since we are interested in segregated states, we must construct stationary solutions for more than one afferent, so must consider stationary solutions in m>1 dimensions with reflecting boundaries at Si = 0. We focus as usual principally on m = 2 afferents, but will also show a stationary solution for m = 3 afferents in order to demonstrate the qualitative similarity between solutions for different numbers of afferents. We construct m-dimensional stationary solutions by additionally erecting reflecting boundaries at large Si, where these boundaries are so far from the bulk of the stationary probability distribution that their presence has negligible impact on the solution. Enclosing the dynamics in a square or a cube guarantees that the stationary solutions cannot possess any net flow and therefore guarantees detailed balance. Using detailed balance, we can then reduce the m = 2 problem from an eigenmatrix problem to an eigenvector problem in the probabilities along one edge of the square. From these edge probabilities, the entire distribution can be reconstructed. Similarly, the m = 3 problem can be reduced from an eigen-rank-3-tensor problem to an eigenmatrix problem in the probabilities on one face of the cube. This problem can itself be somewhat reduced by exploiting the expecting symmetry of the stationary solution in all its variables, so the eigenmatrix problem in, say, n2 variables becomes an eigenvector problem in variables. Nevertheless, constructing m = 3 stationary solutions is very hard, computationally speaking, taking orders of magnitude longer than the construction of m = 2 stationary solutions. We check that these solutions are insensitive to the placement of the upper reflecting boundaries.

Figure 7 shows two examples of stationary probability distributions for m = 2 afferents using our new, standard set of parameters, but with the standard in Figure 7A and the nonstandard in Figure 7B. In both cases, we have used the linear summation model with a threshold Θ = 30 and induction processes based on our tristate switch model of STDP. Superimposed on the contours of constant probability are the flows of the deterministic dynamics. We have used a step size T = 1/250 to generate the stationary solutions. We see two states of high probability corresponding to the two segregated states, with a narrow, low-probability neck connecting them. This neck defines the trajectory along which fluctuations drive a change in state, and its center coincides with the separatrix of the deterministic dynamics. Figure 8 shows an example of the stationary probability distribution for m = 3 afferents, again for the linear summation model and the tristate switch model of STDP, but with Θ = 40 and T = 1/50. This figure shows surfaces of constant probability. Three high probability states are now apparent, corresponding to the three expected segregated states. On each of the three Si = 0 surfaces, we see cross sections of this three-dimensional stationary solution that are qualitatively identical to the two-dimensional solutions in Figure 7. Narrow necks connect all the pairwise combinations of fixed points. Note, in particular, the “hole” in the center of the distribution. This hole reflects the fact that fluctuations that drive a change in state via the growth of any pair of zero-strength afferents are suppressed compared to those involving just one zero-strength afferent, as discussed earlier: fluctuations involving all three afferents arise with much lower probability than fluctuations involving pairs of afferents, one of which is strong.

Figure 7:

Contours of the stationary solutions of the master equations for two afferents of strengths S1 and S2 in the tristate switch model of STDP using the linear summation filter with threshold Θ = 30. Parameters are as in Figure 3, except that we use the standard value in A and the nonstandard value in B. The plasticity step size is T = 1/250 in both cases. In A, the contours correspond, moving out from the maxima on the axes, to probabilities of 10−2, 10−3, 10−4, 10−5, 10−6, 10−7, 10−8, 10−9, and 10−10, while in B, they correspond to 10−5, 10−8, 10−11, 10−14, 10−17, 10−20, 10−23, 10−26, 10−29, and 10−32. In B the two contours for a probability of 10−2 are not shown because they are very close to the axes and indistinguishable from them. The distribution in B is very tightly segregated, reflected in the rapid drop in probability away from the maxima and the extremely low-probability “neck” connecting the two maxima. Also shown are the flows for the deterministic dynamics, indicating that the center of the trajectory along which fluctuations drive a change in state coincides with the separatrix of the deterministic dynamics.

Figure 7:

Contours of the stationary solutions of the master equations for two afferents of strengths S1 and S2 in the tristate switch model of STDP using the linear summation filter with threshold Θ = 30. Parameters are as in Figure 3, except that we use the standard value in A and the nonstandard value in B. The plasticity step size is T = 1/250 in both cases. In A, the contours correspond, moving out from the maxima on the axes, to probabilities of 10−2, 10−3, 10−4, 10−5, 10−6, 10−7, 10−8, 10−9, and 10−10, while in B, they correspond to 10−5, 10−8, 10−11, 10−14, 10−17, 10−20, 10−23, 10−26, 10−29, and 10−32. In B the two contours for a probability of 10−2 are not shown because they are very close to the axes and indistinguishable from them. The distribution in B is very tightly segregated, reflected in the rapid drop in probability away from the maxima and the extremely low-probability “neck” connecting the two maxima. Also shown are the flows for the deterministic dynamics, indicating that the center of the trajectory along which fluctuations drive a change in state coincides with the separatrix of the deterministic dynamics.

Figure 8:

Surfaces of the stationary solution of the master equations for three afferents of strengths S1, S2, and S3 in the tristate switch model of STDP using the linear summation filter with threshold Θ = 40. Parameters are as in Figure 3, and the plasticity step size is T = 1/50. In addition, the third-order correlation coefficient is set to zero. Surfaces of constant probability are displayed, corresponding to probabilities of 4 × 10−3 (the inner cores on each of the three axes), 2 × 10−3 (the middle layers surrounding the three inner cores), and 10−3 (the connected region with the central hole).

Figure 8:

Surfaces of the stationary solution of the master equations for three afferents of strengths S1, S2, and S3 in the tristate switch model of STDP using the linear summation filter with threshold Θ = 40. Parameters are as in Figure 3, and the plasticity step size is T = 1/50. In addition, the third-order correlation coefficient is set to zero. Surfaces of constant probability are displayed, corresponding to probabilities of 4 × 10−3 (the inner cores on each of the three axes), 2 × 10−3 (the middle layers surrounding the three inner cores), and 10−3 (the connected region with the central hole).

The segregated states correspond to maxima of the stationary probability distribution provided that the plasticity step size T is small enough. In order to establish whether there is a bifurcation in the structure of the stationary solution at a critical value of the step size T, it suffices to examine, for example, as a function of S, for different values of T. Results are shown in Figure 9 for the unfiltered tristate switch model of STDP, which is just the linear summation model with Θ = 1. For large T (T = 1/10 in this figure), despite the reflecting boundaries at Si = 0 and the nonzero spontaneous rate at , the stationary probability distribution is nevertheless concentrated around S = 0 and monotonically decreases as S increases. As T decreases (to T = 1/40), the distribution spreads out somewhat and we see indications of a developing, incipient maximum away from S = 0, but there is still only one maximum, at S = 0, and the distribution continues to be monotonically decreasing. For still smaller T (T = 1/80), however, a maximum at nonzero S appears, with the maximum at S = 0 still present and with comparable probability. Finally, for an even smaller T (T = 1/160 in Figure 9), the maximum at S = 0 is small and the maximum for some nonzero S dominates. Thus, the stationary probability distribution for m = 2 afferents exhibits a bifurcation as T is reduced, with the bifurcation in Figure 9 occurring at TC ≈ 1/51. For TTC, the synaptic strength is with high probability, but for TTC, is concentrated around the states with with high probability. In a regime in which TTC or T is only somewhat lower than TC, the dynamics of synaptic strength will be dominated by fluctuations rapidly driving the strength state between both segregated states and . A bifurcation is present for m = 3 afferents too, and we generally expect such bifurcations for any number of afferents. The plasticity step-size parameter T therefore plays the role of a temperature-like parameter, and at the critical point, symmetry breaking between the m afferents occurs, with segregated states emerging as maxima of the stationary probability distribution.

Figure 9:

The stationary probability distribution for two afferents along the cross-section on which the second afferent's strength is zero. Distributions are shown for four different step sizes T as indicated and are each normalized by their maximum value so that they all appear on a similar scale. For larger values of T, has only one maximum, at S = 0, but as T is reduced, a bifurcation occurs, and a second maximum develops corresponding to the emergence of segregated states.

Figure 9:

The stationary probability distribution for two afferents along the cross-section on which the second afferent's strength is zero. Distributions are shown for four different step sizes T as indicated and are each normalized by their maximum value so that they all appear on a similar scale. For larger values of T, has only one maximum, at S = 0, but as T is reduced, a bifurcation occurs, and a second maximum develops corresponding to the emergence of segregated states.

We have seen in section 4 that the presence of an expression filter scales the plasticity step-size parameter T into an effective step-size parameter Teff governed by the threshold Θ and the decay parameter η′. We would therefore expect the location of the bifurcation TC to depend on these filter parameters. For the linear summation model we show, in Figure 10A, the location of the nonzero maximum S of the stationary probability distribution as a function of the step size T for different values of the threshold Θ. The maxima are expressed in units of step size: how many steps from zero are needed to reach a maximum. The sudden jump of the location of a maximum from zero to some nonzero value gives the critical step size TC for each value of Θ. Also shown in this graph is the location of the segregated fixed points of the deterministic dynamics, which is independent of both T and Θ. In units of T, therefore scales inversely with T. For larger step sizes, there is a considerable deviation between the maxima of the stationary master equation and the location of the fixed points of the deterministic dynamics, as we might expect since the deterministic dynamics are derived under the assumption that the step size is small. For smaller step sizes, however, these locations must converge, as we see. In Figure 10B, we plot the critical step size TC as a function of the threshold Θ. We see that for smaller values of Θ, TC scales linearly with Θ, so that the critical step sizes themselves also exhibit scaling for small Θ. However, as Θ increases, TC switches over to asymptotic behavior, giving a maximum possible step size for any linear summation filter dynamics consistent with a nonzero solution of the stationary master equation. Scaling in TC must therefore break down as Θ increases.

Figure 10:

The locations of the critical step sizes TC in the linear summation model for different thresholds Θ. (A) In units of the step size T, the location of the nonzero maximum S of is plotted as a function of T. The sudden onset of a nonzero solution indicates a bifurcation, defining a critical step size, TC. The location of the segregated fixed points of the deterministic dynamics is also shown (line marked “FP”). (B) The dependence of the critical step size TC on the threshold Θ. For smaller Θ, TC scales linearly with Θ, but switches to asymptotic behavior for larger values of Θ.

Figure 10:

The locations of the critical step sizes TC in the linear summation model for different thresholds Θ. (A) In units of the step size T, the location of the nonzero maximum S of is plotted as a function of T. The sudden onset of a nonzero solution indicates a bifurcation, defining a critical step size, TC. The location of the segregated fixed points of the deterministic dynamics is also shown (line marked “FP”). (B) The dependence of the critical step size TC on the threshold Θ. For smaller Θ, TC scales linearly with Θ, but switches to asymptotic behavior for larger values of Θ.

Finally, in Figure 11, we plot on a log-log graph the effective critical step size TeffC against the threshold Θ for both the linear summation model and the low-pass filter model for different values of the decay parameter η′. The effective critical step size is just the critical step size TC divided by the scaling factor induced by the expression process. For the linear summation model this is simply Θ, while for the low-pass filter model, it is given in Equation 4.35. We clearly see the initial invariance of TeffC for small Θ, regardless of the filter model and the value of the decay parameter η′, reflecting the nearly perfect scaling of the critical point when this point is sufficiently below the maximum possible step size consistent with a nonzero maximum of the stationary solution of the master equation. The deviation of the behavior of the low-pass filter model from the linear summation model reflects the additional scaling of the threshold Θ by the decay parameter η′ discussed earlier. The curves for the low-pass filter model for different values of η′ converge as Θ increases and are displaced from that for the linear summation by a factor of 2 in the threshold Θ, reflecting the factor of 2 in Equation 4.35 for large η′ but also large Θ. Similar results to those for the low-pass filter model are obtained for the molecular filter model, so we do not show them.

Figure 11:

The dependence of the effective critical step size, TeffC, on the threshold Θ for the linear summation filter (line marked η′ = 0.00) and for the low-pass filter for different choices of η′ as indicated. The effective critical step size is the critical step size divided by the scaling induced by the filter process. Results are shown on a log-log plot. The initial, almost perfect scaling for both filter models is clear.

Figure 11:

The dependence of the effective critical step size, TeffC, on the threshold Θ for the linear summation filter (line marked η′ = 0.00) and for the low-pass filter for different choices of η′ as indicated. The effective critical step size is the critical step size divided by the scaling induced by the filter process. Results are shown on a log-log plot. The initial, almost perfect scaling for both filter models is clear.

5.2.  Analytical Results for m = 1 Afferent.

So far we have approached the bifurcation by constructing stationary solutions of the master equation in two dimensions. This construction involves solving an eigenvector problem, so it is likely very hard to address analytically. In one dimension, we can write down an explicit solution for the stationary probability distribution that could be amenable to analysis, but of course with only one afferent, competitive dynamics cannot arise. Nevertheless, the key feature of the bifurcation for m>1 afferents is the sudden emergence of probability maxima away from for T < TC. These maxima could correspond, at O(q0), to the fixed hyperplane of nonzero states at (for which competition is absent because the afferents are indistinguishable) or, at higher orders in q, to the m nonzero states (for which competition drives development to only one of these segregated states). It is therefore sufficient to determine whether the stationary probability distribution for m = 1 afferent also exhibits the sudden emergence of a maximum of Pstat(S) away from S = 0. Since the expression process merely modifies an existing bifurcation in the induction process, for simplicity we restrict to a consideration of just the unfiltered tristate switch model of STDP. Further, we consider activity patterns only with q = 0, so that we consider only O(q0) results. This is then equivalent to determining when a fixed hyperplane of states S = S+0>0 appears.

In the presence of the reflecting boundary at S = 0, the stationary probability at S = nT, Pstat(nT), can be written down directly in terms of that at S = 0, Pstat(0), as
5.1
where, for n = 1, the product term is understood to be absent. Notice that if λspont = 0, then ϱ+in(μ, 0) ≡ 0, so that the distribution would collapse to Pstat(nT) = δn,0. From equations 2.5 and 2.6, we may write
5.2
The product of the terms in the numerator is immediate, but we must evaluate the product of the terms in the denominator. We do this by taking the logarithm and making the approximation that the resulting sum can be replaced by an integral:
5.3
We then obtain
5.4
The essential dynamics are in the n-dependent part of the right-hand side of this equation, so we consider only this part and throw away all constant factors. We now also write S = nT, or n = S/T, and take the continuum limit by regarding S as a continuous rather than discrete variable. We must therefore consider whether the function
5.5
for S>0 exhibits a T-dependent maximum in S. Notice that the location of any maximum in this approximation is independent of λspont. Although λspont>0 is vital to prevent the stationary solution from collapsing to a delta function, the location of the bifurcation is insensitive to its exact value.
Differentiating f(S) with respect to S to determine the location of any extrema, they are located at the zeros of the function
5.6
For the particular case in which we take the limit T → 0, g(S) reduces to a form whose zeros can be read off by inspection. From the logarithmic factor, we obtain in this case a zero at S = (λ + μ − λ+)/μ ≡ S+0. This is precisely the location of the fixed point to O(q0) in equation 4.14 obtained from the deterministic dynamics in equation 4.5. Thus, although the form of f(S) was derived by making an approximation, it gives the correct maximum in the limit T → 0.
For T>0, any zeros of g(S) are the solutions of a transcendental equation, so that we cannot express them in closed form. We see that if T is large, then the first term on the right-hand side of equation 5.6 dominates and that there are therefore no zeros for S>0. Hence, g(S) certainly exhibits a T-dependent bifurcation in its positive zeros. We expect the critical value TC to be small, and we also expect S at any maximum to be less than unity. So we make a further approximation by expanding g(S) to second order in both S and T, so that the roots of g(S) are obtained, approximately, from the roots of a quadratic equation. The discriminant of this quadratic equation then determines whether real roots exist, and therefore where the bifurcation occurs. The discriminant is a quartic equation in T, but since we expect TC to be small, we simply linearize the discriminant in T and obtain the critical value of the plasticity step size in the unfiltered, tristate switch model of STDP as
5.7
where we have written Ω = log[λ+/(λ + μ)]. Notice that TC = 0 when λ + μ = λ+, or when S+0 = 0. Thus, despite the several approximations necessary to derive equation 5.7, the critical value of TC vanishes at the point at which the hyperplane of fixed points passes through the origin. The approximations therefore produce results consistent with our expectations about the possible locations of any maxima. Inserting our standard values μ = 50 Hz, λ = 50 Hz and λ+ = λ/0.6 into equation 5.7, we obtain TC ≈ 0.02 = 1/50. This analytical result agrees astonishingly well, given the approximations, with the numerical result for m = 2 afferents read off from Figure 10 for Θ = 1, namely, TC ≈ 1/51.

The fact that we are able to obtain the location of the critical plasticity step size in the unfiltered, tristate switch model of STDP depends fundamentally on the simplicity of the ratio of upward and downward transition rates in equation 5.2. In particular, after taking the logarithm, we can explicitly evaluate the integral in equation 5.3 to obtain a formula for f(S). With a more complicated induction process, or in the presence of an expression filter, so that the induction rates ϱ±in are replaced by the much more complicated expression rates ϱ±ex, we would not in general be able to evaluate the integral and obtain f(S) in closed form. Nevertheless, the scaling of TC observed above can be used to deduce, in our tristate switch STDP model, other critical points in the presence of a filter, at least before the transition to asymptotic behavior occurs. The calculation of TC for the unfiltered induction process is therefore the basic one that underlies all other calculations.

In order to obtain explicit numerical and analytical results in this section, we have had to employ our tristate switch model of STDP for the induction process. However, the hint of the temperature-like role for T owes its origin to the observations concerning the dependence of the quantity 2ΔU2 appearing in the formula for the mean transition time between segregated states in equation 4.31. This hint was then confirmed by the observation of bifurcation behavior, both numerically and analytically, in the stationary solutions of the master equation, for which it was necessary to use an explicit model of the induction process. We would not therefore expect the existence of this bifurcation to depend on the specific induction process considered here, but rather to be a generic feature of models of synaptic plasticity with a discrete plasticity step size, provided that these dynamics are competitive in nature, or have other types of fixed points away from the origin. Furthermore, the precise form of the scaling relations that we derived in equations 4.34 and 4.35 and equivalent but more messy relations for the molecular filter model are completely general, not requiring any knowledge of the induction process.

6.  Discussion

In previous work on our tristate switch model of STDP (Appleby & Elliott, 2005), we have distinguished between plasticity induction and plasticity expression in order to introduce a synaptic filter between them that can control fluctuations in synaptic strength (Elliott, 2008, 2010b, 2011; Elliott & Lagogiannis, 2009). In the presence of a synaptic filter, we found that we could increase the size of the fixed-amplitude step changes in synaptic strength without compromising the long-term stability of states of synaptic strength against fluctuations (Elliott & Lagogiannis, 2009; Elliott, 2011). In this earlier work, we were largely content with a simulation-based approach in order to determine the dependence of the MSI, a measure of stability, on the plasticity step size. The realization that a non-Markovian random walk underlies our tristate switch model of STDP (Elliott, 2010b), however, raised the possibility of an analytical attack on this problem of stability against fluctuations by deploying well-established and powerful techniques from the theory of stochastic processes but adapted to our non-Markovian setting. Moreover, such an approach permits a generality of analysis that has hitherto been impossible, liberating us from the details of our particular (although still favored) tristate switch model of STDP and our earlier, painfully laborious analytical methods, and allowing us to consider general, non-Markovian random walk processes of fixed-amplitude changes in synaptic strength.

Several earlier theoretical studies have used a Fokker-Planck or master equation formalism to move beyond an examination of the mean dynamics of synaptic plasticity and consider fluctuations, although in the context of conventional models of STDP (see, e.g., Kempter, Gerstner, & van Hemmen, 1999; van Rossum et al., 2000; Rubin et al., 2001; Burkitt et al., 2004). Kempter et al. (1999) performed an analysis of an additive model of STDP with inhomogeneous Poisson spiking processes, carefully separating out processes occurring on different timescales. In considering variability in synaptic strength, they restricted to constant pre- and postsynaptic rates and statistically independent pre- and postsynaptic spike trains. The source of variability in their analysis is variability in the spike trains themselves, not other sources, such as those arising in an intrinsically stochastic model of synaptic plasticity, as here. The variance in synaptic strengths in their model grows linearly in time, like a standard diffusion process, and they calculate the ratio of the timescale for strength diffusion to become significant to the timescale for the emergence of structure formation. They argue that this ratio can always be made greater than unity by scaling strength changes, or the learning rate, by a constant, overall factor. Nevertheless, although they can guarantee that structure formation will always occur, it is inevitable that structure will always eventually be washed out by this strength diffusion process. In our approach, we have explicitly rejected modifying the overall scale of strength changes, considering this to be biologically implausible (Elliott, 2008). Moreover, although we have shown in our model that the changes in strength asymptote to a Wiener process with drift when pre- and postsynaptic firing rates are fixed (Elliott, 2010b), we clearly see from the stationary solutions of the master equation (see Figures 7 and 8) that when postsynaptic firing is governed by synaptic transmission, segregated states are never washed out by fluctuations. Rather, fluctuations simply drive the system between different segregated states, and our aim here has been to understand, analytically, the timescale on which this process occurs and to determine how a synaptic filter increases this timescale without our having to reduce the plasticity step size unacceptably, or by adjusting a somewhat artificial, ad hoc “learning rate” parameter.

Typically it has proved rather hard to calculate the diffusion term in the Fokker-Planck equation in many models of STDP (see, e.g., Rubin et al., 2001, and especially Burkitt et al., 2004). The difficulty arises because although synapses on the same neuron may experience different presynaptic spike trains, they necessarily always experience the same postsynaptic spike train. Because STDP is driven by spikes, synaptic strength changes over multiple synapses are therefore not independent processes. Why is it, then, that we have apparently overcome this problem and managed to derive a non-Markovian master equation for the strength of a single synapse (Elliott, 2010b) and, by making a Markovian approximation to enforce renewal, extended it here to include the strengths of multiple synapses in Equation 3.5? The reason is that the structure of our stochastic, tristate switch model of STDP is fundamentally different from more conventional models of STDP. As argued above and demonstrated elsewhere (Elliott & Lagogiannis, 2009), intrinsic fluctuations, due to the stochastic nature of the synaptic switching processes that govern whether fixed amplitude changes in synaptic strength are induced, dominate stimulus-induced fluctuations, at least for the stationary Poisson processes that we have usually considered. Our model requires each individual synapse to possess a tristate plasticity switch and, critically, these switches in different synapses are assumed to be independent. For different synapses experiencing different presynaptic spike trains, the intrinsic fluctuations will tend to decorrelate the changes in synaptic strength due to the common postsynaptic spike train. Furthermore, the extent of this decorrelation improves in the presence of a synaptic filter. With a filter, most spikes, whether pre- or postsynaptic, do not lead to expressed changes in synaptic strength, but possibly only to changes in filter states. In the presence of a filter, further variability in synaptic plasticity is introduced by an amplification of the stochasticity in the induction signals or by introducing yet further sources of intrinsic stochasticity, such as in the molecular filter model with its stochastic decay. Thus, even a pair of synapses experiencing identical pre- and postsynaptic spike trains will typically not undergo correlated changes in synaptic strength because the internal degrees of freedom in both the independent induction mechanisms and the independent expression mechanisms at both synapses will guarantee that the likelihood of simultaneous strength changes at both synapses will be very small. Thus, although correlated strength changes are inevitable in conventional STDP models in which, given a defined pre- and postsynaptic spike train, strength changes are entirely deterministic, in our model of STDP in the presence of a synaptic filter, correlated strength changes are rare. The structure of our stochastic model of STDP and its extension to include synaptic filtering ensures that we may to a very good approximation take the changes in synaptic strength across multiple synapses experiencing the same postsynaptic spike train to be independent. Even without a synaptic filter, for Θ = 1, the agreement between our analytical and numerical results indicates that the degree of decorrelation introduced by the intrinsic fluctuations driven by the stochastic tristate switch model is sufficient to be able to ignore any surviving correlations. Stochastic induction models coupled with synaptic filtering ensure that synapses express changes in synaptic strength only rather slowly, and stochastically, rather than rapidly and deterministically, on the millisecond timescales tacitly assumed by conventional models of STDP. This slowing down and randomizing of plasticity expression events ensures that the rate-based limit is achieved at single synapses even for single plasticity steps and eliminates spike-induced correlations in strength changes across different synapses.

Our principal aim in this letter has been to begin the program of work articulated elsewhere (Elliott, 2010b) and use a Fokker-Planck formalism to compute the mean transition time between meta-stable, segregated states of synaptic strength in stochastic models of synaptic plasticity induction in the presence of a filtering expression mechanism. We have shown in our model that a change in segregation state driven by fluctuations can be cast in terms of a potential barrier problem. We have derived the potential, in equation 4.23, for general induction and expression mechanisms, described by the functions in equations 4.12 and 4.13. This led to a computation of the mean transition time, in equation 4.31, that differs somewhat from the standard, Arrhenius-like form because we do not have to make the usual quadratic approximations to the potential and our problem is bounded by reflecting boundaries. Two approximations and one simplifying assumption were necessary to derive these results. We have approximated the non-Markovian expression processes by Markov processes whose rates are the asymptotic rates of the non-Markovian processes, but this approximation is not a bad one. We also made an adiabatic approximation in averaging the rates of the induction processes over pre- and postsynaptic activity patterns. Such an approximation cannot be avoided in the absence of a detailed model of pre- and postsynaptic firing. Finally, we simplified our analysis by considering only small fluctuations of presynaptic activity patterns around their means, so that the separatrix and the trajectory along which fluctuations drive a change in segregation state could be taken to be a straight line. This simplifying assumption is likely not indispensable since we could in principle derive the equation of the separatrix, but it is an extremely convenient one. Translating the derived mean transition time into an MSI, we found that our analysis captures in detail the qualitative features of our simulation results and that, moreover, our analytical results consistently and uniformly differ from simulation results by only 20%, despite these two approximations.

We have attempted to maintain as much generality as possible by not specifying in detail either the induction process or the expression process unless concrete models are required. Where necessary, we have considered three particular models for the expression of synaptic plasticity. These three expression models have been examined in detail elsewhere (Elliott, 2011), where we focused particularly on the parameter regimes in which plasticity expression times are within reasonable, experimental bounds when the relative drift Δρ in plasticity induction signals is not close to zero, but are large when the relative drift is close to zero. This latter requirement stood as a surrogate for a large mean transition time between segregated states of afferent connectivity in the absence of an explicit calculation of the mean lifetimes of meta-stable states. The first and simplest of these models of plasticity expression that we have considered is the linear summation model. In this model, induction signals are effectively added to a running counter, and expression occurs only when its modulus first exceeds a given threshold. The model controls fluctuations very effectively and is particularly easy to analyze, so it provides considerable insight into the manner in which an expression mechanism can tame fluctuations. For example, the scaling of the plasticity step size T by the linear summation expression model involves a simple factor of Θ. A natural generalization of the linear summation model is to include a deterministic, continuous decay that returns the running counter to zero. Such a model constitutes a low-pass filter, and the decay turns the expression model into an Ornstein-Uhlenbeck process, although in general in the presence of nongaussian and non-Markovian signals. The continuous nature of the induction filter in the low-pass filter model is, however, problematic from a biophysical perspective, so we have also considered a discrete, molecular filter in which the running counter decays to zero in stochastic, unit-sized steps. Both filters with decay are, in fact, broadly similar in their overall behavior, although the molecular filter model can be analyzed exactly while the low-pass filter model can be analyzed only in the diffusion limit. Both decay models permit the plasticity step size T to be increased somewhat more than in the linear summation model, because of the enhanced stability of these plasticity models. Models with decay ensure that periods of stimulus quiescence are not simply periods of synaptic inactivity, because the filter will return to zero in the absence of induction signals. We also find that synaptic plasticity expression can, for all practical purposes, be completely switched off if the decay rate exceeds some relatively sharp threshold level. In this case, the rate of plasticity expression collapses exponentially fast to very close to zero. Expression models with filter decay therefore provide a very simple and natural mechanism for controlling the level of expressed plasticity.

Although our principal aim has been the computation of the mean transition time between segregated states of afferent connectivity, the precise form of the dependence of this time on both the plasticity step size and the expression mechanism leads to two extremely tantalizing observations. The first is that the expression mechanism induces an effective scaling of the actual plasticity step size T into a smaller, effective plasticity step size Teff. In the linear summation expression model, for example, the effective step size is just T/Θ, so the expression threshold scales downward the plasticity step size. Plotting MSIs from simulation against this effective plasticity step size instead of the actual plasticity step size demonstrates this scaling explicitly: all curves become superimposed. Such scaling does not imply, however, that a filtered induction model with step size T is completely equivalent to and can be replaced by a filterless induction process with (actual, not effective) step size Teff. The impact of a filter is to suppress the influence of the slower of the potentiating or depressing induction processes, but this suppression is highly nonlinear, so that for driftless induction processes, both are suppressed (Elliott, 2011). Replacing a process of rate r, say, and step size T by a process of rate rT/Teff and step size Teff does not capture such nonlinearity. This scaling is derived from the mean transition time and therefore reflects only the impact of filtering on the control of fluctuations. In terms of the control of fluctuations, a filtered induction process behaves like a filterless induction process of step size Teff, but more generally these two processes behave differently. The scaling of T into Teff, furthermore, can be only approximate, valid only when T is small, since the scaling relations are ultimately derived from the Fokker-Planck equation, which requires the diffusion limit.

The second intriguing observation is the possibility that the plasticity step size T behaves as a temperature or temperature-like parameter in stochastic models of synaptic plasticity. We might then expect the existence of a critical step size TC at which the dynamics undergo a qualitative change. The quite sharp transitions in the MSIs as a function of T hint at a possible critical value of T, and the fact that T controls the size of fluctuations is at least consistent with this possible interpretation. However, in the mean dynamics from the Fokker-Planck equation, we see no sign of a bifurcation because T appears just as an overall scale. In order to address this question, we were forced to consider the master equation, from which the Fokker-Planck equation is extracted as a small T limit. Doing so, we indeed found, both numerically and analytically, a bifurcation in the structure of the stationary probability distribution at a critical plasticity step size TC, and that this critical value depends on the filtering mechanism and exhibits, when small, the same scaling discussed above.

Thus, we have a very elegant interpretation. The plasticity step size acts precisely like a temperature parameter for induction processes, controlling the size of fluctuations on the mean dynamics from the master equation. When T>TC, the only maximum of the stationary solution of the master equation is at , with all afferents having zero synaptic strength. As T passes downward through TC, a bifurcation occurs, and the stationary solution acquires maxima away from . At TC, spontaneous symmetry breaking occurs. The symmetry between all m afferents is broken, and the system moves toward segregated states in which one afferent is singled out. For , fluctuations dominate, but for TTC, the mean dynamics dominate. The scaling of T into Teff by the filter mechanism acts by cooling the dynamics, reducing the impact of fluctuations, and stabilizing patterns of synaptic connectivity. It is likely that the dynamics undergo a phase transition at the critical step size TC, but we have yet to characterize the process in detail. Certainly this process is strongly reminiscent of ferromagnetism and the Ising model (Huang, 1987). In much earlier, unrelated work, we also showed a phase transition in a model of ocular dominance column formation (Elliott, Howarth, & Shadbolt, 1996), but in that model, the “temperature” was merely a parameter imposed by hand, identical in spirit to that found in a Boltzmann machine (Hinton & Sejnowski, 1983), which is a Hopfield network (Hopfield, 1982) with probabilistic updates. Here, we stress that the plasticity step size T is an intrinsic parameter of stochastic, random-walk models of synaptic plasticity and that we have uncovered, rather than imposed, its role as a temperature or temperature-like parameter in such models.

Because we have attempted to maintain as much generality as possible, most of these conclusions are independent of, or likely independent of, our tristate switch model of STDP. The scaling of T in the control of fluctuations depends on only general induction and expression processes, although it does require the existence of segregated fixed points, since the scaling comes from the mean transition time between them. While it was necessary to employ a concrete model of plasticity induction in order to construct stationary solutions of the master equation, the existence of the bifurcation in these solutions likely does not depend too strongly on the details of the induction model, provided that it can give rise to probability maxima away from the origin. Indeed, it probably does not matter, either, that these maxima correspond to the segregated states of neuronal connectivity that we expect to emerge during neuronal development (Purves & Lichtman, 1985). We have considered such states because of our long-standing interest in developmental synaptic plasticity, but they are also extremely convenient as tools with which to explore issues of long-term stability against fluctuations. The key feature is the existence of nonzero fixed points in synaptic strength, whatever their interpretation, in stochastic models of synaptic plasticity. These fixed points could correspond to developmentally relevant states or states relevant to learning and memory more generally.

If the ubiquity of the bifurcation is as great as we suspect, then the bifurcation coupled with the scaling induced by a synaptic expression filter could have important consequences for the acquisition, stability, and erasing of functionally relevant states of neuronal connectivity. We could consider, for example, a neuron with an intrinsic plasticity step size that is close to the bifurcation point in the absence of a filter. By dynamically adjusting its expression filter, using either the threshold Θ or the decay parameter η′ (or other parameters in other filtering models), the neuron could regulate its own effective plasticity step size, making it less, or more, sensitive to the influence of fluctuations and its patterns of connectivity more, or less, stable. New states of connectivity could be acquired during a less stable phase, and after acquisition, the filter parameters changed rapidly to freeze in structure. Indeed, since each synapse must possess such a filter, these processes could occur in a highly local, synapse-specific manner, freezing some while boiling others. Given the nature of our proposed filtering mechanism, it is easy to conceive that the filter parameters are under modulatory control by other inputs. An alternative scenario is to imagine neurons sprouting and retracting synapses, perhaps keeping overall synaptic strength approximately constant (cf. Peng et al., 2009), (Peng et al., 2009) so that a neuron with more synapses has individually weaker synapses (smaller plasticity step size), while a neuron with fewer synapses has individually stronger synapses (larger plasticity step size). In this case, even without a filter mechanism, the neuron, by regulating the number of its synapses, automatically regulates the stability of its connectivity and its capacity to seek out new patterns of connectivity. Given the growing realization of the importance of ongoing structural plasticity in the brain (for a review, see, Chklovskii, Mel, & Svoboda, 2004; see also Elliott et al., 1996 for a model of anatomical plasticity), (Chklovskii et al., 2004) (Elliott et al., 1996) a network of neurons could “self-anneal” into a functionally relevant state of connectivity and subsequently boil it out and freeze in a new state. Both dynamic filtering and dynamic synaptogenesis could occur in the same neuron, over a large population of neurons with their various filters and synapses in different states, generating extremely rich, diverse possibilities for neuronal functioning. It will be fascinating to explore these ideas in future work.

Appendix:  Timescale for the Approach of Expression Rates to Asymptotic Behavior

We must examine the approach of the non-Markovian expression processes to their asymptotic behaviors in order to assess the approximation of replacing the non-Markovian processes by Markov processes whose rates are the asymptotic rates of the non-Markovian processes. Because we are considering dynamics between expression steps, the postsynaptic firing rate will remain constant provided that the presynaptic firing rates remain constant. For induction processes based on our tristate switch model of STDP, we already know that r±in(t) usually attain their asymptotic limits on timescales shorter than , and in general typically well under 1 second for all but very small pre- and postsynaptic firing rates (Elliott, 2010b). Consider, therefore, Markov induction processes of (constant) rates at any given synapse. Then . With these forms for , we can use equations 2.7 and 2.8 for the linear summation model to determine how quickly the rates of expression steps, r±ex(t), attain their asymptotic limits. Since , we have that
A.1
A.2
The denominators are in fact only a degree Θ and not degree 2Θ polynomial in p since Φ+(p(p) = ρ+ and Φ+(p) + Φ(p) = (p + ρ+ + ρ)/ρ+. The Θ roots of this polynomial can be written as the solutions of Φ+(p) = exp(2πil/Θ), for l = 0, …Θ − 1, from which we obtain the roots as
A.3
We can see that p0 = 0 and Re{pl} < 0 for l>0, so that nonzero asymptotic rates do indeed exist. The rate of approach to these asymptotic limits is governed by the nonzero roots (or root, for Θ = 2) with largest (i.e., least negative) real parts. Since Re{pl} = (ρ+ + ρ)[cos(2πl/Θ) − 1], the relevant roots are p1 and pΘ−1, for Θ>2. Notice that the relative drift in induction rates Δρ affects only the imaginary parts, Im{pl} = (ρ+ + ρρsin(2πl/Θ), so it does not affect the rate of approach to the asymptotic limits. It may seem odd that there is in general an oscillatory component to the rates r±ex(t), but if we set, for example, ρ = 0, then the (now undefective) PDF g+ex(t) is simply that for a gamma process of order Θ with parameter ρ+ and asymptotic rate ρ+/Θ. It is easy to check that even a gamma process has an oscillatory component to the time evolution of its rate, so we should also generally expect such components for expression processes with ρ± ≠ 0. Writing τ−1 = (ρ+ + ρ)[cos(2π/Θ) − 1] ≈ 2π2+ + ρ)/Θ2 for large enough Θ, the time τ ≈ Θ2/[2π2+ + ρ)] sets the scale for the approach to asymptotic expression behavior.

A total average induction rate ρ+ + ρ ≈ 20 Hz is typical for our tristate switch model of STDP, the average being taken over λπ, λp ∈ [0, 100] Hz (Elliott, 2011), so τ ≈ Θ2/400 s. For Θ = 20, τ ≈ 1 s, and for Θ = 40, τ ≈ 4 s. These are timescales characterizing the approach of the expression rates to their asymptotic values. They must be compared to the mean time for expression in equation 2.11. Expression times depend not only on the total induction rate but also on the relative drift in the induction rates, Δρ. With λπ, λp ∈ [0, 100] Hz and λ± = 50 Hz as typical inverse timescales for STDP, the relative induction drift in our STDP model never exceeds 0.5 in modulus but is typically smaller. With this maximal value Δρ = 0.5 and the average value ρ+ + ρ = 20 Hz, the mean expression times are s for Θ = 20 and s for Θ = 40. These are worst-case scenarios for Δρ giving the smallest mean expression times, and in these cases, τ does not exceed . The average (absolute) value of the relative drift is Δρ ≈ 0.17, although this is actually skewed upward by the unbiological regimes in which λπ ≈ 0 Hz and λp ≈ 100 Hz, and vice versa, in the averaging over λπ, λp ∈ [0, 100] Hz. In this case, the mean expression times are s for Θ = 20 and s for Θ = 40, both of which are adequately larger than their corresponding values of τ. For |Δρ| and Θ large enough, from equation 2.11, so we certainly can engineer regimes in which , which requires Δρ>2π2/Θ. Only for does the required value of |Δρ| drop below 0.5, so we can avoid such regimes by not taking Θ too large. In any event, considerations of biomolecular instantiation of the filter X at single synapses also preclude taking Θ too large (Elliott, 2011), so the regimes in which r±ex(t) do not attain their asymptotic limits before the expected time to expression can be ignored. Although these are general estimates for the relaxation of r±ex(t) to their asymptotic limits, these arguments indicate that the approximation of replacing the non-Markovian expression processes by Markov processes with rates is not a bad one, at least for the linear summation model. For the low-pass and molecular filter models, determining the rate of approach of r±ex(t) to asymptotic behavior is considerably harder in full generality, but we expect the same broad conclusion to hold.

Notes

1

To some degree, the existence of meta-stable states is an epistemological question rather than an ontological one.

References

Abarbanel
,
H. D. I.
,
Talathi
,
S. S.
,
Gibb
,
L.
, &
Rabinovich
,
M. I.
(
2005
).
Synaptic plasticity with discrete state synapses
.
Phys. Rev. E
,
72
,
031914
.
Appleby
,
P. A.
, &
Elliott
,
T.
(
2005
).
Synaptic and temporal ensemble interpretation of spike-timing-dependent plasticity
.
Neural Comput.
,
17
,
2316
2336
.
Appleby
,
P. A.
, &
Elliott
,
T.
(
2006
).
Stable competitive dynamics emerge from multispike interactions in a stochastic model of spike-timing-dependent plasticity
.
Neural Comput.
,
18
,
2414
2464
.
Bagal
,
A. A.
,
Kao
,
J. P. Y.
,
Tang
,
C.-M.
, &
Thompson
,
S. M.
(
2005
).
Long-term potentiation of exogenous glutamate responses at single dendritic spines
.
,
102
,
14434
14439
.
Baldassi
,
C.
,
Braunstein
,
A.
,
Brunel
,
N.
, &
Zecchina
,
R.
(
2007
).
Efficient supervised learning in networks with binary synapses
.
,
104
,
11079
11084
.
Burkitt
,
A. N.
,
Meffin
,
H.
, &
Grayden
,
D. B.
(
2004
).
Spike-timing-dependent plasticity: The relationship to rate-based learning for models with weight dynamics determined by a stable fixed point
.
Neural Comput.
,
16
,
885
940
.
Caporale
,
N.
, &
Dan
,
Y.
(
2008
).
Spike timing-dependent plasticity: A Hebbian learning rule
.
Annu. Rev. Neurosci.
,
31
,
25
46
.
Castellani
,
G. C.
,
Bazzani
,
A.
, &
Cooper
,
L. N.
(
2009
).
Towards a microscopic model of bidirectional synaptic plasticity
.
,
106
,
14091
14095
.
Castellani
,
G. C.
,
Quinlan
,
E. M.
,
Cooper
,
L. N.
, &
Shouval
,
H. Z.
(
2001
).
A biophysical model of bidirectional synaptic plasticity: Dependence on AMPA and NMDA receptors
.
,
98
,
12772
12777
.
Chklovskii
,
D. B.
,
Mel
,
B. W.
, &
Svoboda
,
K.
(
2004
).
Cortical rewiring and information storage
.
Nature
,
431
,
782
788
.
Elliott
,
T.
(
2008
).
Temporal dynamics of rate-based plasticity rules in a stochastic model of spike-timing-dependent plasticity
.
Neural Comput.
,
20
,
2253
2307
.
Elliott
,
T.
(
2010a
).
Discrete states of synaptic strength in a stochastic model of spike-timing-dependent plasticity
.
Neural Comput.
,
22
,
244
272
.
Elliott
,
T.
(
2010b
).
A non-Markovian random walk underlies a stochastic model of spike-timing-dependent plasticity
.
Neural Comput.
,
22
,
1180
1230
.
Elliott
,
T.
(
2011
).
The mean time to express synaptic plasticity in integrate- and-express, stochastic models of synaptic plasticity induction
.
Neural Comput.
,
23
(
1
),
124
159
.
Elliott
,
T.
, &
Lagogiannis
,
K. M.
(
2009
).
Taming fluctuations in a stochastic model of spike-timing-dependent plasticity
.
Neural Comput.
,
21
,
3363
3407
.
Elliott
,
T.
,
Howarth
,
C. I.
, &
,
N. R.
(
1996
).
Axonal processes and neural plasticity. I: Ocular dominance columns
.
Cereb. Cortex
,
6
,
781
788
.
Enoki
,
R.
,
Hu
,
Y.-L.
,
Hamilton
,
D.
, &
Fine
,
A.
(
2009
).
Expression of long-term plasticity at individual synapses in hippocampus is graded, bidirectional, and mainly presynaptic: Optical quantal analysis
.
Neuron
,
62
,
242
253
.
Fusi
,
S.
,
Drew
,
P. J.
, &
Abbott
,
L. F.
(
2005
).
Cascade models of synaptically stored memories
.
Neuron
,
45
,
599
611
.
Fusi
,
S.
, &
Senn
,
W.
(
2006
).
Eluding oblivion with smart stochastic selection of synaptic updates
.
Chaos
,
16
,
026112
.
Graupner
,
M.
, &
Brunel
,
N.
(
2007
).
STDP in a bistable synapse model based on CaMKII and associated signaling pathways
.
PLoS Comput. Biol.
,
3
,
e221
.
Gütig
,
R.
,
Aharonov
,
R.
,
Rotter
,
S.
, &
Sompolinsky
,
H.
(
2003
).
Learning input correlations through non-linear temporally asymmetric Hebbian plasticity
.
J. Neurosci.
,
23
,
3697
3714
.
Hinton
,
G. E.
, &
Sejnowski
,
T. J.
(
1983
).
Optimal perceptual inference
. In
Proceedings of the IEEE Computer Science Conference on Computer Vision and Pattern Recognition
(pp.
448
453
).
Los Alamitos, CA
:
IEEE Computer Society Press
.
Hopfield
,
J. J.
(
1982
).
Neural networks and physical systems with emergent collective computational abilities
.
,
79
,
2554
2558
.
Huang
,
K.
(
1987
).
Statistical Mechanics
.
Hoboken, NJ
:
Wiley
.
Kempter
,
R.
,
Gerstner
,
W.
, &
van Hemmen
,
J. L.
(
1999
).
Hebbian learning and spiking neurons
.
Phys. Rev. E
,
59
,
4498
4514
.
Leibold
,
C.
, &
Kempter
,
R.
(
2008
).
Sparseness constrains the prolongation of memory lifetime via synaptic metaplasticity
.
Cereb. Cortex
,
18
,
67
77
.
Montgomery
,
J. M.
, &
,
D. V.
(
2002
).
State-dependent heterogeneity in synaptic depression between pyramidal cell pairs
.
Neuron
,
33
,
765
777
.
Montgomery
,
J. M.
, &
,
D. V.
(
2004
).
Discrete synaptic states define a major mechanism of synapse plasticity
.
Trends Neurosci.
,
27
,
744
750
.
O'Connor
,
D. H.
,
Wittenberg
,
G. M.
, &
Wang
,
S. S.-H.
(
2005a
).
Dissection of bidirectional synaptic plasticity into saturable unidirectional process
.
J. Neurophysiol.
,
94
,
1565
1573
.
O'Connor
,
D. H.
,
Wittenberg
,
G. M.
, &
Wang
,
S. S.-H.
(
2005b
).
Graded bidirectional synaptic plasticity is composed of switch-like unitary events
.
,
102
,
9679
9684
.
Peng
,
Y.-R.
,
He
,
S.
,
Marie
,
H.
,
Zeng
,
S.-Y.
,
Ma
,
J.
,
Tan
,
Z.-.J.
, et al
. (
2009
).
Coordinated changes in dendritic arborization and synaptic strength during neural circuit development
.
Neuron
,
61
,
71
84
.
Petersen
,
C. C. H.
,
Malenka
,
R. C.
,
Nicoll
,
R. A.
, &
Hopfield
,
J. J.
(
1998
).
All-or-none potentiation at CA3-CA1 synapses
.
,
95
,
4732
4737
.
Purves
,
D.
, &
Lichtman
,
J. W.
(
1985
).
Principles of neural development
.
Sunderland, MA
:
Sinauer
.
Romani
,
S.
,
Amit
,
D. J.
, &
Amit
,
Y.
(
2008
).
Optimizing one-shot learning with binary synapses
.
Neural Comput.
,
20
,
1928
1950
.
Rubin
,
J. E.
,
Gerkin
,
R. C.
,
Bi
,
G. Q.
, &
Chow
,
C. C.
(
2005
).
Calcium time course as a signal for spike-timing-dependent plasticity
.
J. Neurophysiol.
,
93
,
2600
2613
.
Rubin
,
J.
,
Lee
,
D. D.
, &
Sompolinsky
,
H.
(
2001
).
Equilibrium properties of temporally asymmetric Hebbian plasticity
.
Phys. Rev. Lett.
,
86
,
364
367
.
Satel
,
J.
,
Trappenberg
,
T.
, &
Fine
,
A.
(
2009
).
Are binary synapses superior to graded weight representations in stochastic attractor networks
?
Cogn. Neurodyn.
,
3
,
243
250
.
Senn
,
W.
(
2002
).
Beyond spike timing: The role of nonlinear plasticity and unreliable synapses
.
Biol. Cybern.
,
87
,
344
355
.
Senn
,
W.
, &
Fusi
,
S.
(
2005
).
Convergence of stochastic learning in perceptrons with binary synapses
.
Phys. Rev. E
,
71
,
061907
.
Senn
,
W.
,
Markram
,
H.
, &
Tsodyks
,
M.
(
2001
).
An algorithm for modifying neurotransmitter release probability based on pre- and postsynaptic spike timing
.
Neural Comput.
,
13
,
35
67
.
Shouval
,
H. Z.
,
Bear
,
M. F.
, &
Cooper
,
L. N.
(
2002
).
A unified model of NMDA receptor-dependent bidirectional synaptic plasticity
.
,
99
,
10831
10836
.
Sobczyk
,
A.
, &
Svoboda
,
K.
(
2007
).
Activity-dependent plasticity of the NMDA-receptor fractional Ca2+ current
.
Neuron
,
53
,
17
24
.
Song
,
S.
,
Miller
,
K. D.
, &
Abbott
,
L. F.
(
2000
).
Competitive Hebbian learning through spike-timing-dependent synaptic plasticity
.
Nature Neurosci.
,
3
,
919
926
.
Swindale
,
N. V.
(
1996
).
The development of topography in the visual cortex: A review of models
.
Network: Comput. Neural. Syst.
,
7
,
161
247
.
van Ooyen
,
A.
(
2001
).
Competition in the development of nerve connections: a review of models
.
Network: Comput. Neural. Syst.
,
12
,
R1
R47
.
van Rossum
,
M. C. W.
,
Bi
,
G. Q.
, &
Turrigiano
,
G. G.
(
2000
).
Stable Hebbian learning from spike timing–dependent plasticity
.
J. Neurosci.
,
20
,
8812
8821
.
Willshaw
,
D. J.
,
Duneman
,
O. P.
, &
Longuet-Higgins
,
H.
(
1969
).
Nonholographic associative memory
.
Nature
,
222
,
960
962
.
Wittenberg
,
G. M.
, &
Wang
,
S. S.-H.
(
2006
).
Malleability of spike-timing-dependent plasticity at the CA3-CA1 synapse
.
J. Neurosci.
,
26
,
6610
6617
.
Yasuda
,
R.
,
Sabatini
,
B. L.
, &
Svoboda
,
K.
(
2003
).
Plasticity of calcium channels in dendritic spines
.
Nature Neurosci.
,
6
,
948
955
.