Synaptic motility and functional stability in the whisker cortex

The high motility of synaptic weights raises the question of how the brain can retain its functionality in the face of constant synaptic remodeling. Here we used the whisker system of rats and mice to study the interplay between synaptic plasticity (motility) and the transmission of sensory signals downstream. Rats and mice probe their surroundings by rhythmically moving their whiskers back and forth. The azimuthal position of a whisker can be estimated from the activity of whisking neurons that respond selectively to a preferred phase along the whisking cycle. These preferred phases are widely distributed on the ring. However, simple models for the transmission of the whisking signal downstream predict a distribution of preferred phases that is an order of magnitude narrower than empirically observed. Here, we suggest that synaptic plasticity in the form of spike-timing-dependent plasticity (STDP) may provide a solution to this conundrum. This hypothesis is addressed in the framework of a modeling study that investigated the STDP dynamics in a population of synapses that propagates the whisking signal downstream. The findings showed that for a wide range of parameters, STDP dynamics do not relax to a fixed point. As a result, the preferred phases of downstream neurons drift in time at a non-uniform velocity which in turn, induces a non-uniform distribution of the preferred phases of the downstream population. This demonstrates how functionality, in terms of the distribution of preferred phases, can be retained not simply despite, but because of the constant synaptic motility. Our analysis leads to several key empirical predictions to test this hypothesis.

The preferred phases of whisking neurons are widely distributed on the ring in a non-uniform manner, Fig. 1(g).The distribution of preferred phases in specific brain regions can be approximated by a von Mises (circular normal) distribution, Pr(φ) = e κ cos(φ−ψ) 2πI 0 (κ) where I 0 (κ) is the modified Bessel function of order 0.
The parameters ψ and κ quantify the mean and width (or sharpness) of the distribution.Typically, the value of κ is around 1, whereas the mean, ψ, depends on the brain region [8,13].
The distribution of preferred phases poses a challenge to our understanding of the propagation of sensory signals in the central nervous system.To delve into the source of this conundrum, we consider a simple model for the transmission of whisking signals from layer 4 downstream to layer 2/3 (L2/3) in the whisker cortex.In this simplistic model, we ignore the effect of recurrent connectivity within L2/3, as well as the possible contribution of direct thalamic input Fig. 1(g).Excitatory neurons in layer 4 of the whisker cortex almost do not respond to whisking, whereas about 30% of the inhibitory neurons in layer 4 (L4I) are whisking neurons [22,23].Thus, the main source of the whisking signal to layer 2/3 originates from L4I neurons.
We model the activity of the downstream L2/3 neuron by a delayed linear response to its inputs (see [24]).This higher level of abstraction (compared to the frequently used linear non-linear Poisson model [25][26][27][28][29][30][31]) is applied to facilitate analysis, and is expressed as where I ex denotes an excitatory drive that is independent of the whisking phase, w k is the synaptic weight of the kth L4I neuron and d > 0 is the delay.Parameter ρ k denotes the activity of neuron k of the population of L4I neurons that serve as feed-forward input to the downstream L2/3 neuron.We further assume that the preferred phases of the upstream population, {φ k } N k=1 , are i.i.d.according to Eq. ( 1) with mean ψ L4I and width κ L4I .Thus, the preferred phase of the downstream neuron is determined by a weighted average of the preferred phases of its inputs, weighted by their synaptic strengths.
However, synaptic weights are highly volatile [33][34][35][36].Fluctuations of about 50% in the synaptic weights were shown to occur over a period of several days [35,36].Additionally, activity-independent fluctuations were reported to be on the same order of magnitude [36].If synaptic plasticity is purely activity-independent, the synaptic weights, {w k } N k=1 , will be random and independent of the preferred phases of the upstream L4I neurons.Thus, the preferred phase of a downstream L2/3 neuron is determined by random pulling of the preferred phases of the N L4I neurons that serve as its input.As a result, the width of the distribution of the preferred phases of L2/3 whisking neurons will decay to zero as the number of pooled L4I neurons, N , grows, Fig. 1(f).Even pooling responses of as few as N = 20 L4I neurons results in κ L2/3 ≈ 10, which is an order of magnitude larger than typically observed [37].Raising the question: what mechanism can generate this distribution?
Recently, we studied the transmission of whisking signals from excitatory neurons in the thalamus downstream to layer L4I neurons.We showed that activity-dependent plasticity in the form of spike-timing-dependent plasticity (STDP, as detailed in section II B) can provide a mechanism that generates a non-trivial distribution of preferred phases governed by the STDP rule [37].However, this work only considered the plasticity of excitatory synapses and capitalized on the strong positive feedback inherent to excitatory STDP [38].In contrast, inhibitory STDP is typically characterized by a negative feedback mechanism [24,39].Thus, it remains unclear whether STDP can account for the distribution of the preferred phases of L2/3 neurons.
Here we examined the hypothesis that STDP may underlie the unexplained distribution of preferred phases in L2/3.We analyzed under which conditions inhibitory STDP of L4I-to-L2/3 synapses can generate a non-trivial distribution of preferred phases, and investigated how the STDP rule shapes the resultant distribution.Below, we first define the network architecture and the response statistics of L4I neurons to whisking, in section II A. The STDP model is defined in section II B. Next, in section II C we study the STDP dynamics of a single synapse in the limit of weak coupling.This will serve as a 'free theory' that will later shed light on the STDP dynamics of a population of synapses.Then, in section II D we analyze the STDP dynamics in the special case of an isotropic L4I population, i.e., with a uniform distribution of preferred phases, κ L4I = 0. Subsequently, we turn to investigate the effect of κ L4I > 0 and study how the STDP rule shapes the resultant distribution of phases in L2/3.Finally, we summarize our results, discuss the limitations of this study, and propose several empirical predictions to test our hypothesis.

