Closed-Form Treatment of the Interactions between Neuronal Activity and Timing-Dependent Plasticity in Networks of Linear Neurons

Network activity and network connectivity mutually influence each other. Especially for fast processes, like spike-timing-dependent plasticity (STDP), which depends on the interaction of few (two) signals, the question arises how these interactions are continuously altering the behavior and structure of the network. To address this question a time-continuous treatment of plasticity is required. However, this is - even in simple recurrent network structures - currently not possible. Thus, here we develop for a linear differential Hebbian learning system a method by which we can analytically investigate the dynamics and stability of the connections in recurrent networks. We use noisy periodic external input signals, which through the recurrent connections lead to complex actual ongoing inputs and observe that large stable ranges emerge in these networks without boundaries or weight-normalization. Somewhat counter-intuitively, we find that about 40% of these cases are obtained with a long-term potentiation-dominated STDP curve. Noise can reduce stability in some cases, but generally this does not occur. Instead stable domains are often enlarged. This study is a first step toward a better understanding of the ongoing interactions between activity and plasticity in recurrent networks using STDP. The results suggest that stability of (sub-)networks should generically be present also in larger structures.

All these observations pose a problem, because STDP can generically lead to the situation that individual pre-and post-synaptic spike pairs influence the synaptic strength of their connection and thereby in a recurrent way immediately also the activity in the network. The question arises, thus, whether such fast synaptic changes would carry any functional significance or whether they will be averaged out by the ongoing network activity? At the moment this issue is left open because resolving it would require measuring or calculating the step-by-step changes of synaptic connectivity together with their activity for a large number of neurons and ideally for every spike pair.
Experimentally this is at the moment still impossible and so far there are also no theoretical tools available (apart from simulations) to address the issue of ongoing mutual interactions between activity and plasticity in networks. For this, analytical methods would be particularly helpful with which it would be possible to predict the dynamics of synaptic connections under STDP.
As discussed above, this is typically investigated with mean-field methods (van Rossum et al., 2000;Gütig et al., 2003) that have average values in their focus. However, here we would like to examine the fine temporal dynamics of each weight as, fundamentally, the IntroductIon At the network level recent theoretical advances have made it possible to analytically calculate the activity patterns for certain types of simplified neuron models as long as one keeps the synaptic connectivity fixed (Memmesheimer and Timme, 2006a,b). Concerning plasticity, the situation is reversed. When fixing the activity pattern or rather, when considering long-term activity averages, it is for many plasticity rules possible to calculate the overall development of synaptic weights. Older results exist which provide specific solutions for estimating synaptic strengths in certain types of networks and learning rules, all of which however need to constrain structure or dynamics of the system in different ways (Hopfield, 1982;Miller and MacKay, 1994;Roberts, 2000;van Rossum et al., 2000;Kempter et al., 2001;Burkitt et al., 2007;Kolodziejski et al., 2008) and many of them do not resemble spike-timing-dependent plasticity (STDP, Magee and Johnston, 1997;Markram et al., 1997). Concerning STDP, in large networks one finds that over longer times the distribution of synaptic efficiencies usually develops into either an unimodal (van Rossum et al., 2000;Gütig et al., 2003) or a bimodal (Song et al., 2000;Izhikevich et al., 2004) shape. These results mainly rely on simulations (Song et al., 2000;van Rossum et al., 2000) or mean-field approximations (van Rossum et al., 2000;Gütig et al., 2003). Additionally, in order to keep the weights beyond a given value, soft (van Rossum et al., 2000;Morrison et al., 2007)

stdP-lIke PlastIcIty rule
The general one-neuron system used as the basis of this study consists of N synapses with strength  i that receive input from neurons i with continuous values x i . Each input produces an excitatory post-synaptic potential (EPSP) which is modeled by filter functions h i (see Figure 1A for an example). The output of the neuron is, thus: where μ (set to 0.01 throughout the article) is the plasticity rate and F [] and G [] are linear functionals. However, although the analytical solution that we will present later on holds for the generally formalized plasticity rule (Eq. 2), we will only consider F = 1 (where 1 is the identity) and G = d/dt. This is called differential Hebbian learning (Kosco, 1986;Klopf, 1988) and allows for the learning of temporal sequences of input events (Porr and Wörgötter, 2003). It resembles STDP (Roberts, 1999;Wörgötter and Porr, 2004). Another important setting for F and G is conventional Hebbian learning with F = G = 1. To avoid that weight changes will follow spurious random correlations one generally assumes that learning is a slow process, where inputs change much faster than weights, with d i / i << d(x i * h i )/ x i * h i , μ → 0. In spite of this separation of time scales, such a simplification makes it still possible to investigate the dynamics of each weight separately. The separation simplifies Eq. 2 and we neglect all temporal derivatives of  i on the right hand side which leads to    as G [] is linear. If we take  i as the i-th component of a vector , we write its most general form where  denotes the transposition of matrix . In the results section we will show that following closedform solution is possible t as such is quite complex. Nevertheless, we will also show that for our quasi-static assumption μ → 0 the approxi- where I is the identity matrix, is arguable.

