Multiplexing rhythmic information by spike timing dependent plasticity.

Rhythmic activity has been associated with a wide range of cognitive processes including the encoding of sensory information, navigation, the transfer of information and others. Rhythmic activity in the brain has also been suggested to be used for multiplexing information. Multiplexing is the ability to transmit more than one signal via the same channel. Here we focus on frequency division multiplexing, in which different signals are transmitted in different frequency bands. Recent work showed that spike-timing-dependent plasticity (STDP) can facilitate the transfer of rhythmic activity downstream the information processing pathway. However, STDP has also been known to generate strong winner-take-all like competition between subgroups of correlated synaptic inputs. This competition between different rhythmicity channels, induced by STDP, may prevent the multiplexing of information. Thus, raising doubts whether STDP is consistent with the idea of multiplexing. This study explores whether STDP can facilitate the multiplexing of information across multiple frequency channels, and if so, under what conditions. We address this question in a modelling study, investigating the STDP dynamics of two populations synapsing downstream onto the same neuron in a feed-forward manner. Each population was assumed to exhibit rhythmic activity, albeit in a different frequency band. Our theory reveals that the winner-take-all like competitions between the two populations is limited, in the sense that different rhythmic populations will not necessarily fully suppress each other. Furthermore, we found that for a wide range of parameters, the network converged to a solution in which the downstream neuron responded to both rhythms. Yet, the synaptic weights themselves did not converge to a fixed point, rather remained dynamic. These findings imply that STDP can support the multiplexing of rhythmic information, and demonstrate how functionality (multiplexing of information) can be retained in the face of continuous remodeling of all the synaptic weights. The constraints on the types of STDP rules that can support multiplexing provide a natural test for our theory.

Introduction

Neuronal oscillations have been described and studied for more than a century [1-14]. Rhythmic activity in the central nervous system has been associated with: attention, learning, encoding of external stimuli, consolidation of memory and motor output [8–10, 12, 14–23]. Rhythmic activity has also been suggested to support multiplexing in the central nervous system [24, 25]. Multiplexing is the ability to transmit two or more different signals via the same channel. The two main forms of multiplexing are: (i) Time division multiplexing, in which different time slots are allocated for the transmission of the different signals, e.g., distributing information for different clients by the same server. (ii) Frequency division multiplexing, in which different signals are transmitted via different frequency bands, e.g., the allocation of different frequency bands for different radio stations.

Various forms of multiplexing have been proposed to be utilized by the brain [24-34]. Caruso et al. [32] suggested that time division multiplexing is used in the auditory system to represent different objects. They demonstrated that some neurons shift their response from one object to another in time (similar to fluctuating focus of attention) in a manner that correlated with behaviour. Frequency division multiplexing has been suggested by Teng & Peoppel [30] (they also suggested other methods as well) to be utilized by the auditory system for encoding different features of the same auditory object in the theta and gamma channels; thus, frequency division multiplexing may also be used for binding. Here, we focus on frequency division multiplexing. Our aim is to study a mechanism that can shape synaptic connectivity to facilitate frequency division multiplexing.

The transfer of even a single oscillatory signal downstream is not necessarily trivial and requires some mechanism to prevent destructive interference and maintain the rhythmic component [35]. Recently, it has been suggested that synaptic plasticity, and especially spike-timing-dependent plasticity (STDP), can provide such a mechanism. STDP can be thought of as a generalization of Hebb’s rule that neurons that “fire together wire together” [36] to the temporal domain. In STDP, the amount of potentiation and depression depends on the temporal relation between the pre- and post-synaptic firing [2, 37–44]. Luz and Shamir [35] analyzed the characteristics of the STDP rule that will enable the transfer of a single frequency channel downstream.

Multiplexing requires the transfer of more than one frequency channel. However, STDP has been shown to generate a winner-take-all like competition between subgroups of correlated pools of neurons [45-49]. Consequently, one may expect that the transfer of one frequency channel will suppress the other; thus, raising serious doubts whether STDP is consistent with multiplexing in the brain. Can STDP develop the capacity for transmitting rhythmic activity in more than one frequency band spontaneously, and facilitate multiplexing of information?

We address this question here in the framework of a modelling study. Below, we define the network architecture and the STDP learning rule. We then derive a mean-field approximation for the STDP dynamics in the limit of slow learning rate for a threshold-linear Poisson downstream neuron model. Analysing the STDP dynamics yields constraints on the STDP rules that enable multiplexing. Next, we test the generalisation of our understanding beyond the simplified analytical toy model using numerical simulations. Finally, we summarize our results and discuss how STDP can yield robustness of function in the face of constant synaptic remodelling.

Results

The pre-synaptic populations

We model a system of two excitatory populations of N neurons, each responding to a different feature of an external stimulus, Fig 1. The external stimulus is characterized by two feature variables, D1 and D2, to which populations 1 and 2 respond, respectively. The response of each population is further assumed to be rhythmic, albeit in a different frequency, representing the different features of the stimulus.

Fig 1

Network architecture.

A schematic description of the network architecture showing two pre-synaptic populations, each oscillating at a different frequency. The output of these pre-synaptic neurons, serves as a feed-forward input to a single post-synaptic neuron.

The spiking activity of neuron k ∈ {1, …N} in population η ∈ {1, 2}, (where are the spike times) is a doubly stochastic process. Given the ‘intensity’ D of feature η ∈ {1, 2} of the external stimulus, the spiking activity, ρ(t), follows an independent inhomogeneous Poisson process statistics with a mean rate (mean over the Poisson distribution given the intensity variables D1 and D2) that is given by: where ν is the angular frequency of oscillations for neurons in population η, γ is the modulation to the mean ratio of the firing rate, and ϕ is the preferred phase of the kth neuron from population η. The preferred phases are assumed to be evenly spaced on the ring. Thus, it is convenient to think of the neurons in each population as organized on a ring according to their preferred phases of firing.

As the intensity parameters, D, represent features of the external stimulus they fluctuate on a timescale which is typically longer than the characteristic timescale of the neural response. For simplicity we assumed that D1 and D2 are independent random variables with identical distributions: where 〈…〉 denotes averaging with respect to the neuronal noise and stimulus statistics. The essence of multiplexing is to enable the transmission of different information channels; hence, the assumption of independence of D1 and D2 represents fluctuations of different features. This assumption also drives the winner-take-all competition between the two populations. We further assume that the stimulus changes at a slower timescale than the neural responses; thus, for example typical values for the timescale for the stimulus are on the order of 1s, neural response will follow the stimulus within an order of 10ms and spike at about 10Hz.

As we are interested in studying multiplexing, we assume that the two populations are synapsing in a feed-forward manner onto the same downstream neuron.

The downstream neuron model

Spike time correlations are the driving force of STDP dynamics [39, 45, 46, 50, 51]. Correlated pairs are more likely to affect the firing of the downstream (post-synaptic) neuron, and as a result, to modify their synaptic connections [45, 46, 49]. To analyze the STDP dynamics we need a simplified model for the post-synaptic firing that will enable us to compute the pre-post cross-correlations, and in particular, their dependence on the synaptic weights.

Following past studies [35, 39, 48, 51–53], the post-synaptic neuron is modeled as a linear Poisson neuron with a characteristic delay d > 0. The mean firing rate of the post-synaptic neuron at time t, r(t), is given by where w is the synaptic weight of the kth neuron of population η.

