Hierarchical clusters in neuronal populations with plasticity

TL;DR

This paper investigates frequency clustering in Hodgkin-Huxley neuron networks with spike timing-dependent plasticity, revealing how neural populations split into synchronized groups and how this may relate to slow oscillations in neural signals.

Contribution

It introduces a phenomenological model for multi-clustering dynamics considering synaptic plasticity and identifies plasticity rules that lead to such clustering phenomena.

Findings

Neural populations can split into multiple synchronized clusters.

Amplitude of mean field exhibits low-frequency modulations.

Plasticity rules can induce multi-clustering behavior.

Abstract

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.

Equations

$$\begin{array}{rl}C{\dot{V}}_i& =I_i-g_{Na}{m}_i^3{h}_i(V_i-E_{Na})-g_k{n}_i^4(V_i-E_K)-g_L(V_i-E_L)\\ & -\frac{(V_i-E_r)}{N}\sum _{j=1}^N{\kappa }_{ij}{s}_j,\\ {\dot{m}}_i& ={\alpha }_m(V_i)(1-m_i)-{\beta }_m(V_i)m_i,\\ {\dot{h}}_i& ={\alpha }_h(V_i)(1-h_i)-{\beta }_h(V_i)h_i,\\ {\dot{n}}_i& ={\alpha }_n(V_i)(1-n_i)-{\beta }_n(V_i)n_i,\\ {\dot{s}}_i& =\frac{5(1-s_i)}{1+e^{(}\frac{-V_i+3}{8})}-s_i.\end{array}$$

$${\alpha }_m(V)$$

$${\beta }_m(V)$$

$${\alpha }_h(V)$$

$${\beta }_h(V)$$

$${\alpha }_n(V)$$

$${\beta }_n(V)$$

$${\kappa }_{ij}\to \{\begin{array}{ll}0,& \text{ if }{\kappa }_{ij}+\delta W(\Delta t_{ij})<0\\ {\kappa }_{ij}+\delta W(\Delta t_{ij}),& \text{ if }0\le {\kappa }_{ij}+\delta W(\Delta t_{ij})\le {\kappa }_{\max }\\ {\kappa }_{\max },& \text{ if }{\kappa }_{ij}+\delta W(\Delta t_{ij})>{\kappa }_{\max }\end{array}$$

$$W(\Delta t_{ij})=c_p{e}^{-\frac{|\Delta t_{ij}|}{{\tau }_p}}-c_d{e}^{-\frac{|\Delta t_{ij}|}{{\tau }_d}}$$

$$I_i^{input}(t)=I(V_r-V_i(t))\sum _{{\tau }_{i,k}<t}\alpha (t-{\tau }_{i,k})e^{-\alpha (t-{\tau }_{i,k})})$$

$${\dot{\phi }}_1$$

$${\dot{\phi }}_2$$

$$\dot{\phi }=\omega -F(\phi )$$

$$\Delta T=\left[T\frac{{\phi }_1}{2\pi }-T\frac{{\phi }_2}{2\pi }\right]\mathrm{mod}T=\left[T\frac{\phi }{2\pi }\right]\mathrm{mod}T\approx \frac{\phi \mathrm{mod}2\pi }{\bar{\omega }}.$$

$$\frac{\delta }{T}(W(T-\Delta T)+W(\Delta T))=\frac{\delta \bar{\omega }}{2\pi }G(\phi ),$$

$$G(\phi ):=W\left(\frac{2\pi -(\phi \mathrm{mod}2\pi )}{\bar{\omega }}\right)+W\left(\frac{\phi \mathrm{mod}2\pi }{\bar{\omega }}\right).$$

$$\dot{\phi }$$

$$\dot{\sigma }$$

Full text

Hierarchical clusters in neuronal populations with plasticity

Vera Röhr1, Rico Berner2,3, Ewandson L. Lameu4,5,

Oleksandr V. Popovych6,7, Serhiy Yanchuk3

1 Neurotechnology Group, Technische Universität Berlin; 2 Institut für Theoretische Physik, Technische Universität Berlin; 3 Institut für Mathematik, Technische Universität Berlin; 4 National Institute for Space Research (INPE), Brazil; 5 Institut für Physik, Humboldt-Universität zu Berlin 6 Institute of Neuroscience and Medicine - Brain $\&$ Behaviour (INM-7), Research Centre Jülich, Germany; 7Institute for Systems Neuroscience, Medical Faculty, Heinrich-Heine University Düsseldorf, Germany

