Frequency cluster formation and slow oscillations in neural populations with plasticity

We report the phenomenon of frequency clustering in a network of Hodgkin-Huxley neurons with spike timing-dependent plasticity. The clustering leads to a splitting of a neural population into a few groups synchronized at different frequencies. In this regime, the amplitude of the mean field undergoes low-frequency modulations, which may contribute to the mechanism of the emergence of slow oscillations of neural activity observed in spectral power of local field potentials or electroencephalographic signals at high frequencies. In addition to numerical simulations of such multi-clusters, we investigate the mechanisms of the observed phenomena using the simplest case of two clusters. In particular, we propose a phenomenological model which describes the dynamics of two clusters taking into account the adaptation of coupling weights. We also determine the set of plasticity functions (update rules), which lead to multi-clustering.


Introduction
Clustering of the dynamics and coupling is observed at several scales of the brain structure and function. For example, in the data measured by the functional magnetic resonance imaging (fMRI), the brain networks form functional clusters that can be seen in the matrices of the functional and effective connectivities for task-based and task-free (resting state) paradigms [1][2][3][4][5][6]. Clustering has also been observed for dynamic functional connectivity, where the time courses of the connectivity exhibit a few discrete states with well pronounced clusters [7]. Disruption of such clustered states of connectivity may be associated with some brain disorders [8,9]. It is therefore important to investigate the emergence of clustering in neural populations that we address in this study.
Neural networks are able to adapt their structure depending on the activity of the nodes or external stimuli [10]. One of the possible mechanisms of such an adaptation, which may lead a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 also reported on clustering in the neural populations with plasticity [35,42,43]. These clusters have been observed for different models that ranged from simple phase oscillators to the models of spiking and bursting neurons and demonstrate stability with respect to heterogeneity of the interacting neurons and random perturbations [35,42,43]. In this paper we provide a simple phenomenological model and explain a mechanism governed by synaptic plasticity of the stabilization of such clusters in a neural population.
The structure of the paper is as follows. In the first section we present the model. The next section shows numerically observed multi-clusters. The detailed mechanisms of the stability of frequency clusters is explained afterwards using the simplest case of two clusters. Then we propose a phenomenological model, which describes the dynamics of two clusters taking into account the adaptation of the weights. The model is shown to reflect not only qualitative, but also some basic quantitative properties of the two-cluster formation. We also determine the set of plasticity functions (update rules), which lead to the clustering.

Model
The network of N HH neurons is described by the following system [24,25,44,45] C _ V i ¼ I i À g Na m 3 i h i ðV i À E Na Þ À g k n 4 i ðV i À E K Þ À g L ðV i À E L Þ À _ n i ¼ a n ðV i Þð1 À n i Þ À b n ðV i Þn i ; Here V i is the potential of the i-th neuron with the corresponding equilibrium potentials E Na = 50mV, E K = −77mV, and E l = −54.4mV. C = 1μF/cm 2 . Our choice of E r = 20mV corresponds to the excitatory neurons. m, h, and n are gating variables, and their dynamics depend on opening and closing rates The parameters are g Na = 120mS/cm 2 , g K = 36mS/cm 2 , and g l = 0.3mS/cm 2 . The constant current I i is set to 9μA/cm 2 so that the individual neurons are identical and oscillatory.
The synaptic input current from j-th neuron is scaled by the synaptic strength κ ij , which changes due to plasticity. The adaptation of κ ij occurs discontinuously whenever one of the neurons i or j spikes. More specifically, the discontinuous change is given by the following plasticity function where Δt ij = t i − t j is the spike time difference between the postsynaptic and presynaptic neurons; δ > 0 is a small parameter determining the size of the single update; κ max > 0 is the maximal coupling; and the plasticity function [14][15][16]21] is with positive parameters c p , τ p , c d , and τ d . We also assume no autapses and set κ ii = 0. Example of the considered plasticity function W used in our simulations is shown in Fig 1. This is a symmetric function, which corresponds to a potentiation of the coupling weights of the neurons with highly correlated firing. As we will discuss at the end of the results section, there is a family of plasticity functions of similar form that allow for the frequency clustering.