Temporal correlations & order parameters

The utility of the linear neuron model is that the pre-post correlations are given as a linear combination of the correlations of the pre-synaptic populations. The cross-correlation between pre-synaptic neurons at time difference Δt is given by: where δ = 1 when η = ξ and 0 otherwise, is the Kronecker delta function.

The correlation between the jth neuron in population 1 and the post-synaptic neuron can therefore be written as in which δ(t) is the Dirac delta function, and the correlations are determined by global order parameters and , where is the mean synaptic weight and is its first Fourier component. For N ≫ 1 these parameters are defined as follows and The phase ψ is determined by the condition that is real non-negative. Note that the coupling between the two populations is only expressed through the last term of the correlation function, Eq (5).

The STDP rule

Following [46, 50, 51] we model the synaptic modification Δw following either a pre- or post-synaptic spike as: The STDP rule, Eq (8), is written as a sum of two processes: potentiation (+, increase in the synaptic weight) and depression (-, decrease). We further assume a separation of variables by writing each process as a product of a weight dependent function, f±(w), and a temporal kernel, K±(Δt). The term Δt = t − t is the time difference between pre- and post-synaptic spiking. Here we assumed, for simplicity, that all pairs of pre and post spike times contribute additively to the learning process via Eq (8). Note, however, that the temporal kernels of the STDP rule, K±(Δt) have a finite support. Here we normalized the kernels, ∫K±(Δt)dΔt = 1. The parameter λ is the learning rate. It is assumed that the learning process is slower than the neuronal spiking activity and the timescale of changes in the external stimulus. Thus, the synaptic weights are relatively fixed on timescales characterizing changes in the external stimulus and the neural response. Here, we used the synaptic weight dependent functions of the form of [46]: where α > 0 is the relative strength of depression and μ ∈ [0, 1] controls the non-linearity of the learning rule. The functions f(w)± ensure that the synaptic weights are confined to the region w ∈ [0, 1]. Gütig and colleagues [46] showed that the relevant parameter regime for the emergence of a non-trivial structure is α > 1 and small μ. Gütig and colleagues have also showed that the limit of μ = 0, termed the additive model enhances the competitive nature of STDP dynamics, whereas the limit of μ = 1, termed the linear model, greatly suppresses the competitive nature.

Empirical studies reported a large repertoire of temporal kernels for STDP rules [37, 38, 40–42, 44, 54–56]. Here we focus on two families of STDP rules: 1. A temporally asymmetric kernel [37, 41, 44, 55]. 2. A temporally symmetric kernel [38, 42, 44, 56].

For the temporally asymmetric kernel we use the exponential model, Fig 2a: where Δt = t − t, Θ(x) is the Heaviside function, and τ± is the characteristic timescale of the potentiation (+) or depression (−). We assume that τ− > τ+ as typically reported.

Fig 2

The STDP rules.

The temporal kernels, K(Δt) = K+(Δt) − K−(Δt), of the STDP rules are shown as a function of pre-post spike time difference, Δt = t − t, for (a) The asymmetric learning rule, Eq (11) with τ− = 50ms. (b) The symmetric learning rule, Eq (12) with τ+ = 20ms. Different values of τ+ in (a) and τ− in (b) are depicted by color as shown in the legend.

For the temporally symmetric learning rule we use a difference of Gaussians model, Fig 2b: where τ± is the temporal width. In this case, the order of firing is not important; only the absolute time difference. We further assume, in both models, that τ+ < τ−, as is typically reported.

STDP dynamics in the limit of slow learning

Due to noisy neuronal activities, the learning dynamics is stochastic. However, in the limit of a slow learning rate, λ → 0, the fluctuations become negligible and one can obtain deterministic dynamic equations for the (mean) synaptic weights (see [50] for a detailed derivation) where η = 1, 2 and In Methods we derive the dynamics of the global order parameters using the correlation structure induced by the rhythmic activity, Eq (5). Note that in the dynamics of the order parameters, Eq (26), does not appear explicitly in the dynamics of , and vice-versa. This results from the linearity of the post synaptic neuron model we chose, Eq (3).

The homogeneous solution, Winner-take-all and multiplexing

Fig 3 shows three example results of simulating the STDP dynamics in the limit of slow learning, Eqs (25) and (26). Panels a & b show the dynamics of the synaptic weights of both populations, color coded by their preferred phase of firing. Panel c depicts the spectrum of the downstream (post-synaptic) neuron firing. Initially, as the synaptic weights are random, the input to the downstream neuron has almost no rhythmic component. In the example of Fig 3a–3c the synaptic weights of both populations converge to a homogeneous solution, in which all the synaptic weights from input neurons of different preferred phases and different populations are the same. As a result, the input to the downstream neuron has no rhythmic component and its activity shows no peak at any non-trivial frequency. The homogeneous solution is expected to be stable for large μ [46].

Fig 3

Three examples: Homogeneous solution, winner-take-all and multiplexing.

Simulation results of the STDP dynamics in the limit of slow learning and a linear Poisson downstream neuron, Eq (13), are shown for three example cases. (a)-(c) The homogeneous solution, using μ = 0.1 and α = 1.05. (d)-(f) Winner-take-all competition, using μ = 0.001 and α = 1.1. (g)-(i) Multiplexing, using μ = 0.01 and α = 1.05. Panels a, b, d, e, g and h show the synaptic weights as a function of time. Each trace depicts the dynamics of a single synaptic weight. The synapses (traces) are differentiated by color according to the preferred phases of the corresponding pre-synaptic neurons, see legend. Panels c, f and i show the spectrogram of the downstream neuron. The horizontal dashed white lines depict the location of the rhythmic channels. In addition to the rhythmic channel (f) and rhythmic channels (i) one can identify sidebands in the spectrograms (f & i). Note that: 1) These additional bands are attenuated by about 100dB relative to the rhythmic channels. 2) They reflect the ongoing STDP dynamics; hence, they do not appear in c, and freezing the STDP abolishes these additional bands (results not shown). In all these examples the asymmetric learning rule, Eq (11), was used. The initial conditions of the synaptic weights were random. The synaptic weights at time zero were independent and identically distributed uniformly on [0, 1]. The parameters that were used in these examples are: ν1 = 11Hz, ν2 = 14Hz, N = 120, λ = 0.001, γ = 1, D = 10Hz, σ = 0.8, d = 10ms, τ+ = 20ms and τ− = 50ms.

In the example of Fig 3d–3f rhythmic activity is transferred. As can be seen from Fig 3e, the homogeneous solution is not stable and the STDP dynamics causes the development of preference in the synaptic weights of population 2 to certain phases over the others. This allows the transmission of the rhythmic signal downstream. However, the STDP dynamics induces a winner-take-all competition and population 2 fully suppress population 1; hence, only one rhythmic signal is transmitted downstream and multiplexing is not enabled. The example of Fig 3g–3i depicts the desired scenario of multiplexing. In this case the homogeneous solution is unstable and both populations develop a phase preference; thus, enabling the transmission of rhythmic activity for both signals downstream.

