Frequency cluster formation and slow oscillations in neural populations with plasticity

Summary

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.

Keywords:

Population dynamics – functional magnetic resonance imaging – Neurons – Synapses – Neural networks – Action potentials – Synaptic plasticity – Neuronal plasticity

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–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 to persistent changes in neural connections and relate to learning and memory, is synaptic plasticity [11]. The efficacy of synapses to transmit the electrical potential between neurons may increase or decrease depending on the mutual neural activity, which results in short- or long-term potentiation or, respectively, depression of synapses [12, 13]. An example is spike timing-dependent plasticity (STDP) which describes the synaptic weight change as a function of the difference of spiking times between pre- and post-synaptic neurons [14–18].

One of the famous plasticity rules, the Hebbian rule, assumes that the modifications of the synaptic weights are driven by correlations in the firing activity of pre- and post-synaptic neurons. More specifically, it assumes that those connections are potentiated, for which one neuron contributes to the firing of another [11]. Nevertheless, in many publications, the Hebbian rule is considered in a more narrow sense of a closeness between the spiking times: the smaller the distances between the spikes are the higher is the potentiation of the corresponding synapse [19, 20]. In this work, we are dealing with spike-based learning rules rather than rate-based.

Previous studies of neural networks with STDP showed that such networks can evolve and create various coupling structures. For instance, the weights can exhibit stable localized spatial structures, that can be interpreted as receptive fields [21]. These structures can be either unidirectionally of bidirectionally coupled, depending on the plasticity rule or external input properties. The STDP mechanism plays an important role in temporal coding of information by spikes [14, 21]. On the one hand, a synchronized firing in neural ensembles with STDP can be stabilized through potentiation of synaptic coupling by stimulation-induced transient synchronization of neurons [22–25]. On the other hand, a desynchronized state can lead to a depression of synaptic weights [22, 23]. Thus, neural networks with plasticity are prone to a co-existence of different stable dynamical and structural states, which may be realized by choosing appropriate initial conditions or stimulation procedures.

Human brain networks demonstrate different degrees of modularity, sometimes with hierarchical features [26–30]. Recently, a hierarchical clustering was observed in phenomenological models of adaptive networks of phase oscillators [31–33]. As a result of an adaptation, the network evolved into groups of strongly-connected clusters, while the coupling between the groups was depressed. The stability analysis of such clusters reveals [33, 34] that the preferred stable cluster configuration corresponds to significantly different sizes of the clusters. The dynamics within each cluster are frequency-synchronized, while the frequencies between clusters differ. Thus, self-organized emergence of clusters leads to the emergence of different collective frequencies in the system. The multi-stability of such clusters was also observed in ensembles of Morris-Lecar bursting neurons with STDP in [35].

In this paper we report on the phenomenon of clustering with respect to connectivity and frequencies in a network of adaptively coupled Hodgkin-Huxley (HH) neurons. The spike timing-dependent adaptation is considered to be symmetric as experimentally found for hippocampal synapses [36] and can also be derived from asymmetric STDP for an “effective time window” [37]. Then the observed clusters are bidirectionally coupled [35]. Splitting of a neural population to a few clusters synchronized at different frequencies could lead to a slow waxing and waning of the amplitude of the mean field, where the clusters transiently gather together and move apart as the time evolves [35]. The frequency of such a modulation of the mean neural activity could be much smaller than the firing rate of individual neurons and depends on the differences between the clusters’ frequencies. The emergence of synchronized clusters could explain the origin of the low-frequency modulation of the spectral power of macroscopic brain signals like local field potentials (LFP) or electroencephalographic (EEG) signals in higher frequency bands, which also correlates with slow oscillations of the blood oxygen level-dependent (BOLD) signal measured by fMRI [38–41]. Several other modeling studies have 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.

Materials and methods

Model

The network of N HH neurons is described by the following system [24, 25, 44, 45]

The parameters are g_Na = 120mS/cm², g_K = 36mS/cm², and g_l = 0.3mS/cm². The constant current I_i is set to 9μA/cm² 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

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.

Results

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.

Fig 2(A) shows the initial coupling weights, and Fig 2(C) 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 for two clusters. The system in Fig 3 possesses the same parameters as in Fig 2, and the difference is just another realization of random initial conditions. In contrast to the synchronized state, the final state shown in Fig 3(F) 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 smaller cluster in Fig 3(C) and 3(F) is formed through the merging of transient clusters and finally establishes at t ≈ 17000 (dashed curve in 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 τ i , k + 1 - τ i , k = Δ τ i , k ∼ N ( 14 m s , 4 m s ). α is set to 24/〈Δτ_i,k〉.

The numerical simulations Fig 6 show that the clustering is still observed under the influence of random input I i i n p u t ( t ) of intensity I. More specifically, for sufficiently weak perturbations with I < 0.01, all three clusters survive (Fig 6(A)). With increasing the amplitude I, the clusters start to decouple. The smaller clusters are affected first (Fig 6(B)–6(D)), they start desynchronizing at I = 0.01. The biggest cluster keeps shrinking in size while I is increased and finally for I = 0.07 the whole network decouples (Fig 6(E)).

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 Ω_s and Ω_b. If their sizes are different, then Ω_s ≠ Ω_b, and the clusters arrive in-phase periodically with the ‘beating’ frequency ΔΩ = Ω_b − Ω_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 ΔΩ, 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 ΔΩ, such an episode is short, and the inter-cluster coupling remains small, which keeps the clusters oscillating at different frequencies in a stable manner.

Fig 8(A)–8(C) displays an example of two-cluster stable state with N_s = 8. 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).

Fig 8(G) and 8(H) show how the spiking frequency of the clusters change over time. 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 ∑ i = 1 N 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 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 ΔΩ 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–49]

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. 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⁻¹(ω/σ_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 ω 1 = ω ¯ + ω / 2, ω 2 = ω ¯ - ω / 2 (here we introduced the mean frequency ω ¯), the spiking period of the both clusters can be approximated as T ≈ 2 π / ω ¯ 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

Since the update of σ is proportional to the obtained function, and taking into account the smallness of δ, this update can be written as σ ˙ = ε 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 σ ˙ is discontinuous at the boundaries σ = 0 and σ = σ_max, i.e. σ ˙ = max { 0 , ε G ( φ ) } for σ = 0 and σ ˙ = min { 0 , ε 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

Properties of the model

Phase space of system (10)–(11) is two dimensional with (φ, σ) ∈ S¹ × [0, σ_max]. The nullclines are given by G(φ) = 0 for σ ˙ = 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 σ ∈ (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. 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 ) = 2 π t - t k t k + 1 - t k + 2 π k for t ∈ [ t k , t k + 1 ), where {t₁, …, t_n, …} are spiking times with t_k < t_k+1 [47]. Correspondingly, the phase difference is φ_HH(t) = φ₁(t) − φ₂(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 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–41].

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 two-dimensional model.