Numerical observation of synchrony and frequency clustering
In order to investigate the dynamics of network (1), we initialize the neurons and the coupling randomly and integrate the system numerically. For the parameter values τ p = 2, τ d = 5, c p = 2, c d = 1.6, and κ max = 1.5 we observe two phenomena: complete synchronization and the emergence of frequency clusters hierarchical in size, see Figs 2 and 3, respectively.   illustrates the spike times of the neurons at the beginning of the simulation. One can observe that while the neurons start with an incoherent spiking, they enhance the coherence already after a few spikes due to the interaction between them as well as the plasticity. The plasticity potentiates the connections of neurons that fire together. A complete synchronization is established and the coupling weights increase to κ max after a transient (Fig 2(B) and 2(D)). In a completely synchronized state, the individual neurons spike simultaneously, hence, the spike time differences Δt ij = 0.
The emergence of frequency clusters is shown in Fig 3  consist of two groups of synchronized neurons. These cluster states also manifest themselves as two groups of strongly coupled elements in the coupling matrix κ (Fig 3(C)). The coupling weights between the neurons from the different groups is very small or zero.
We observe that the largest cluster is formed rather quickly as time evolves, whereas the formation of the small cluster takes much more time. This is illustrated in Fig 4, where the time courses of the mean coupling within each of the two clusters are shown. The average coupling within the big cluster reaches its maximum fast (at t � 1000, solid curve in Fig 4), whereas the  Fig 4).
For the states with more clusters, each new formed cluster is significantly smaller than the previous one, see Fig 5, where three clusters are shown. The spiking period of the cluster appears to be proportional to its size: the bigger the cluster the larger is the period. Simulation of the cases with even more clusters becomes computationally expensive due to large transients.
Clustering with independent random input. To investigate the robustness of our findings, we added an α-train as additional independent random input to the membrane potential V i of every neuron: The Eq (4) models a postsynaptic potential (PSP) that is received by the neuron at certain random times τ i,k . The inter-spike interval is Gaussian distributed

Two clusters in more detail
In this section we numerically show that depending on the relative size of the two clusters, such two-cluster states can be either dynamically stable or transient leading to complete synchronization. In order to investigate the cluster stability, we initialize the system in a two cluster state with the number of neurons N s in the small cluster and N b = N − N s in the big cluster. The total number of neurons is set to N = 50. The inter-cluster couplings are set to zero initially while the intra-cluster couplings equals κ max . All neurons in the same cluster are initialized with the same initial conditions, so the clusters are fully synchronized at t = 0. Fig 7(A) shows frequency difference of two uncoupled clusters as a function of the size of the small cluster. The frequency difference demonstrates an almost linear dependence on the cluster size and decays as the size of the smaller cluster increases. Moreover, we also observe that clusters with sufficiently different sizes are stable while the clusters of similar sizes, in the considered case with N s > 8, are transient, merge into a single cluster and eventually lead to a stable completely synchronous state, see Fig 8. Although the threshold of how different the clusters should be in order to be stable is certainly model dependent, the synchronization of similar clusters is a general property. The merging of two clusters can be explained qualitatively as follows. Initially uncoupled clusters evolve each with their natural frequencies O s and O b . If their sizes are different, then O s 6 ¼ O b , and the clusters arrive in-phase periodically with the 'beating' frequency ΔO = O b − O s . As soon as such an in-phase episode occurs, the interspike intervals Δt ij between any two neurons from different clusters become small and, hence, due to the plasticity rule W (see Eq (2) and Fig 1) the inter-cluster coupling weights increase. Moreover, the duration of such an in-phase episode depends on the frequency difference between the clusters. As a result, for a small frequency difference ΔO, the time interval where the clusters are practically in-phase is sufficiently long in order to potentiate the coupling weights to their maximum value. This unites the two clusters into one. In contrast, for large ΔO, such an episode is short, and the inter-cluster coupling remains small, which keeps the clusters oscillating at different frequencies in a stable manner.   Starting from the two-cluster state, after t = 1500ms, the coupling between the clusters increases, see Fig 8(B) due to the "in-phase episode" when the clusters are synchronous. Afterwards, however, the inter-cluster coupling weights return to their initial configuration (Fig 8(C)), since the spike time differences for neurons from different clusters are again far enough apart to cause the depression of the inter-cluster synapses. Such a process is repeated every time the clusters meet and is typical for the stable cluster states. A typical case of transient clusters is presented in Fig  8(D)-8(F) for N s = 9. The inter-cluster coupling is again potentiated when the clusters meet, but it does not decrease again, and the clusters merge in a single cluster of a fully coupled and synchronized regime (Fig 8(F)). The transient time that could be elapsed until the cluster fusion depends on the cluster size as illustrated in Fig 7(B). During the in-phase episode, the cluster with the higher natural spiking rate slows down significantly, while the slower cluster (with larger number of neurons N b ) speeds up a little. For a stable cluster state the cluster frequencies again deviate from each other (Fig 8(G)), whereas all neurons fire with the same frequency when the clusters unite into one (Fig 8(H)). We found this phenomenon for different numbers of neurons and different κ max . Increasing κ max increases the initial period difference, but the behavior in general stays the same. Fig 9 shows the dynamics of the mean synaptic activity SðtÞ ¼ 1 N P N i¼1 s i ðtÞ of the network in the case of two stable clusters, which models the dynamics of LFP. During the in-phase episodes of the two clusters, S(t) has a higher amplitude, because both clusters spike synchronously. The maximum amplitude is generated by maximum synchronization in the network. The low amplitude of S(t), on the other hand, corresponds to the time intervals when the clusters are out of phase. In the latter case, the mean synaptic activity shows two peaks, the higher Frequency cluster formation and slow oscillations in neural populations with plasticity peak is generated by the larger cluster and the lower by the smaller one, see Fig 9(B). For the considered case, the synchronized oscillations of individual neurons in the clusters take place at a time scale of several milliseconds (period *15 ms, Fig 9(B)), see also Fig 8(G) and 8(H). The neurons are tonically spiking. The frequency difference ΔO between clusters is, however, of the order of sub-Hz, because the corresponding cluster frequencies are close to each other (Fig 8(G) and 8(H)). Then the modulation of S(t) is observed at a much slower timescale of a few seconds, which is of two orders of magnitude slower than the intrinsic neural firing, see Fig 9(A), as observed in empirical data of the brain activity [40,41].
In the following section, a phenomenological model is introduced in order to further investigate the dynamics of two clusters.