The above examples differ in the parameters that define their respective STDP rules. Below we aim obtain insight as to the conditions that allow multiplexing. Eq (23) describes high dimensional coupled non-linear dynamics for the synaptic weights. Studying the development of rhythmic activity is, thus, not a trivial task. To this end, we take an indirect approach. We study the stability of the homogeneous solution, for all i, in which the rhythmic activity does not propagate downstream, Fig 3a–3c. Specifically, we investigate the conditions in which the homogeneous solution is unstable, and the STDP dynamics can evolve to a solution that has the capacity to transmit rhythmic information in both channels downstream. A complete derivation of the stability analysis can be found in Methods.

The homogeneous solution

The symmetry of the STDP dynamics, Eq (23), with respect to rotation guarantees the existence of a uniform solution where ∀j ∈ {1, …N} and with η = 1, 2. Solving the fixed point equation for the homogeneous solution yields where Due to the scaling of X± with N, α is not expected to be far from 1. From symmetry :

Substituting the homogeneous solution into the post-synaptic firing rate equation, Eq (3) yields Thus, in the homogeneous solution, the post-synaptic neuron will fire at a constant rate in time and the rhythmic information will not be relayed downstream.

Stability of the homogeneous solution

Performing standard stability analysis, we consider small fluctuations around the homogeneous fixed point, w = w* + δw, and expand to first order in the fluctuations: where is the stability matrix. Analysis of the stability matrix yields four prominent eigenvalues, see Methods. The first, , represents fluctuations in the uniform direction, in which all the synapses are either potentiating or depressing together, is always stable Fig 4a. Furthermore, provides a stabilizing term in other modes of fluctuation. We distinguish two regimes according to the relative strength of depression, α. For α < α, , and the uniform solution is expected to remain stable. For α > α, , and structure may emerge for sufficiently low values of μ, Fig 4a. Thus, in order for the STDP dynamics to develop structure, whether winner-take-all of multiplexing, the relative strength of depression, α, must be sufficiently high and μ has to be low.

Fig 4

The uniform and winner-take-all eigenvalues.

(a) The uniform eigenvalue, , is shown as a function of μ, for different values of α/αc, depicted by color. (b) The competitive eigenvalue of the winner-take-all mode, λWTA, is shown as a function of μ, for different values of α/αc, depicted by color. (Inset) Enlarged section of the figure, showing the eigenvalue corresponding to the choice of parameters of Fig 8a–8d, depicted by a blue asterisk. (c) The maximal value of λWTA in the interval μ ∈ [0, 1] is shown as a function of α. Note that for α ≥ αc, is obtained on the boundary μ = 0. The eigenvalues were computed for the asymmetric STDP rule, Eq (11). Unless stated otherwise, the following parameters were used: σ = 0.6, D = 10Hz, N = 120, τ− = 50ms, τ+ = 20ms and d = 10ms.

Fig 8

STDP dynamics with a conductance based downstream neuron.

Results of two numerical simulation of STDP dynamics with a conductance based downstream neuron are presented: (a)-(d) using a downstream neuron with a linear f-I curve, and (e)-(h) using a downstream neuron with a non-linear f-I curve, see Details of numerical simulations in Methods. (a), (b), (e) and (f) The synaptic weights are shown as a function of time for population 1 (in a and e) and population 2 (in b and f). The different traces show the dynamics of different synapses colored by the preferred phase of their pre-synaptic neuron, see legend. (c) and (g) The dynamics of the order parameters: the mean, , and first Fourier component, , are shown as a function of time for both populations, see legend. (d) and (h) The dynamics of the phases, ψ1 and ψ2, is shown as a function of time in red and pink, respectively. Here we used the temporally asymmetric STDP rule, Eq (11). Additional parameters are: N1 = N2 = 120, , , α = 1.1 and μ = 0.011. The learning rate λ in the non-linear case is 5 times larger than in the linear case. For further details see Methods.

The second eigenvalue, λWTA, represents a winner-take-all like mode of fluctuations, in which synapses in one population suppress the other. Clearly, a winner-take-all competition will prevent multiplexing. Typically, for sufficiently high values of α and low values of μ the winner-take-all competitive mode will become unstable; hence, suppressing the option of multiplexing, Fig 4b and 4c (see Methods for complete analysis).

The last two prominent eigenvalues, (η = 1, 2), are due to the rhythmic modes, see Eq (36) in Methods. Thus, , represents fluctuations in the synaptic weights of population η, in which a phase preference is developed, enabling the transmission of rhythmic activity in (angular) frequency ν. Examining the rhythmic eigenvalues reveals that they are composed of a sum of two terms, see Methods. The first term is similar to λWTA and can become positive (unstable) for sufficiently large values of α and low values of μ, and is largely independent of the temporal structure of the STDP rule. The second term depends on the rhythmic activity and the temporal structure of the STDP rule. It is this second term that can enable the rhythmic modes to develop while the competitive winner-take-all mode is suppressed.

Figs 5 and 6 show the dependence of the rhythmic eigenvalues on the various parameters that govern the STDP dynamics for the temporally asymmetric, Eq (11), and the temporally symmetric, Eq (12), learning rules, respectively. The temporal structure of the STDP rule, Eqs (11) and (12), as well as the delay, d, and the modulation depth, γ, determine the frequency dependence of , whereas increasing α and decreasing μ shifts the curve up, increasing the range of unstable frequencies.

Fig 5

The rhythmic eigenvalue, , for the asymmetric learning rule, Eq (11).

(a)-(d) The rhythmic eigenvalue, , is shown as a function of frequency, , for different values of: the delay, d, in (a), relative strength of depression, α/αc in (b), μ in (c), and of the potentiation time constant, τ+, in (d)—as depicted by color. The black (11Hz) and purple (14Hz) stars in (d) show the eigenvalues with the parameters used in the simulations as shown in Fig 8a–8d. (e) and (f) The phase diagram of the system showing regions of different types of solutions in the plane of [τ+, τ−] in (e) and the plane of [μ, α] in (f), as determined by the signs of and λWTA. The abbreviations to the right of the panels are: NR—non rhythmic, R—rhythmic (multiplexing) and WTA—winner-take-all. Unless stated otherwise, the parameters used in this figure are: γ = 1, σ = 0.6, D = 10Hz, N = 120, τ− = 50ms, τ+ = 20ms, μ = 0.01, α = 1.1 and d = 10ms. In (d) μ = 0.011.

Fig 6

The rhythmic eigenvalue, , for the symmetric learning rule, Eq (12).

(a)-(d) The rhythmic eigenvalue, , is shown as a function of frequency, , for different values of: the delay, d, in (a), relative strength of depression, α/αc in (b), μ in (c), and of the potentiation time constant, τ+, in (d)—as depicted by color. The blue (11Hz) and red (14Hz) stars in (d) show the eigenvalues with the parameters used in the simulations as shown in S4a–S4d Fig. (e) and (f) The phase diagram of the system showing regions of different types of solutions in the plane of [τ+, τ−] in (e) and the plane of [μ, α] in (f), as determined by the signs of and λWTA. The abbreviations to the right of the panels are: NR—non rhythmic, R—rhythmic (multiplexing) and WTA—winner-take-all. Unless stated otherwise, the parameters used in these figures are: γ = 1, σ = 0.6, D = 10Hz, N = 120, τ− = 50ms, τ+ = 20ms, μ = 0.01, α = 1.1 and d = 10ms. In (d) α = 1.05 and μ = 0.011, in (e) μ = 0.005 was taken.