Abstract

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.

Author summary

Synaptic plasticity is one of the key mechanisms that allow for neural networks to adapt their structure. Depending on the mutual neural activity, the efficacy of synapses may change, which results in short- or long-term potentiation or depression of a synapse. In this paper we investigate the structural changes that are caused by the spike timing-dependent plasticity (STDP), where the synaptic weights are adapted depending on the difference of spiking times between the pre- and postsynaptic neuron. The adaptation is considered to be symmetric as experimentally found for hippocampal synapses [1] and can also be derived from asymmetric STDP under certain conditions [2]. In an adaptive network of Hodgkin-Huxley neurons we observe the emergence of clusters of neurons that are synchronized at different frequencies in different clusters. Our study shows that such a self-organized cluster formation is robust against changes of plasticity function. While the spiking frequency of each synchronized cluster appear to be on the timescale of individual neurons, the amplitude of the mean field of such cluster system can evolve at a few orders of magnitude slower timescales. The reported slow modulations of the mean neural activity in a neural population with plasticity can explain the emergence of slow cortical oscillations detected in the empirical data of spectral power of local field potentials (LFP) and electroencephalographic (EEG) signals, which correlate with spontaneous fluctuations of the blood oxygen level-dependent (BOLD) signals measured by functional magnetic resonance imaging (fMRI) [3, 4].

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 [10, 5, 6, 7, 8, 9]. 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 [11]. Disruption of such clustered states of connectivity may be associated with some brain disorders [12, 13]. 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 [14]. 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 [15]. 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 [16, 17]. 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 [18, 19, 20, 21, 22].

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 [15]. 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 [23, 24]. 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 [25]. 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 [18, 25]. 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 [26, 27, 28, 29]. On the other hand, a desynchronized state can lead to a depression of synaptic weights [26, 27]. 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 [30, 31, 32, 33, 34]. Recently, a hierarchical clustering was observed in phenomenological models of adaptive networks of phase oscillators [35, 36, 37]. 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 [37, 38] 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 [39].

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 [1] and can also be derived from asymmetric STDP for an ”effective time window” [2]. Then the observed clusters are bidirectionally coupled [39]. 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 [39]. 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 [3, 4, 40, 41]. Several other modeling studies have also reported on clustering in the neural populations with plasticity [42, 43, 39]. 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 [42, 43, 39]. 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 [44, 45, 28, 29]

[TABLE]

Here $V_{i}$ is the potential of the $i$-th neuron with the corresponding equilibrium potentials $E_{Na}=50$mV, $E_{K}=-77$mV, and $E_{l}=-54.4$mV. $C=1\mu F/cm^{2}$. Our choice of $E_{r}=20$mV corresponds to the excitatory neurons. $m$, $h$, and $n$ are gating variables, and their dynamics depend on opening and closing rates

[TABLE]

