Multistability in a star network of Kuramoto-type oscillators with synaptic plasticity

We analyze multistability in a star-type network of phase oscillators with coupling weights governed by phase-difference-dependent plasticity. It is shown that a network with N leaves can evolve into \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2^N$$\end{document}2N various asymptotic states, characterized by different values of the coupling strength between the hub and the leaves. Starting from the simple case of two coupled oscillators, we develop an analytical approach based on two small parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon$$\end{document}ε and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document}μ, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon$$\end{document}ε is the ratio of the time scales of the phase variables and synaptic weights, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document}μ defines the sharpness of the plasticity boundary function. The limit \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu \rightarrow 0$$\end{document}μ→0 corresponds to a hard boundary. The analytical results obtained on the model of two oscillators are generalized for multi-leaf star networks. Multistability with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2^N$$\end{document}2N various asymptotic states is numerically demonstrated for one-, two-, three- and nine-leaf star-type networks.


Model
We consider a star network with N + 1 nodes described by Kuramoto-type phase oscillators: Here θ 0 is the phase of the central node (hub) and ω 0 is its natural frequency (spiking rate). The variables θ j and the parameters ω j for j = 1, . . . , N are the phases and the natural frequencies of the leaves, respectively. The hub is forced by all leaves, and each leaf is forced by the hub. The parameter A k indicates the synaptic weight of the directional link from the kth leaf to the hub, and B j is the synaptic weight of the directional link from the hub to the jth leaf. We modify the synaptic weights A j and B j in dependence on the phase difference between the hub and the jth leaf. For PDDP, we adopt STDP rules typical of excitatory synapses 6,15 . When the leaf phase is ahead of the hub phase ( ϕ j < 0 ), the synaptic weight A j increases, and the synaptic weight B j decreases. In the opposite case, ϕ j > 0 , the weight A j decreases, and the weight B j increases. More specifically, we modify synaptic weights using the time-continuous form of the PDDP rule introduced in Ref. 10 : Here ε ≪ 1 is a small but not vanishing parameter that takes into account the slow change in synaptic weights relative to the fast phase dynamics. The learning windows over which post-(pre-) synaptic spikes will cause synaptic potentiation (depression) are indicated as τ + ( τ − ). Following experimental evidences 6-8 , we assume τ − > τ + . In numerical experiments below, we choose 10 : τ − = 0.3 and τ + = 0.15 . The function F(x) is introduced to satisfy the boundary conditions, that is, to keep the synapses from achieving unrealistically large values or becoming inhibitory, namely 0 ≤ A j ≤ α and 0 ≤ B j ≤ α . Various types of boundary functions are considered in the literature. A soft boundary is associated with the linear function F(x) = x as in Refs. 12,13 , and a hard boundary is determined by the Heaviside step function H(x), F(x) = H(x) as in Refs. 14,15 . Interpolation between soft and hard boundaries can be achieved by a power function 16 with a sharpness parameter µ varying in the interval 0 < µ ≤ 1 . For µ = 1 this function corresponds to the soft boundary, and for µ → 0 it is the case of the hard boundary. In this paper, to analyze the transition to a hard boundary, we will mainly use the sigmoid function Like function (4), this function vanishes at x = 0 and monotonically increases with an increase of x. The parameter µ determines the sharpness of the transition from zero to one: F(x) ≈ x/µ for x ≪ µ and F(x) ≈ 1 for x ≫ µ . The hard boundary corresponds to the limit µ → 0 . Note that the function (5) is analytic at x = 0 for any µ , while the function (4) is non-analytic at x = 0 for 0 < µ < 1.
From the Eq. (1) we can derive equations for phase differences (2): Thus, the dynamics of a star-type plastic network with N leaves is determined by a closed system of 3N differential Eqs. (3) and (6). In this system, 2N Eq. (3) describe the slow adaptation of synaptic weights, and N Eq. (6) determine the fast variation of phase differences.