The above analysis provides intuition as to the required conditions for multiplexing. Essentially, one expects that if the competitive eigenvalue is stable, λWTA < 0, and both rhythmic eigenvalues are unstable, , then the synaptic weights will evolve towards a state that will allow the transfer of both rhythms downstream. This is summarized in the phase diagram of the system, which depicts the different types of behaviour as a function of the parameters that characterize the STDP rule, panels e & f of Figs 5 and 6. Note, however, that the computation of these phase diagrams is based on local analysis of the homogeneous solution (the stability of λWTA, and instability of ). To study the non-local behaviour, a numerical investigation is required.

Fig 7 depicts the results of a numerical simulation of the STDP dynamics with λWTA ≈ −0.003 < 0, , and . Fig 7a and 7b show the dynamics of the synaptic weights. Initially, the synaptic weights were homogeneous (up to a small noise component, see caption) and no rhythm was transmitted downstream. However, through a process of spontaneous symmetry breaking both populations developed a phase preference; thus, enabling multiplexing.

Fig 7

Multiplexing as a limit cycle solution.

Simulation results of the STDP dynamics in the limit of slow learning and a linear Poisson downstream neuron, Eq (13). (a) and (b) The synaptic weights are shown as a function of time, for populations 1 and 2 in (a) and (b), respectively. Each trace depicts the dynamics of a single synaptic weight. The synapses (traces) are differentiated by color according to the preferred phases of the pre-synaptic neurons, see legend. (c) The dynamics of the order parameters: mean, , and the magnitude of the first Fourier component, , are shown for populations 1 and 2, see legend. (d) The dynamics of the phases, ψ1 and ψ2, is shown in red and pink, respectively, as a function of time. The synaptic weights at tie zero were random, stochastically independent with identical uniform distribution on [0.45, 0.55]. The parameters used in this simulation are: N = 120, , , λ = 0.001, γ = 1, σ = 0.6, D = 10Hz, N = 120, τ− = 50ms, τ+ = 20ms, μ = 0.01, α = 1.05 and d = 10ms.

Examining the dynamics of the order parameters, one can see how the rhythmic components, and , evolve in time, rising from zero towards a fixed point value, Fig 7c. In contrast with the order parameters, and , the synaptic weights themselves do not reach a fixed point, rather they remain dynamic. How can the order parameters remain fixed while the entire synaptic population is constantly changing? The solution to this puzzle is provided by examining the dynamics of ψ, see Eq (7), the phases of the rhythmic inputs, Fig 7d. As can be seen from the figure, the phases, ψ1 and ψ2, drift on the circle with constant velocities. Thus, the STDP dynamics converge to a limit cycle solution, in which the synaptic weights profile remains fixed—relative to to its phase, w(ϕ, t) = w(ϕ − ψ(t)), while the phases, themselves drift in time, ψ(t) = ψ(0) + vt. Qualitatively similar results can be obtained for the temporally symmetric STDP rule, Eq (12), see S3 Fig in Supporting information.

Conductance based downstream neuron model

The above analysis relies on the choice of a linear neuron model, Eq (3), which resulted in the lack of explicit interaction term between the two rhythms and , see Eq (26) in Methods. Non-linearity in the response of the downstream neuron to its inputs will generate interaction between the two rhythms. This interaction may increase the competition between the two rhythms that will prevent multiplexing. How robust are our results with respect to non-linear response of the downstream neuron? This issue is addressed below by studying the STDP dynamics in a conductance based Hodgkin-Huxley type model for the downstream neuron.

We used the conductance based model of Shriki et al. [57] for the downstream neuron. This choice was motivated by the ability to control the degree of non-linearity of the neuron’s response to its inputs. Often the response of a neuron is quantified using an f-I curve, which maps the frequency (f) of the neuronal spiking response to a certain level of injected current (I). In the Shriki model [57], a strong transient potassium A-current yields a threshold linear f-I curve to a good approximation. Thus, by adjusting the strength of the A-current we can control the ‘linearity’ of the f-I curve of the downstream neuron.

Fig 8 presents the results of simulating the temporally asymmetric STDP dynamics with a conductance based downstream neuron. In Fig 8a–8d the downstream neuron is characterized by a strong A-current. In this regime, the f-I curve of the downstream neuron is well approximated by a threshold-linear function (see Fig 1 in [57]). Consequently, it is reasonable to expect that our analytical results will hold, in the limit of a slow learning rate. Indeed, even though the initial conditions of all the synapses are uniform, the uniform solution loses its stability, and a structure that shows phase preference emerges, Fig 8a and 8b. After about half an hour of simulation time, the STDP dynamics of each sub-population converges to an approximately periodic solution. The order parameters and appear to converge to a fixed point, Fig 8c, while the phases ψ1 and ψ2 continue to drift with a relatively fixed velocity, Fig 8d. For the specific choice of parameters used in this simulation, the competitive winner-take-all eigenvalue is stable, see inset Fig 4b, whereas the rhythmic eigenvalue is unstable, see stars in Fig 5d.

Fig 8e–8h depict the STDP dynamics, for a non-linear downstream neuron. To this end we used the Shriki model [57] with no A-current. As can be seen from the figure, the system converges to a dynamical solution that shows some degree of similarity with the linear neuron (Fig 8a–8d). Specifically, the system relaxes to a dynamical solution that enables the transmission of both rhythms, i.e., multiplexing. However, in the non-linear case the order parameters and do not converge to a fixed point but fluctuate around some mean value. Qualitatively similar results can be obtained for the temporally symmetric STDP rule, Eq (12), see S4 Fig in Supporting information.

Discussion

Rhythmic activity in the brain has fascinated and puzzled neuroscientists for more than a century. Nevertheless, the utility of rhythmic activity remains enigmatic. One explanation frequently put forward is that of the multiplexing of information. Our work provides some measure of support for this hypothesis from the theory of unsupervised learning.

We studied the computational implications of a microscopic learning rule, namely STDP, in the absence of a reward or a teacher signal. Previous work showed that STDP generates a strong winner-take-all like competition between subgroups of correlated neurons, thus effectively eliminating the possibility of multiplexing [46]. Our work demonstrates that rhythmic activity does not necessarily generate competition between different rhythmic signals. Moreover, we found that under a wide range of parameters STDP dynamics will develop spontaneously the capacity for multiplexing.

Not every learning rule, i.e., choice of parameters that describes the STDP update rule, will support multiplexing. This observation provides a natural test for our theory. Clearly, if multiplexing has evolved via a process of STDP, then the STDP rule must exhibit instability with respect to the rhythmic modes and stability against fluctuations in the winner-take-all direction. These constraints serve as basic predictions of our theory.

