Rhythmogenesis evolves as a consequence of long-term plasticity of inhibitory synapses

Brain rhythms are widely believed to reflect numerous cognitive processes. Changes in rhythmicity have been associated with pathological states. However, the mechanism underlying these rhythms remains unknown. Here, we present a theoretical analysis of the evolvement of rhythm generating capabilities in neuronal circuits. We tested the hypothesis that brain rhythms can be acquired via an intrinsic unsupervised learning process of activity dependent plasticity. Specifically, we focused on spike timing dependent plasticity (STDP) of inhibitory synapses. We detail how rhythmicity can develop via STDP under certain conditions that serve as a natural prediction of the hypothesis. We show how global features of the STDP rule govern and stabilize the resultant rhythmic activity. Finally, we demonstrate how rhythmicity is retained even in the face of synaptic variability. This study suggests a role for inhibitory plasticity that is beyond homeostatic processes.


Results
The neuronal response model. We explored STDP dynamics in a model of two neuronal populations with reciprocal inhibition. The spiking activity of individual neurons in each population was modelled as an inhomogeneous Poisson process with a mean firing rate that obey the following dynamics: where N i is the number of neurons in population i = 1, 2, r ix is the firing rate of neuron x in population i that receives external excitatory input I i . For simplicity we take I 1 = I 2 ≡ I. Throughout this paper, the function g(x) will be taken to be a threshold linear function of its input, = = + ⌊ ⌋ g x x x ( ) for x > 0 and 0 otherwise (see also 79 ). The term a ix represents the adaptation variable of neuron x in population i, and parameter A denotes the adaptation strength. J ix,jy ≥ 0 is the strength of the inhibitory coupling from neuron y in population j to neuron x in population i.
Parameter τ m is the membrane time constant and τ a is the time constant of the adaptation. It is assumed that adaptation is a slower process than the neural response to its input, τ a > τ m . Thus, the neuronal firing rate follows changes in its input with a time scale of τ m and then adapts its rate in response to a constant input with a time scale of τ a by decreasing its firing rate by a factor of 1 + A. We further assume, for simplicity, that the populations are relatively homogeneous. Thus, we omit the sub-indices x and y from Eqs (1)(2)(3)(4). r i represents the mean activity in population i, and J ij the mean synaptic weight from a pre-synaptic neuron in population j to a post-synaptic neuron in population i, see Eqs (21)(22)(23)(24) in Methods. In the limit of slow adaptation, τ τ ≡ → / 0 m a  , a complete analytical solution is possible; see the phase diagram section and the limit cycle calculations in Methods. Unless noted otherwise the results are given in the → 0  limit and time is measured in units of the adaptation time constant, τ a . solution. In this region the system relaxes to a limit cycle (i.e., a periodic oscillatory solution, denoted by L.C. on the phase diagram) of anti-phase oscillation, Fig. 1B. In this case, the limit cycle solution has two phases. During phase 1, population 1 is dominant and active, r 1 > 0, whereas population 2 is quiescent, r 2 = 0. However, this state is not stable. Due to the adaptation, the activity of population 1 will decrease until population 2 is released from its inhibition and will further suppress population 1. During phase 2, population 2 is dominant and population 1 is quiescent. In the limit of slow adaptation,  → 0, a complete solution for the limit cycle can be derived; see the limit cycle solution section in Methods.
We denote by T i the dominance time of population i, and by T = T 1 + T 2 the period of the oscillations; see Fig. 1B. Along the diagonal of the phase diagram, = =Ĵ J J 12 21 , the dominance times are equal, T 1 = T 2 = T/2, and the oscillation period monotonically increases from zero on the boundary of the stable Fusion solution, = J 1, to infinity on the boundary of the Rival solutions, = = = + J J J A 1 12 21 , Fig. 1C. The dominance time of population 1, T 1 , diverges to infinity on the boundary of Rival 1 state, , and similarly T 2 , diverges on the boundary of Rival 2 state; see Eqs (40 and 43) in the limit cycle solution section in Methods. Thus, the basic features of the oscillatory solution can be understood from the geometry of the phase diagram.