Phenomenological model
Model derivation. In this section we introduce a reduced qualitative model for the coupling and phase difference of two clusters. The model is based on the assumption that oscillators are synchronized identically within each cluster and the coupling between the clusters is weak. As a result, the interaction between oscillatory clusters can be described in the framework of two coupled phase oscillators that are interacting via their phase differences [46][47][48][49] where ω 1 and ω 2 are the natural frequencies of the individual clusters, F 1 and F 2 are effective interaction functions. For the phase difference φ = φ 1 − φ 2 , this system reads where ω = ω 1 − ω 2 is the difference of the natural frequencies, and F(φ) = F 1 (φ) − F 2 (−φ). Since the clusters are synchronized for a sufficiently small frequency mismatch ω, the periodic interaction function F(φ) must satisfy F(0) = 0 and F 0 (0) > 0. The latter means that there is a stable equilibrium φ = 0 for small ω. Aiming at a qualitative insight, we further simplify the model by assuming that F(φ) = σ sin(φ + α), where sin φ can be viewed as a first Fourier harmonic of the interaction function and σ as an effective coupling weight. The parameter α = sin −1 (ω/σ max ) is a constant phase shift assuring that the phase difference of the synchronized cluster is zero. In fact, for small ω, this parameter is also small and it does not play important role in the qualitative behavior of the model apart from a small shift of the synchronized state to φ = 0.
Another component of the model is the plasticity-driven changes of the coupling σ. In order to derive the equation for σ, we consider the STDP update in the case of a periodic motion of the clusters. We assume that the coupling σ is proportional to an averaged coupling between the clusters. This is a natural assumption in the case of weakly coupled systems. Let us find out how the update of the intercluster coupling depends on the phase difference φ. For a given phase difference φ and the frequencies , the spiking period of the both clusters can be approximated as T � 2p= � o up to small terms of order ω, and the distance ΔT between the spikes of two clusters Since the spike time differences ΔT and T − ΔT occur recursively, see Fig 10(B), the updates per unit time sum to the function where Since the update of σ is proportional to the obtained function, and taking into account the smallness of δ, this update can be written as _ s ¼ εGðφÞ, where ε is a small parameter of the coupling adaptation that controls the scale separation between the fast dynamics of the clusters and the slow dynamics of the coupling.
Additionally, the coupling strength σ(t) should be bounded to the interval [0, σ max ] by imposing cut-off conditions. More specifically, the derivative _ s is discontinuous at the boundaries σ = 0 and σ = σ max , i.e. _ s ¼ max f0; εGðφÞg for σ = 0 and _ s ¼ min f0; εGðφÞg for σ = σ max . The considered cut-off corresponds to "hard" bound conditions [50]. Another possibility would be "soft" or "multiplicative" bounds [51], when the update is proportional to the distance to the boundary. We consider here the hard bound, since it corresponds to the hard bound of the STDP rule for HH system. The final phenomenological model reads as follows with frequency mismatch ω > 0 and α = sin −1 (ω/σ max ). Properties of the model. Phase space of system (10)-(11) is two dimensional with (φ, σ) 2 S 1 × [0, σ max ]. The nullclines are given by G(φ) = 0 for _ s ¼ 0 and σ = ω/sin(φ + α) for _ φ ¼ 0 in the internal points of the phase space. For the parameter values as in Fig 11, the φnullcline corresponds to the two lines φ = φ � � 0.23 and φ = −φ � , while the σ-nullcline to a U-shaped nonlinear curve (grey lines in Fig 11).
There is just one fixed point (φ � , σ � = ω/sin φ � ) of saddle-type within the region σ 2 (0, σ max ). This point is given by the intersection of the nullclines. Fig 11 shows this fixed point and its stable and unstable separatrixes (black lines). An additional fixed point as well as periodic attractor emerge in system (10)-(11) due to the non-smoothness at the boundaries. More specifically, three situations are observed: (I): One globally stable fixed point S = (0, σ max ) which corresponds to the fusion of the two clusters into one. The coupling σ = σ max and the phase difference is zero at the fixed point, see Fig 11. All orbits are approaching this stable fixed point with time. This corresponding phase portrait is shown in Fig 11(A).