Khamechian and colleagues [28] suggested a model in which visual information from the ventral and dorsal pathways is transmitted to the prefrontal cortex by means of two well separated rhythms. Specifically, they proposed that the ventral pathway conveys information to the prefrontal cortex in the gamma band (40-70 Hz), whereas information from the dorsal pathway is relayed in the high gamma band (180-220 Hz). As the frequency bands are so segregated, does this necessarily imply two different STDP rules with different temporal characteristics are used, one for each stream? No. Interestingly, Khamechian and colleagues showed that for the dorsal stream to transmit information efficiently, the distribution of preferred phases must change from uniform to a non-uniform distribution. This is consistent with the homogeneous solution, implying that the rhythmic eigenvalue of the higher gamma is stable (as well as the competitive eigenvalue). Alternatively, the preferred phases of dorsal neurons can fluctuate from trial to trial, yielding a homogeneous correlation structure, on the timescale of plasticity. In this case, STDP will not develop a rhythmic solution regardless of the plasticity rule. Thus, to pursue our theory empirically, one has first to characterize the distribution of preferred phases of the upstream population, and especially, the stability of the relative phases over time. Secondly, one has to establish multiplexing: is there a subgroup of neurons that receives and responds-to both streams of information? Or is there a segregation of signals also at the level of the downstream population? Third, one has to characterize the STDP rule and compute its eigenvalues.

In our work we have made several simplifying assumptions to facilitate the analysis:

We limited the discussion to multiplexing of only two rhythms. Is our theory general? Can STDP support multiplexing of more than two rhythms?

We assumed the two populations are symmetric. This reduced the number of free parameters (e.g., N, D, σ and so on). Will asymmetry generate a winner-take-all competition in which the ‘stronger’ population suppresses the other; thus, preventing multiplexing?

We further assumed symmetry within each population. Specifically, we assumed that the preferred phases of the neurons in each population are evenly spaced on the ring. However, if the preferred phases of the neurons are not distributed uniformly, then even a homogeneous solution will transmit the rhythm downstream. In that case, can multiplexing emerge in a trivial manner irrespective of the STDP rule?

STDP can support multiplexing of more than one rhythm, as illustrated in S1 Text in the section Supporting information that shows the multiplexing of three rhythms. Can STDP support the multiplexing of 10, 100, 1000 rhythms? Is there a capacity limit? We believe that the nature of rhythmic activity in the brain limits the number of channels that can be efficiently multiplexed. This is due to the fact the rhythmic activity in the brain is relatively wide band and efficient multiplexing requires the different signals to be well segregated in frequency. Consequently, we believe that if frequency division multiplexing is used in the brain, then the number of multiplexed signals is small.

Question number 2 is addressed in S2 Text where we show that fine tuning of a symmetry between the two populations is not required in order to obtain multiplexing. Regarding question number 3, while it is true that a non uniform phase distribution will make the transmission of rhythmic activity easier, it will not abolish the winner-take-all like competition between the two signals. Furthermore, by selectively potentiating certain phases and depressing others, STDP has the ability to amplify the transmitted rhythmic component. This issue and the question of transmitting information about the phase of the rhythmic activity are beyond the scope of the current work and will be addressed elsewhere.

In Luz & Shamir [35], due to the underlying U(1) symmetry, the system converged to a limit cycle solution. For a single population, the order parameter, , will drift on the ring |w| = const with a constant velocity. Here, in the linear neuron model, in the limit of a slow learning one expects the system to converge to the product space of two limit cycles. As there is no reason to expect that the ratio of drift velocities of the two populations will be rational, the order parameters , will, most likely, cover the torus uniformly. Nevertheless, the reduced dynamics of each population will exhibit a limit cycle. This intuition relies on the lack of interaction between and in the dynamics of the order parameters, Eq (26). Essentially, the interaction between the two populations is mediated solely via their mean component, . However, introducing non-linearity to the response properties of the post-synaptic neuron will induce an interaction between the modes. Similarly, for any finite learning rate, λ ≠ 0, the rhythmic modes will not be orthogonal and consequently will be correlated.

A post-synaptic neuron with a non-linear f-I curve and a finite learning rate is expected to induce an interaction between the two populations. Consequently, for a finite learning rate the order parameters , will not be confined to a torus and will fluctuate around (in contrast with on) the ring. Traces for this behaviour can be seen by the fact that drift velocity in the numerical simulations is not constant (compare Fig 8d and 8h) and the global order parameters and do not converge to a fixed point, but remain to fluctuate around some mean value (compare Fig 8c and 8g). Thus, the system converges to a strange attractor around the torus. Nevertheless, in this example, the post-synaptic neuron responds to both rhythms; hence, multiplexing does not require a linear neuron.

Synaptic weights in the central nervous system are highly volatile and demonstrate high turnover rates as well as considerable size changes that correlate with the synaptic weight [58-65]. How can the brain retain functionality in the face of these considerable changes in synaptic connectivity? Our work demonstrates how functionality, in terms of retaining the ability to transmit downstream rhythmic information in several channels, can be retained even when the entire synaptic population is modified throughout its entire dynamic range. Here, robustness of function is ensured by the dynamics of the global order parameters.

Methods

STDP dynamics of the order parameters

Using the correlation structure, Eqs (5) and (14) yields where and are the Fourier transforms of the STDP kernels Note that for our specific choice of kernels, , by construction.

The dynamics of the synaptic weights can be written in terms of the order parameters, and (see Eqs (6) and (7)). In the continuum limit, Eq (13) becomes where Integrating Eq (23) over ϕ yields the dynamics of the order parameters where Note, that in Eq (26), does not appear explicitly in the dynamics of .

Analysis of the stability matrix

Below we analyze the stability matrix with respect to fluctuations around the homogeneous solution, M, Eq (19). Using Eqs (13) and (20), the fluctuations can be written as with Without loss of generality, taking η = 1, ξ = 2 yields In the homogeneous fixed point, Eq (15): where

Studying Eq (30), the stability matrix, M, has four prominent eigenvalues. Two are in the subspace of the uniform directions of the two populations, and two are in directions of the of the rhythmic modes. As the uniform modes of fluctuations, , span an invariant subspace of the stability matrix, , we can study the restricted matrix, , defined by:

The matrix has two eigenvectors, and , and the corresponding eigenvalues are

The first eigenvector represents the uniform mode of fluctuations and its eigenvalue is always negative, Fig 4a; hence, the homogeneous solution is always stable with respect to uniform fluctuations. Furthermore, serves as a stabilizing term in other modes of fluctuations. We distinguish two regimes: α < α and α > α. For α < α, , and the uniform solution is expected to remain stable. For α > α, , and structure may emerge for sufficiently low values of μ, Fig 4a.

The second eigenvalue represents a winner-take-all like mode of fluctuations, in which synapses in one population suppress the other, and will prevent multiplexing. For α > α, limλWTA = 2(α−1). For the temporally asymmetric learning rule, Eq (11), X− = 0; consequently α > 1. In the temporally symmetric difference of Gaussians STDP model, Eq (12), α > 1 if and only if τ+ < τ−, which is the typical case. For α < α in the limit of small μ the divergence of stabilizes fluctuations in this mode. In this case (α < α), λWTA reaches its maximum at an intermediate value of μ ∈ (0, 1), Fig 4b. For a small range of α < α, α ≈ α this maximum can be positive, see Fig 4c. Fig 4b depicts λWTA as a function of μ for different values of α, shown by color. Note that λWTA depends on the temporal structure of the STDP rule solely via the value of α and X−; however, its sign is independent of X−. As can seen in the figure, for α > α > 1, λWTA is a decreasing function of μ and the homogeneous solution loses stability in the competitive winner-take-all direction in the limit of small μ. Note that the competitive eigenvalue is not identical in the asymmetric and symmetric rules due to the fact that . However, as , λWTA behaves qualitatively the same in both types of STDP rules.