The correlation function.
One key factor that shapes STDP dynamics is the pre-post cross-correlation function. Because neuronal activities follow independent inhomogeneous Poisson processes statistics, the cross-correlation of different neurons is given by the product of their mean firing rates. Specifically, we are be interested in the temporal average of these correlations (see below). For a periodic solution we define Figure 1D shows the temporal average cross-correlation, Γ ij (Δ), for finite  (green and blue) and in the limit of  → 0 in black. Note that the main difference is the slight deviation in the oscillation period due to finite , which is more important at low T. A detailed derivation of the cross-correlation functions appears in Methods. To analyze the STDP dynamics it is convenient to use the following quantities: 21 12 as shown in Fig. 1E,F, respectively, as a function of the time difference, Δ, for T = 2 and different values of T 1 (differentiated by color). In general, Γ ± (Δ) are periodic functions of the time difference, Δ, with a period of T. Γ + (Δ) is a positive even function of the time difference, Δ, that is symmetric with respect to T/2, whereas Γ − (Δ) is an odd function of Δ that is anti-symmetric with respect to T/2. Importantly, on the diagonal of the phase diagram, from symmetry, one obtains that Γ − (Δ) = 0.
The STDP rule. The above analysis was carried out for fixed values of the synaptic weights, assuming that the time scales in which the synaptic weights change are much longer than the characteristic times of the neuronal population dynamics, τ m and τ a (see e.g. 24,29,32,36 ). Next we consider the effect of STDP. We assume that initially the synaptic weights are relatively weak (i.e., near the origin of the phase diagram in the Fusion state) and examine how activity dependent plasticity shapes its evolution. Hence, the STDP dynamics can be thought of as a flow on the phase diagram. We are interested in understanding how the features of the STDP rule shape this flow. In particular, we aim to elucidate when this flow leads the system into the limit cycle region. Following Luz and Shamir 36 we write the STDP rule as the sum of two processes, potentiation and depression, where ΔJ is the synaptic weight difference associated with pre and post spikes with a time difference of Δt = t post − t pre . The functions K ± (t) are the temporal kernels for the potentiation (+) and depression (−) of the STDP rule, respectively, and α is the relative strength of the depression. Parameter λ is the learning rate. We assume that the learning process occurs on a slower time scale than the adaptation. Specifically, here we focus on the family of temporally a-symmetric exponential learning rules: where Θ(x) is the Heaviside step function, and τ ± denote the characteristic time scales of the potentiation (+) and depression (−) branches of the rule. The parameter H = ±1 governs the nature of the learning rule, with H = 1 for a "Hebbian" rule (i.e., potentiating at the causal branch, when the post fires after pre, Δt > 0), and H = −1 for the "Anti-Hebbian" STDP rule. Below we analyze the mean field approximation in the limit of λ → 0.  post/pre . In the limit of slow learning, λ → 0, the right hand side of Eq. (10) can be replaced by its temporal mean, yielding (see also 24,29,32,36 ), ij ij In regions of the phase diagram where a stable fixed point solution exists, i.e., = ⁎ r t r ( ) i i , the correlation function is given by the product of the time independent means, Γ = ⁎ ⁎ t r r ( ) 1 2 , and one obtains that =   J J 12 21 . As the firing rates are non-negative and the temporal kernels of the potentiation and depression, K ± , have an integral of one, the sign of  J is determined by 1 − α. As a corollary, the synaptic weights will flow towards the region of the limit cycle solution from initial conditions close to the origin in the phase diagram if α < 1. This result holds for any choice of temporal structure for the STDP rule. In particular it is independent of the Hebbianity (the value of H) of the STDP rule. Note that a similar condition (α < 1) was assumed to be the biologically relevant choice for inhibitory plasticity in 35 . Thus, initial conditions of weak synaptic coefficients (J ij close to the origin) will flow towards the region of the limit cycle solution and will enter it near the diagonal, J 21 = J 12 .
Order parameters of the STDP dynamics. In the region of the limit cycle the STDP dynamics do not necessarily flow in parallel to the identity line, but rather depend on the specific limit cycle solution and on the temporal structure of the STDP rule. It is convenient to formulate the STDP dynamics in terms of the mean and relative synaptic weights, For Γ ± see Eqs (6 and 7) and Fig. 1E,F. On the diagonal, J 12 = J 21 , due to the symmetry of the limit cycle solution Γ 12 (t) = Γ 21 (t), and as a result = −  J 0. The mean correlation, Γ + , on the other hand, is a positive even function of time with a period of T. Near the boundary of stable Fusion, the oscillation frequency diverges, T → 0. In this limit (for  → 0) the limit cycle solution for the neuronal responses will approach a square wave solution (with 50% duty cycle on J 12 = J 21 ) transitioning between 0 and 2I/(2 + A) in anti-phase. The mean correlation function, Γ + (Δ), will approach a triangular wave starting at 0 for Δ = 0 and peaking at 2I 2 /(2 + A) 2 for Δ = T/2. Consequently, for T → 0, the integral on the right hand side of Eq. (14) will be dominated by the DC component of Γ + , yielding (2 ) 2 in this limit. Hence, the same condition that allows the STDP dynamics to enter the limit cycle region from the Fusion region will also cause it flow in the positive J + direction after entering the Limit cycle region.
STDP dynamics along the diagonal. Equation (14) provides two non-linear equations for J + and for J − that are coupled in a non trivial manner via the dependence of the correlations on the synaptic weights. However, on the diagonal of the phase diagram the situation is simplified: since = −  J 0 the problem is reduced to a one dimensional flow. To analyze the dynamics of J + on the diagonal it is convenient to write it as the sum of two terms: ,pot/dep / Figure 2A,B show +  J ,pot and +  J ,dep , respectively, on the diagonal as a function of the oscillation period, T (note that T is a function of Ĵ , see e.g. The dynamics of J + along the diagonal are determined by the weighted sum of both +  J ,pot and α − +  J ,dep . +  J will be positive for α < 1 for small T -near the crossing from the Fusion region. For τ + < τ − and 1 > α > α c (τ + , τ − ) (see Methods), +  J will change its sign at T*; thus, the fixed point (note = −  J 0 on the diagonal) at T* will be stable along the J + direction. This scenario is illustrated in Fig. 2D that shows +  J on the diagonal as a function of T (for different values of A, depicted by color). Interestingly, for this choice of exponential kernels for the STDP rule, the fixed point does not depend on the adaptation strength, A. The oscillation period at the fixed point, T*, is zero for α = 1 and diverges as α approaches a critical value α c (τ + , τ − ), Fig. 2E,F, see subsection Calculation of α c in Methods. For fixed α ≤ 1 and τ + , T* is minimal for τ − → ∞, increases monotonically as τ − decreases and will diverge for a critical value τ −,c < τ + such that α c (τ + , τ − ) = α. For τ − < τ −,c (and α ≥ 1) there will be no fixed point along the diagonal and the STDP dynamics along the diagonal will flow outside of the limit cycle region.
STDP dynamics away from the diagonal. The stability of the STDP fixed point requires stability in the J − direction as well. On the diagonal = −  J 0. A small perturbation in the direction of J − will affect J − dynamics via the cross-correlation term Γ − (Δ), Eq. (14). The cross-correlations depend on the synaptic weight via the dominance times, T 1 and T 2 . Hence, for a small perturbation around the diagonal, The geometry of the phase diagram ( Fig. 1A) reveals that increasing (decreasing) J − results in advancing towards the Rival 1 (Rival 2) region, and consequently increasing T 1 (T 2 ) and (decreasing) Consequently, a fixed point (on the diagonal) that is stable in the J − direction for Hebbian plasticity will be unstable for Anti-Hebbian plasticity and vice versa. Figure 4 shows the flow induced by the STDP on the phase diagram for the (A) Hebbian and (B) Anti-Hebbian learning rules. As can be seen, the Anti-Hebbian learning rule is unable to converge to a state that allows oscillatory activity. In contrast, the Hebbian STDP generates symmetric (T 1 = T 2 ) anti-phase oscillatory activity in which the oscillation period is determined and controlled by the relative strength of the depression, α. This specific learning rule provides robustness with respect to the strength of adaptation, A. Fluctuations in A do not affect the period of the oscillation.

Discussion
We examined whether rhythmic activity can emerge via an unsupervised learning process of STDP. Our main result is that under a wide range of parameters, rhythmicity can develop via STDP. Specifically, we found that to develop the capacity for rhythmic activity, the STDP rule must obey the following conditions (i) a bias towards potentiation, α < 1, will lead the system into the oscillatory region of the phase diagram, (ii) a longer characteristic time for depression than for potentiation, τ − > τ + , will enable the existence of a fixed point on the diagonal that can be governed by the exact value of α, and (iii) the stability of the fixed point in the orthogonal direction is governed by the 'Hebbianity' of the plasticity rule.
Using the framework of a simplified toy-model, Magnasco and colleagues studied the computational implications of neuronal plasticity in recurrent networks 80 . It was claimed that Anti-Hebbian plasticity rules drive the network into a state of near criticality. This raises the question of why we did not find traces of criticality in our model. The explanation has to do with the differences between our models. The most significant is the temporal structure of the learning rule, which together with the neuronal cross-correlations is the driving force of STDP dynamics. In their work, Magnasco and colleagues used an instantaneous plasticity rule. Consequently, only correlations at a zero time difference, which are inherently symmetric, contributed to their learning dynamics. As a result, only the symmetric part of their synaptic connectivity pattern was affected by the dynamics. In contrast, in our model both J + and J − , which denote the symmetric and anti-symmetric parts of the connectivity patternrespectively, evolve with time. When one allows for a non-trivial STDP rule, much richer dynamical behaviors can develop 60 . Nevertheless, it is important to emphasize we do not claim that that Hebbian but not Anti-Hebbian plasticity will induce rhythmogenesis. We found that due to inherent symmetry if the Hebbian STDP fails to yield rhythmogenesis then the Anti-Hebbian can, and vice-versa.
Control of rhythmic activity. STDP may also provide a mechanism for selecting and stabilizing oscillations; for example, the oscillation frequency can be governed and manipulated by the relative strength of the depression, α, or changes in the time constants of the STDP rule, τ ± , see Fig. 2F. Disruption of the STDP rule may result in changes to the learned oscillation frequency.
Simplifying assumptions. The analysis of STDP dynamics in recurrent networks is challenging. To facilitate the analysis we used the framework of a simplified model for the neuronal responses and made several simplifying assumptions. We assumed a separation of three time scales τ τ λ −   m a 1 . The separation of the neuronal time constant from that of the adaptation enabled us to obtain an analytic expression for the temporal correlations that drive the STDP dynamics. The assumption that long term synaptic plasticity occurs on a longer time scale allowed us to consider STDP dynamics as a flow on the phase diagram. Numerous studies have employed phase diagram description to depict the possible dynamical states of the network as a function of various parameters. Our approach to STDP dynamics adds another layer to this description. Figure 5 shows numerical solutions for the STDP dynamics. The vector field depicts our analytic solution to the STDP dynamics using exact expressions for the cross-correlations (see calculation of the cross-correlation function in Methods) in the limit of  = 0. The red, green and blue traces show the results of numerically simulating the STDP dynamic with a neuronal model with a small but finite  = .
0 001. As can be seen from the figure, the numerical results with  = . ); mainly no population is ever fully suppressed. This affects the cross-correlation function, which in turn will modify the flow along the phase diagram. As a result, the STDP dynamics will converge to a different fixed point. Nevertheless, the STDP dynamics still converges to a state of anti-phase oscillations, Fig. 5C. Thus, although quantitatively the results are different, similar qualitative behavior is obtained. The limit of small  enabled us to obtain complete analytical expressions for the cross-correlations.
The interplay between short and long term plasticity processes deserves consideration. Oscillations would not be possible in this model without short term plasticity; here, adaptation. Thus, short term plasticity plays a major role in shaping the temporal structure of the neuronal cross-correlations, Γ ij (t) that drive the STDP dynamics, which in turn, may or may not converge to a state that allows this oscillatory behavior. It is interesting to note that short term plasticity, specifically the value of A, affects and shapes the phase diagram. Decreasing the value of A to zero, for example, will shrink the region of oscillatory activity to zero and rhythmic activity will no longer be possible.
The reflection of the flow on the phase diagram with respect to the diagonal when reflecting the STDP rule with respect to time stems from the inherent symmetry of the cross-correlation function which drives the dynamics (Γ ij (Δ) = Γ ji (−Δ)); hence, it is general and holds regardless of the choice of model. Certain other assumptions can easily be relaxed. For example, we assumed symmetry between the two competing populations. However, using the (threshold) linearity of our model one can easily rescale the neuronal responses to allow for different inputs and adaptation strengths. On the other hand, the independence of the fixed point, T*, on the adaptation strength, A, is specific to this model and for the choice of an exponentially decaying STDP rule.
Choice of network architecture. A central assumption in this study was the choice of (a reciprocal inhibition) architecture. Because it is not possible to analyze a system with an undefined architecture some choice had to be made. The specific choice of architecture was made to obtain a model that could be fully analyzed. However, the choice of architecture (including the short-term-plasticity mechanism) shapes the phase diagram, allows for the different regions of dynamical solutions (fixed points, In/Out of/Anti -phase oscillations, etc.) and determines the cross-correlations. Additionally, under certain conditions, propagation delays may also have a major effect on the computational outcome of the STDP dynamics 81 . For example, The dotted line in Fig. 5A shows the numerical results of simulating the STDP dynamics in a neuronal model that also incorporates within population inhibition, see section neuronal dynamics with local inhibition term in Methods. The STDP dynamics with local inhibition follows a different flow on the phase diagram and converges to a different point. Nevertheless, the STDP dynamics converges to a state of anti-phase oscillations, Fig. 5E. Incorporating a local inhibition term will not only modify the flow on the phase diagram, but will also change the phase diagram itself. Consequently, the effect of the network architecture on STDP dynamics should not be underestimated. Because this effect is highly non-linear, one cannot generalize these results to other architectures in a straightforward manner. Nevertheless, the approach delineated here; namely, studying the induced flow on the phase diagram of the system, can be applied to other models in the limit of a slow learning rate.
Robustness to synaptic variability. Yet another central simplifying assumption we made throughout our analysis was that the population-mean firing rates and the mean synaptic weights were representative of single neuron firing rates and single synaptic weights. This allowed us to explore a model with reduced dimensionality (the phase diagram was analyzed in 2 dimensions instead of 2N 1 N 2 dimensions) and only to study the dynamics of the global order parameters such as the effective (or mean) couplings between the two populations. However, it is not at all obvious that the mean synaptic weight is indeed a good representative of the synaptic weight distribution, since for instance, the neuronal populations may break down into sub-clusters. Figure 5A shows that even when individual synapses are free to potentiate and depress independently, the mean weights follow the predicted flow on the phase diagram. Moreover, the mean synaptic weights remain representative of the individual weights that are distributed around them, Fig. 5B. Nevertheless, the distribution around the mean weights is not trivial. Consequently, different neurons may receive different levels of inhibition and miss-tuning of the oscillation frequency might occur. Figure 5C shows that the firing rates of different neurons (vertically shifted) in the two populations (differentiated by color) are identical in spite of the synaptic weights distribution. Moreover, even though synaptic weight distribution is different in the three examples shown (differentiated by color), the oscillation period is almost identical. Thus, functionality, in terms of obtaining a specific oscillation frequency, is retained even in the face of synaptic variability. What is the source of this remarkable outcome? We believe that this results from the fact that the STDP dynamics (e.g. Eq. (11)) only depend on the synaptic weights via the cross-correlations, which in turn, are determined by the oscillation period and dominance times. Thus, the fixed point of the STDP dynamics itself is determined by the oscillation period due to the activity dependence of the plasticity rule. On the other hand, to obtain this rhythmic activity it was also essential to have an architecture with two distinct inhibitory populations.