A. Network model and order parameters
We study STDP dynamics in a purely feed-forward architecture Fig. 1(g).Each downstream L2/3 neuron receives input from N L4I whisking neurons and a constant excitatory input.The excitatory input is assumed to be non-rhythmic and independent of the whisking phase.To facilitate the analysis, the downstream neurons are assumed to be delayed linear neurons following Eq.( 2).
The spiking activity of the L4I whisking neurons is modeled by independent inhomogeneous Poisson processes with mean instantaneous firing rates (intensity), given by where ) is the spike train of neuron i with t i,k denoting the time of the kth spike of the L4I neuron i with preferred phase φ i .Parameter D is the mean firing rate during whisking (in the air, averaged over one cycle), and γ is the modulation depth.The phase along the whisking cycle is taken to be θ(t) = νt with ν = ν/(2π) = 1/T w denoting the whisking frequency.The angular brackets • • • denote averaging with respect to the statistics of the noisy neuronal responses.Under these assumptions, the firing rate of the downstream neuron will also oscillate at the whisking frequency, ν, with where from Eqs. ( 2) and (3) the mean firing rate, amplitude modulation and preferred phase of the downstream neuron are determined by global order parameters that characterize the synaptic weights profile: where w is the mean synaptic weight, and we iψ is its 'population vector' with R ∋ w ≥ 0. Thus yielding: The STDP model STDP can be thought of as a microscopic unsupervised learning rule in which the synaptic weight is modified according to the temporal relation of pre-and post-spike times.A wide range of STDP rules has been observed empirically [40][41][42][43][44][45][46][47][48].For example, Bi and Poo reported a temporally asymmetric STPD rule in which an excitatory synapse was potentiated (strengthened) when the firing was causal, i.e., the pre-synaptic neuron fired about 20ms before post, and depressed (weakened) when the firing was acausal [41].The STDP of inhibitory synapses has also been reported.Woodin and colleagues described a temporally symmetric rule [48], whereas Haass and colleagues reported (in a different brain region) a temporally asymmetric STDP rule for inhibitory synapses [49].
We write the STDP rule as a sum of two processes, potentiation, and depression (see also [24,28,37,[50][51][52][53]: where ∆w is the change in the synaptic weight w, λ is the learning rate, and ∆t = t post − t pre is the time difference between the pre-and post-synaptic spike times.The first term on the right-hand side of Eq. ( 8) represents the potentiation (+) and the second term the depression (−).For mathematical convenience, we assumed a separation of variables, and thus write each process (potentiation and depression) as a product of a weight-dependent function f ± (w) and a temporal kernel K ± (∆t).
As in Gütig et al [50], we used the following choice for the synaptic dependence for the STDP rule where µ ∈ [0, 1] controls the non-linearity of the learning rule.
In our numerical analysis, we used Gaussian kernels for the temporal dependence of the STDP rule: where τ ± and T ± are the widths and centers of the temporal kernels, respectively.Specifically, we focus on two families of learning rules.One, as in Woodin et al [48], is a temporally symmetric difference of Gaussians 'Mexican hat' learning rule, in which T + = T − = 0. We shall term the upright Mexican hat rule, τ + < τ − , Hebbian, and the inverted Mexican hat rule, τ + > τ − , anti-Hebbian.
In the limit of slow learning, λ → 0, STDP dynamics effectively averages the noise in the neuronal activity.The fluctuations in the synaptic weights become negligible and the stochastic dynamics of the synaptic weights can then be replaced with deterministic dynamics for their means (see [28,37,50,53,55,56], yielding where and is the temporal average of the cross-correlation between the jth L4I pre-synaptic neuron and the L2/3 postsynaptic neuron, averaged over one period of the whisking cycle, T w see section IV A.

C. STDP dynamics of a single synapse
We begin by analyzing the simple case in which only a single L4I-to-L2/3 synapse is plastic.Although highly artificial, this case constitutes a 'free theory' that will prove pivotal later.In the limit of large N , the activity of the downstream neuron will be almost unaffected by the activity of the single plastic neuron.In this case, the firing rate of the downstream neuron will oscillate at the whisking frequency, ν, with a fixed preferred phase, φ post , determined by its non-plastic inputs, following Eqs.( 4) and (7).
In the limit of weak coupling, the cross-correlation between the plastic pre-synaptic neuron and the postsynaptic cell can be approximated by the product of their mean responses, yielding for large N ) where φ = φ pre − φ post = φ pre − (π + ψ + dν) is the phase difference between the pre-and post-post synaptic neurons.Substituting Eq. ( 15) into Eqs.( 12) and ( 13) one obtains that the synaptic weight converges to a single fixed point where η = γpreγpost 2 , and K± and K± e iΩ± are the Fourier transforms of the STDP kernels: Note that in our specific choice of kernels, K± = 1, by construction.
Part of the utility of this simplified scenario is that it enables a complete understanding of how different parameters affect the resultant fixed-point, Eq. ( 16).In the limit of weak coupling, the STDP dynamics of a single inhibitory synapse is similar to that of an excitatory synapse, which was analyzed in [28].In brief, the synapse is potentiated or depressed depending on the phase difference between the pre-and the post-synaptic neurons, φ.The parameter µ governs the smoothness of the transition between potentiated and depressed synapses.In the additive rule, i.e., µ = 0, the synapses are either potentiated to 1 or depressed to 0 with a sharp transition.As the value of µ increases, the transition becomes smoother, Figs.2(a)-2(b).
The temporal structure of the STDP rule determines which phases will be potentiated and which will be depressed.It is convenient to define the width of the profile as the range of values of φ such that w * (φ) ≥ 1/2, and critical values of φ by the condition with the definitions Thus, the width of the synaptic weights profile is always π [57].For a temporally symmetric learning rule Ω ± = 0, yielding a symmetric profile, φ post − φ c = ±π/2 .A Hebbian learning rule will potentiate synapses with a phase difference that is small in absolute values, |φ| < π/2.Whereas, an anti-Hebbian temporally symmetric rule will potentiate synapses with |φ| > π/2, Fig. 2(a).
For the temporally asymmetric delta rule, in the limit of τ ± → 0, K± = 1 and

D. STDP in an isotropic model
We now turn to study the STDP dynamics of a population of N L4I-to-L2/3 inhibitory synapses, Fig. 1(g).First, we consider the case of an isotropic population, κ L4I = 0, taking the preferred phases of the L4I neurons to be evenly spaced on the ring, φ k = 2πk/N .The biological scenario of κ L4I > 0 is addressed in section II E.

The uniform fixed point
In the isotropic case, κ L4I = 0, a uniform solution, where w k = w * , ∀k ∈ {1 . . .N } always exists.In the uniform solution w = w * and w = 0. Consequently, γ post = 0 and the whisking signal is not transmitted downstream.Thus, the uniform solution prevents the transmission of the whisking signal.Below we study the stability of the uniform solution in the limit of large N ; see section IV C for details.In this limit, there are two types of uniform fixed-point solutions: In type 1, potentiation and depression cancel each other out via the weight-dependence of the STDP rule, . In type 2, the mean inhibitory input to the downstream neuron balances the excitatory input, I ex − Dw * 2 = 0. Along the uniform direction, the stable fixed point is given by min{w * 1 , w * 2 }.To further study the stability of the fixed points we expand the dynamics to first order in fluctuations δw(φ j ) = w(φ j ) − w * , yielding where M is the stability matrix.
The stability matrix has two prominent eigenvalues.One, m u , is in the uniform direction, w.The other, m w , is in the 'whisking' direction, w (with degeneracy 2).The eigenvalues of the stability matrix around w * 1 in this limit were calculated and yielded Thus, the uniform fixed point w * 1 is stable with respect to fluctuations in the uniform direction, m u < 0, for I ex > D/2.In the type 1 fixed point, m u provides a stabilizing term against fluctuations in the whisking direction, m w , Eq. ( 26).This stabilizing term can become arbitrarily small by decreasing µ; thus enabling fluctuations in the whisking direction to develop and destabilize the homogeneous fixed point.For I ex > D/2 there is a critical value of µ, such that for any µ < µ crit , the uniform solution is unstable in the whisking direction.Since µ ≥ 0, a critical value exists if and only if cos(α 0 ) > 0.
An analysis of the stability of type 2 (balanced) uniform solution, w * 2 , in the limit of large N , indicates that the two prominent eigenvalues are given by Type 2 fixed-point is stable with respect to uniform fluctuations for I ex < D/2.In contrast to the type 1 fixedpoint, here m u does not provide a stabilizing term against fluctuations in other directions -as expected from the balanced solution; see [24].Similar to type 1 fixed-point, the balanced fixed-point tends to lose stability in the whisking direction as cos(α 0 ) increases.For both types of uniform fixed points, in the limit of the additive learning rule, µ → 0, the uniform fixed point is unstable for cos(α 0 ) > 0. Thus, the temporal structure of the learning rule governs the stability of the uniform fixed point via a single parameter, α 0 , as defined in Eq. (21).
For the temporally asymmetric rule, we consider the limit of τ ± → 0, in which the STDP kernels, K ± (t) converge to a Dirac delta.In this case K cos ).Thus, for example, in a Hebbian temporally asymmetric rule T − < 0 and T + > 0, for Thus overall, the uniform solution does not allow for the transmission of the whisking signal downstream.However, for a wide range of parameters, this solution is unstable.In this case, one expects that a non-uniform synaptic weights profile will develop and facilitate the transmission of the whisking signal downstream.

The limit cycle solution
Figure 3 shows the STDP dynamics when the uniform fixed-point is unstable.As shown in the figure, the synaptic weights do not converge to a non-uniform fixed point.Rather, the dynamics converge to a limit cycle, in which all the synaptic weights vary across their entire dynamic range, Figs.3(a)-3(c).Nevertheless, the whisking signal is transmitted downstream, and its power, D post γ post , remains steady.This results from the fact that the global order parameters w and w (see Eqs. ( 5) and ( 6)) converge to a stable fixed-point, Fig. 3(d).In contrast, the phase of the downstream neuron, φ post = ψ(t) + νd − π drifts in time with a constant drift velocity v drift , Fig. 3(e).
Can a non-uniform fixed-point exist?Let us assume that the synaptic weights converged to a non-uniform fixed-point.In this case, the downstream neuron would fire rhythmically at a fixed rate, modulation depth and phase according to Eqs. ( 4) and (7).Hence, the system is under the terms of section STDP for a single synapse, as detailed in section II C. As a result, the synaptic weights must obey the 'free theory', Eq. ( 16), that determines the order parameters: w, w, and ψ.The order parameters of the 'free theory' dictate the activity of the downstream neuron, via D post , γ post , and φ post (see Eq. ( 7)), which, in turn, must be consistent with the initially assumed activity.This self-consistency argument is illustrated in Figs.3(g)-3(i).Assuming the system converged to a non-uniform fixed point and that w.l.g.φ (ass) post = 0, the synaptic weight profile, given by Eq. ( 16), is depicted by the black line.The phase of the rhythmic input to the downstream neuron, ψ + π, is shown by the dashed red line.Consequently, the preferred phase of the downstream neuron is given by φ (s.c.) post = ψ + νd + π (green dashed line).In Fig. 3(g In contrast, in Figs.3(h)-3(i) the self-consistency condition is not satisfied.In Fig. 3(h Thus, a non-uniform fixed-point (for κ L4I = 0) is a special case in which v drift = 0.In the generic case v drift = 0. We find that, for small µ, the drift velocity can be approximated by (see section IV D 3): compare the solid line, Eq. ( 30), and open squares (numerical) in Fig. 3(f).Thus, the drift velocity is determined by the STDP rule; however, only via K and α 0 .
In particular, the sign of v drift is determined solely by α 0 in accordance with the self-consistency argument above.
E. STDP dynamics for κL4I > 0 Figure 4(a) shows the STDP dynamics in the case of an un-isotropic upstream population, κ L4I > 0. As in the isotropic case, for a wide range of parameters, the STDP dynamics converge to a limit cycle.In contrast with the isotropic case, the drift velocity is not constant and depends on the phase along the whisking cycle, Figs.4(b)-4(c).Consequently, the time the downstream neuron spends at each phase is not constant and is proportional to the inverse of the drift velocity in that phase, Fig. 4(f).Thus, STDP dynamics induces a distribution over time for the preferred phase of single L2/3 neurons, which in turn is translated into a distribution over the downstream population, Figs.4(d)-4(e).
Figure 4(h) shows the color-coded distribution of preferred phases of the downstream population, for different values of κ L4I .In the limit of an isotropic upstream population, κ L4I = 0, the drift velocity is constant and the induced distribution of preferred phases in layer 2/3 will be uniform.For small κ L4I , the peak of the distribution of the preferred phase in L2/3 will shift by either π/2+dν or −π/2 + dν, relative to L4I, depending on the STDP rule.As κ L4I grows, the distribution of preferred phases in the upstream layer (L4I) converges to a delta function, the STDP dynamics will converge to a fixed-point and distribution of preferred phases in L2/3 will converge to a delta function at dν, as shown in inset Fig .4(h).

III. DISCUSSION
STDP can be viewed as a microscopic (unsupervised) learning rule.Typically, one pictures learning as a process that converges to an optimum.Thus, STDP dynamics is expected to relax to a fixed point or to a manifold attractor.However, empirical findings show that synaptic weights do not converge to a fixed point.Rather, synaptic weights remain volatile and vary by ∼ 50% over a period of several days [35,36].Moreover, activityindependent plasticity on a similar order of magnitude has also been reported [36].
Thus, while activity-dependent plasticity (e.g.STDP) is considered to be an organizing force that pulls the system towards an attractor, activity-independent plasticity is thought to introduce noise into the learning dynamics.This noise will cause the synaptic weights to fluctuate around a point attractor or induce a random walk on a manifold attractor.Consequently, it has been argued that activity-independent plasticity can generate 'drifting representations' [58][59][60][61][62][63][64].That is, the neural representation of an external stimulus will not be fixed in time.
An example of a drifting representation was described by Driscoll and colleagues [58].They showed that while the distribution of preferred stimuli across a population remained stable, the preferred stimuli of individual neurons varied over a period of several days.
Here we showed that activity-dependent plasticity can also provide a mechanism that generates drifting repre- sentations.Moreover, this mechanism induces a distribution over time for the preferred stimuli of single neurons.This distribution is then translated into a distribution over the ensemble of neurons.
Our theory leads to several empirical predictions.The first and the most natural prediction is that preferred phases are dynamic and have a consistent drift velocity that is tightly related to the distribution of the preferred phases.In other words, drift velocity is expected to be minimal (maximal) at the most (least) common preferred phase.
Second, according to our theory, the distribution of preferred phases in the downstream population is governed by the STDP rule.A wide distribution, κ L2/3 ≈ 1, implies that the fixed point solution of the STDP dynamics is unstable.Furthermore, the STDP rule also determines whether the mean phase of the downstream population will be advanced or delayed relative to the upstream.These observations yield quantitatively testable predictions on the temporal structure of the STDP rule of L4I-to-L2/3 synapses.
Third, according to our theory, the distribution of the downstream layer is shifted and broadened due to STDP, when compared to the random pooling model, Fig. 1(f).Consequently 'turning off' STDP will cause the mean phase at the downstream population to shift towards that of the upstream, and the distribution is expected to become narrower.We do not yet know how to manipulate the STDP rule.However, STDP dynamics are driven by the product of the STDP rule and the crosscorrelations of the neural activity.The cross-correlations can be modified by altering the sensory input to the system, e.g., by inducing artificial whisking.Such manipulations have been performed in the past, illustrating activity-dependent plasticity in the whisker system [65].Thus, inducing artificial whisking at high frequencies is expected to weaken the STDP; see e.g.Fig. 2(d), and consequently cause the distribution to shift its center and become narrower.
Drifting representation has implications to our understanding of the neural code.Consider a linear readout that estimates the phase along the whisking cycle from the responses of L2/3 neurons.A linear estimator, for example, can be a weighted sum of the neural responses.An optimal linear estimator is a linear estimator with a choice of weights that maximizes the signal-to-noise ratio.Typically, for a linear readout the signal is embedded in the co-variation of the stimulus (whisking phase) and the neural responses [66,67].However, if the neural representation of the stimulus is constantly changing, so is the 'signal'.As a result, the weights of the optimal linear estimator will have to adapt constantly.
One alternative to consider is that the readout may not be optimal (at least not optimized to estimate the whisking phase) [67].As the distribution of preferred phases is not uniform, a linear readout that pools the neural responses with random weights will carry information on the whisking phase.Even though the accuracy of a random pooling readout is inferior to that of the optimal linear estimator, it is robust to drifting representation.
In our work, we made several simplifying assumptions.We ignored the possible contribution of recurrent connectivity within L2/3, which may 1) affect the pre-post cross-correlations, and 2) be plastic itself.In addition, the pre-post correlations were governed solely by rhythmic activity; however, a strong collective signal such as touch may also affect the STDP dynamics.Studying these effects was beyond the scope of the current work and will be addressed elsewhere.
Nevertheless, learning in the central nervous system has been addressed empirically on two separate levels.On the microscopic level, the STDP rule has been investigated in preparations that lacked the functionality of the system.On the macroscopic level, functionality and behavior were studied in preparations that do not enable the investigation of the synaptic learning rules.The theory developed here bridges this gap, and draws direct links between these two levels that can thus serve as testable predictions for our hypothesis.

IV. SUPPLEMENTARY MATERIALS A. The cross-correlation
The cross-correlation between the pre-and postsynaptic firing is the driving force of STDP dynamics.In our model, the pre-post correlations are governed by the correlations within the upstream population.The crosscorrelation between neurons j and k of the pre-synaptic population obey Using Eq. ( 2), the cross-correlation between the jth L4I neuron and the downstream L2/3 neuron is given by where the order parameters are defined in Eqs. ( 5) and (6).In the limit of N ≫ 1 the cross-correlation can be written as follows where D post , γ post , and φ post are given in Eq. (7).

C. Stability of the homogeneous fixed-points
In the homogeneous case, κ = 0, the STDP dynamics in the limit of large N , Eq. ( 36), has two solutions, Eqs.(22) and (23).
To examine the stability of the first solution, w * 1 , we consider small perturbations around the fixed point w where yielding with the notation ∆f = f − − f + and Recalling that ∆f (w * 1 ) = 0, the eigenvalue corresponding to uniform fluctuations is given by Note that in this case the uniform eigenvalue does not depend on the learning rule.
The eigenvalue in the whisking direction is where K and α 0 were defined in Eq. ( 21).The corresponding eigenvalues around the balanced fixed point, Eq. ( 23), obey and

D. Derivation of drift velocity in the additive model
The phase of the synaptic weights profile, ψ, is given by Eq. ( 6) via the condition 0 ≤ w ∈ R, namely: where we take the continuum limit.Note that Pr(φ) is the probability density of the preferred phases in the upstream (L4I) layer, Eq. ( 1).
Taking the temporal derivative of Eq. ( 45) yields the dynamic equation for the phase, ψ, To utilize Eq. ( 46) we need a better understanding of the shape of the moving profile of synaptic weights, w(φ).Below, we first study the synaptic weights profile in the limit of small µ, and derive a semi-phenomenological model for the moving profile.Next, we use a selfconsistent argument to determine the parameters that characterize this profile.Then, we use Eq. ( 46) to compute the drift velocity.
1.The semi-phenomenological model for small µ Fig. 3(h) shows an example of a moving profile with a positive drift velocity and small µ.The synaptic weights profile consists of four distinct regions.Two are saturated regions, in which the synapses are either fully depressed or fully potentiated.The other two are fronts; a preceding front and a receding front.
In the saturated regions, w(φ) = 0 or w(φ) = 1, and due to the weight dependence of the STDP rule, f ± (w), ẇ(φ) = 0.As a result, the saturated regions do not contribute to the drift velocity; see Eq. ( 46).In the interior of the fronts w(φ) ∈ (0, 1), and in the limit of the additive model, µ → 0, f + (w) = f − (w) = 1.Consequently, in the fronts: The edges of the fronts, φ c 1/2 , are given by the condition of the vanishing temporal derivative, ẇ(φ) = 0, in Eq. (47), yielding (compare to Eq. ( 20)) To compute the drift velocity we use the following semi-phenomenological model for the synaptic weights profile (see Fig. ( 5)): Note that Eq. ( 49) must be taken with a grain of salt, as φ is an angular variable.Figures 3(h)-3(i) illustrates the agreement between our semi-phenomenological model and a moving profile obtained from simulating the STDP dynamics.The parameter ∆ denotes the width of the front.
The synaptic weights profile, in the limit of small κ, is fully defined by φ c1 , φ c2 and ∆.The profile determines ψ via Eq.( 6), which in turn must be consistent with φ c1 and φ c2 as determined by Eq. ( 48).This self-consistent constraint determines the width of the fronts, ∆.
Below we determine the width of the front, ∆, and use it to compute the drift velocity.First, in the homogeneous case, κ = 0.Then, we extend our results to a nonuniform distribution of preferred phases for 0 < κ ≪ 1, in a perturbative manner.

Uniform distribution
In the case of a uniform distribution, κ = 0, instead using Eq. ( 45) to determine ψ, we demand that ψ is the 'center of mass', ψ com , of the synaptic weights profile: Using Eq. ( 49), the denominator of Eq. ( 50) The numerator of Eq. ( 50) gives yielding from the self-consistency constrain, Eq. ( 50), Using the fact that the saturated regions do not contribute to the drift velocity, Eq. ( 46) yields substituting the value of ∆, Eq. ( 53), yields Eq. ( 30): drift does not depend on ψ.The superscript 0 denotes that this is the zero-order in κ.Below we now extend the calculation of the drift velocity in a perturbative manner in κ.

Nonuniform distribution in the limit of small κ
To study the leading non-trivial order of the drift velocity, v drift , in κ, we write where v (n) drift and ∆ (n) (n ∈ N ∪ {0}) are independent of κ.The drift velocity, v drift , can be estimated from Eq. ( 46), once the synaptic weights profile w(φ) is established, given the width of the front, ∆.
where for 0 < κ ≪ 1, we approximated the probability density of preferred phases, Eq. (1), by In Eqs. ( 59) and ( 60) we explicitly denote the dependence of the synaptic weight profile on the width of the front.Note that as Eq. ( 60) is O(κ), it is sufficient to take only the zero order in κ for the width.

E. Numerical simulations & pre-synaptic phase distributions
Scripts of the numerical simulations were written in Matlab.Figure 2 was produced with Mathematica.Unless stated otherwise, the numerical results presented in this paper were obtained by solving Eq. ( 36) with the Euler method.The cftool with the Nonlinear Least Squares method and a confidence level of 95% for all parameters was used for fitting the data presented in Figs.4(d)-4(e).
Unless stated otherwise, we simulated the STDP dynamics in the mean field limit without quenched disorder.To this end, the preferred phase, φ k , of the kth neuron in a population of N pre-synaptic L4I neurons was set by the condition φ k −π Pr(ϕ)dϕ = k/N .

FIG. 1 :
FIG.1: Transmission of the whisking signal.(a) In the lemniscal pathway[10,11,14] the whisking signal is relayed downstream via the trigeminal ganglion (TG) to the trigeminal nuclei (TN) of the brainstem down to the ventral posterior medial (VPM) nucleus of the thalamus.From the thalamus, the whisking signal in the lemniscal pathway is transmitted downstream to layer 4 of the cortex (L4I -mainly inhibitory neurons), and then down to layer 2/3 (L2/3)[12][13][14][15].The polar histograms represent the distributions of preferred phases along the pathway, based on[7,13,16,18] -which are presented purely for purposes of illustration.Note, that the whisking signal is also transmitted via the paralemniscal pathway[10,18,32].(b)-(e) encoding of the whisking signal by a whisking neuron -illustration.(b) The angular position of a whisker, β, is shown as a function of time.The angle is often modeled as β(t) = β mid (t) + ∆β(t) cos(φ(t)), where β mid (t) and ∆β(t) are the midpoint and the whisking amplitude, respectively, that change slowly in time.(c) The whisking phase φ as a function of time is φ(t) = (νt) mod2π , where ν is the angular frequency of the whisking.(d) Raster plot.The phases in which a single model whisking neuron spiked are marked (black dots) for different whisking cycles (on the y-axis)).The preferred phase of the model neuron at φ = π/2.(e) Tuning curve.The normalized mean firing rate of the model neuron in (d) is shown as a function of the phase along the whisking cycle.(f) The expected width in the downstream layer, κ L2/3 , in the random pooling model, is depicted as a function of the number of pooled L4I neurons, N , for random pooling.The width of the distribution, κ L2/3 , was estimated from 10000 repetitions of drawing N preferred phases of the upstream population with κ L4I = 1.(g) Feed-forward architecture of the L4I-L2/3 section of the pathway.