(II):
Coexistence of the stable fixed point S = (0, σ max ) and a stable periodic orbit, see Fig  11(B). As in the case (I), the fixed point corresponds to the merging of two clusters. The periodic orbit corresponds to two simultaneously existing clusters. The clusters possess different frequencies and, as a result, the phase difference is not bounded and rotate along the circular direction φ. Part of the periodic orbit is located on the boundary σ = 0, i.e. vanishing inter-cluster coupling. The coupling σ(t) increases between −φ � and φ � and decreases otherwise. In fact, one can parameterise the coupling σ by the phase φ on the periodic attractor. In the case when σ(φ � ) < σ � , the solution returns to the boundary σ = 0, moves along it till the orbit reaches the point (−φ � , 0), and the periodic motion repeats.
(III): When max φ G(φ) < 0 then there exists globally stable periodic solution φ = ωt + φ 0 , σ = 0. In such a case, the fixed point on the boundary disappears. Formally, this corresponds to an uncoupling between the clusters. However, in the original HH system, this parameter regime corresponds to complete uncoupling of all oscillators because of the depression of all synapses.
In fact, the parameter boundary between the cases (I) and (II) is determined by the condition σ(φ � ) = σ � , which can be interpreted geometrically as hitting the point (−φ � , 0) by the stable manifold of the saddle equilibrium point, see Fig 11(C). In this special case, the saddle equilibrium attracts the whole set of points from the phase space that is below the stable manifolds, see white area in Fig 11(C). In case (II), the separation between the basins of attraction of the fixed point and the periodic orbit are given by the saddle equilibrium and its stable manifolds. A sufficient condition for the case (III) is given by c d � c p and τ d � τ p . Under these conditions G(φ) � 0 for all φ.
Summarizing, the case (II) corresponds to the situation when clusters are stable and do not merge into one. For this, initial conditions must belong to the basin of attraction of the periodic solution (Fig 11(B), blue domain). The analysis of the phenomenological model indicates that the cluster case always coexists with stable complete synchronization.
Comparison of the model and cluster dynamics in HH network. In order to compare dynamics of the phenomenological model (10)-(11) and the original system (1)-(2), we ran a series of simulations of the HH network for parameter values that allow for a stable two-cluster solution. The phases of the clusters are calculated as φ 1;2 ðtÞ ¼ 2p tÀ t k t kþ1 À t k þ 2pk for t 2 ½t k ; t kþ1 Þ, where {t 1 , . . ., t n , . . .} are spiking times with t k < t k+1 [47]. Correspondingly, the phase difference is φ HH (t) = φ 1 (t) − φ 2 (t). The coupling measure σ HH is given by the mean inter-cluster coupling.
Extracting the quantities σ HH and φ HH from the numerically computed solutions of HH system (1)-(3) we obtain a two-dimensional projection of the solution to the plane (φ HH , σ HH ), see Fig 12. The discontinuities in the orbits are related to the discrete STDP updates. Additionally, since the phases φ 1,2 (t) can be firstly accessed after the both clusters fired, some of the area of the phase diagram (see white area in Fig 12) was not accessible. This "empty" area corresponds to anti-phase initial conditions, which are very sensitive, and, after each cluster fires, they appear immediately either in the red or blue area. Nevertheless, the behavior has the same qualitative features as in the phenomenological model, compare Figs 11 and 12. Criteria for the emergence of clusters. Model (10)- (11) can be used to describe plasticity functions, which lead to multiple clusters. For this, we investigate numerically the condition σ (φ � ) = σ � . More specifically, system (10)- (11) was initialized at the point (−φ � , 0) and numerically integrated forward in time. If σ(φ � ) < σ � , the two clusters are stable and do not merge. This procedure can be repeated for different parameter values.
In order to restrict the set of plasticity parameters, we fix τ p = 2 and τ d = 5 and vary c p and c d . The results of the simulation are shown in Fig 13(A). The white, black and grey parameter areas correspond to the appearance of stable periodic solution of (10)-(11) (case (II)), globally stable fixed point (case (I)) and the case (III), respectively.
In order to compare the parameter regions obtained for the phenomenological model ( Fig  13(A)) with those for the original HH system, we ran numerical simulations of system (1)- (3) with N = 50 neurons and N s = 7 neurons in the small cluster. Starting from the two-cluster state, we monitor the dynamics of the clusters. Fig 13(B) shows the results: the white region corresponds to the case when the clusters survive and stay apart after the simulation time 3000 ms, black-when the clusters merge into one synchronous group, and grey-when the clusters split into uncoupled neurons. This behaviour stays qualitatively the same for different cluster sizes. However, depending on the frequency difference between the clusters, the set of parameters allowing stable cluster states may change its size.
Comparison of the results for the phenomenological system and the HH system in the Fig 13(A) and 13(B) shows that the phenomenological model provides a reasonable approximation.