Methods
Phase diagram and limit cycle calculations. The fixed points of the dynamics. From Eqs (1-4) we obtain the dynamics of the mean firing rates in each population  where δ ≡ − ⁎ x x x , yielding the four eigenvalues for the stability matrix: , The sum of the pair of eigenvalues λ + ± , 1 2 and their product is  + + > A J ; hence, these eigenvalues are always stable. On the other hand, for the pair of eigenvalues λ − ± , 1 2 the sum is  + − + J (1 ), which is negative if and only if inhibition is sufficiently weak, < + J 1  (in that case their product will also be positive, assuming  is small). Thus, the Fusion state loses its stability when reciprocal inhibition becomes sufficiently strong,  > + J 1 .
The limit cycle solution. In the region of the phase diagram where no stable fixed point exists the network dynamics relaxes to anti-phase oscillations. Below we provide a detailed solution for the limit cycle in the limit of . The limit cycle is solved using the anti-phase oscillations ansatz. First the neuronal dynamics is solved for each phase, where the dynamics are linear. This provides a piecewise solution with several parameters to be determined. Then we apply two sets of constraints: periodicity and transition.
Assuming the anti-phase oscillations ansatz we separate the cycle into two phases. During phase-1 population 1 is dominant and fully suppresses population 2, for times t ∈ (0, T 1 ). In the limit of slow adaptation, → 0  , dynamics during phase-1 are given by: where we measure time in units of τ a . Eqs (29)(30)(31)(32) can be easily solved, yielding Similarly, during phase-2, when population 2 is dominant and fully suppresses population 1, = ′ + ∈ t t T 1 T T T ( , ) 1 1 2 Continuity of the adaptation variables, a i , dictates that, for example, the initial conditions of Eq. (36), a 2 (T 1 ), will be given from Eq.
We now need to determine four parameters: a 1 (0), a 2 (0), T 1 and T 2 . These parameters are determined by two sets of constraints. One is periodicity, namely i i 1 2 yielding, The second set of constraints is given by the transition conditions. Specifically, the transition time from phase-1 to phase-2 at T 1 is not arbitrary; rather, T 1 is a special point in time in which population 2 is released from being fully suppressed, such that the net input to population 2 changes its sign from negative to positive; thus, , which obeys J 12 J 21 → 1; hence, the limit of the zero oscillation period is obtained on the boundary of stable Fusion (note that these calculations were done for → 0  ).
, dominance times are equal, T 1 = T 2 = T/2, Consequently, the oscillation period, T, increases monotonically along the diagonal of the phase-diagram from zero at the transition to Fusion ( = J 1) to infinity at the transition to the Rival states ( = + J A 1 ).

