The Effect of STDP Temporal Kernel Structure on the Learning Dynamics of Single Excitatory and Inhibitory Synapses

Spike-Timing Dependent Plasticity (STDP) is characterized by a wide range of temporal kernels. However, much of the theoretical work has focused on a specific kernel – the “temporally asymmetric Hebbian” learning rules. Previous studies linked excitatory STDP to positive feedback that can account for the emergence of response selectivity. Inhibitory plasticity was associated with negative feedback that can balance the excitatory and inhibitory inputs. Here we study the possible computational role of the temporal structure of the STDP. We represent the STDP as a superposition of two processes: potentiation and depression. This allows us to model a wide range of experimentally observed STDP kernels, from Hebbian to anti-Hebbian, by varying a single parameter. We investigate STDP dynamics of a single excitatory or inhibitory synapse in purely feed-forward architecture. We derive a mean-field-Fokker-Planck dynamics for the synaptic weight and analyze the effect of STDP structure on the fixed points of the mean field dynamics. We find a phase transition along the Hebbian to anti-Hebbian parameter from a phase that is characterized by a unimodal distribution of the synaptic weight, in which the STDP dynamics is governed by negative feedback, to a phase with positive feedback characterized by a bimodal distribution. The critical point of this transition depends on general properties of the STDP dynamics and not on the fine details. Namely, the dynamics is affected by the pre-post correlations only via a single number that quantifies its overlap with the STDP kernel. We find that by manipulating the STDP temporal kernel, negative feedback can be induced in excitatory synapses and positive feedback in inhibitory. Moreover, there is an exact symmetry between inhibitory and excitatory plasticity, i.e., for every STDP rule of inhibitory synapse there exists an STDP rule for excitatory synapse, such that their dynamics is identical.