1 FIG. 3 :
FIG. 3: STDP dynamics for κ = 0. (a) and (b) The synaptic weights profile, w(φ, t), is shown as a function of the preferred phases of the upstream neurons, φ, at different times (differentiated by color).(c) Trace depicts synaptic weights are shown as a function of time.The synaptic weights are differentiated by color according to the preferred phase of the upstream neuron.(d) and (e) Dynamics of the order parameters w (red) & w (black) are shown in d, and ψ in e. (f) Drift velocity is shown as a function of T − .The black line follows Eq. (30), and red squares depict the estimation of the drift velocity from the numerical simulations.(g)-(i) Synaptic weights profile and the self-consistency argument.The black lines depict the profile assuming zero-drift, i.e., according to Eq. (30).The vertical dashed grey lines (aligning with dashed black lines) show φ c1,2 .The input to the downstream neuron, ψ 0 + π, and the phase of the downstream neuron, ψ 0 + π + νd, assuming zero drift (i.e., the black profile), are shown by the dashed red and green lines, respectively.The parabolic approximation of the semi-phenomenological model, Eq. (49), is depicted by the dashed magenta lines.For comparison, the synaptic weights profile obtained by numerical simulation of the STDP dynamics is shown in blue as a function of φ − φ post .The parameters used to generate the figure are as follows.In (a) -(e), the symmetric anti-Hebbian rule was used with random uniformly distributed initial conditions for synaptic weights in the interval [0.3, 0.7], and λ = 10 −3 , D = 10hz, I = 10hz, τ − = 20ms, τ + = 50ms, T ± = 0ms, µ = 10 −4 , ν = 7hz, d = 5ms and γ = 1.In (f) -(i) the asymmetric Delta rule was used with λ = 0.01, D = 10hz, I = 6hz, T + = 36ms, µ = 0, ν = 20hz, d = 12ms, γ = 1, and T − = 37, 42, 32ms in (g), (h), and (i), respectively.Initial conditions were w(φ, t = 0) = 0.5 + 0.3 cos(φ).In all simulations, N = 150 was used.