Methods
Two oscillators under PDDP. The case of two oscillators allows for analytical treatment. The analytical results obtained in this section are used to predict the asymptotic behavior of multi-leaf star networks. In the case of two oscillators, N = 1 and the system's state is characterized by three dynamical variables A 1 , B 1 and ϕ 1 . Since the variable ϕ 1 is fast with respect to the variables A 1 and B 1 , we can eliminate it from the system's equations. The slow variables A 1 , B 1 are governed by the Eq. (3) with j = 1 , and the Eq. (6) for the fast phase variable ϕ 1 reads where is the sum of the synaptic weights, and is the frequency detuning between the oscillators. Without loss of generality, we assume ω 0 > ω 1 , so the parameter is positive. Due to the small parameter ε , we can neglect the variations of the slow variables 38 A 1 and B 1 on the fast time scale of the variable ϕ 1 . For the constant K, the solution of Eq. (7) depends on the ratio K/� . When this ratio is greater than one, the phases of the oscillators are locked, and when it is less than one the oscillators are not synchronized. Below we consider these two cases separately. (i) Case K > � . Under this condition, the Eq. (7) has two fixed points, one of which is stable: It determines the phase difference of the locked oscillators in dependence of the variable K. Substituting this phase difference in Eq. (3) and taking into account that ϕ * 1 (K) ∈ (0, π) , we obtain the following approximate system to describe the slow dynamics of the synaptic weights A 1 and B 1 : The condition F(0) = 0 implies that this system has a fixed point Linear stability analysis of the system (11) shows that the fixed point (12) is stable due to F ′ (0) > 0 . This fixed point describes an asymptotic synchronized state of the system in which the oscillators become unidirectionally coupled. The directional link from the slower oscillator (with the natural frequency ω 1 ) to the faster oscillator (with the natural frequency ω 0 ) becomes disconnected, while the synapse associated with the directional link from a faster oscillator to the slower oscillator, reaches the maximum value α . The instantaneous frequencies of the synchronized oscillators become equal to the natural frequency of the faster oscillator: θ 1 =θ 0 = ω 0 . Since Eq. (11) are valid for A 1 + B 1 > � , the necessary condition for the existence of the synchronized state (12) is It is important to note that the above synchronized state also exists in the original (non-reduced) system of Eqs. (3) and (7) without the assumption that the parameter ε is small. Indeed, in the three-dimensional state space (A 1 , B 1 , ϕ 1 ) of this system there is always a stable fixed point (A * 1 , B * 1 , ϕ * 1 ) = (0, α, arcsin(�/α)) if the condition (13) is met. Thus, the synchronized state with disabled synaptic connection A 1 = 0 and maximum synaptic connection B 1 = α is exact for any values of the parameters ε and µ.
(ii) Case K < � . Under this condition, the oscillators are not synchronized and the phase difference ϕ 1 shows fast oscillations. Again, we neglect the slow variations of A 1 and B 1 on the fast time scale. For constant K, the oscillation period can be obtained from the Eq. (7) in analytical form: Now an approximate system of equations for the slow variables A 1 and B 1 can be obtained by averaging the Eq. (3) over the period of fast oscillations 39 : where (7) ϕ 1 = � − K sin(ϕ 1 ),  (8), (10), (14) and (16) constitute a complete system of reduced differential equations for describing the dynamics of the system in the entire plane ( A 1 , B 1 ) of slow variables. Unlike the original system of three differential Eqs. (3) and (7), here the fast variable ϕ 1 is eliminated. The line A 1 + B 1 = in the plane of slow variables divides the areas of synchronized and unsynchronized motion of oscillators (see the red dashed line in Fig. 1a). The region A 1 + B 1 > � corresponds to synchronized motion and is described by Eq. (11), while the region A 1 + B 1 < � represents unsynchronized motion and is defined by Eq. (15).
To justify the validity of the reduced system Eqs. (11) and (15), in Fig. 1 we compare the phase portrait in the plane of variables ( A 1 , B 1 ) obtained from the reduced system ( Fig. 1a) with the projection of the solution of the original system of three differential Eqs. (3) and (7) to the same plane ( A 1 , B 1 ) (Fig. 1b). For ε = 0.01 , we see excellent agreement between these two results. The red star is the fixed point (12) described above, which corresponds to the synchronized state of oscillators with a maximum synaptic weight B 1 and a vanishing synaptic weight A 1 . We see that there is another stable fixed point located in the unsynchronized region below the red dashed line marked by green square. This means that the system is bistable with coexisting synchronized and unsynchronized steady states. The basins of attraction of stable fixed points are separated by the stable separatrix of the saddle point, which is marked with a red circle. The saddle is in the unsynchronized region, just below the red dashed line. Strictly speaking, the term "fixed point" in the unsynchronized region is correct only within the framework of the reduced system of Eq. (15). In the original system of Eqs. (3) and (7), the variables A 1 and B 1 show high-frequency oscillations around the usynchronized steady state with a small amplitude proportional to ε . The period of these oscillations is determined by the Eq. (14). At ε → 0 the amplitude of these oscillations tends to zero. Numerical simulation of the original system of Eqs. (3) and (7) showed that convergence to asymptotic states does not depend on the initial conditions of the fast variable ϕ 1 for almost all initial conditions of the slow variables (A 1 , B 1 ) , except for the case when they are in the close vicinity of the stable separatrix of the saddle point.
Dependence of unsynchronized steady state on the sharpness parameter µ . The reduced system of Eq. (15) is convenient for analyzing the dependence of an unsynchronized steady state on the sharpness parameter µ of the boundary function. Especially in the limit of the hard boundary, µ → 0 , the results can be obtained in analytical form. For the power boundary function (4), the stationary solutions of Eq. (15) at any µ can be demonstrated graphically. By equating the right hand sides of Eq. (15) to zero, we obtain two equations, The stationary value of the variable K = A 1 + B 1 satisfies the equation K = f 1 (K) + f 2 (K) . Figure 2a shows the graphical solution of this equation for different values of the parameter µ . We see that the stationary point K * obtained as intersection of the curve f 1 (K) + f 2 (K) with the identity line moves towards zero as µ is decreased. .
(17) www.nature.com/scientificreports/ Since the synaptic weights A 1 and B 1 are always positive, the stationary values of A * 1 and B * 1 also move towards zero with decreasing µ . This is demonstrated in Fig. 2b, were the dependence of the stationary values A * 1 and B * 1 on the parameter µ is depicted. Thus in the limit of the hard boundary µ → 0 the fixed point responsible for the unsynchronized state of the oscillators approaches the origin of the phase space, . This means that in the limit µ → 0 , the oscillatorsbecomes completely disconnected. To verify whether this property is universal, in Fig. 2c we plot the dependence of the stationary values of A * 1 and B * 1 on the parameter µ for the case of the sigmoid boundary function (5). We see that here the stationary values of A * 1 and B * 1 also tend to zero for µ → 0. Now we generalize the above results for any continuous boundary function F(x), which for 0 ≤ x ≤ α has the following properties: (i) F(x) is a monotonically increasing function, i.e., F ′ (x) > 0 ; (ii) in the interval 0 ≤ x ≤ µ ≪ α , the function increases sharply from F(0) = 0 to F(µ) 1 , and in the interval µ ≤ x ≤ α , it slowly increases to a value F(α) ≈ 1 . In the examples above, we saw that for small µ the fixed point corresponding to the unsynchronized state was placed in the region K ≪ . Here we assume that this condition is valid for any boundary function with the above properties and will show that the obtained results confirm this assumption.
The assumption K ≪ allows us to simplify the Eq. (15) and treat the problem analytically. Under this assumption, the Eqs. (14) and (16) can be approximated as T(K) ≈ 2π/� and G ± (K, ϕ 1 ) ≈ exp (ϕ 1 /τ ± )/� , respectively. Substituting these expressions into Eq. (15) and evaluating the corresponding integrals, we get We see that the equations for A 1 and B 1 are identical and independent of each other. This means that their stationary solutions coincide A * 1 = B * 1 and it is enough to analyze only one of the two equations, say, the Eq. (18a). The stationary solution of Eq. (18a) satisfies: where is a parameter depending on PDDP learning windows τ + and τ − . At any τ + < τ − , this parameter is less than one. Specifically, for τ + = 0.15 and τ − = 0.3 used in this paper, the value of this parameter is q = 0.5 . For small µ , the function F(x) is close to 1 for all x except for a small region 0 ≤ x ≤ µ ≪ α , where it sharply increases from F(0) = 0 to F(µ) 1 . From this property it follows that for q < 1 the Eq. (19) can be fulfilled only at  (22) and (23), respectively, which are valid for small values of the parameter µ. www.nature.com/scientificreports/ where F (−1) denotes the inverse of the function F. This solution is stable since F ′ (A * 1 ) > 0 and F ′ (α − A * 1 ) > 0 . Due to the assumptions made for the function F(x), the value of F (−1) (q) tends to zero for µ → 0 . Therefore, for any boundary function with the above properties, the fixed point (A * 1 , B * 1 ) of the unsynchronized state approaches the origin of the phase plane when µ → 0 . Choosing a sufficiently small value of µ , we can always satisfy the condition K ≪ , which was used to derive the Eq. (18).
The general Eq. (21) can be used to estimate unsynchronized stationary states for specific boundary functions (4) and (5) presented above. In the case of the power boundary function F(x) = x µ , we get According to this law, the fixed point very quickly approaches the coordinate origin as µ decreases. In the case of a sigmoid boundary function F(x) = tanh(x/µ) , the fixed point approaches the origin linearly with decreasing µ: The dependencies (22) and (23) are shown by dashed curves in Fig. 2b,c, respectively. For small µ they are in good agreement with the results obtained by numerical solution of Eq. (15).
Numerical examples in Fig. 3(I) and (II) columns show the dynamics of the system state variables for two different initial conditions that converge to different asymptotic modes: the synchronized state (I) and the unsynchronized state (II). The results are presented for the sigmoid boundary function (5) with the parameter µ = 0.2 . Panels (a) and (b) show the dynamics of the slow variables A 1 (t) and B 1 (t) obtained from the original system of Eqs. (3) and (7), as well as the reduced system of Eqs. (11) and (15). For ε = 0.001 , these results are indistinguishable in the figure. In panel (a), the synaptic weight A 1 asymptotically vanishes ( A 1 → 0 ), and the synaptic weight B 1 reaches the maximum value ( B 1 → α = 1 ). The phases of the oscillators remain locked during this process, as can be seen from the dynamics of the phase difference ϕ 1 (t) , shown in panel (c). In panel (b), both synaptic weights approach small values close to the value of the parameter µ = 0.2 . Here the dynamics of The phases are locked when the sum K = A 1 + B 1 of the synaptic weights is greater than the frequency detuning = 0.5 . When K becomes smaller than , the phase www.nature.com/scientificreports/ difference experiences fast oscillations. Panels (e) and (f) demonstrate the accuracy of the reduced Eqs. (11) and (15) depending on the parameter ε . The dynamics of the absolute value of the deviation δK(t) = K red (t) − K(t) in semi-log plot is shown for different values of the parameter ε : 0.1 (blue curves), 0.01 (red curves), and 0.001 (orange curves). Here K(t) is the sum of the synaptic weights A 1 (t) + B 1 (t) obtained from the original system of Eqs. (3) and (7), and K red (t) is the same variable obtained from the reduced system of Eqs. (11) and (15). We see that the deviation |δK| decreases in proportion to the ε parameter. Figure 3g shows the dynamics of the variable K(t) obtained from the original system of Eqs. (3) and (7), for fixed ε = 0.001 and different values of the parameter µ : 0.2 (red curve), 0.1 (blue curve), 0.01 (green curve), and 0.001 (black curve). These results are presented for the sigmoid boundary function F(x) = tanh(x/µ) . The lower yellow curve shows the result for the Heaviside step boundary function, F(x) = H(x) . The initial conditions are the same as in the (II) column, which correspond to the unsynchronized asymptotic state of the system. As expected from the above analytical results [see Eq. (23)], the asymptotic values of the variable K decrease in proportion to the µ parameter. Note that Eq. (23) was derived in the limit ε → 0 and predicts stationary asymptotic values of K. Here we observe small amplitude oscillations in the asymptotic dynamics due to the finite value of ε . The oscillations with a small amplitude of order ε remain also in the case of the Heaviside step boundary function (see the lower yellow curve in the figure).
To conclude this section, we formulate the main results concerning the asymptotics of two coupled oscillators with PDDP. Below, we will be interested in the case of small parameters ε and µ , since this case allows to predict the number of possible states in a star network with an arbitrary number of leaves. In this case the system of two oscillators can demonstrate bistability with coexisting synchronized and unsynchronized attractors. The synchronized attractor is characterized by a unidirectional coupling, in which directional link from a slower oscillator to a faster oscillator is disconnected, A 1 = 0 , and the synapse associated with the directional link from a faster oscillator to a slower oscillator is maximal, B 1 = α . This attractor exists only when the maximum value α of the synaptic weight is greater than the frequency mismatch between the oscillators. The unsynchronized attractor is characterized by low values of both synaptic weights A 1 and B 1 . In the limit µ → 0 and ε → 0 , the oscillators in the unsynchronized state are completely disconnected, A 1 = B 1 = 0 . In the next section, we will use these results to predict the asymptotic behavior of multi-leaf star networks.