Conclusion
Our results show that adaptive neural networks are able to generate self-consistently dynamics with different frequency bands. In our case, each cluster corresponds to a strongly connected Frequency cluster formation and slow oscillations in neural populations with plasticity component with a fixed frequency. Due to a sufficiently large difference of the cluster sizes and frequencies, the inter-cluster interactions are depreciated, while the intra-cluster interactions are potentiated. In this study, we describe the mechanisms behind the formation and stabilization of these clusters. In particular, we explain why the significant difference between the cluster sizes is important for the decoupling of the clusters. From a larger perspective, the decoupling of the clusters in our case is analogous to the decoupling of timescales in systems with multiple timescales.
Furthermore, we present a two-dimensional phenomenological model which allows for a detailed study of the clustering mechanisms. Despite of the approximations made by the derivation, the model coincides surprisingly well with the adaptive Hodgkin-Huxley network. Using the phenomenological model, we find parameter regions of the plasticity function, where stable frequency clustering can be observed.
Clustering behavior also emerges at the brain scale, where synchronized communities of brain regions constituting large distributed functional networks can intermittently be formed and dissolved [52,53]. Such clustering dynamics can shape the structured spontaneous brain activity at rest as measured by fMRI. In this study, we show that slow oscillations based on the modulation of synchronized neural activity can already be formed at the resolution level of a single neural population if adaptive synapses are taken into account. These modulations of the amplitude of the mean field can be generated in a stable manner [Fig 9], see also Ref. [35]. The mechanism relies on fluctuations of the extent of synchronization of tonically firing neurons. This is caused by the splitting of the neural population into clusters and the corresponding cluster dynamics. It might contribute to the emergence of slow brain rhythms of electrical (LFP, EEG) and metabolic (BOLD) brain activity reported by [38][39][40][41]. Frequency cluster formation and slow oscillations in neural populations with plasticity However, other mechanisms for generating slow oscillations are possible. The papers [54,55] discussed the emergence of slow oscillatory activity (< 1Hz) that can be observed in vivo in the cortex during slow-wave sleep, under anesthesia or in vitro in neural populations. The suggested mechanism relies on the corresponding modulation of the firing of individual neurons, and the slow oscillation at the population level was proposed to be the result of very slow bursting of individual neurons that synchronize across the neural population. In contrast, the present work shows that the slow oscillations of the population mean field can also emerge when the firing of individual neurons is not affected. The neurons may tonically fire at high frequencies. The amplitude of the population mean field then oscillates at much lower frequencies due to the slow modulation caused by the cluster dynamics.
Additionally, we would like to mention that the observed frequency clustering resembles phenomenologically the weak chimera states [56,57] where clusters with different frequencies are formed in symmetrically coupled oscillators without adaptation. However the properties and mechanisms of the appearance of such clusters are different from those presented here, which are essentially based on the slow adaptation.
To conclude, we observe self-organised emergence of clusters in neural networks with STDP. The clustering splits the neural population into groups synchronised at different frequencies, which determine the dynamics of the clusters. These cluster dynamics might play a role in low frequency oscillations during the resting state and can be described by a twodimensional model.