4 FIG. 4 :
FIG. 4: STDP for κ > 0. (a) Synaptic weight dynamics.Each trace depicts the temporal evolution of a single synaptic weight, differentiated by color according to its preferred phase.(b)-(c) Dynamics of the order parameters: w (red) and w (black) are shown in b, and ψ is shown in c. (d)-(e) Polar histograms of the distributions of preferred phases of the L4I neurons (in d) with κ L4I = 0.6 and ψ L4I = 0.25π and the L2/3 neuron (in e) with κ L2/3 ≈ 1.2 and ψ L2/3 ≈ 2.3 rad.(f) Histogram depicting the distribution of preferred phases of the downstream, L2/3, neuron over time.The characteristics of the distribution are identical to those in panel (e).The red dots show the inverse of the drift velocity of the L2/3 neuron, normalized such that the non-normalized data are multiplied by the ratio of the maximum value of the histogram to the maximum value of the non-normalized data.Each red dot is the average of 200 data points.(h) The distribution of the preferred phases of the downstream, L2/3, neurons, ψ L2/3 , is shown in color code as a function of κ L4I .In each column, the un-normalized distribution: Pr(ψ)/max{Pr(ψ)}, is depicted.The inset shows the distribution width downstream, in layer 2/3, κ L2/3 , as a function of the width of the distribution upstream, κ L4I (blue).The identity line is shown (dashed black) for comparison.Here we used the following parameters: N = 150, ν = ν/(2π) = 10hz, D L4I = 10hz, I ex = 8hz, and d = 14ms.The temporally symmetric STDP rule, Eq. (11), was used with τ − = 20ms, τ + = 50ms, T ± = 0ms, µ = 0.001, and λ = 10 −3 .Initial conditions were random with a uniform distribution in the interval [0.3, 0.7].All the data presented here were sampled from simulations with an integer number of cycles.For better visibility, the data presented in (a)-(c) shows a non-integer number of cycles.

FIG. 5 :
FIG.5: Synaptic weights profile.The semi-phenomenological model for the moving profile of synaptic weights as a function of φ, as described in Eq.(49).