Results
We show that a star network with N leaves described by Kuramoto-type phase oscillators Eq. (1), with synaptic weights satisfying the PDDP rule (3), can evolve into 2 N various stable configurations. This result was obtained under the assumption of small parameters ε and µ . We use the analytic results obtained for two oscillators to predict the asymptotic behavior of an arbitrary star network. A generalization for an N-leaf star network is obtained by gradually enlarging the number of leaves. We prove our prediction numerically for the cases of two, three, and nine leaves.
For simplicity, we assume that there are no oscillators with coinciding natural frequencies. Then, without loss of generality, let us number the natural frequencies of the leaves in ascending order: We will consider all possible options for choosing the natural frequency ω 0 of the hub, which can fall into different frequency intervals of the leaves, i.e., it can satisfy N + 1 various inequalities: ω 0 < ω 1 , ω j < ω 0 < ω j+1 for j = 1, . . . , N − 1 , and ω N < ω 0 . Finally, we assume that the limiting value α of synaptic weights is greater than the frequency mismatches � j = |ω 0 − ω j | between the hub and all leaves: The last inequalities are a generalization of the necessary condition (13) for the existence of synchronized states.
Two-leaf star network. We begin our analysis with a two-leaf star network. We assume that each leaf in the network can be either synchronized or unsynchronized with the hub. We also assume that the fastest oscillator always wins in the synchronization process, so that unidirectional links from the fastest oscillator to the slower oscillators are established with the maximum synaptic weights. Since the fastest oscillator dominates the slower oscillators, no more than one synchronized link can go to the hub. Using these assumptions, we constructed seven possible configurations, shown in Fig. 4a-g. Solid arrows correspond to the synchronized states of the oscillators, in which the connection is unidirectional with the maximum synaptic weight. Dotted arrows correspond to unsynchronized states between oscillators, in which synaptic weights in both directions are small. We assume that they vanish in the limit ε → 0 and µ → 0 . The presented configurations can be realized only when certain frequency inequalities are met. These inequalities are written above each configuration. The inequalities were obtained using the Eqs. In configuration (a), the coupling weights A 1 and A 2 are zeros, and both leaves are synchronized with the hub. In this case, the phases of the oscillators satisfy the equations: www.nature.com/scientificreports/ The hub is not influenced by leaves, so the instantaneous frequencies of the leaves must be equal to the natural frequency of the hub, θ 1 =θ 2 = ω 0 , to ensure synchronization. This leads to two equations for the phase differences as defined by Eq. (2) [ ϕ 1,2 = θ 0 − θ 1,2 ]: sin ϕ 1,2 = ω 0 − ω 1,2 /B 1,2 . From Eq. (3a) we conclude that the phase differences ϕ 1,2 must be in the interval (−π, 0) in order to satisfy the conditions A 1 = A 2 = 0 . Then from the stationary solution of the Eq. (3b), we get B 1 = B 2 = α , and from the above equations for the phase differences, we obtain two inequalities: ω 1 < ω 0 and ω 2 < ω 0 . In Fig. 4a-g we have presented only the second inequality corresponding to this configuration, since the first inequality holds due to the assumption (24). As a second illustrative example, consider configuration (c). Here, as in the previous example, all oscillators are synchronized, but now A 1 and B 2 are zeros, and the Eq. (1) for the phases can be written as: The second leaf is not influenced by the hub and it oscillates freely with its natural frequency ω 2 . To ensure synchronization between all oscillators, we must require θ 0 =θ 1 = ω 2 . As a result, we get two equations for the phase differences ϕ 1,2 : sin (ϕ 1 ) = (ω 2 − ω 1 )/B 1 and sin (ϕ 2 ) = (ω 0 − ω 2 )/A 2 . It follows from the first equation that the phase difference ϕ 1 is in the interval (0, π) , since ω 2 > ω 1 . This provides A 1 = 0 and B 1 = α as stationary solutions to the Eqs. (3a) and (3b). To satisfy the condition B 2 = 0 , we must require that the phase difference ϕ 2 is in the interval (−π , 0) . Then from the stationary solution of the Eq. (3b) we get A 2 = α , and from the above equations for the phase difference ϕ 2 we obtain the inequality ω 0 < ω 2 , which gives the necessary condition for this configuration.  3) and (6). The results presented in three rows correspond to three different possibilities of the hub frequency to fall into different frequency intervals of the leaves: ω 0 < ω 1 (top row), ω 1 < ω 0 < ω 2 (middle row), and ω 2 < ω 0 (bottom row). The top row of (h) (from left to right) match configurations (g,b,e,a). For the middle and bottom rows, the corresponding configurations are, respectively, as follows: (g,d,e,c) and (g,d,f,c). All frequencies are taken from the array (0.55, 0.85, 1). The parameter ε = 0.001 , and the parameter µ = 0.01 corresponds to the sigmoid boundary function (5). The values of the parameters α , τ + and τ − are the same as in Fig. 1 www.nature.com/scientificreports/ Our third example is configuration (e). Here, the first leaf is synchronized with the hub and the second leaf is not synchronized. Considering that A 1 = 0 and the weights A 2 and B 2 are small, the Eq. (1) for the phases take the form: Here O(A 2 ) and O(B 2 ) are small terms proportional to the weights A 2 and B 2 respectively, which we assume to be zero for ε → 0 and µ → 0 . If we neglect these terms, then we will see that the hub and the second leaf freely oscillate with their natural frequencies. The synchronization condition of the first leaf with the hub, θ 1 = ω 0 , leads to the equation sin (ϕ 1 ) = (ω 0 − ω 1 )/B 1 . From the Eq. (3a) it follows that the phase difference ϕ 1 must be in the interval (−π , 0) in order to satisfy the condition A 1 = 0 . Then from the stationary solution of the Eq. (3b), we get B 1 = α , and from the above equation for the phase difference ϕ 1 , we obtain the inequality ω 1 < ω 0 , which is the necessary condition for the realization of configuration (e).
For further classification of asymptotic states, it is convenient to introduce the following notations. We encode each asymptotic configuration of an N-leaf network as an array of N characters (s 1 s 2 . . . s N ) , where s j defines the state of the jth leaf. Each element s j of the array is represented by one of three possible symbols: "0" (zero), " 1 H " (one with subscript H), and " 1 L " (one with subscript L). The symbol "0" denotes an unsynchronized state of the leaf with small (vanishing for ε → 0 and µ → 0 ) synaptic weights. The symbols " 1 H " and " 1 L " represent the synchronized state with maximum coupling strength directed to the hub and leaf, respectively. The corresponding codes are written under each configuration in Fig. 4.
All asymptotic configurations shown in Fig. 4a-g were numerically validated for different sets of hub and leaves frequencies. Using these configurations, we can predict all possible asymptotic states of a two-leaf network for any given value of the hub frequency when it falls into different frequency intervals of the leaves. When the hub frequency is less than the frequencies of both leaves, ω 0 < ω 1 , then, depending on the initial conditions, the configurations (c), (d), (f) and (g) can be asymptotically attained. When the hub frequency falls within the interval between leaf frequencies, ω 1 < ω 0 < ω 2 , the admissible asymptotic configurations are (c), (d), (e), and (g). Finally, when the hub frequency is greater than the frequencies of both leaves, ω 2 < ω 0 , we identify (a), (b), (e), and (g) as admissible configurations. Thus, at any value of the hub frequency, a two-leaf star network can asymptotically reach four different configurations, characterized by different values of synaptic weights.
In Fig. 4h, we confirm the above statement with a specific numerical example. Taking different initial conditions, we integrated the system of Eqs. (3) and (6) over a long period of time until the values of the weights A j and B j became saturated, and then displayed these values in color code. The results presented in three rows correspond to three different possibilities of the hub frequency to fall into different frequency intervals of the leaves: ω 0 < ω 1 (top row), ω 1 < ω 0 < ω 2 (middle row), and ω 2 < ω 0 (bottom row). In all cases, four different asymptotic states were found, as predicted above. On each row, we order the asymptotic states so that we get similar sequences of configuration codes. If we remove the subscripts at " 1 H " and " 1 L " in these codes, we get an identical sequences of configuration codes for each row: (0 0) , (0 1) , (1 0) , and (1 1) . This is a binary sequence of natural numbers from 0 to 3, which gives us a convenient way to index the various asymptotic configurations for any given hub frequency. We will use this way of numbering different asymptotic configurations when considering star networks with a larger number of leaves.
Three-leaf star network. The classification of various asymptotic states in a three-leaf star network can be done in the same way as for a two-leaf star network. Using arguments similar to those given above, we have constructed fifteen possible asymptotic configurations of the three-leaf star network, which are shown in Fig. 5. They are designated by letters from (a) to (o). The three characters in brackets below each configuration represent its code. The required frequency inequality is indicated above each configuration. Figure 5 allows predicting all possible asymptotic states of a three-leaf star network for any given value of the hub frequency when it falls into different frequency intervals of the leaves. In all cases, there are eight different admissible asymptotic states presented in Table 1. The first column lists four different possible inequalities for the hub frequency, and the next eight columns show the codes of the corresponding allowed configurations. For example, when the hub frequency falls within the interval between the frequencies of the second and third leaves, ω 2 < ω 0 < ω 3 , the codes listed from left to right in the third row correspond to configurations (o), (h), (i), (g), (m), (f), (j), and (e) in Fig. 5. We order the configurations in the same way as in the case of a two-leaf star network. After removing the subscripts at " 1 H " and " 1 L " in the configuration codes, we get identical binary sequence of natural numbers from 0 to 7 in all rows of the Table 1. We use this numbers to index the corresponding asymptotic configurations for any given hub frequency.
All asymptotic configurations presented in Table 1 were numerically validated for different sets of hub and leaves frequencies. An example of numerical testing of asymptotic states of a tree-leaf star network, when the hub frequency is in the interval between the frequencies of the second and third leaves, ω 2 < ω 0 < ω 3 , is shown in Fig. 6a. By integrating the system of Eqs. (3) and (6) with many different initial conditions over a long period of time, we obtained eight different asymptotic configurations, as predicted in the third column of the Table 1. The probability distribution of the various asymptotic configurations, numbered as described above, is shown in Fig. 6b. This graph was constructed using 1000 randomly selected initial conditions. The initial values of the slow variables A j and B j were randomly and independently taken from uniform distribution [0, α] . Each of the 1000 numerical experiments resulted in one of eight predicted asymptotic configurations, and no other asymptotic  (Fig. 5e) is most likely, and configuration (0 0 0) with all unsynchronized oscillators (Fig. 5o) is the least likely. N-leaf star network. The results obtained for two-and three-leaf star networks can be generalized for a star network with an arbitrary number of leaves. Assuming that each leaf in a N-leaf star network can be either synchronized or unsynchronized with the hub, we can construct 2 N different asymptotic configurations for any given hub frequency. Let us denote the code of the nth asymptotic configuration of an N-leaf star network as C (k) n , where n is a natural number that varies from 0 to 2 N − 1 , and the superscript (k) identifies a specific frequency range of the hub. We assign k = 1 when the hub frequency is less than the lowest frequency of the leaves ω 0 < ω 1 . The values k = 2, . . . , N correspond to the inequalities ω k−1 < ω 0 < ω k . Finally, we assign  Table 1. Theoretically predicted asymptotic configurations of a three-leaf star network for four different possibilities of the hub frequency to fall into different frequency intervals of the leaves. The first column lists the possible inequalities for the hub frequency, and the next eight columns show the codes of the corresponding allowed configurations. www.nature.com/scientificreports/ k = N + 1 , when the hub frequency is greater than the highest frequency of the leaves, ω N < ω 0 . Summarizing the results of sections dedicated to two and three leaf networks, we found that the code C n . To demonstrate these rules, let's write the code for the third asymptotic configuration of a three-leaf star network when the hub frequency is in the range ω 1 < ω 0 < ω 2 . In this case k = 2 , N = 3 and n = 3 . The number n in binary is (0 1 1) . Assigning the rightmost digit 1 in the positions p ≥ 2 to the subscript H, and the remaining digits 1 to the subscript L, we obtain the code C (2) 3 = (0 1 L 1 H ) , which is the same as shown in the second row of the fifth column of Table 1.
As an example of applying the above rules to a more complex case, we present the code of the 25th asymptotic configuration of a five-leaf star network when the hub frequency is in the range ω 4 < ω 0 < ω 5 : To verify the validity of the predicted asymptotic states in a star network with a large number of leaves, we will consider the dynamics of a 2N-dimensional state vector in the complete phase space of slow variables. By integrating the system of Eqs. (3) and (6) with different initial conditions, we will check whether the state of the network approaches the predicted asymptotic states. The state vector of the nth predicted asymptotic configuration, determined by its code C (k) n = (s 1 s 2 . . . s N ) is constructed as follows. In the first half of the vector components, we assign A * j = α if the corresponding character s j in the code is " 1 H ", and we assign A * j = 0 if the character s j is " 1 L " or "0". In the second half of the vector components, we assign B * j = α if the character s j is " 1 L ", and we assign B * j = 0 if the character s j is " 1 H " or "0". For example, the state vector of the 25th predicted asymptotic configuration in a five-leaf star network defined by the code Eq. (29) is: In Fig. 7, we confirm the validity of the predicted asymptotic states for a nine-leaf star network when the hub frequency is in the range ω 8 < ω 0 < ω 9 . Panel (a) shows the results obtained with the sigmoid boundary function (5) at µ = 0.01 , and panel (b) corresponds to the Heaviside step boundary function. To numerically check the existence of the predicted asymptotic configurations in the system of Eqs. (3) and (6), we prepared 512 different initial states R(0) so that they are at the same small distance |R(0) − R (9) n | = 0.05 from different    (6) with different initial conditions. The leaves natural frequencies are (ω 1 , ω 2 , ω 3 ) = (0.55, 0.7, 1). The hub natural frequency ω 0 = 0.85 satisfies the inequality ω 2 < ω 0 < ω 3 . The resulting distributions are labeled with codes that agree with the theoretically predicted asymptotic configurations shown in the third column of Table 1. (b) Probability distribution of eight different asymptotic configurations for a three-leaf star network, when the hub frequency is in the interval between the frequencies of the second and third leaves, ω 2 < ω 0 < ω 3 . In panel (a), the configurations are numbered according to the decimal representation of the binary codes shown above each pattern in panel (a) (see main text for details). The distribution is constructed using 1000 randomly selected initial conditions for Eqs. (3) and (6). The initial values of the weights A j and B j were randomly and independently taken from the uniform distribution [0, α].  (9) n , n = 0, . . . , 511 . These initialdistances are depicted asblue squares on the top horizontal line. The evolution of the distances obtained by integrating the Eqs. (3) and (6) for the initial conditions chosen in this way is presented by two snapshots. The yellow dots show the values of the corresponding distances |R(t) − R (9) n | at time t = 300 , and the red circles at time t = 76,000 . We see that for both boundary functions, the distances from the predicted states decrease with time, which means that these states are stable asymptotic network configurations. However, the distances (except for the case n = 511 ) do not vanish even after a fairly long time. This is because we integrate the Eqs. (3) and (6) for finite values of the parameters µ and ε , while our theoretical prediction of asymptotic states is based on assumption µ → 0 and ε → 0 . Exceptional configuration n = 511 corresponds to the state with all synchronized oscillators. The distance from this state tends to zero, since this state is an exact solution of the Eqs. (3) and (6) even for finite µ and ε . All other predicted asymptotic states are approximate solutions of the system of Eqs. (3) and (6). Note that the Heaviside boundary function results in better match of the asymptotic solutions of the Eqs. (3) and (6) with the predicted ones than the sigmoid function; the final distances |R(t) − R (9) n | at time t = 76,000 in Fig. 7b are smaller than in Fig. 7a. This is due to the fact that in the case of the sigmoid boundary function, both parameters µ and ε are finite, and in the case of the Heaviside boundary function only ε is finite.