Relation to spike-timing-dependent plasticity
In order to show that F = 1 and G = d/dt resembles STDP we use a simplified version of the system, namely a system with N = 2 (see Figure 1C), however, where only one of the synapses is plastic ( 1 ); the other stays fixed ( 2 = 1). In order to get the correct STDP protocol, we now present two pulses to this system, one at t = 0 to the input connected with the plastic synapse, hence a pre-synaptic spike, and another at t = T to the other input. As we are using a linear model and as the corresponding weight  2 is set to a fixed value of 1, a pulse at the input is also a pulse at the output, hence creating the pendant of a post-synaptic spike. The equivalence of relating differential Hebbian plasticity with STDP arises essentially from this and its implications are discussed in detail in Saudargiene et al. (2004) and Tamosiunaite et al. (2007). Recently also Clopath et al. (2010) related voltage-based Hebbian learning to STDP.
For such a network structure, it is possible to analytically calculate the shape of the STDP curve. We start with Eq. 3 where we focus on the first entry of (t):  1 (t). The second synapse stays fixed:  2 (t) =  2 .
of the method developed in this paper. The next structure is a small network with two neurons and only one plastic synapse (see Figure 2A). This network structure is equivalent to just one neuron with a plastic recurrent connection of total delay d. This is due to our assumption of linearity. The last structure we consider is a network with three neurons and two plastic synapses (see Figure 2C). Here the equivalence is to a single neuron with two plastic recurrent connections. Note, these two structures could be considered as part of a larger recurrent network, where -after some intermittent stages -activity arrives back at its origin (with unchanging synapses in between). Hence, recurrent structure like those could be seen as simple network building blocks and in the end we shortly discuss larger network structures (i.e., a ring of neurons) that have similar stability properties as those building blocks. For the input we will use input spike distributions with first moment P (expectation value) that comes in three different conditions. The firing will be without any noise and thus, periodic, with Gaussian noise using a mean value of P and a standard deviation of σ = 5 ms and, additionally, we will also apply a purely Poisson distributed firing with rate P (hence no periodic input).
As our model is linear, it can not as such produce spikes on its own. However, we could at any point introduce a threshold so that each time a pulse exceeds this threshold a spike would travel along the axon to the next neuron via its synapses. In our first order approximation we do not need this threshold, however, we could apply it afterward.

InPut Pulse tIMIng and aMPlItude
The external input to the system is a sequence of δ-function spikes, which are being transformed into PSPs by convolving the normalized filter function h with the input sequence. As the model is linear, outputs are graded pulses, which result from the weighted linear summation of all inputs at the soma of the neuron (Eq. 1). As we do not use spike thresholds, these outputs directly provide the recurrent part of the input to the system. We need to calculate the complete input time function, notable the timing and amplitude of all pulses that enter a neuron. We investigate recurrent systems with an external input of periodicity P with delays d i , i = 1,…,R (in this study we only consider R = 1 and 2 but the analysis holds for R ∈ N). To get a better intuition for recurrences in networks like those described above (Figures 2A,C) we can compare each recurrence to a modulo operation. We find all pulse timings pt by iterating the map pt c+1 = mod (pt c + d i , P) until pt pt c N s + = = 1 0 . The total number of relevant pulse times N s = P/gcd(P,d 1 ,…,d R ) depends on the greatest common divisor (gcd) of P and all R delays d i . For instance, all relevant pulse times for P = 75 and d = 60 are pt c = {0, 60, 45, 30, 15}. This gives us the timings of the pulses. In order to calculate the amplitude G of each pulse of a system with R recurrent connections, we need to go further and solve this linear system of equations: G = (I − L) −1 λ where L is a matrix that handles the delayed recurrences and λ the external input; I is the identity matrix. The details of the derivation are provided in Section "Pulse Timings and Amplitude" in the Appendix. Having calculated the pulse amplitudes we are now able to calculate the weight change.
First, we calculate the influence of a single pulse on the weight  1 . To this end we set x 2 (t) = 0 for all t and x 1 (t) = δ(t) which is a spike at t = 0. It simplifies the convolution to a temporal shift in the filter function h h t t t t h d This leaves us with a first order differential equation with following solution: Filters usually have the property that they decay to 0 after a while which turns the exponential function into 1 and results in no weight change. Thus, a single pulse or a rate produced by a single input does not have any influence on the weight when not combined with another pulse. For these pulse-pulse correlations we have to set the other input to have a pulse at t = T: x 2 (t) = δ(t − T). If we use a simplification (with all the details in Kolodziejski et al., 2008) the result of the differential Eq. 5 for very long times writes where we used here and throughout the manuscript Hence, different values of the pulse timing T lead to different ∆ values, which, when plotted in Figure 1B against T, resemble a fully symmetrical STDP curve (for more details see Porr and Wörgötter, 2003;Saudargiene et al., 2004;Tamosiunaite et al., 2007;Kolodziejski et al., 2008). The shape and symmetry of the STDP curve depends purely on the filters' shape, thus on the post-synaptic potentials' (PSPs) shapes (see Figure 1A). Thus, with different filters h for x 1 and for x 2 one could get either LTP or long-term depression (LTD)dominated curves (see Figure 2c in Porr and Wörgötter, 2003).
We used h(t) with parameters a = 0.009, β = 0.0099, and σ = 0.029. The maximum of h(t) is at t max = (log(β) − log(a))/(β − a). The width of the STDP curve measured at biological synapses is usually about T = ±50 ms (Bi and Poo, 1998;Caporale and Dan, 2008). Thus, as our filter ceases to 10 −4 at around t ≈ 1000, we define 1000 time steps as 50 ms (which corresponds to a = 0.18 ms −1 and β = 0.198 ms −1 ) to which we will refer throughout the text.
Not only the shape of the STDP curve is determined by the filter function (Eq. 8) but also the shape of the pulses. As our model is linear, each pulse could be thought of as a weighted delta function (pseudo spike) convoluted with the filter function.