Introduction
Spike timing dependent plasticity (STDP) is a generalization of the celebrated Hebb postulate that ''neurons that fire together wire together'' to the temporal domain, according to the temporal order of the presynaptic and postsynaptic spike times. A temporally asymmetric Hebbian (TAH) plasticity rule has been reported in experimental STDP studies of excitatory synapses [1][2][3], in which an excitatory synapse undergoes long-term potentiation when presynaptic firing precedes the postsynaptic firing and long-term depression is induced when the temporal firing order is reversed, e.g., Figure 1A.
Many theoretical studies [4][5][6][7][8][9] that followed these experiments used an exponentially decaying function to represent the temporal structure of the STDP. Throughout this paper we term this STDP pattern the ''standard exponential TAH''. Gütig and colleagues [7] also provided a convenient mathematical description for the dependence of STDP on the synaptic weight in the standard exponential TAH STDP rule: where w[ 0,1 ½ is the dynamic parameter that describes the synaptic strength; Dw is the modification of w following pre (2) or post (+) synaptic firing; Dt is the time difference between the presynaptic and postsynaptic firing; l is the learning rate; t is the temporal decay constant and m[ 0,1 ½ and aw0 are dimensionless parameters of the model that characterize the weight dependent component of the STDP rule. This representation introduces a convenient separation of variables, in which the synaptic update is given as a product of two functions. One function is the temporal kernel of the STDP rule, i.e. K Dt ð Þ, and the other is the weight dependent STDP component, i.e. f + w ð Þ. For convenience, throughout this paper we shall adopt the notation of Gütig and colleagues for the weight dependence of the STDP rule, f z={ w ð Þ, equations (3) - (4). This function, f z={ w ð Þ, is characterized by two parameters: the relative strength of depression -a, and the degree of non-linearity in w of the learning rule -m. Note, that other choices for f z={ w ð Þ have also been used in the past [5], [10], [11].
Properties of the ''standard exponential TAH'' As previously shown [6,7], the standard exponential TAH model can generate positive feedback that induces bi-stability in the learning dynamics of an excitatory synapse. For a qualitative intuition into this phenomenon, consider the case of a weightindependent STDP rule, also termed the additive model, i.e., m~0. If the synaptic weight is sufficiently strong, there is a relatively high probability that a presynaptic spike will be followed by a postsynaptic spike. Hence, causal events (i.e., Dtw0 post firing after pre) are more likely to occur than a-causal events (with Dtv0). Because the STDP rule of the standard exponential TAH model implies LTP for Dtw0 there is a greater likelihood for LTP than for LTD. Thus, a ''strong'' synapse will tend to become stronger. On the other hand, if the synaptic weight is sufficiently weak, then pre and post firing will be approximately uncorrelated. As a result, the stochastic learning dynamics will randomly sample the area under the STDP temporal kernel. Here we need to consider two types of parameter settings. If the area under the causal branch in equation (1) is larger than the area under the acausal branch, av1, the net effect is LTP for weak synapses as well. Thus, in this case, all synapses will potentiate until they reach their upper saturation bound at 1. Hence, the regime of av1, in this case, is not interesting. On the other hand, if the area under the a-causal branch is larger than the area under the causal branch, aw1, random sampling of the STDP temporal kernel by the stochastic learning dynamics (in the limit of weakly correlated pre-post firing, mentioned above) will result in LTD. Thus, in the interesting regime, aw1, a ''weak'' synapse will tend to become weaker; thus producing the positive feedback mechanism that can generate bi-stability.
It was further shown [7] that this positive feedback can be weakened by introducing the weight dependent STDP component via the non-linearity parameter m in equations (3) and (4). Setting mw0 decreases the potentiation close to the upper saturation bound and decreases the depression close to the lower saturation bound; thus, for sufficiently large values of m the learning dynamics will lose its bi-stability.
Experimental studies have found that the temporally asymmetric Hebbian rule is not limited to excitatory synapses and has been reported in inhibitory synapses as well [12]. Similar reasoning shows that in the case of inhibitory synapses the standard exponential TAH induces negative feedback to the STDP dynamics. It was shown [13] that this negative feedback acts as a homeostatic mechanism that can balance feed-forward inhibitory and excitatory inputs. Interestingly, Vogels and colleagues [14] studied a temporally symmetric STDP rule for inhibitory synapses, and reported that this type of plasticity rule also results in negative feedback that can balance the feed-forward excitation. This raises the question whether inhibitory plasticity always results in a negative feedback regardless of the temporal structure of the STDP rule? On the other hand, theoretical studies have shown that the inherent positive feedback of excitatory STDP causes the learned excitatory weights to be sensitive to the correlation structure of the pre-synaptic excitatory population -for different choices of STDP rules [5,7,10]. Does STDP dynamics of excitatory synapses always characterized by a positive feedback?

Outline
Although theoretical research has emphasized the standard exponential model, empirical findings report a wide range of temporal kernels for both excitatory and inhibitory STDP; e.g., [1,12,[15][16][17][18][19], (see also the comprehensive reviews by Caporale and Dan [20] and Vogels and colleagues [21]). Here we study the effect of the temporal structure of the STDP kernel on the resultant synaptic weight for both excitatory and inhibitory synapses. This is done in the framework of learning of a single synapse in a purely feed-forward architecture, as depicted in Figure 2. First, we suggest a useful STDP model that qualitatively captures these diverse empirical findings. Below we define our STDP model. This model serves to study a large family of STDP learning rules. We derive a mean field Fokker-Planck approximation to the learning dynamics and show that it is governed by two  (7) and (8) with the ''standard exponential TAH'' as a reference. Each plot (normalized to a maximal value of 1 in the LTP branch) qualitatively corresponds to some experimental data. In all plots, the blue curve represents the potentiation branch K z , the red curve represents the depression branch {K { and the dashed black curve represents the superposition/ sum of K z {K { . For simplicity, all plots were drawn with the same t~20ms. (A) The ''standard exponential TAH'' [1,18]. (B) +K + Dt; 1,t ð Þ Alternate approximation to the standard exponential TAH [1,18]. (C) +K + Dt; {1,t ð ÞTemporally asymmetric Anti-Hebbian STDP [15]. (D) +K + Dt; 0:75,t ð ÞTAH variation [12,19]. (E) +K + Dt; 0,t ð Þ Temporally symmetric Hebbian STDP [16,17]. (F) +K + Dt; {0:75,t ð ÞVariation to a temporally asymmetric Anti-Hebbian STDP [19] doi:10.1371/journal.pone.0101109.g001 global constants that characterize the STDP temporal kernel. Stability analysis of the mean-field solution reveals that the STDP temporal kernels can be classified into two distinct types: Class-I, which is always mono-stable, and Class-II that can bifurcate to bistability. Finally, we discuss the symmetry between inhibitory and excitatory STDP dynamics.

Generalization of the STDP rule
In order to analyze various families of STDP temporal kernels found in experimental studies [1,12,[15][16][17][18][19] we represent the STDP as the sum of two independent processes: one for potentiation and the other for depression. The synaptic update rule that we use throughout this paper is given by: Note that the main distinction between equation (5) and (1), is that here, equation (5), the z={ signs denote potentiation and depression, respectively; while in equation (1) the z={ signs denote causal/a-causal branch. Thus, in our model for every Dt the synapse is affected by both potentiation and depression; whereas, according to the model of equation (1)

Skew-Normal kernel
Here we used the Skew-Normal distribution function to fit the temporal kernels of the STDP rule, K + Dt ð Þ. Note that the specific choice of the Skew-Normal distribution is arbitrary and is not critical for the analysis below. Other types of functions may serve as well. The ''Skew-Normal distribution'' is defined by: where j is the temporal shift, t is the temporal decay constant, and Q is a dimensionless constant that affects the skewness of the curve and erf x ð Þ is the Gaussian error function. It is also useful to reduce the number of parameters that define the STDP temporal kernel. Thus, we define: where h[ {1,1 ½ is a single continuous dimensionless parameter of the model that characterizes the STDP temporal kernel and t is the time constant of the exponential decay of the potentiation branch. The mapping of j!h 1{h 2 À Á ensures that the temporal shift parameter, j, will be zero for h~{1,0,1. In order to obtain temporally symmetric Mexican hat STDP rule for h~0 one needs to demand t dep wt pot , where t pot~t and t dep~t 1z0:5 1{h 2 À Á 2 . We also required t dep~tpot for h~+1, in order to be compatible with several previous studies. This reduction in parameters was chosen in order to capture the qualitative characteristics of various experimental data; however, other choices are also possible. Figure 1B-F illustrates how one can shift continuously from a temporally asymmetric Hebbian kernel (h~1, Figure 1B) to a temporally asymmetric anti-Hebbian kernel (h~{1, Figure 1F). Figure 1A shows the temporal kernel of the standard exponential TAH model, compare with h~1, Figure 1B.

''Mean field'' Fokker-Planck approximation
We study the STDP dynamics of a single feed-forward synapse to a postsynaptic cell receiving other feed-forward inputs through synapses that are not plastic. We assume that all inputs to the cell obey Poisson process statistics with constant mean firing rate, r pre ; that the presynaptic firing of the studied synapse is uncorrelated with all other inputs to the postsynaptic neuron; and that the synaptic coupling of a single synapse is sufficiently weak. The STDP dynamics is governed by two factors: the STDP rule and the pre-post correlations. To define the dynamics one needs to describe how the pre-post correlations depend on the dynamical variable, w. Under the above conditions one may assume that the contribution of a single pre-synaptic neuron that is uncorrelated with the rest of the feed-forward input to the post-synaptic neuron will be small. Thus, it is reasonable to approximate the pre-post correlation function (see Methods -equation (26)) up to a first order in the synaptic strength w (e.g., [8,[22][23][24]), yielding: where r pre=post is the instantaneous firing of the pre/post synaptic cell represented by a train of delta functions at the neuron's spike times (see Methods), r pre =r post is the pre/post synaptic mean firing rate; and the function c D ð Þ describes the change in the conditional mean firing rate of the postsynaptic neuron at time Dzt following a presynaptic spike at time t. Note that we use upper case C to represent the full pre-post correlations, C~Sr pre r post T, whereas c denotes the first order term in the synaptic weight, w, of these correlations.
In the limit of a slow learning rate, l?0, one obtains the meanfield Fokker-Planck approximation to the stochastic STDP dynamic (see Methods -equation (27)), and using the linear approximation of the pre-post correlations, equation (9), yields: ÞdD denotes the mean over time (using equation (9) with c D ð Þ~0 for Dv0). In our choice of parameterization, K z={ are set to have the same integral; i.e., K z~K{~ K K. The difference between the strength of potentiation and depression of the STDP rule is controlled by the parameter a (equation (4)). Substituting expressions (3) & (4) into equation (10) yields: where X z={ :cK z={ K z={ are constants that govern the mean-field dynamics. A fixed point solution, w Ã , of the meanfield Fokker-Planck dynamics, S _ w wT w ð Þ~0, satisfies: Numerical simulations -the steady state of STDP learning We performed a series of numerical simulations to test the approximation of the analytical result of the mean-field approximation at the limit of vanishing learning rate, using a conductance based integrate-and-fire postsynaptic neuron with Poisson feedforward inputs (see Methods for details; a complete software package generating all the numerical results in this manuscript can be downloaded as File S1). We simulate a single postsynaptic neuron receiving feed-forward input from a population of N E~1 20 excitatory neurons and N I~4 0 inhibitory neurons firing independently according to a homogeneous Poisson process with rate r pre~1 0spikes=s. All synapses except one (either excitatory or inhibitory) were set at a constant strength (of 0.5). The initial conditions for the plastic synapse were as specified bellow.
We first estimated the spike triggered average (STA) firing rate of a single presynaptic neuron triggered on postsynaptic firing, in order to approximate the function c D ð Þ, equation (9). Figure 3 shows the STAs of excitatory (A) and inhibitory (B) synapses for varying levels of synaptic weights (color coded), as were estimated numerically (dots). The dashed lines show smooth curve fits to the STA. The specific temporal structure of these curves depends upon particular details of the neuronal model. Nevertheless, the linear dependence on the synaptic weight is generic for weak synapses; thus, in line with the assumed linearity of the model, equation (9). The STA shows the conditional mean firing rate of the presynaptic neuron, given that the post fired at time t~0. In the limit of weak coupling, w?0, pre and post firing are statistically independent and the conditional mean equals the mean firing rate of the pre, r pre~1 0spikes=s. For an excitatory synapse, as the synaptic weight is increased the probability of a post spike following pre will also increase. Consequently, so will the likelihood of finding a pre spike during a certain time interval preceding a post spike. Hence, the STA of an excitatory synapse is expected to show higher amplitude for stronger synapse (as shown in A). Correspondingly, the STA of an inhibitory synapse is expected to show a more negative amplitude for stronger synapse (as shown in B).
To fit the STA with an analytic function we used c D ð Þã both the inhibitory and excitatory cases. This ad hoc approximation serves to enable the numerical integration that calculates the constants X z={ that govern equation (11) for the mean field approximation to the STDP dynamics.
All the richness of physiological details that characterize the response of the post-synaptic neuron affect the STDP dynamics only via the two constants X z and X { . These two constants X z={ denote the overlap between the temporal structure of the pre-post correlations, c D ð Þ, and the temporal kernel, K z={ , of the potentiation/depression kernel, respectively, Figure 4. Consequently, as c D ð Þ, is positive for excitatory synapses and negative for inhibitory synapses -so are the constants X z={ . In addition, as the correlations in our model are causal, c D ð Þ~0 for tv0, the constant X z (X { ) is expected to decay to zero when the STDP kernel K z (K { ) vanishes from the causal branch, h?1 (h?{1). For the specific choice of parameters in our simulations, DX z D obtains its maximal value at h~1. However, one may imagine other choice of parameters in which DX z D will obtain its maximal value at h[ 0,1 ð Þ. Note, from Figure 4, that the crossing of the X z and X { curves, is coincidentally almost the same for both synapse types, and is obtained at &{0:2. The significance of this point is discussed below.
Fixed point solutions for the STDP dynamics Figure 5 shows w Ã as a function of a for different values of m (color coded, note that a and m are the parameters that characterize the synaptic weight dependence of the STDP rule, equations (3) and (4)). The panels depict different STDP setups that differ in terms of the temporal kernels as well as the type of synapse (excitatory/inhibitory). These two factors affect the mean field equations via X z={ . The dashed lines show the solution to the fixed point equation, equation (12), using the numerically calculated X z={ . The fixed points were also estimated numerically by directly simulating the STDP dynamics in a conductance based integrate and fire neuron (circles and error bars).
For the estimation of the steady state value of the synaptic weight, the simulations were set to run for 5 hours of simulation time, which, according to manual offline analyses of convergence time scales, is much more than twice the time required for the system to converge and fluctuate around its steady state. The circles and error bars depict the mean 6 standard deviation of the synaptic weight, as estimated from the last 2.5 hours of the simulation (weights were recorded at a 1 Hz sample rate). Note the high agreement between the fixed point solution (w Ã ) of equation (12), and the asymptotic synaptic weight (w 0 ) as estimated by the numerical simulation (regression coefficient of 1+5 : 10 {4 with R 2 w0:999 when performing a regression test on the entire set {w 0 ,w Ã } presented in each of the panels).
The panels of Figures 5A and 5B compare the standard exponential TAH rule of equation (2), in A, and our current STDP model with h~1 in B, for a representative set of parameters {(a i ,m j )} applied to the examined synapse (middle column for inhibitory synapse and right column for excitatory synapse). Note that some lines may overlap each other near the boundaries: w 0~0 ,1. As is evident from the figures, the results of the two models coincide. In particular, the Hebbian STDP dynamics of inhibitory synapses is characterized by a one to one function w Ã of a and there is no bi-stability, as previously reported [13,14]. On the other hand, the Hebbian STDP of excitatory synapses is characterized by bi-stable solutions at low levels of m below a certain critical value, see e.g., [6,7]. Thus, the current model with h~1 coincides with previous results.
The panels of Figure 5F show the results of a temporally asymmetric Anti-Hebbian STDP with h~{1. In striking contrast to the Hebbian STDP, in this case, inhibitory plasticity is characterized by bi-stability whereas, the excitatory plasticity is characterized by mono-stability.
The panels of Figures 5C and 5E explore two other types of asymmetric rules (Hebbian and Anti-Hebbian respectively). These results show similar behavior as 5B and 5F in terms of the classification of STDP kernels discussed in the next section.
The panels of Figure 5D show the results of the symmetric STDP with h~0 -note, that the dynamics of inhibitory synapse under the symmetric STDP rule, is characterized by a one to one function w Ã of a corresponding to negative feedback, as previously reported [14].

Stability of the fixed point solution
The stability of the fixed point solution w Ã to equation (12) is determined by the sign of the partial derivative of the dynamical equation, equation (11), with respect to the synaptic weight: On the other hand, examination of Figure 5 suggests that the stability of the fixed point is governed by the sign of Lw Ã =La. Taking the logarithm and the derivative with respect to a of both sides of equation (12), one obtains: sign where the last equality holds as aw0. At the fixed point, substituting equation (12) into equation (13) one obtains:   (11) for the excitatory and inhibitory synapses of the neuronal model used in our numerical simulations, as a function of h. These values were calculated using numerical integration (see File S1) with K z={ Dt; h,t ð Þas defined by equations (7) and (8), with t~20msec as set throughout the simulations, and with the fitted formula for c D ð Þ. doi:10.1371/journal.pone.0101109.g004 The Effect of STDP Temporal Kernel PLOS ONE | www.plosone.org  (12) (dotted lines), is compared to the asymptotic synaptic weight (w 0 ) (circles), of a single synapse learning dynamics for various learning rules as defined by equation (5). Each of the panels in the middle column (for inhibitory synapse) and in the right column (for excitatory synapse) explores the weight dependent STDP component, f + w ð Þ of equations (3) and (4), for representative set of m (shown by different colors as depicted in the legend) as a function of a. The different rows correspond to different STDP kernels, K + Dt ð Þ as shown by the panels in the left column. The circles and error bars represent the mean and standard deviation of the synaptic weight (w 0 ), calculated over the trailing 50% of each learning dynamics simulation (see Methods). The mean field constants {X z ,X { } were numerically calculated using the c D ð Þ constants estimated as in Figure 3. The dotted lines were computed by equation (12) that was calculated for 10,000 sequential values of w in 0,1 ½ . To this end, we replaced m~0 with m~10 {6 in order to use equation (12) to plot the dashed red line. Initial conditions for the simulations: for the majority of the simulations we have simply used w~0:5 as initial condition for the plastic synaptic weight. In order to show the bi-stable solutions in panels (A2, B2, F1), for m~0,0:01,0:02 and a~1:01,1:03, . . . 1:19 we ran two simulations one with initial condition w~0 and another with initial condition w~1. (A0-F0) are the STDP kernels (as in Figure 1 Yielding: Hence, for w Ã X z w{1, the fixed point solution in Figure 5 is stable in segments with negative slope, and unstable in segments with positive slope. Note that in our simulation setup DX z Dv1 (cf. Figure 4); thus, the condition w Ã X z w{1 holds for all values of h in our case.
Revisiting the different scenarios depicted in Figure 5, we note the existence of two qualitatively different behaviors; namely, one that can only show mono-stability (A, C, and F) and the other has the potential for bi-stability (in panels B, D, and E). We use this behavior to classify the different STDP temporal kernels that are parameterized by the single variable h. We shall term ''class-I temporal kernels'' the temporal kernels such that w Ã is monostable for all aw0; m[ 0,1 ½ . We shall term ''class-II temporal kernels'' the temporal kernels such that w Ã is bi-stable for some m[ 0,m c ½ Þ and some a m ð Þ. Note that this classification depends on the type of synapse (which via c D ð Þ together with h determine X z={ ). In addition, we note the existence of a special solution at w Ã~1 =2 that is invariant to m, and enables us to obtain a simple condition for this classification. In class-I kernels the derivative Lw Ã =La at w Ã~1 =2 is always negative, whereas in class-II models there is a critical value of m below which the derivative changes its sign.
The ''m-invariant'' solution and the critical m Vm, the solution of the fixed point equation, equation (12), at w Ã~ŵ w:1=2 is m-invariant. For a given STDP temporal kernel (h), i.e. a given set of {X z ,X { } (see Figure 4; and note that X z={ are also determined by the pre-post correlation structure via c D ð Þ), the solution ofŵ w~1=2 is obtained with:â Substituting the m-invariant solution, equation (18), into equation (15), yields Thus, the condition for instability of the m-invariant solution is: Thus, for m c v0 the m-invariant solution,ŵ w, is stable for all values of m[ 0,1 ½ and the STDP rule is class-I for that synapse. On the other hand, if m c w0 the STDP rule is class-II. This classification depends solely on the values of {X z ,X { }. In our simulation setup DX z={ Dv1 (see Figure 4), thus the classification of the parameter combinations is simply determined by the sign of (X z {X { ); i.e. the manifold that is determined by the condition {X z~X{ } separates the parameter space (that characterizes the STDP rule and the synapse) between class-I and class-II.
Bimodal distribution nearâ a Figure 6 depicts (using numerical simulations with set of class-II parameters) the bifurcation plots for the learning dynamics for inhibitory (A, B) and excitatory (C, D) synapses. For inhibitory synapses the anti-Hebbian (h~{1) plasticity rules were chosen, and for the excitatory synapses, the Hebbian (h~1). The panels show the resultant distribution of the synaptic weight color-coded after 216101 of 5 hours of simulations for 21 values of the bifurcating argument (either m or a) along the abscissa. In order to calculate the synaptic weight distribution for the set of parameters without the bias of initial conditions, 101 simulations were performed with different initial weight values evenly spaced from 0 to 1. The rationale for running the simulations for 5 hours each was to make sure that the learning dynamics had reached a steady state regime and the synaptic weight fluctuated around it for the entire trailing 2.5 simulation hours. During these trailing 2.5 simulation hours, the synaptic weights were recorded at a 1 Hz sample rate. For the estimation of the weight distribution, all the samples from the 101 simulations (differing only by their initial conditions) were used with 40 evenly spread bins between 0 and 1.
As expected from the analysis, there was a bifurcation along the m dimension (top panels), in which above m c the distribution was uni-modal whereas below m c the distribution was bi-modal. Along the a dimension (bottom panels) the distribution resembled the theoretical (dashed) curves of Figure 5 (without the unstable segment of Lw Ã =Law0).

Symmetry and phase transition along h
The high degree of similarity between the simulation results for inhibitory and excitatory synapses ( Figure 5) stems from the fact that they obey the same mean-field equation (11), albeit with a different set of parameters. Thus, an excitatory synapse, w exc , with a specific choice of parameters {a,m,X z ,X { } obeys the exact same mean-field equation as (1{w inh ), where w inh is an inhibitory synapse with the transformed set of parameters 1 a and a somewhat different learning constant (note that X z={ are positive for excitatory synapses and negative for inhibitory ones, see Figure 4). This symmetry is illustrated for different STDP temporal kernels in Figure 7, where the mean field fixed point, w Ã , is plotted as a function of a for different values of h (color coded) at m&0. The different h were chosen around h c which is defined by the condition X z~X{ (see Figure 4) to display the phase transition from class-I to class-II along this parameter. Coincidentally, in our simulations and the chosen model (equations (7) and (8)), this specific h c was almost the same for excitatory and inhibitory synapses; i.e. for both synapses h c w{0:22 and h c v{0:21(see Figure 4). Under these conditions, for an excitatory synapse, hƒ{0:22 defines the class-I kernels, and h §{0:21 the class-II, whereas for an inhibitory synapse, hƒ{0:22 defines the class-II kernels, and h §{0:21 the class-I.

Discussion
The computational role of the temporal kernel of STDP has been studied in the past. Câteau and Fukai [8] provided a robust Fokker-Planck derivation and analyzed the effects of the structure of the STDP temporal kernel. However, their analysis focused on excitatory synapses and the additive learning rule (m~0). Previous studies have linked the Hebbian STDP of inhibition with negative feedback which acts as a homeostatic mechanism that balances the excitatory input to the postsynaptic cell [13,14]. Positive feedback and bi-stability of STDP dynamics have been reported only for excitation, and linked to sensitivity to the input correlation structure [6,7]. Here it was shown that the STDP of both excitation and inhibition can produce either positive or negative feedback depending on the parameters of the STDP model. Thus, for example, it was reported that both a temporally asymmetric Hebbian STDP (h~1) and a temporally symmetric learning rule (h~0) for inhibitory synapses generate negative feedback [13,14]. These reports are in-line with our finding that the critical h for transition from negative to positive feedback for inhibition is negative (h c &{0:2).
In general, STDP dynamics of single synapses was classified here into two distinct types. With class-I temporal kernels, the  dynamics is characterized by a negative feedback and has a single stable fixed point. In contrast, class-II temporal kernels are characterized by a sub-parameter regime in which the system is bistable (has positive feedback), and another sub-parameter regime with negative feedback. However, the mechanism that generates the negative feedback, (i.e., the stabilizing mechanism) in the two classes is different in nature. Whereas in class-I the negative feedback is governed by the convolution of the pre-post correlations with the temporal kernel, (i.e. the mean field constants X z={ , similar to the homeostatic mechanism in [13]), in class-II, the stabilizing mechanism is the non-linear weight dependent STDP component, f + w ð Þ. Hence, there is no reason a-priori to assume that the negative feedback in class-II should act as a homeostatic mechanism.
We found that there is no qualitative difference between the STDP of excitatory and inhibitory synapses and that both can exhibit class-I and class-II dynamics. Moreover, there is an exact symmetry between the excitatory and inhibitory STDP under a specific mapping of the parameters {a,m,X z ,X { }. This symmetry results from the fact that the mean-field dynamics depend solely on the global mean field constants X z={ . It is important to note that although neural dynamics is rich and diverse, due to the separation of time scales in our problem, the STDP dynamics only depends on these fine details via the global mean field constants X z={ .
Certain extensions to our work can be easily implemented into our model without altering the formalism. For example, empirical studies report different time constants for depression and potentiation, e.g. [1]. However, although in our simulations we used identical time constants at DhD~1, for DhDv1 the depression time constant is larger than the potentiation time constant in our simulations. Moreover, our analytical theory depends on the time constants only via X z={ . Consequently, changing time constants or any other manipulation to the temporal kernel can be incorporated into our mean-field theory by modifying X z={ . Similarly, assuming separation of time-scales between short term and long term plasticity, the effect of short term plasticity can be incorporated by modifying X z={ accordingly.
STDP has also been reported to vary with the dendritic location, e.g. [18,25]. For a single synapse this effect can also be modeled by a modification of the parameters X z={ . However, the importance of the dendritic dependence of STDP may reside in the interaction with other plastic synapses along different locations on the dendrite. Network dynamics of a 'population' of plastic synapses is beyond the scope of the current paper and will be addressed elsewhere.
In our model we assumed that the contribution of different "STDP events" (i.e., pre-post spike pairs) to the plastic synapse are summed linearly over all pairs of pre and post spikes, see e.g. equation (21). However, empirical findings suggest that this assumption is a mere simplification, and that STDP depends on pairing frequency as well as triplets of spike time and bursts of activity, e.g. [3,[26][27][28][29][30]. The computational implications of these and other non-linear interaction of spike pairs in the learning rule, as well as the incorporation of non-trivial temporal structure into the correlations of the pre-synaptic inputs to the cell are beyond scope of the current paper.
Empirical studies have reported a high variability of STDP temporal kernels over different brain regions, locations on the dendrite and experimental conditions, e.g., [1,12,15,[17][18][19]. Here we represented the STDP rule as the sum of two separate processes, one for potentiation and one for depression with an additional parameter, h, that allows us to continuously modify the temporal kernel and qualitatively obtain a wide spectrum of reported data. Representation of STDP by two processes has been suggested in the past. Graupner and Brunel [31], for example, proposed a model for synaptic plasticity in which the two processes (long term potentiation and depression) are controlled by calcium level. Thus, in their model the control parameter is a dynamical variable that may alter the plasticity rule in response to varying conditions. In our work, however, we did not model the dynamics of h. Moreover, we assumed that h remains constant during timescales that are relevant for synaptic plasticity. It is, nevertheless, tempting to speculate on a metaplasticity process [32,33] in which the temporal structure of the STDP rule is not hard wired and can be controlled and modified by the central nervous system. Thus, in addition to controlling the learning rate, l, or the relative strength of potentiation-depression, a, a metaplasticity rule may affect the learning process by modifying the degree of 'Hebbianitty', h. Such a hypothesis, if true, may account for the wide range of STDP kernels reported in the experimental literature. How can such a hypothesis be probed? One option for addressing this issue is to try and characterize h during different time points and study its dynamics. One would expect to find that h (for excitatory synapses) decreases with time in cases where the neural network has been reported to becomes less sensitive to its input statistics, for example during developmental changes.

''Mean field'' Fokker-Planck approach for the learning dynamics
From the synaptic update rule, equation (5), changes in the synaptic weight, w, at time t, result from either pre or post synaptic firing at time t, affecting both the depression and potentiation branches (functions) of the adaptation rule. Thus: where X stands for Excitation or Inhibition, N X is the number of synapses, t ½ z :max t,0 ð Þ is the dimensionless time value (in seconds), and t i j n o j are the spike times of synapse i. For the temporal characteristic of the a-shape response we chose to use t E~tI~5 ms, and for the conductance coefficient g 0 X our constant is scaled by N X as elaborated below.
In order to estimate the postsynaptic membrane potential in equation (28), the software performs the integration of the synaptic and leak currents using the Euler method with a Dt~1ms step size. The rationale for using such a low resolution step size and its verification are discussed below.
Modeling presynaptic activity. Throughout the simulations in this work, presynaptic activities were modeled by an independent homogeneous Poisson processes, with stationary mean firing rate r pre~1 0 spikes=s. To this end, each of the inputs was approximated by a Bernoulli process generating binary vectors defined over discrete time bins of Dt~1ms. These vectors were then filtered using a discrete convolution a-shaped kernel (as defined above) with a limited length of 10t X (after which this kernel function is zero for all practical purposes). In all simulations we used: N E~1 20,N I~4 0.
Conductance constants. In order to be compatible with previous studies; e.g., [7,13], and to have simulations that are executed with a robust and generic software package accompanying this manuscript as File S1, we scaled the synaptic conductance inversely to the number of synaptic inputs in our simulations. We used the following scaling formula g 0 X~g R X S X , with: g R E~3 0nS, S E~1 000=N E , g R I~5 0nS and S I~4 00=N I , where N E ,N I are the number of excitatory and inhibitory presynaptic inputs, respectively.
The learning rate. The simulations of the STDP process were carried out to obtain the asymptotic weight distribution of the plastic synapse. Convergence to the asymptotic region was accelerated by manipulating the learning rate constant l of equation (1). The software code was designed to support a given vector of l for each minute of the simulation. Specifically we used the following formula to generate this vector: Postsynaptic spike time accuracy vs. simulation step size resolution. Figure 5 shows the remarkable match between the fixed point solution (w Ã ) of equation (12), and the asymptotic synaptic weight (w 0 ) of the simulations; the regression coefficient on the entire set {w 0 ,w Ã } in all the panels is 1+5 : 10 {4 with R 2 w0:999 when using an integration step of size 1ms. Tests of this kind were performed on simulations using integration steps ranging from 0:1ms to 1ms in two calculation modes (see below), and it was found that higher resolution provides a better match to the analytical solution. However, the key feature that contributes to this high degree of similarity between the analysis and the simulations (more than an order of magnitude for the error term 1{R 2 ) was the definition of the spike times of the postsynaptic cell rather than a 106 decrease of the integration step size. The spike times of an integrate and fire neuron are defined as the times in which its membrane potential crossed the firing threshold, t Ã . However, in the numerical simulations we used discrete times, nDt bin f g N n~0 . In previous work we define the time of the post-synaptic firing by the last discrete time preceding the threshold-crossing time to: t post~n Dt such that nDtƒt Ã v(nz1)Dt. This choice may change the causal order of pre-post firing (from pre before post to simultaneous firing) at time intervals of the time-bin. Consequently, it will affect the STDP rule -mainly when kernels that are discontinuous at zero are used. Here we defined the spike time of the post-synaptic neuron to be: t post~( nz1=2)Dt such that nDtƒt Ã v(nz1)Dt (i.e., shifted by half a time-bin from previous definition); thus, this manipulation retains the causality of firing.

Supporting Information
File S1 This package (1Syn-STDP4PLOS.zip) is a Matlab set of scripts and utilities that includes all the numerical simulations that were used to produce the figures in this manuscript. It also contains all the scripts that generated the figures. The scripts in the main folder are divided into two categories. The files that begin with ''Bat'' execute the numerical simulations, and the ones that begin with ''Plot'' generate the figures. All the supporting numerical utilities are stored in the sub folder ''CommonLib''. (ZIP) Author Contributions