The rhythmic modes are eigenvectors of the stability matrix with eigenvalues The first two terms in the right hand side of Eq (36) contain the stabilizing term , and their dependence on α and μ is similar to that of λWTA; compare with Eq (35). The last term depends on the real part of the Fourier transform of the delayed STDP rule at the specific frequency of oscillations, ν. This last term can destabilize the system in a direction that will enable the propagation of rhythmic activity downstream while keeping the competitive WTA mode stable depending on the interplay between the rhythmic activity and the temporal structure of the STDP rule. For , is an increasing function of the modulation to the mean ratio γ. If in addition α > α > 1 then is an increasing function of σ as well.

In the low frequency limit, , and depends on the characteristic timescales of τ+, τ−, and d only via α. For large frequencies . In this limit the STDP dynamics loses its sensitivity to the rhythmic activity. Consequently, the resultant modulation of the synaptic weights profile, , will become negligible; hence, effectively rhythmic information will not propagate downstream even if the rhythmic eigenvalue is unstable [35]. Thus, the intermediate frequency region is the most relevant for multiplexing.

In the case of the temporally symmetric kernel, Eq (12), the value of is given by where . Fig 6a shows the rhythmic eigenvalue, , as a function of the oscillation frequency, , for different values of the delay, d as depicted by color. Since typically, τ+ < τ−, for finite ν, . Consequently, will be dominated by the potentiation term, , except for the very low frequency range of . Typical values for the delay, d, are 1-10ms, whereas typical values for the characteristic timescales for the STDP, τ±, are tens of ms. As a result, the specific value of the delay, d, does not affect the rhythmic eigenvalue much, and the system becomes unstable in the rhythmic direction for .

Increasing the relative strength of the depression, α, strengthens the stabilizing term ; however, Δf(w*) scales approximately linearly with α such that the rhythmic eigenvalue is elevated, causing the frequency range in which to widen, Fig 6b. Similarly, for α > α > 1, increasing μ strengthens and reduces the frequency range in which , Fig 6c. Decreasing the characteristic timescale of potentiation, τ+ increases the frequency region with an unstable rhythmic eigenvalue; however, when τ+ becomes comparable to the delay, d, the oscillatory nature of in ν is revealed, Fig 6d.

In the case of the temporally asymmetric kernel, Eq (11), the value of is given by with . The main difference between the temporally symmetric and the asymmetric rules is that due to the discontinuity of the asymmetric STDP kernel, decay algebraically rather than exponentially fast with ν. As a result, the phase , plays a more central role in controlling the stability of the rhythmic eigenvalue. As above, since typically, τ+ < τ−, then . Fig 5a shows the rhythmic eigenvalue, , for different values of d. The dashed vertical lines depict the frequencies at which the potentiation term changes its sign, . As can be seen from the figure, for this choice of parameters the upper cutoff of the central frequency range in which the rhythmic eigenvalue is unstable is dominated by ν*, which is governed by the delay.

The effects of parameters α and μ show similar trends as for the symmetric STDP rule. Specifically, increasing μ or decreasing α, in general, shrinks the region in which fluctuations in the rhythmic direction are unstable, Fig 5b and 5c. Increasing the characteristic timescale of potentiation, τ+ beyond that of the depression, makes the depression term, , more dominant, Fig 5d. In this case the lower frequency cutoff will be dominated by the change of sign in the depression term; i.e., by the angular frequency ν*, such that .

Details of the numerical simulations

The conductance based model of the downstream neuron

We used the conductance based model of Shriki et al. [57]. The model is fully defined and studied in [57]. Here, we briefly describe the model dynamics, as used in our numerical simulations. The dynamics of the membrane potential is given by: where the leak current is given by I = g(V−E). I and I are the sodium and potassium currents, respectively, and are given by and . The relaxation equations of the the gating variables x = h, n are dx/dt = (x∞ − x)/τ. The time independent functions x∞ = h∞, n∞, m∞ and τ are: x∞ = α/(α + β) and τ = 0.1/(α + β), with α = − 0.1(V + 30)/(exp(− 0.1(V + 30)) − 1), β = 4 exp(− (V + 55)/18), α = 0.07 exp(− (V + 44)/20), β = 1/(exp(− 0.1(V + 14)) + 1), α = − 0.01(V + 34)/(exp(− 0.1(V + 34)) − 1) and β = 0.125 exp(− (V + 44)/80).

The A-current, I, that linearizes the f-I curve is given by , where a∞ = 1/(exp(−(V + 50)/20) + 1) and db/dt = (b∞−b)/τ. The time independent function of b∞ is b∞ = 1/(exp((V + 80)/6) + 1) with the voltage independent time constant τ.

The term g is the total conductance of the pre-synaptic population η and can be written as follows Here, N is the number of neurons in population η, is the synaptic weight of the jth neuron from population η and [y]+ = max(0, y). The s spike of the jth neuron is denoted by . We used with S = 1000/N, τ = 5ms and (see [35, 46]).

The membrane capacity is c = 0.1μF/cm2. The sodium, potassium and leak conductances are , and g = 0.05mS/cm2. For the conductance based downstream neuron model with a linear f-I curve we used the following parameters: g = 20mS/cm2 and τ = 20ms. For the conductance based downstream neuron with a non-linear f-I curve we took g = 0. The reversal potentials of the ionic and synaptic currents are E = 55mV, , , E = Eexc = 0mV.

Modeling pre-synaptic activity

Pre-synaptic activities were modeled by independent inhomogeneous Poisson processes, with time dependent mean firing rate given by Eq (1), with γ = 1. Every second of simulation time D1 and D2 were independently sampled from a uniform distribution with a minimum of 7Hz and a maximum of 13Hz, D = (7+ U(0, 6))Hz. Each pre-synaptic neuron, spiked according to an approximated Bernoulli process, with a probability of p ≈ rΔt, where r is the mean firing rate (Eq (1)) and Δt = 1ms. The number of pre-synaptic neurons in each population was N = 120.

STDP

The learning rate of the simulations with a linear f-I curve is λ = 0.01, Fig 8a–8d and S4a–S4d Fig. In the non-linear cases the learning rate is λ = 0.05, Fig 8e–8h and S4e–S4h Fig. In both the linear and non-linear cases we used μ = 0.011. In order to update the synaptic weights we relied on the separation of time scales between the synaptic dynamics of Eq (40); hence, the synaptic weights were updated every 1s of simulation.

Asymmetric learning rule: The ratio of depression to potentiation is α = 1.1 and the characteristic decay times were chosen to be τ+ = 20ms and τ− = 50ms.

Symmetric learning rule: Here, we chose α = 1.05, furthermore, based on our analysis, we chose the ratio of decay times to be ∼10, τ+ = 5ms, τ− = 50ms.

Initial conditions for all neurons were uniform; i.e., w(ϕ, t = 0) = 0.5, η = 1, 2.

Multiplexing three signals.

(PDF)

Click here for additional data file.

Multiplexing asymmetric signals.

(PDF)

Click here for additional data file.

Multiplexing three signals for the temporally asymmetric STDP rule.