The parameters are $g_{Na}=120$mS/cm2, $g_{K}=36$mS/cm2, and $g_{l}=0.3$mS/cm2. The constant current $I_{i}$ is set to $9\mu A/\mathrm{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 $\kappa_{ij}$, which changes due to plasticity. The adaptation of $\kappa_{ij}$ occurs discontinuously whenever one of the neurons $i$ or $j$ spikes. More specifically, the discontinuous change is given by the following plasticity function

[TABLE]

where $\Delta t_{ij}=t_{i}-t_{j}$ is the spike time difference between the postsynaptic and presynaptic neurons; $\delta>0$ is a small parameter determining the size of the single update; $\kappa_{\max}>0$ is the maximal coupling; and the plasticity function [18, 19, 25, 20] is

[TABLE]

with positive parameters $c_{p},\tau_{p},c_{d},$ and $\tau_{d}$. We also assume no autapses and set $\kappa_{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.

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 $\tau_{p}=2,\tau_{d}=5,c_{p}=2,c_{d}=1.6$, and $\kappa_{\max}=1.5$ we observe two phenomena: complete synchronization and the emergence of frequency clusters hierarchical in size, see Figs. 2 and 3, respectively.

Figure 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 $\kappa_{\max}$ after a transient (Fig. 2(B),(D)). In a completely synchronized state, the individual neurons spike simultaneously, hence, the spike time differences $\Delta 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 $\kappa$ (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\approx 1000$, solid curve in Fig. 4), whereas the smaller cluster in Fig. 3(C, F) is formed through the merging of transient clusters and finally establishes at $t\approx 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 $\alpha$-train as additional independent random input to the membrane potential $V_{i}$ of every neuron:

[TABLE]

The Eq. (4) models a postsynaptic potential (PSP) that is received by the neuron at certain random times $\tau_{i,k}$. The inter-spike interval is Gaussian distributed $\tau_{i,k+1}-\tau_{i,k}=\Delta\tau_{i,k}\sim\mathcal{N}(14ms,4ms)$. $\alpha$ is set to ${24}/{\langle\Delta\tau_{i,k}\rangle}.$

The numerical simulations Fig. 6 show that the clustering is still observed under the influence of random input $I_{i}^{input}(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-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 $\kappa_{\max}$. All neurons in the same cluster are initialized with the same initial conditions, so the clusters are fully synchronized at $t=0$.

Figure 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 $\Omega_{s}$ and $\Omega_{b}$. If their sizes are different, then $\Omega_{s}\neq\Omega_{b}$, and the clusters arrive in-phase periodically with the ’beating’ frequency $\Delta\Omega=\Omega_{b}-\Omega_{s}$. As soon as such an in-phase episode occurs, the interspike intervals $\Delta 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 $\Delta\Omega$, 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 $\Delta\Omega$, such an episode is short, and the inter-cluster coupling remains small, which keeps the clusters oscillating at different frequencies in a stable manner.

Figure 8(A)-(C) displays an example of two-cluster stable state with $N_{s}=8$. Starting from the two-cluster state, after $t=1500$ms, 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)-(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).

Figures 8(G, 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 $\kappa_{\max}$. Increasing $\kappa_{\max}$ increases the initial period difference, but the behavior in general stays the same.

Figure 9 shows the dynamics of the mean synaptic activity $S(t)=\frac{1}{N}\sum_{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 $\sim 15$ ms, Fig. 9(B)), see also Fig. 8(G), (H). The neurons are tonically spiking. The frequency difference $\Delta\Omega$ between clusters is, however, of the order of sub-Hz, because the corresponding cluster frequencies are close to each other (Fig. 8(G), (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 [3, 4].

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]

[TABLE]

where $\omega_{1}$ and $\omega_{2}$ are the natural frequencies of the individual clusters, $F_{1}$ and $F_{2}$ are effective interaction functions. For the phase difference $\varphi=\varphi_{1}-\varphi_{2}$, this system reads

[TABLE]

where $\omega=\omega_{1}-\omega_{2}$ is the difference of the natural frequencies, and $F(\varphi)=F_{1}(\varphi)-F_{2}(-\varphi)$.

Since the clusters are synchronized for a sufficiently small frequency mismatch $\omega$, the periodic interaction function $F(\varphi)$ must satisfy $F(0)=0$ and $F^{\prime}(0)>0$. The latter means that there is a stable equilibrium $\varphi=0$ for small $\omega$. Aiming at a qualitative insight, we further simplify the model by assuming that $F(\varphi)=\sigma\sin(\varphi+\alpha)$, where $\sin\varphi$ can be viewed as a first Fourier harmonic of the interaction function and $\sigma$ as an effective coupling weight. The parameter $\alpha=\sin^{-1}(\omega/\sigma_{\max})$ is a constant phase shift assuring that the phase difference of the synchronized cluster is zero. In fact, for small $\omega$, 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 $\varphi=0$.

Another component of the model is the plasticity-driven changes of the coupling $\sigma$. In order to derive the equation for $\sigma$, we consider the STDP update in the case of a periodic motion of the clusters. We assume that the coupling $\sigma$ 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 $\varphi$. For a given phase difference $\varphi$ and the frequencies $\omega_{1}=\bar{\omega}+\omega/2$, $\omega_{2}=\bar{\omega}-\omega/2$ (here we introduced the mean frequency $\bar{\omega}$), the spiking period of the both clusters can be approximated as $T\approx 2\pi/\bar{\omega}$ up to small terms of order $\omega$, and the distance $\Delta T$ between the spikes of two clusters

[TABLE]

Since the spike time differences $\Delta T$ and $T-\Delta T$ occur recursively, see Fig. 10(B), the updates per unit time sum to the function

[TABLE]

where

[TABLE]

Since the update of $\sigma$ is proportional to the obtained function, and taking into account the smallness of $\delta$, this update can be written as $\dot{\sigma}=\varepsilon G(\varphi)$, where $\varepsilon$ 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 $\sigma(t)$ should be bounded to the interval $[0,\sigma_{\max}]$ by imposing cut-off conditions. More specifically, the derivative $\dot{\sigma}$ is discontinuous at the boundaries $\sigma=0$ and $\sigma=\sigma_{\max}$, i.e. $\dot{\sigma}=\max\{0,\varepsilon G(\varphi)\}$ for $\sigma=0$ and $\dot{\sigma}=\min\{0,\varepsilon G(\varphi)\}$ for $\sigma=\sigma_{\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

[TABLE]

with frequency mismatch $\omega>0$ and $\alpha=\sin^{-1}(\omega/\sigma_{\max})$.

Properties of the model

Phase space of system (10)-(11) is two dimensional with $(\varphi,\sigma)\in S^{1}\times[0,\sigma_{\max}]$. The nullclines are given by $G(\varphi)=0$ for $\dot{\sigma}=0$ and $\sigma={\omega}/{\sin(\varphi+\alpha)}$ for $\dot{\varphi}=0$ in the internal points of the phase space. For the parameter values as in Fig. 11, the $\varphi$-nullcline corresponds to the two lines $\varphi=\varphi^{*}\approx 0.23$ and $\varphi=-\varphi^{*}$, while the $\sigma$-nullcline to a U-shaped nonlinear curve (grey lines in Fig. 11).

There is just one fixed point $(\varphi^{*},\sigma^{*}=\omega/\sin\varphi^{*})$ of saddle-type within the region $\sigma\in(0,\sigma_{\max})$. This point is given by the intersection of the nullclines. Figure 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,\sigma_{\max})$ which corresponds to the fusion of the two clusters into one. The coupling $\sigma=\sigma_{\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,\sigma_{\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 $\varphi$. Part of the periodic orbit is located on the boundary $\sigma=0$, i.e. vanishing inter-cluster coupling. The coupling $\sigma(t)$ increases between $-\varphi^{*}$ and $\varphi^{*}$ and decreases otherwise. In fact, one can parameterise the coupling $\sigma$ by the phase $\varphi$ on the periodic attractor. In the case when $\sigma(\varphi^{*})<\sigma^{*}$, the solution returns to the boundary $\sigma=0$, moves along it till the orbit reaches the point $(-\varphi^{*},0)$, and the periodic motion repeats.

(III): When $\max_{\varphi}G(\varphi)<0$ then there exists globally stable periodic solution $\varphi=\omega t+\varphi_{0}$, $\sigma=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 $\sigma(\varphi^{*})=\sigma^{*}$, which can be interpreted geometrically as hitting the point $(-\varphi^{*},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}\geq c_{p}$ and $\tau_{d}\geq\tau_{p}$. Under these conditions $G(\varphi)\leq 0$ for all $\varphi$.

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 $\varphi_{1,2}(t)=2\pi\frac{t-t_{k}}{t_{{k+1}}-t_{k}}+2\pi k\text{ for }t\in[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 $\varphi_{HH}(t)=\varphi_{1}(t)-\varphi_{2}(t)$. The coupling measure $\sigma_{HH}$ is given by the mean inter-cluster coupling.

Extracting the quantities $\sigma_{HH}$ and $\varphi_{HH}$ from the numerically computed solutions of HH system (1)-(3) we obtain a two-dimensional projection of the solution to the plane $(\varphi_{HH},\sigma_{HH})$, see Fig. 12. The discontinuities in the orbits are related to the discrete STDP updates. Additionally, since the phases $\varphi_{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 $\sigma(\varphi^{*})=\sigma^{*}$. More specifically, system (10)-(11) was initialized at the point $(-\varphi^{*},0)$ and numerically integrated forward in time. If $\sigma(\varphi^{*})<\sigma^{*}$, 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 $\tau_{p}=2$ and $\tau_{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. Figure 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 Figs. 13(A,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. [39]. 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 [3, 40, 4, 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.

Acknowledgments

VR, SY and RB acknowledge the financial support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Project 384950143, Project 411803875 and Project 308748074. OVP acknowledges the support by the European Union’s Horizon 2020 Research and Innovation Programme under Grant Agreement 826421 (VirtualBrainCloud). ELL acknowledges the financial support by the São Paulo Research Foundation (processes FAPESP 2016/23398-8 and 2017/13502-5).