network structures and InPuts
Although our method would allow investigating general network structures, we will concentrate on only three very simple, but fundamental ones.
The first structure consists of a single neuron with two plastic synapses (see Figure 1D). Although quite simple, the analytical solution of this feed-forward structure demonstrates the power Now, as we know the complete analytical solution of Eq. 3, we investigate approximations and their errors in order to judge their usefulness for further considerations. To estimate the approximation errors we can use δ functions as inputs to the system. Furthermore, we assume that all h i = h are equal. The pulses are again modeled as δ functions δ(t − t i ) for times t i which simplifies the convolution to a temporal shift in the filter function where t i and t j are the pulse timings of neuron x i and x j respectively. We will use the filter shown in Figure 1A, given by Eq. 8.
The different approximation errors are exemplified in Figure 1E. For this we are using a single pulse pair at two synapses (see Figure 1C) for which we calculate the final synaptic strength    = →∞ lim ( ) τ τ (Eq. 9). This has been performed for differential Hebbian learning, but we remark that the error is identical for Hebbian learning. For this setup, weight changes are computed in three ways: without any approximations, yielding   (Eqs. 9 and 11); using the truncated solution only, yielding   E k ,( ) ; and using the truncated solution while also expanding the exponential function, yielding   S k ,( ) . Thus, we use   and compare it to approximations   ⋅,( ) , k calculating the error as:  This is plotted in Figure 1E for different approximations against the plasticity rate μ on a log-log scale where we set the maximal value of the input filter h to 1. As approximations E,(k) and S,(k) become very similar for k > 2, only four curves are shown. We observe that the behavior of the difference-error ∆ ·,(k) follows the order of the approximation used. The error for the linear expansion approximation S,(k = 2) is slightly higher than that from its corresponding truncation approximation E,(k = 2). However, using a plasticity rate of μ = 0.001 already results in a differenceerror value of 10 −8 compared to 10 −2 when using the same value for μ and the maximum of h. Therefore, one can in most applications use even the simplest possible linear approximation S,(k = 2) to calculate the change in synaptic strength.
As this calculation has been based on two pulses at two synapses only, we need to ask how the error develops when using N synapses and complex pulse trains. For this we first consider pulse trains, which are grouped "vertically" into groups with each input firing at most once. Filters of pulses within a group will overlap but we assume that grouping is possible such that adjacent groups are spaced with a temporal distance sufficient to prevent overlap between filter responses of temporally adjacent groups. Thus we calculate  in the same way as above leading to: When using such a temporal tiling,  B ⋅,( ) k depends only on the pulse timing matrix T with elements T ij = t j − t i . Then, we get the synaptic strength after M groups by calculating the product over all groups This solution is easy to compute. Because a product of matrices results in a summation of matrix elements, the error does not increase exponentially but only linearly in M. Because of this it follows that even after 10000 pulses the error is still of an order of only 10 −4 given the example above (see Figure 1E).
Finally we estimate how the error behaves when filters overlap. This mainly happens during bursts of pulses with temporarily high spiking frequencies, which are, in general, rare events. However, using the solution which assumes independent temporal intervals instead of the time-continuous calculation (Eq. 4), only includes results analytIcal solutIon of weIght dynaMIcs Here we are going to develop a closed-form for the matrix B(t) which governs the temporal development (Eq. 4) of the weights (t) according to Eq. 3. This solution is not trivial as the matrix A(t) is also a function of time. This problem is often found in quantum mechanics and the main problem is that matrices usually do not commute. A solution exists, however, it includes an infinite series, called the Magnus series (see Magnus, 1954 for more details), with where (t 0 ) are the synaptic strengths at time t 0 , hence before plasticity, and V(t) is the solution of following equation Here the braces η ξ η ξ ξ ξ , [ Thus, Eq. 9 combined with Eq. 11, gives us analytically the time development of all weights connected to a neuron under Hebbian plasticity in the limit of small plasticity rates μ. With this we are able to calculate without simulations in principle directly the synaptic strengths of N synapses given N different pulse trains. This property will be useful to analytically calculate the fixed-point values of the weights of our small network structures.

Approximations for the analytical solution
Before we apply the solution to our network structures, we transform it into a computable form and provide error estimates. As the commutators in the infinite series in Eq. 11 are generally non-zero we truncate the series and neglect iterations above degree (k). We write the truncated solution as B E k For two synapses this is solved directly later on, most often, however this needs to be calculated by expanding the exponential function. We denote this approximation, where we neglect terms of both the exponential as well as the Magnus series greater than of order (k), with an S instead of an E, i.e., S,(k), thus B S k ( / . ( )) ). ( ( = + = = ⋅ ≤ I Σ 2 1 1 V Notice that in the limit k → ∞ the approximation transforms into the general solution (Eq. 9). This solution is computable for arbitrary input patterns.
October 2010 | Volume 4 | Article 134 | 5 Kolodziejski et al. Stability in networks with STDP number of periods P: (n) = (n + ·P) for  ≥ 1 and n large enough. Here in this study we concentrate on the analytical solution for  = 1. Therefore, in order to calculate those fixed-point weights, the following system of equations needs to be solved: with n large enough so that the weights have already converged. ˆ( , ) B B = + n P n is defined to include all relevant pulse correlations; it is, thus, similar to  B. The roots of this system represent the fixedpoints which are stable if all eigenvalues of B are negative.
In the next step we will use the E,(k = 2) simplification for B which is I + μ ∫ U A(τ)dτ with U being an interval that includes all relevant correlations. A detailed definition is given in Section "Fixed-Point Analysis" in the Appendix. With this we will arrive at an analytical solution of the weight dynamics of the recurrent structures. If we include this simplification into Eq. 12 we get Now we have to write down matrix A to get its integral. With R recurrences and a single input, A is of dimension R + 1. The first R rows and columns describe the in-and outputs to the plastic synapses and the last row and the last column the in-and output to the constant input synapse. We define its input as K, which is a sum of δ functions always at the beginning of each period. The inputs to the plastic synapses are given by the pulse amplitudes G (Eq. A.2), however, delayed by different delay times d i . Putting everything together, we have to solve following integral: an additional error of order k = 2 due to the linearity of the filter functions h. The error after matrix multiplication results in the square of the lowest term of the Magnus series (Eq. 11).
Thus, the easily computable group decomposition  B ⋅,( ) k will yield accurate enough results even for long, non-bursting pulse trains.

fIxed-PoInt analysIs
Now we are able to calculate the temporal development of the weight dynamics in our networks. Most interesting are those developments that eventually lead to a stable, thus non-diverging, weight dynamics. Here, we develop the equations with which configurations of parameters (i.e., periodicity and delays of the recurrences) can be identified that lead to stable weights.
As mentioned, some of those (P, d i ) configurations converge to stable weight configurations, others do not. In order to reveal whether a configuration is stable, hence whether the plastic weights do not diverge, the nullclines (   = 0) need to have stable crossing points in the open region ]−1, +1[ R with R being the number of recurrent connections (±1 would lead to diverging activity since a constant input is added periodically for a theoretically infinite number of times). A crossing point outside this interval is never stable because of the recurrent connections involved [(n + 1) = (n) with  > 1 always diverges].
For the calculations we place all input pulses on a discrete grid, but calculate correlations by continuous analytical integration. Hence, the weight change in the discretized version of Eq. 4 is (n) = B(n,m)·(m) where B(n,m) denotes B(t) with lower and upper boundaries m and n respectively. Weights are ( periodically or asymptotically) stable if they take the same value after a certain 10 -20

Figure 1 | Setup (A-D) and verification (e) of the approximation of eq. 4. (A)
Two pulses x i and x j at times t i and t j respectively are convolved with filter h. (B) Depending on the time difference T, differential Hebbian plasticity leads to a STDP-like weight change curve ∆ (Porr and Wörgötter, 2003;Kolodziejski et al., 2008). The maximum is at t max = (log(β) − log(a))/(β − a). (C) A simple network with only a single plastic synapse with which the relation of differential Hebbian plasticity to spike-timing-dependent plasticity is shown. (D) The used feed-forward network with two plastic synapses. (e) Here we show the degree of consistency between our general solution and the proposed approximations.
To this end we plot the difference ∆ ⋅  (see Approximations for the Analytical Solution) between the approximation and the exact solution of Eq. 3 for one input pulse pair against the plasticity rate μ on a log-log scale. Here, a filter function h with a = 0.1, β = 0.2, σ = 0.25, and max t h(t) = 1 is used (Eq. 8). The temporal difference T = t i − t j , between the two input pulses was varied over the length of the used filter functions (here between 1 and 100 steps) and error bars representing the standard deviation are given.
. τ τ τ Here we use two input neurons (N = 2), which receive inputs at t = 0 and T, respectively. Overall, this results in following solution for the weight after a pulse pair ( l im ( )) Having found the analytical solution we now determine dynamics of the temporal development by calculating the eigenvalues λ This shows clearly that the weights oscillate around 0 and that the oscillations of the weights are damped if the real parts of the eigenvalues λ i are smaller than 0. This constraint holds for almost all T values as ℜ( )= λ i cos µν  T ( ) − = 1 0 is only true for a multiple of the number 2 π. However, 0 ≤ < µν π  T holds since we assume μ to be small. Furthermore, µν  T = 0 is trivial and does not lead to any weight change at all. Hence, the weights will continuously shrink to 0.

Recurrent networks
Next, we investigate the stability of the two neuron model with a single recurrence (see Figure 2A) for which we need the previously provided solution for Eq. 13. In small panels to the right of Figure 2B we show three different temporal developments of the plastic weight. The top panel depicts a divergent weight (P = 100 ms and d = 85 ms), the middle panel a convergent weight development (P = 100 ms and d = 60 ms) and the bottom panel an oscillating weight (P = 59 ms and d = 94 ms). In the main panel of Figure 2B the weight values of all configurations for 0 < P < 100 ms and 0 < d < 100 ms (which is 0-2.0 times the characteristic STDP length of 50 ms) are plotted. In case of convergence the final weight is depicted, however, if the weight diverges, no point is plotted at all (white regions). For oscillating weights we plot the mean value. Figure 2B reveals that if (P, d) is stable, then (nP, nd) with n ≥ 1 is also stable (lines through origin). However, weights will not be the same for (P, d) and (nP, nd). The stability statement only holds for configurations in which the weight has not reached a value whose absolute value is ≥1. This immediately leads to divergence and is the reason why the lines through the origin not always start at the origin.
We skip the detailed calculations (including the definition for s ) here, given in Section "Fixed-Point Analysis" in the Appendix, and directly write the result where only correlations in the interval U have been considered and the pulse amplitudes G (Eq. A.2) were applied: Obviously, the number of fixed points is not restricted to one and there exist many configuration with more than one fixedpoint. However, we will show below that if a configuration is stable, than it does not matter how many fixed points in the interval ]−1, 1[ exist.
Here we would like to stress again that we are not interested in average quantities (no mean-field approach) but in the fine temporal dynamics of each single weight and additionally that a nondiverging weight is also a necessary condition for a non-diverging overall activity as we do not use any kind of boundaries. Now, we have all tools ready to analyze the network structures (Figures 1C and 2A,C) starting without additional noise, thereafter with noise and finally by also considering asymmetrical STDP rules.

analysIs of the network structures wIthout noIse
First, we will investigate the structures, which we introduced in the beginning of this study, in the absence of noise. The feed-forward structure always displays oscillation of weights and activity whereas the recurrent structures show more complex behavior. Depending on the network parameters, the weights either diverge, converge to a fixed-point or oscillate.

Feed-forward single neuron model
Here we look in more detail at symmetrical differential Hebbian plasticity with two plastic synapses (see Figure 1D) where the simplification E,(k = 2) is analytically fully solvable. We have also based the error analysis provided in Figure 1E on these calculations. For the approximation E,(k = 2) the matrix results in October 2010 | Volume 4 | Article 134 | 7 Kolodziejski et al.

Stability in networks with STDP
three different values of P: 0.5, 1.0, and 1.5 times the characteristic STDP time window (we leave out the 2.0 times as we do not get anything novel at higher periodicity values according to Figure 2B). Starting with P = 25 ms ( Figure 2D) it is apparent that only a few inhibitory connections evolved (only about 10%). However, when moving to higher P values (Figures 2E,F), more and more inhibitory connections show up. If we had plotted (not shown) the ratio of excitatory connections, we would find that at about P = 100 ms the ratio of excitatory connections reaches 0.5 and for all smaller values the ratio is always >0.5, thus more excitatory connections develop. We note that all properties mentioned in this paragraph remain essentially the same, when adding noise to the input. Figures 2E,F reveal that the main regions (upper-right, middle) move to higher (d 1 , d 2 ) configurations and enlarge at the same time. In addition, new regions originate at smaller (d 1 , d 2 ) configurations.
The total number of fixed points N 0 scaled by the number of the possible configurations, of which there are P 2 , also grows with increasing periodicity (see Figure 3, black dots). In the next section we will compare the number of fixed points to two different noisy input protocols.

analysIs of the network structures wIth noIse
Above we used deterministic timing at the input to our system. As this is not realistic, we will now use two different types of noise. The first one is Gaussian jitter around the given periodicity value Additionally, Figure 2B shows that excitatory and inhibitory regions always alternate and that the two regions for d < P (below diagonal) are almost of the same size as the many regions for d > P (above diagonal). Furthermore, if the first delayed pulse is closer to its input pulse (i.e., d < P/2), than the final weight is inhibitory. By contrast, if the first delayed pulse is closer to the consecutive input pulse (i.e., P/2 < d < P), than the final weight is excitatory. In the excitatory case, the pulse at the plastic synapse is followed by another pulse. This resembles a causal relationship which should lead to a weight growth. The opposite is true for the inhibitory case.
We also conducted a perturbation analysis in which we perturbed the weight for each configuration with increasing amplitude until the weight could not converge back to its fixed point. It turns out that the crucial measure for stability is whether the absolute value of the sum of the fixed-point value and the perturbation amplitude exceeds 1. This is true for almost all configurations with stable weights which is why we do not need to plot the results of the perturbation analysis.
Another property of the system is that the final weights are identical for (P, d + nP) configurations for all n ≥ 1 (see Figure 2B, in particular the black dashed line at P = 40 ms). Thus, if the weight for e.g., P = 40 ms and d = 25 ms converges to , so does the weight for P = 40 ms and d = 65 ms. We use this property to simplify the convergence plots of the network structure with three neurons or two quasi recurrences (Figures 2D-F) and only depict fixed points for delay values up to the periodicity P. We show the fixed points for or even losing their stability. However, in a few trials they still happen to survive which leads to the regions with many rarely stable configurations.
It is clearly visible that the number of stable configurations for the structure with a single recurrence substantially decreases. This is, however, not the case for the double recurrence structure. There, the number of stable configurations increases (see red border in Figures 4G,H) because of the compact structure. However, this effect starts above P = 30 ms which is revealed by Figure 3. In this figure we plot the number of stable (non-divergent) configurations N 0 normalized by the maximum possible configurations (P 2 ) against the periodicity P. Only between about 30 and 60 ms noise increases the number of stable configurations by about 20%. Above 60 ms the smaller non-compact areas that lose stability dominate over the border areas that gained stability. Above a certain periodicity the smaller areas become large enough so that noise becomes advantageous again and border regions become stable. Therefore, the effect that noise slightly increases stability reoccurs periodically. Figure 3 also shows clearly that Poisson statistics leads to slightly more stable configurations than adding Gaussian noise. Most important, however, is the fact that noise does not deteriorate the system, i.e., the number of stable configurations stays almost unchanged.
Similar to the feed-forward structure, weights sometimes change from positive to negative values (or vice versa), hence they switch from being excitatory to being inhibitory which is not known from biology. However, this could be easily solved if we exchange our linear synapse with two non-linear synapses (one that is linear in the positive regime and 0 elsewhere and the other vice versa), so that such a push-pull mechanism covers the full linear range. This is the conventional way of addressing this issue (Pollen and Ronner, 1982;Ferster, 1988;Palmer et al., 1991;Heeger, 1993;Wörgötter et al., 1998). Note, however, in general we found that without noise switching is rather rare in our recurrent structures anyways. Even with additional noise weights tend to switch only if the final weight is close to 0.

analysIs of the network structures wIth asyMMetrIc stdP rules
In biological systems, the STDP curve is not symmetrical (Bi and Poo, 1998;Caporale and Dan, 2008). In order to achieve an asymmetric differential Hebbian plasticity curve, we leave the positive part which describes the strengthening of the weight (LTP) as it is and shrink or stretch the negative part which describes the weakening (LTD). This is done by multiplying negative weight changes with a factor r.

Feed-forward structure
For the feed-forward structure we calculated the temporal development analytically (see Section "Feed-Forward Structure" in the Appendix) and find that the general behavior is not changed by the asymmetry. This becomes visible from the eigenvalues of the system: Both oscillation and damping still exist, however, damping constant and frequency of the oscillation is scaled by the asymmetry r. P with a standard deviation of 5 ms. The second one is Poisson firing with rate P. Due to the stochastic nature of the signals, analytical results cannot be obtained anymore and we have to rely on simulations.

Feed-forward structure
The general behavior of the feed-forward structure does not change in the presence of either noise. Although the oscillations are not "perfect" anymore, weights are still oscillating and eventually they decay to 0.

Recurrent structures
For the recurrent structures the situation is different as noise is able to change an initially unstable configuration to now exhibit a convergent weight development and vice versa. In Figure 4 we show the stability of the weights for different configurations. Here, red stands for highly stable weights (i.e., all trials lead to convergence) and blue for rarely stable ones (i.e., only a few trials lead to convergence). Stability in between where half of the trials lead to convergence is very rare, so we included those configurations to the rarely stable ones. In order to compare the noise results with those for no noise, we include the no noise stable regions in gray in Figure 4. Additionally, the color of each configuration that is stable both with and without noise is lighter (i.e., light red and light blue). We do not plot the actual final weight values as those are not largely different from the no noise ones.
For the two neuron structure, a few of the branches in Figures 4A,B completely disappeared, others moved slightly. The two main branches for d < T (below diagonal) remain highly stable, however, only in the "core". The more diffuse, non-compact regions, i.e., regions with detached stable configurations become either completely unstable or at least rarely stable. Noise affects the input periodicity, thus configurations that are stable for a large range of different P values survive when noise has been added. By contrast, configurations within non-compact regions are always deflected toward unstable configurations, thus, compromising  Figure 4) with respect to the periodicity P. Additionally, N 0 is normalized by the maximum possible configurations (P 2 ) and we plot the results without noise (black), with Poisson distributed firing with rate P (blue) and Gaussian distributed periodicity with mean P and a variance of 5 ms (red).
October 2010 | Volume 4 | Article 134 | 9 Kolodziejski et al. Stability in networks with STDP asymmetric STDP curve yields more stable configurations than a symmetric one. Here, a dominant LTP part is advantageous over a dominant LTD part. This situation changes when we introduce noise. An STDP curve with a dominant LTD part results now in more stable configurations than a LTP-dominant STDP curve. However, almost all asymmetric STDP curves lead to many more stable configurations as compared to a symmetric STDP curve.
Here we note again that making the input spiking obeying the Poisson statistics is advantageous over adding Gaussian noise (if comparing absolute numbers). Next, we consider the double recurrence structure ( Figure  5B) where the number of stable configurations were also counted up to P = 50 ms. Here, the no noise result looks more similar to the results with noise; two peaks, one for an LTP-dominant and one for an LTD-dominant STDP curve, are visible. However, adding noise to the system again changes the relative heights of the peaks: the STDP curve with more LTD becomes advantageous over the STDP curve with more LTP. Similar to the single

Recurrent structures
Now we come to the recurrent structures. Here, similar to the introduction of noise, the situation is different. With an additional asymmetry the number of stable configurations changes. In the following, we depict the number of stable configuration achieved with a certain asymmetry r = 2 a (logarithmic representation) by N a . Thus, the number of stable configuration for the symmetric STDP curve is still N 0 . Then, to better compare the different results, we normalize N a to the symmetric number N 0 plotting N a /N 0 in Figure 5. It is worthwhile to note that the total (not the relative) number of stable configurations without noise is reduced for any value of a for the single recurrence structure and stays almost unchanged for the double recurrence structure (see Figure 3).
First, we look in more detail at the single recurrence structure (Figure 5A). We plot the normalized number N a /N 0 of stable configurations up to a periodicity of P = 50 ms. The respective N 0 values can be found in the caption of Figure 5. The no noise results have a peak for positive and for negative asymmetries, thus an   at the level of PSPs, which are here modeled by ways of the filter function h. This way closed-form solutions can be found. Spiking could be introduced by ways of a threshold. Introducing such a threshold essentially leads to the reduction of activity (eliminating all small linearly summed signals) but the general recurrence is not altered. The first network investigated was a simple feed-forward structure (see Figure 1C) which displays oscillating weights which decay to 0. This observation is interesting and not expected as this network and its inputs represent a fully symmetrical situation. Hence, oscillation should occur, one would think, but without damping. Asymmetrical STDP curves do not alter this observation, because the asymmetry of the STDP curve only affects the frequency of the weight oscillation.
The focus of the current study, however, was to investigate recurrent network structures without averaging the activity. The two small recurrent networks studied here can be seen as fundamental network building blocks. This is due to the fact that the direct recurrent connections used here could also be replaced by a network, the output of which providing the recurrent signals. This seems to be similar to Roberts (2004), Burkitt et al. (2007), and Gilson et al. (2009), however, in those studies only average activities are considered. Roberts (2004) does not treat plasticity, Burkitt et al. (2007) and Gilson et al. (2009) do not have self-connections, and Roberts (2004) and Burkitt et al. (2007) do not have delays included. Without delays, selfconnections would be cumbersome to implement in a general way (i.e., for different delays a different number of neurons need to be interposed). However, Burkitt et al. (2007) and also Gilson et al. (2009) could easily extend their models to include self-connections.
As a side remark, almost all of the above mentioned studies only consider pair-wise STDP correlations although a method for three-pair couplings exists (Pfister and Gerstner, 2006). By contrast, differential Hebbian plasticity treats multicorrelations naturally.
In the much cited study of Bi and Poo (1998) a clear asymmetry is visible in the STDP curve (Figure 7 of Bi and Poo, 1998): the part which describes the strengthening of synaptic efficiency, LTP, is more pronounced than the LTD part. In general, asymmetrical recurrence structure, there are more stable configurations when using Poisson noise, however, with Gaussian noise the boost in stability for asymmetrical as compared to the symmetrical STDP curve is larger. If we compare the number of stable configurations between LTP-dominant and LTD-dominant STDP curves we find that without noise about 55.0% LTP-dominant configurations of the total number of configurations exists. With Gaussian noise this ratio decreases to 41.7% and with Poisson statistics to 44.0%.
The dynamics of the structures examined here (and the reasoning should extend to all possible structures) would allow only for diverging weights if we set a → ±∞. As this would correspond to an STDP curve with pure LTP or LTD respectively, the dynamics pushes the weight only in one direction (either positive or negative) which leads to no fixed point (see Eqs. 14 and 15). This is also indicated in Figure 5 where the number of configurations leading to convergence is severely reduced even for a = ±4.

dIscussIon
Real neurons often display rich, non-stationary, firing patterns by which all synaptic weights will be affected. The so far existing solutions which describe Hebbian learning, on the other hand, constrain the temporal dynamics of the system or limit plasticity to a subset of synapses. With the solution presented here we can calculate weight changes for the first time without these restrictions. This is a valuable step forward in our understanding of synaptic dynamics of general Hebbian plasticity as well as of STDP (which is a special case of the general Hebbian plasticity rule) in different networks. Specifically, we have presented the time-continuous solution for the synaptic change of general Hebbian plasticity (Eqs 9 and 11), its approximation for general spiking or continuous inputs as well as a specific solution for non-bursting pulse trains. Of practical importance is the fact that the error of the computable approximations remains small even for long pulse trains.
We think that there is in general no intuitive access, which would allow us to understand the temporal development of activity and plasticity in networks with recurrent connections. Even small structures driven by periodic activity, such as the ones investigated here, display unexpected and counterintuitive behavior. The systems investigated are linear and summation takes place stable regions. The difference between Gaussian and Poisson statistic is small, however, we found that Poisson statistics, thus the statistics most often used to describe neuronal spiking (Bair et al., 1994;Shadlen and Newsome, 1998;Dayan and Abbott, 2001), yields more stability. More importantly, noise does not deteriorate but even improves the system. The benefits of noise for neuronal signal processing has been discussed in several other studies (Fellous et al., 2003;Deco and Romo, 2008). The fact that synaptic stability benefits from noise, especially in the more complex double-recurrent structures, adds to this discussion. In our networks noise helps to reach out to nearby stable regions to exploit their stability. The temporal development of multi-synapse systems and the conditions of stability are still not well understood. Some convergence conditions have been found (see for example Hopfield, 1982;Miller and MacKay, 1994;Roberts, 2000;van Rossum et al., 2000;Kempter et al., 2001;Burkitt et al., 2007;Kolodziejski et al., 2008), however, in general the synaptic strengths of such networks will diverge or oscillate. This is undesired, because network stability is important for the formation of, e.g., stable memories or receptive fields. Using the time-continuous solution for linear Hebbian plasticity described here serves therefore as a starting point to better understand mechanisms, structures and conditions for which stable network configurations will emerge. The rich dynamics, which govern many closed-loop adaptive, network based physical systems can, thus, now be better understood and predicted, which might have substantial future influence for the guided design of network controlled systems.
An important next step to improve our understanding of the interaction of plasticity and activity is to increase the complexity of the investigated network structures. Raising the number of recurrences is straightforward as we have already presented a general method to solve such networks. It will be more challenging to include plastic synapses at more than one neuron. Figure 6A shows as a small outlook such a possible network structure where we used N = 5 "building blocks" (see Figure 2C) each receiving constant periodic input (i.e., pulses), activity from itself (selfconnection) and from the previous neuron (ring-connection). In this example we set all delays to the same value (d i = d and d is = d) and the input has periodicity P (see Figures 6C,D). Each neuron i receives its first input spike at time t = (i − 1)·P/N. Both the self-connection and the ring-connection are eligible to changes and Figures 6C,D show that, similar to our fundamental building blocks (see Figures 1D and 2A,C), configurations which lead to converging weights exist. The temporal development depicted in Figure 6B exemplifies the typical behavior found in those networks: due to its symmetry all self-connections and all ring-connections converge each to practically the same weight value. Many possible network topologies could visioned using those or similar building blocks. Hence, investigating those structures in more detail goes beyond the scope of this paper.
Also more complex input protocols must be investigated. The current study suggest that dynamically stable recurrent substructures can exist within larger networks. Similar results exist for the generation of stable activity patterns in (non-plastic) networks (Memmesheimer and Timme, 2006a,b). These studies may allow us to gradually approach a better understanding of dynamically STDP curves are found in many studies (see Caporale and Dan, 2008 for a review) where sometimes the situation is reversed (LTD more pronounced than LTP). Many theoretical studies of large-scale network development use STDP curves with stronger LTD part, because in many cases it has been found that otherwise weights will diverge (Song et al., 2000;Izhikevich et al., 2004). The current study, on the other hand, suggests that collective divergent behavior could also be prevented by more specific network designs. Here we were able to show that stability without boundaries is possible. We observe that -at least in those small linear networks -many fixed points exist for different types of STDP curves. Notably there are almost always more fixed points for asymmetrical as compared to the symmetrical STDP curve. Furthermore, the notion that an LTD dominance would be beneficial for stability does not seem to be unequivocally correct. Even in the presence of noise large numbers of fixed points are found when using LTP-dominated STDP curves.
The studies of Izhikevich et al. (2004), Burkitt et al. (2007), Morrison et al. (2007), Lubenov and Siapas (2008), Gilson et al. (2009), andClopath et al. (2010) have considered fully as well as randomly connected networks with external drive, where they need to apply soft or even hard boundaries to stabilize weights. By contrast here we have looked at small recurrent networks without averaging and without boundaries. Thus, these approaches cannot be directly compared. Our results, however, appear promising. It may well be that networks with specific topologies can be constructed (or generated by self-organization), where stability exists using STDP for a large number of input patters without boundaries or weight-normalization.
Additionally, there exists a Hebbian rule that is not using boundaries either. It has been shown that the so called BCM rule (Bienenstock et al., 1982) is a rate description-based rule which is analogous to the STDP rule (Izhikevich and Desai, 2003;Pfister and Gerstner, 2006). However, to get weights stabilized, BCM is using an additional variable that pushes the weights so that the output activity is close to a set value. This rule, besides that it relies on a rate model, uses the activity to limit the weights in an indirect way, similar to other homeostasis mechanisms (Turrigiano and Nelson, 2004). Thus, it is possible that the BCM rule among others makes the weights diverge or rather diffuse which, nevertheless, would lead to a state with a finite and stable mean activity. Such a behavior does not emerge in our model as a single diverging weight will eventually reach infinity resulting in an infinite activity. Besides, our focus here is on the fine temporal dynamics of the weights.
Furthermore, we observe that noise has an influence on the recurrent structures (see Figures 2A,C). For some input periodicities P the number of stable configurations decreases, for other settings it grows. Growth is found especially for the double-recurrent structure. A reduced number of stable configurations is mainly visible in areas with many initially stable isolated points from where such a point can be easily "nudged out" to an unstable domain by noise. Such a nudging effect, however, is also the reason why the number of stable configurations can increase. Our results show that stable regions either enlarge at their boundaries or within dense but not fully covered areas. Here initially unstable data points are "nudged in" to adjacent changing modular networks going beyond the currently still dominating unifying approaches, which mostly only consider feedforward structures or randomly connected networks with average activities. This will require intensive studies but we hope that the current contribution can provide useful analytical tools and first insights into these very difficult problems.

aPPendIx Pulse tIMIngs and aMPlItude
In the main text we only calculated the timings of the pulses. In order to determine the amplitude of each pulse of a system with R recurrent connections, we need to go further and solve a linear system of equations. Doing this we need to assume that the weights eventually converge, therefore being constant. In more detail, in each period of length P there are at most P pulses, each having different amplitude. Within each period each pulse gets multiplied by certain weights (which we assume to be constant) and is then delayed. Additionally, the input provides the neuron with another pulse of amplitude 1. Hence, if we write the pulse amplitudes in a vector G and put the multiplicative factors (i.e., the weights) into a matrix L always at the delayed position, we find the pulse amplitude G of the next period according to G LG  The solution results in equations that involve high polynomial terms and are, thus, complicated to simplify. An example for the pulse timings with P = 10, d 1 = 4 and d 2 = 6 (here in arbitrary units) can be found below: where U is defined to be a meaningful range. It is, however, simpler to restrict the pulse timings to a meaningful range S (defined later) to arrive at a comprehensive description of the equations. In this case, the integral needs to range from 0 to ∞ and we have to start with all possible pulses for the input in this range. The input K t t nP n s ( ) ( ) = − = ∞ Σ δ (concerning s see below) to the constant synapse is a sum of δ functions always at the beginning of a period. As already mentioned in the main part, the δ function δ(t − t i ) for pulse time t i simplifies the convolution to a temporal shift in the filter function h t t h d h t t i i Further on, the inputs x i to the plastic synapses are given by the pulse amplitudes G, however, delayed by different delay times d i . Thus, we have to exchange the inputs x i with the full representation of all pulses of the whole temporal development from time s · P (at which the pulse amplitudes G already converged) to infinity: x n P d k Next, we neglect all terms but one on the pre-synaptic site (index n) as spiking is repetitive anyways. Additionally, we reduce the range of the sum from infinity to a meaningful range S that is determined by the time window W of the STDP curve under consideration (e.g., 50 ms in the main part): S s s = − [ , ] with s W P =  / where  is the ceiling function.

fIxed-PoInt analysIs
In order to get the fixed points, we start with Eq. 14  This solution was used to calculate the difference ∆ ·,(k) for different values of T in Figure 1E.