Simulation results of the STDP dynamics in the limit of slow learning and a linear Poisson downstream neuron, Eq (13), with three input signals. (a), (b) and (c) The synaptic weights are shown as a function of time, for populations 1,2 and 3, respectively. Each trace depicts the dynamics of a single synaptic weight. The synapses (traces) are differentiated by color according to the preferred phases of thier pre-synaptic neurons, see legend. The initial conditions of the synaptic weights were random with uniform distribution on the interval [0, 1]. The parameters used in this simulation are: N = 120, λ = 0.001, , , , D = 10Hz, σ = 0.81, γ = 0.9. We simulated the temporally asymmetric STDP rule, Eq (11), with τ− = 50ms, τ+ = 5ms, μ = 0.01, α = 1.01. The delay of the downstream neuron was d = 10ms.

(TIF)

Click here for additional data file.

Multiplexing asymmetric signals for the temporally asymmetric STDP rule.

Simulation results of the STDP dynamics in the limit of slow learning and a linear Poisson downstream neuron, Eq (13), for two asymmetric signals. (a) and (b) The synaptic weights are shown as a function of time, for populations 1 and 2, respectively. Each trace depicts the dynamics of a single synaptic weight. The synapses (traces) are differentiated by color according to the preferred phases of their pre-synaptic neurons, see legend. The initial conditions of the synaptic weights were random with uniform distribution on the interval [0, 1]. The parameters used in this simulation are N = 120, λ = 0.001, , , A1 = A2 = 10Hz, D1 = 15Hz, D2 = 10Hz. We simulated the temporally asymmetric STDP rule, Eq (11), with: τ− = 50ms, τ+ = 20ms, μ = 0.01, α = 1.05. The delay of the downstream neuron was d = 10ms.

(TIF)

Click here for additional data file.

Multiplexing as a limit cycle solution, symmetric learning rule.

Simulation results of the STDP dynamics in the limit of slow learning and a linear Poisson downstream neuron, Eq (13). (a) and (b) The synaptic weights are shown as a function of time, for populations 1 and 2 in (a) and (b), respectively. Each trace depicts the dynamics of a single synaptic weight. The synapses (traces) are differentiated by color according to the preferred phases of their pre-synaptic neurons, see legend. (c) The dynamics of the order parameters: mean, , and the magnitude of the first Fourier component, , are shown for populations 1 and 2, see legend. (d) The dynamics of the phases, ψ1 and ψ2, is shown in red and pink, respectively, as a function of time. The parameters used in this simulation are: N = 120, , , λ = 0.001, γ = 1, σ = 0.6, D = 10Hz, N = 120, τ− = 50ms, τ+ = 5ms, μ = 0.001, α = 1.05 and d = 10ms.

(TIF)

Click here for additional data file.

STDP dynamics with conductance based downstream neuron for the temporally symmetric STDP rule.

Results of two numerical simulation of STDP dynamics with a conductance based downstream neuron are presented: (a)-(d) using a downstream neuron with a linear f-I curve, and (e)-(h) using a downstream neuron with a non-linear f-I curve, see Details of numerical simulations Methods. (a), (b), (e) and (f) The synaptic weights are shown as a function of time for population 1 (in a and e) and population 2 (in b and f). The different traces show the dynamics of different synapses colored by the preferred phase of the pre-synaptic neuron, see legend. (c) and (g) The dynamics of the order parameters: the mean, , and first Fourier component, , are shown as a function of time for both populations, see legend. (d) and (h) The dynamics of the phases, ψ1 and ψ2, is shown as a function of time in red and pink, respectively. Here we used the temporally symmetric STDP rule, Eq (12). Additional parameters are: , , α = 1.05 and μ = 0.011. The learning rate λ in the non-linear case is 5 times larger than in the linear case. Further details of the numerical simulations appear in Methods.

(TIF)

Click here for additional data file.

25 Feb 2020

Dear Mr. Sherf,

Thank you very much for submitting your manuscript "Multiplexing rhythmic information by spike timing dependent plasticity" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

Note that both reviewers express confusion about precisely what is meant by multiplexing in this work, so clearly more attention is needed to clarifying and illustrating this point.  Similarly, there is agreement that the calculations could be presented more effectively, with the likely best option being to move the core calculations to Methods and make use of supplementary material opportunities to provide more details and explanations.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Jonathan Rubin

Associate Editor

PLOS Computational Biology

Lyle Graham

Deputy Editor

PLOS Computational Biology

***********************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: Multiplexing is a method by which multiple signals are combined into one transmittable signal. There are various methods for multiplexing. For example, multiplexing can be achieved using alternating time slots, a method known as time-division multiplexing. In this study, the authors demonstrate that the spike-timing-dependent-plasticity, rather than supporting a Winner-Take-All dynamics, can underlie multiplexing. I think that they mean something like time-division multiplexing but I could be wrong.

I find the results of this paper both interesting and novel and I strongly recommend that it would be published in PCB. However, I think that the presentation of the results can be improved, and that additional analysis would be useful, as discussed below.

Major comments:

1. Three models are described: a linear Poisson model, a conductance-based model whose f-I curve is approximately threshold linear and a conductance-based model whose f-I curve substantially deviates from threshold linearity. For the linear model, it is demonstrated that for a range of parameters, the homogeneous fixed point is unstable in the rhythmic direction. For the two conductance-based models, numeric simulations demonstrate dynamic instability of the synaptic weights and their convergence to something that looks like a limit cycle.

a) In order for the reader to be able to compare the linear and conductance-based models, the authors should plot the dynamics of the weights in the linear model. As discussed in lines 239-240, the stability analysis is local, and while providing useful insight it would be useful to see the results of a simulation of the (deterministic approximation) of the linear model.

b) In order to demonstrate multiplexing, the authors should show (using their existing simulations) the temporal dynamics of the total input to the output neuron. This will also clarify the kind of multiplexing that operates in this problem, and what information about the two incoming signals is maintained over time.

c) Most of the paper is devoted to the stability analysis of the fixed-point. This part is pretty technical and obscures the results. I strongly recommend that this analysis will be moved to the Methods section, leaving the space in the Results section to explaining only the outline of the stability analysis and the intuition behind the instability of the homogeneous fixed point and the emergence of phase oscillations. It is also recommended to add one or two parameter-space plots, which depict the dependence of the qualitatively different behaviors of the linear model (estimated from the fixed-point analysis) on the different parameters.

Additional comments:

1. The abstract is too technical for a journal like PCB. It should be simplified and clarified.