Discussion
We analyzed the multistability in a star network of Kuramoto-type phase oscillators in the presence of phasedifference-depended plasticity. A star-type network is the simplest network model that captures many of the typical properties of real-world networks; it can be considered an essential building block of neural networks [29][30][31] . The relative simplicity of the network allowed us to accurately estimate the number of possible asymptotic configurations which can be attained during the plastic evolution of the network. We have shown that, depending on the initial conditions, an N-leaf star plastic network can evolve into 2 N various stable configurations characterized by different values of the coupling strength between the hub and the leaves. This result is independent of the relative magnitude of the natural frequency of the hub in relation to the natural frequencies of the leaves. We classified the various asymptotic states of the network by introducing configuration codes. The code is an array of N symbols that define the asymptotic state of each leaf. Each leaf can be synchronized or unsynchronized with the hub.
Results for an arbitrary N-leaf star network were obtained by generalizing the results derived from the analysis of networks with a small number of leaves. We started our analysis from a simple model consisting of two oscillators. Using the small parameter ε , which determines the ratio of the characteristic time scales of phase variables and synaptic weights, we obtained a reduced system of two differential equations for slowly varying synaptic weights. These equations were treated analytically. We have shown that this system is bistable with coexisting synchronized and unsynchronized attractors. The synchronized attractor is characterized by a unidirectional coupling, in which the directional link from a slower oscillator to a faster oscillator is disconnected, and the synapse associated with the directional link from a faster oscillator to a slower oscillator is maximal. The unsynchronized attractor is characterized by low values of both synaptic weights. In the limit of the hard boundary, when the sharpness parameter µ of the PDDP boundary function tends to 0, the oscillators in the unsynchronized state are completely disconnected. Using these results, we predicted all possible asymptotic states for a two-and three-leaf star networks and confirmed our prediction with numerical examples. Next, we generalized our prediction for an arbitrary N-leaf star network. In the limit of two small parameters ε → 0 and µ → 0 , we developed an algorithm for constructing 2 N different asymptotic configurations. As an example, we demonstrated the applicability of this algorithm to a nine-leaf star network. Using our algorithm, we predicted  (6) for a nine-leaf star network with 512 different initial conditions R(0) , each of which is close to the state R (9) n of a particular predicted asymptotic configuration with number n = 0, . . . , 511 . Panels (a) and (b) correspond to the sigmoid boundary function with µ = 0.01 and the Heaviside step boundary function, respectively. The frequencies (ω 1 , . . . , ω 8 , ω 0 , ω 9 ) , written in ascending order, are equidistantly distributed in the interval [0.6, 1]. The states R(0) are chosen so that the initial distances |R(0) − R (9) n | shown in blue squares are the same for all configurations. The yellow dots show the values of the corresponding distances |R(t) − R (9) n | at time t = 300 , and the red circles at time t = 76,000 . Parameter values: ε = 0.001 , τ + = 0.15 , τ − = 0.3 , and α = 1. www.nature.com/scientificreports/ 512 different asymptotic configurations and numerically proved their existence in a network model with small values of the parameters ε and µ. Multistability of neuronal networks is relevant to a number of applications. For instance, as shown numerically in Refs. 9,40-42 , plasticity-mediated multistability enables desynchronizing and/or decoupling stimulus patterns to move neuronal networks from dynamic states characterized by abnormally strong synchrony and correspondingly increased synaptic weights to dynamic states with reduced synchrony and reduced synaptic strengths. Coordinated Reset (CR) stimulation 43 is a spatio-temporally patterned desynchronization stimulation protocol which was computationally designed to induce cumulative and long-lasting therapeutic desynchronizing effects by shifting neuronal networks from abnormally synchronized attractors with strong synaptic weights to desynchronized attractors with reduced synaptic weights 9,44 . Cumulative and long-lasting therapeutic and desynchronizing effects were observed in Parkinson's disease (PD) patients 45 and Parkinsonian monkeys 46,47 treated with CR-deep brain stimulation as well as in PD patients treated with vibrotactile CR stimulation delivered to the fingertips 48,49 . By the same token, acoustic CR stimulation induced cumulative and long-lasting reduction of tinnitus loudness and tinnitus annoyance along with a reduction of tinnitus-related abnormal cortical synchrony 50,51 .
Our analysis strategy, extending analytical results obtained in small systems to larger networks, enables thorough analysis of plastic mechanisms of high-dimensional dynamical systems, as e.g. also performed by applying two-neuron loop analysis in the context of recurrent networks of oscillatory neurons with propagation delays 11 . As a further plan of this study, it is interesting to check the universality of the results obtained. The question is whether the number of different asymptotic states in a star network remains the same if we replace phase oscillators with realistic neuron models and use the more precise STDP rule. A computationally attractive model for this problem would be a star network of integrate-and-fire neurons, which enables event-driven simulations 20 . It is also interesting to extend this study to more complex network topologies. The existence of desynchronized states with low coupling weights allows us to assume that our model will lead to the appearance of multi-cluster states in complex network topologies, as was observed in Refs. 52,53 . The effect of noise is also of great interest as it may induce new states in the network 27,54 . License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.