Calculation of the cross-correlation function.
Calculation of the (temporally averaged) crosscorrelation function, Eq. (5), is done using the analytical solution for the neuronal responses in the limit of slow adaptation, → 0  . These correlations arise from co-fluctuation of the firing rates of the neurons and affect the STDP dynamics via their overlap with the STDP rule; thus, the relevant timescales are determined by the temporal structure of the STDP rule, Eq. (9). When the system relaxes to a fixed point solution, 1, 2), the cross-correlations are constant in time, Thus, correlations will be zero in the Rival states; hence, there will be no STDP. In the Fusion state the cross-correlations will be symmetric, Γ 12 (t) = Γ 21 (t). As a result, the STDP dynamics for J 12 and J 21 will be identical and the flow will be in the uniform direction, parallel to the diagonal line.
At the Limit cycle we use the analytical solution, Eqs (33)(34)(35)(36)(37)(38)(39)(40), to calculate the cross-correlations in a straightforward manner. For Δ ∈ T T [0, min { , }] 1 2 we obtain Along the diagonal, on the edge of the stable Fusion state region, T → 0, the cross-correlation will resemble a triangular chainsaw function (in the → 0  limit) with period T and peak 2I 2 /(2 + A) 2 . Consequently, as T goes to zero, the overlap between the cross-correlation function and the STDP rule will be governed by the DC component, yielding The above expressions for the cross-correlations were given in terms of the dominance times, {T i } instead of the effective couplings J ij . The translation to the synaptic weights from the dominance times is possible by Eq. (43). However, because we were interested in studying the ability to learn and stabilize a specific oscillatory activity, it was more convenient to think about the dynamics in terms of the dominance times. Similarly, to consider stability with respect to the J − direction we utilized the derivative of Γ − = Γ 21 − Γ 12 with respect to T − = T 1 − T 2 . On the diagonal, = ≡ T T T A T 3 [1 ] Calculation of α c . On the diagonal T 1 = T 2 = T/2, in the limit of slow oscillations, T → ∞, one obtains Hence, if α is less than a critical value α c = N(τ + )/N(τ − ), then +  J will always be positive (along the diagonal). On the other hand, if α is larger than α c then +  J will be negative for sufficiently large T, and a fixed point will exist if α < 1.
Neuronal dynamics with a local inhibition term. In Fig. 5