2. The authors discuss too briefly the possibility of multiplexing in the brain and the relevance of their work to specific experiments. I recommend discussing these issues in more details in the Introduction and Discussion sections. Specifically, I wonder if the results presented in the current manuscript can be related to recently-described olfactory (https://www.nature.com/articles/nrn3668 ) and auditory (https://www.nature.com/articles/s41467-018-05121-8) multiplexing.

3. Correct me if I am wrong but if I understand it correctly, there are three time scales in the problem: t1=1/nu, the time scale of the oscillations, t2, the time-scale of the changes in the intensity D and t3, the time-scale of synaptic changes (this is in addition to the time-scale of the STDP rule). The assumption is that t1<4. What would happen if the number of inputs N associated with each of the signals is not equal for the two signals, or alternatively, if the magnitudes of the two signals are not equal on average?

5. What would happen if the number of signals is larger than 2? In other words, what is the "capacity" of the neuron for multiplexing? To clarify, for both this and the previous point, an analytical result is not requested and it would be sufficient to provide a result that is based on a small number of numeric simulations of the linear model for a particular set of parameters.

6. Demonstrating that multiplexing is not restricted to the linear case using a non-linear conductance-based model is important, but does not add much insight. Indeed, only a single paragraph is devoted to this case in the Results section. Given the similarity between Figs. 6 and 7, I ask the authors to consider moving Fig. 7 from the main text to the Supporting Information section.

Minot comments:

1. None of the authors seem to be associated with affiliation 3.

2. Five lines before the end of the Abstract: typo, "and did".

3. Line 24: was � has been

4. Line 36: "a a"

5. Line 68-69: "The essence of multiplexing…" This is a very important sentence that should appear earlier, in the Abstract and Introduction sections.

Reviewer #2: Review is uploaded as an attachment. Just in case a text version here:

The authors discuss how the synaptic weights from two the presynaptic populations with different frequencies are changed under the STDP. They employ an intricate Fokker-Plank formalism to solve the evolution of weights under the assumption of the time-scale separation between neuronal and synaptic dynamics. Authors derive boundaries for the stability of the rhythmic and winner-take-all eigenvalues to set values of parameters for the numerical simulations of the conductance-based neuronal model for numerical simulations.

\\subsection*{Major}

The main problem with the paper is that the results are not yet strongly interpretable in neuroscience terms. The authors claim that they show how STDP \\emph{``enables the transmission of rhythmic activity downstream''}; however, they demonstrate only that the weights are not converging to any stable fixed point. The continuous change of the synaptic strength seems to signify that the mini-network is possibly not learning anything.

If authors would show some cases where this continuous weight dynamics would indeed serve to learn something, it will make the neuroscience case much stronger.

Maybe if the phases of the neurons for different frequencies would be differently distributed (more or less uniformly, possibly randomly, but correlated)? How strong the results depend on the uniform distribution of the phases?

There could be, for example, some information transfer between the population phase and firing of the output neuron. Or maybe under a structured input there could be a non-trivial stable weights configuration.

It would be very useful if the authors could explain how the dynamics they observe can be interpreted as multiplexing of information (also adding on how they would define it in this case).

Another major problem with the paper is that the computations are notoriously hard to follow. Considering that they constitute the substantial part of the paper, I would suggest explaining some transitions more in the supplementary, adding the intuitive explanation, when necessary repeating the previous papers. I believe that after studying in detail the earlier publications of the senior author (Luz and Shamir, 2016, 2014, 2012), I could get a better understanding of the mathematics involved. However, the paper should be reasonably self-contained to be generally comprehensible without a long investigation of the other manuscripts.

\\subsection*{ Intermediate}

\\item Why the Eq. 1 represents a mean rate? Average over what? I imagine that it is an equation describing the inhomogeneous rate.

\\item I think the brackets in Eq. 2b are misplaced ($\\delta_{\\eta, \\xi}$ should be outside of the brackets)

\\item The $\\delta(\\Delta t)$ is not defined in Eq.~4. It cannot be the Dirac-$\\delta$ (as there should be a value for every $\\Delta t$), so I am slightly puzzled.

\\item the figure with the STDP kernels would be helpful

\\item The Fig. 1 is very beautiful but can come much earlier in the text. It also has the same basic frequency for red and blue populations, though the claim is that they are different. And overall, the red neurons spike-trains are identical copies of the blue ones. It is easy to fix, and it would then be more consistent.

\\item For the understanding of the results concerning the stability of rhythmic mode would be interesting to add the figures for the synaptic dynamics (similar to 5 and 6) in a qualitatively different parameter regime.

\\subsection*{Minor}

\\begin{enumerate}

\\item Fig 2,3,4 are very hard to read, as the lines are very fine, and the numbers are tiny.

\\item In the caption of (c) it is written $\\mu \\in (0,1)$ but you probably wanted to write $\\mu \\in [0,1)$.

\\item Authors name and surname are in the wrong order

\\end{enumerate}

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biology data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, PLOS recommends that you deposit laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions, please see

Submitted filename: PLOC_CB_Shamir_review.pdf

Click here for additional data file.

18 Apr 2020

Submitted filename: Response to reviewers.pdf

Click here for additional data file.

20 May 2020

Dear Mr. Sherf,

Thank you very much for submitting your manuscript "Multiplexing rhythmic information by spike timing dependent plasticity" for consideration at PLOS Computational Biology. As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. The reviewers appreciated the attention to an important topic. Based on the reviews, we are likely to accept this manuscript for publication, providing that you modify the manuscript according to the review recommendations.

Please prepare and submit your revised manuscript within 30 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to all review comments, and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Thank you again for your submission to our journal. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Jonathan Rubin

Associate Editor

PLOS Computational Biology

Lyle Graham

Deputy Editor

PLOS Computational Biology

***********************

A link appears below if there are any accompanying review attachments. If you believe any reviews to be missing, please contact ploscompbiol@plos.org immediately:

[LINK]

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: The authors have successfully addressed all my comments.

Reviewer #2: I can congratulate the authors; the manuscript has improved in readability significantly. However, on this path, there are still some things to do that can make reading more enjoyable.

I agree with reviewer 1, the scientific merits of the presented work will allow it to be published in PLOS CB. Nevertheless, I suggest investing time in editing still and making the scientific results accessible to the reader.

I appreciate adding the Fig.2, panel a) probably should have Delta K on the x-axis. The figure also could be larger; there is no need for panels to be smaller than legends. For both panels, not indicated in the legend tau could be listed in the caption (tau_ for (a), tau_+ for (b))

Some things might be formulated more precisely, allowing the reader to understand what is happening. For example, what does “limited” mean in the context of: “Our theory reveals that the winner-take-all like competitions between the two populations is limited. ”

It will also help to explicitly state what “homogeneous solution” means (153) (that all the weight from input neurons of different preferred phases and different populations are the same).

Figure 3. helped a lot with the understanding of the different options; it is central for the presentation of the result and could be made larger. It could also become more informative by indicating that frequencies with non-zero powers observed in the panels c), f), and i) correspond to the frequencies of the two populations (maybe indicate particular frequencies with slim lines?).

What are the third non-zero frequency arising in the panels (f) and (i)? Does its presence in panel i) allows reconstructing the “lost” input frequency. In the caption: “Synapses (traces) are differentiated by color according to their preferred phase.” The preferred phase belongs not to the synapse but the presynaptic neuron.

Phase diagrams have a very small description font. Possibly the descriptions can be united, and the font enlarged.

I did not find an explanation of what is the function delta(x) (without indices as seen in Eq. 4 and 5). It seems to be an indicator function of 0.

Instead of Spikes/sec - Hz. On one occasion, misprinted “mses” instead of sec. Generally, as a unit, “s” typically stands for seconds.

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biology data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, PLOS recommends that you deposit laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see

26 May 2020

Submitted filename: Response to reviewers.pdf

Click here for additional data file.

29 May 2020

Dear Mr. Sherf,

We are pleased to inform you that your manuscript 'Multiplexing rhythmic information by spike timing dependent plasticity' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Jonathan Rubin

Associate Editor

PLOS Computational Biology

Lyle Graham

Deputy Editor

PLOS Computational Biology

***********************************************************

19 Jun 2020

PCOMPBIOL-D-19-02160R2

Multiplexing rhythmic information by spike timing dependent plasticity

Dear Dr Sherf,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Sarah Hammond

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol