Winner-take-all in a phase oscillator system with adaptation

We consider a system of generalized phase oscillators with a central element and radial connections. In contrast to conventional phase oscillators of the Kuramoto type, the dynamic variables in our system include not only the phase of each oscillator but also the natural frequency of the central oscillator, and the connection strengths from the peripheral oscillators to the central oscillator. With appropriate parameter values the system demonstrates winner-take-all behavior in terms of the competition between peripheral oscillators for the synchronization with the central oscillator. Conditions for the winner-take-all regime are derived for stationary and non-stationary types of system dynamics. Bifurcation analysis of the transition from stationary to non-stationary winner-take-all dynamics is presented. A new bifurcation type called a Saddle Node on Invariant Torus (SNIT) bifurcation was observed and is described in detail. Computer simulations of the system allow an optimal choice of parameters for winner-take-all implementation.

We use generalized phase oscillators as the elements of the WTA system. Phase oscillators of the Kuramoto type 52 have been widely used for describing dynamics of Josephson-junction arrays, neutrino flavor oscillations, semiconductor laser arrays, coupled magnetic systems, and neural networks. Reviews of the mathematical theory of phase oscillator systems and their applications can be found in the publications [53][54][55][56] . Systems with phase oscillators and radial connection architectures have been studied in the papers [57][58][59][60][61][62][63] . Conventionally, a phase oscillator is described by a single variable, its phase. The natural frequencies of the oscillators and coupling strengths are the parameters of the system. Generalized phase oscillators differ from original Kuramoto oscillators by the transformation of some parameters of the system into dynamical variables. In the case of the system considered here the variables include the natural frequency of the CO and the strengths of connections from POs to the CO. The natural frequency of the CO is adapted in the direction of its current value. The connection strength of a PO is adapted as a function of the similarity between its phase and the phase of the CO. To obtain the WTA regime the connection strengths from POs to the CO should be positive and the connection strengths from the CO to POs should be negative.
The advantage of our WTA phase oscillator system is the possibility of implementing it in hardware as a laser optical device or as a Josephson junction array. Note that there are no connections between POs, so the whole number of connections in the system is 2n.
WTA systems of generalized phase oscillators with a central element have previously been applied to attention modeling and visual search 64,65 . Here we present for the first time a rigorous mathematical analysis of their dynamics.
We develop a new mathematical theory of WTA oscillatory networks. In the case of POs with identical natural frequencies we derive conditions when only one PO wins the competition. This winner is synchronized with the CO, while all other POs are in antiphase to the CO. These WTA dynamics correspond to the existence of a stable equilibrium in the phase space of the system.
Using perturbation theory and bifurcation analysis, we demonstrate that in the case of non-identical natural frequencies of POs the system can demonstrate both stable and oscillatory versions of the WTA regime depending on the relations between the parameters of the system. In the stable case, the winning PO works nearly inphase with the CO, while the absolute values of the difference between the phases of other POs and the CO exceed π/2. In the oscillatory version of WTA there is a single winner which works nearly inphase with the CO, while the phases of other POs oscillate far from the phase of the CO or run indefinitely in the positive or negative direction. We show that the appearance of the oscillatory WTA regime is due to a Saddle Node on Invariant Torus (SNIT) bifurcation which is a generalization of the well-known Saddle Node on Invariant Circle (SNIC) bifurcation (see 66,67 ). Note that the SNIT bifurcation of two cycles on a 2-dimensional torus was studied among others by C. Baesens et al. 68,69 . The same type of bifurcations has been addressed in considerable detail for the 3-dimensional case, though without the reinjection, by A. Chenciner [70][71][72] and C. Baesens and R. S. MacKay 73 .
Finally, using massive simulations of the model we clarify how the parameters of the system affect the results of WTA.

Model Formulation
The system that we consider contains a central oscillator (CO) and a set of n peripheral oscillators (POs). The CO interacts with POs through feedforward and feedback connections. There are no lateral connections between POs. The dynamics of the system are described by the following equations where θ θ θ ω … … ∈ × × +    a a ( , , , , , , ) n n n n 0 1 0 1 1 are the variables, ω i , i = 1, …, n, α, β, γ, b, c are the parameters, In (1)-(4) θ i are the current phases of the oscillators, ω i are the natural frequencies of the oscillators, a i and b are the connection strengths between the oscillators. We will also associate a i with the amplitudes of oscillations of POs. Positive values of connection strengths correspond to synchronizing interaction and negative values correspond to desynchronizing interaction. We always suppose that α, β, γ, c > 0.
The meaning of equations (1)-(4) will be explained a bit below. Right now let us introduce the restrictions on the functions used in these equations. The functions f, g, h are assumed to be 2π-periodic and satisfy the following conditions 0, ( ) 0, Scientific REPORTS | (2018) 8:416 | DOI: 10.1038/s41598-017-18666-3 1, ( ) 0, (0) ( ) 0 (7) Thus the functions f and g are odd and the function h is even. Periodicity and oddness of the functions f and g imply the conditions f(0) = f(π) = g(0) = g(π) = 0. We assume that the functions f(x), g(x) do not have other zeros in the interval (0,π). We also require that h(x) is monotonic on [0, π] (and [−π, 0]). Equations (1)-(4) can be considered as a generalization of standard Kuramoto equations for phase oscillators 52 by introducing in addition to two phase equations an equation for the adaptation of the natural frequency of the CO ω 0 and the equations for the adaptation of the amplitudes a i . The meaning of equation (3) According to this equation the natural frequency of the CO is adapted in the direction of the current frequency. The parameter α controls the speed of adaptation. The meaning of equation (4) is that there is resonant increase of the amplitude of oscillations of the i-th PO to the level c + γ if this PO works inphase with the CO, otherwise this amplitude decreases to a low level c. The parameter β controls the speed of amplitude adaptation.
System (1)-(4) is similar to the systems invented in 65 with the aim to organize a competition between POs for the synchronization with the CO in such a way that in a typical case only one PO can win the competition. We will show that this is possible if b < 0. Other types of dynamics including multistable and chaotic regimes are also possible depending on the values of the parameters.
Two types of dynamics can be associated with the WTA regime in system (1)-(4). In the stationary case we say that the i-th PO wins the competition if it is the only PO that asymptotically has the amplitude of oscillations equal to c + γ, while the amplitudes of other POs are equal to c. In the non-stationary case the amplitudes of POs are not constants anymore. We say that the i-th PO wins the competition if its amplitude is asymptotically concentrated in a region above the value c + γ − δ, while the amplitudes of other POs vary in a region below the value c + δ, where δ is a small number. Our aim is to find conditions when stationary and non-stationary WTA regimes take place. For this reason we investigate the dynamics of system (1)-(4) in two cases, corresponding to identical and non-identical POs.
Some results about the dynamics of (1)-(4) can be obtained for the general form of the functions (5)-(7). More advanced and complicated results demand a specification of the functions f, g, h. In this case the following types of the functions will be used: with the parameters ν > 1,  σ 1, μ π ∈ (0, ) (see Fig. 1). The function f(x) has a single maximum in the interval (0, π) and a symmetric minimum in (−π, 0) such that Thus, increasing the value of the parameter ν we can continuously move the location of the maximum (and minimum) point from π/2 (and −π/2) to zero ( Fig. 1(a)). The parameter ν also specifies the slope of the function f at the zero point: One can check that f′(π) = 0 for ν > 1, that is the last condition of (5) is satisfied. The parameter σ controls the "width" of the function h. Increased values of σ make sharper the peak at h(0) ( Fig. 1(b)). Below we will show how the values of the parameters ν and σ influence WTA results.

Using phase differences
i i i Identical peripheral oscillators. Consider a symmetric case of equal natural frequencies of POs i Then system (12)- (14) obtains the form System (16) has permutation symmetry S n : the permutation of any two pairs (ϕ i , a i ) and (ϕ j , a j ) does not change the system. We denote , , , 0, , k k n k Using this notation we can formulate the following statement.
Since the functions f, g do not intersect the abscissa in (0, π), system (16) has 2 n fixed points P k corresponding to different values of k.
For fixed k there are C n k symmetric points P k whose multiplication is conditioned by permutation symmetry. The point P 1 corresponds to inphase synchronization of the CO with exactly one PO, while other POs are in antiphase to the CO. This case can be considered as a WTA procedure when one PO wins the competition for the synchronization with the CO. This PO is called the winner, while other POs are called losers.
The following statement describes stability of fixed points presented in Proposition 1.

Proposition 2.
System (12)- (14) with identical POs can have either the stable equilibrium P n (full synchronization) for b > 0 or n stable fixed points P 1 (WTA procedure) together with one stable point P 0 for b < 0. Other 2 n − n − 2 fixed points P k , k = 2, …, n − 1, are unstable points (saddles) for any values of the parameters. The point P n is stable if n ≥ 2 and Proposition 2 is proved in Appendix. For functions (8)-(10) with f′(0) = ν and g′(0) = 1 the conditions of stability for the point P 1 (19) become simpler If b < 0 and the number of POs n is large enough then the system does not have any stable point P 1 .
The system has the Andronov-Hopf (AH) bifurcation at the point P k when k(c + γ)f′(0) + nb = 0 and it has the pitchfork (PF) bifurcation at the point P k when α(c + γ)f′(0) = 0. The AH bifurcation implies the existence of stable limit cycles around the former fixed points P 0 , P 1 , and P n , in other cases the limit cycles are of the saddle type. The PF bifurcation changes the dimension of the stable and unstable invariant manifolds for P 1 , …, P n−1 , but it can not make these points stable.
System (16) of identical oscillators has the invariant manifold for arbitrary values of other parameters. This manifold corresponds to full synchronization. The dynamics on this manifold are described by 3D system (16) for n = 1. This system has two fixed points (0, ω, c + γ) and (π, ω, c) that correspond to inphase and antiphase synchronization of the CO with the PO. The first point is stable if (c + γ)f′(0) > −bg′(0), α(c + γ)f′(0) > 0. These conditions can be satisfied for both positive and negative values of b. The second point is stable in two directions when b < 0 and neutral in the third direction (one eigenvalue of the point is zero). According to the permutation symmetry S n system (16) has also the invariant manifolds that correspond to the regime when two oscillators (i-th and j-th POs) form a synchronous cluster. In a similar way the system can have a cluster with k synchronous POs for any k = 2, …, n.

Remark 1.
The results presented in Proposition 2 correlate with our previous results 60 about a star-like phase oscillator model with constant amplitudes of oscillations a i (t) = a = const. It has been shown that in this model different stable regimes Φ k coexist for different k = 1, … n. The additional condition f′(π) = 0 insures that only three types of stable regimes are possible: global synchronization Φ n , WTA regime Φ 1 , and no-winner regime Φ 0 (all other points Φ k are saddles).

Non-identical peripheral oscillators. Stationary solution.
Equations for connection strengths (14) are linear non-homogeneous differential equations that have the solutions a t c a t c e e e h d ( ) The assumption 0 ≤ h(x) ≤ 1 implies that The following proposition presents the boundaries for amplitude variables. Proposition 3. For any ε > 0 there is a large enough moment of time t 1 such that for any t > t 1 i This inequality is valid for any initial conditions and parameter values. Let us start from the symmetric case when all natural frequencies are equal, ω i = ω. In this case the system has n different symmetric fixed points P 1 (see the previous subsection) which are stable for some parameter values (Proposition 2). To distinguish between these points, the following notation is helpful: The perturbation of natural frequencies leads to the displacement of the points P 1,l . It also breaks the S n symmetry of these points (the coordinates of the points are non-symmetric in the general case), but it cannot change the stability of the points because they are hyperbolic. Denote Q l a perturbed location of the point P 1,l when ω l = ω + Δ l , l = 1, …, n, Δ l are relatively small: |Δ l | ≤ |b|. The point Q l describes the WTA case when the l-th PO is the winner.
The following statement describes the stability of different WTA points Q l in the case when the functions f, g, h are specified as (8)-(10) with large ν and μ < π/2.
has n fixed points Q l that correspond to the synchronization of the l-th PO with the CO. In this case ω 0 ≈ ω l , that is the CO obtains the frequency near the frequency of the "winning" PO. Moreover, these oscillators work nearly inphase, while loser POs are radically incoherent with the CO.
The coordinates of Q l can be approximately represented as l l l n l l

Non-identical peripheral oscillators. Non-stationary solution.
If the natural frequencies of POs are distributed in a large range, the CO has not enough strength to phase-lock all the POs. A stationary solution that describes WTA type of dynamics becomes impossible. Nevertheless, the main feature of the WTA regime when one PO works nearly inphase with the CO, while other POs are out of phase can survive even in this case. We consider here the solutions for which the variables ϕ l , ω 0 , a 1 , …, a n are close to constants: the variable ϕ l is near zero, the variable ω 0 is near ω l , the variable a l varies slightly below the level c + γ, and the variables a i , i ≠ l, vary slightly above the level c. These assumptions together with (8,10) lead to the following approximate equations for the phases ϕ i : For a given l system (23) is a system on the torus −  l n 1 . Note that the torus −  l n 1 is not an invariant set of system (12)- (14), but it is located very close to the true (n − 1)-dimensional invariant manifold  − l n 1 (which has a more intricate structure). The manifold −  l n 1 is stable in the transversal n + 2 dimensions and local dynamics inside this manifold are equivalent to the dynamics of system (23) on l n 1 its geometry is nontrivial in the whole phase space × × n n T R R even in the case of identical oscillators. It becomes more and more complex with the complication of the natural frequencies distribution. This is why a strict mathematical proof of the existence of such manifold is not an easy task. We presume that such proof could be obtained using the normal hyperbolicity theory. Here we only note that the assumption about the existence of such manifolds  l n 1 − with dynamics (23) is helpful for explaining all bifurcation transitions when the natural frequencies of oscillators are changed. The conclusions obtained by this reasoning are fully confirmed by numerical computations.
Without loss of generality let us order all natural frequencies of POs as We have shown (see Methods) that besides stable points Q l there exist other saddle fixed points S l,j , j ≠ l, with only one different coordinate ϕ i . The local saddle-node bifurcation of these two points is also a global saddle-node on invariant cycle (SNIC) bifurcation along the mentioned coordinate ϕ i (which belongs to the manifold −  l n 1 ). Figure 2(a-c) shows disappearance of the stationary WTA regime Q l and appearance of the non-stationary WTA regime LC l . Figure 3(a-c) illustrates the following proposition:

Proposition 5. SNIC bifurcations occur as a result of merging between the stable point Q l and the saddle point S l when
A stable limit cycle LC l appears along the 1-dimensional manifold of the saddle S l with the unbounded coordinate ϕ n (or ϕ 1 in the second case) and bounded phase coordinates ϕ i , i = 1, …, n − 1 (or ϕ i , i = 2, …, n, respectively). If ω 1 − ω n = |b|, two SNIC bifurcations occur simultaneously with the points Q 1 and Q n giving birth to two stable limit cycles LC 1 and LC n .
The periodic trajectory LC 1 rotates in the positive direction along ϕ n with the average frequency and the rotation goes in the negative direction along ϕ n . The amplitudes of POs whose phases run on LC 1 are also bounded (according to Proposition 3), a 1 ≈ c + γ, a n oscillates in a small neighborhood above the value c and all other amplitudes are approximately equal to a i = c. Another SNIC bifurcation at the point Q l occurs when the distance ω l − ω n (or the distance ω 1 − ω l ) is equal to |b| while the distance ω 1 − ω l (respectively, ω l − ω n ) is smaller than |b|. SNIC bifurcations happen with all the points Q s , s < l, (or with Q s , s > l, for ω 1 − ω l = |b|) before it happens with Q l because ω s − ω n ≥ ω l − ω n = |b|.
Besides limit cycles, the non-stationary WTA regime can be associated with more complex solutions if the number of POs n ≥ 3 and there are several (more then one) POs for which the variables ϕ i (i ≠ l) are not bounded but constantly run in the positive or negative direction in phase space. The dynamics on the high-dimensional toroidal manifold l n 1  − , n ≥ 3, is adequately described by system (23). More complex WTA regimes (with some number of PO phases running relative to the phase of the CO) appear as a generalization of the SNIC bifurcation. We call this bifurcation the Saddle-Node (fold) bifurcation on Invariant Torus (SNIT). Similarly to the SNIC bifurcation of two points on the 1D torus (circle) 1  , the SNIT bifurcation is a saddle-node (fold) bifurcation of the stable and saddle 1D cycles on the 2D torus 2  (which leads to full transitivity of the torus), or the bifurcation of stable and saddle 2D tori in 3D phase space  3 , or (in the general case) the bifurcation of the stable and saddle (m − 1)-dimensional tori on m-dimensional torus  . m These bifurcations are described by the following preposition (see Methods for more details): The dynamics on the invariant torus LT l m can be very complex starting from the dimension m = 2. A lot of useful information on this subject can be found in the classical works 68,69,74 (especially for the case  2 ). SNIT transitions of LC l and SL l on 2  are schematically presented in Figs 2(c-e) and 3(c-e) (two simultaneous SNIT bifurcations on two parallel 2D tori), The SNIT bifurcation of 2D tori LT l 2 and ST l 2 on  3 is presented in Fig. 3(e,f). In the two dimensional case the SNIT bifurcation occurs in the literature under alternative names, like a boundary of partial mode-locking (if one takes a Poincaré map, then it is the boundary of an Arnold tongue) or a saddle-node periodic orbit with the global reinjection. The principal parts of the structure of a similar bifurcation but without reinjection for m = 3 were described in detail in 71,73 . As far as we know, there is no mentioning of SNIT bifurcations in the literature for the cases m ≥ 4 and for m = 3 with the global reinjection. It is presented here for the firs time.
A diagram of POs distribution on the plane (a i ,ϕ i ) is schematically shown in Fig. 4. All amplitudes are located inside the area a i ∈ [c, c + γ], the winning PO moves somewhere near the point (c + γ, 0) (the area is shown in blue color). Loser POs are distributed in the red area. Some of them form a cluster near the point (c,π), while others run around in the red ring area close to the circle a i = c.

Remark 2.
The results obtained in this subsection can be generalized for the functions f, g, h more general or slightly different from those defined in (8)- (10). For example, similar results can be obtained for an arbitrary periodic function g(x) (not necessary odd) that satisfies the second and third conditions in (6) In this case bifurcation conditions depend on the values of x max and x min of the function g(x) in addition to the values of ω i and b. The results similar to those above can also be obtained without assuming oddness of the function f(x). The only critical assumptions are high values for the slope f′(0) and zero slope f′(π). Evenness of the function h(x) is also not necessary, but other conditions (7) must be fulfilled.
No-winner solutions. According to Proposition 2, system (12)- (14) with identical POs has not only WTA fixed points P 1 but also the fixed point P 0 = (ϕ 1 , …, ϕ n , ω 0 , a 1 , …, a n ) = (π, …, π, ω, c, … c) that corresponds to a regime, where all POs are in antiphase to the CO (that is there are no winners). In contrast to the points P 1 (which can be asymptotically stable when b < 0) the point P 0 is stable in 2n directions and it is neutral along the last direction (see Lemma 1 in Methods). The neutral direction of P 0 is along the straight line that can be described as one parametric set of initial conditions for the natural frequency of the CO ~L a a ac i n ( ) { ( , , , , , , ) : According to Proposition 1, the fixed point P 0 is isolated, there are no other fixed points in a small neighborhood of this point. However, considering specific functions (8)-(10) for the system we have: dϕ i /dt ≈ 0, dω 0 /dt ≈ 0, dh/dt = 0, i = 1, … n, when the phase point belongs to  ω L( ) 0 in a small neighborhood of P 0 . As a result, computer simulations of the system "perceive" the points of ω  L( ) 0 in the neighborhood of P 0 as fixed points that are neutral along L and stable in other directions.
Consider the system with functions (8)-(10) and perturbed frequencies ω i = ω + Δ i , i = 1, …, n. The perturbed point P 0 (ω 1 , …, ω n ) moves smoothly in phase space from the location P 0 (ω, …, ω) = P 0 with the perturbation of natural frequencies of POs. The perturbation of the frequencies does not change amplitude coordinates a i = c for a wide region of natural frequencies. Our computer simulations show that the perturbed point P 0 (ω 1 , …, ω n ) remains stable in 2n directions and neutral in a single direction. The neutral line L of the point P 0 (ω 1 , …, ω n ) moves with this point and rotates in phase space depending on the frequency distribution. The coordinate ω 0 of the perturbed point is approximately equal to (ω max + ω min )/2, i.e. only the fastest and slowest POs determine the average frequency of the CO in the no-winner situation. There are a few scenarios of disappearance of the point P 0 (ω 1 , …, ω n ) which are based on a saddle-node bifurcation of this point with another saddle point (or points) depending on the system dimension and frequency distribution. In any case, the region where the mentioned fixed point exists belongs to the region which is bounded by the condition ω max − ω min < 2|b|. According to (26), the point P 0 has n eigenvalues A perturbation of ω i leads to the (nonuniform) decrease of the corresponding eigenvalues |λ i (P 0 (ω 1 , …, ω n ))|, i = 2, …, n + 1, which also indicates a restriction put on the attraction basin of P 0 (ω 1 , …, ω n ).

Examples of dynamics.
In this subsection we give some examples of different types of dynamics of system (12)-(14), pointing out the bifurcations that split these examples.
The following examples will be considered: 1) Identical natural frequencies of POs: ω i = ω, i = 1, …, n. According to Proposition 1, system (12)-(14) has n S N -symmetrical stable fixed points P 1,i with coordinates (17). 2) Non-identical natural frequencies of POs that satisfy the inequalities In this case system (12)- (14) has n non-symmetrical stable fixed points Q l , l = 1, …, n, with coordinates 3) The natural frequencies of POs satisfy the equality In this case two simultaneous saddle-node bifurcations occur with two stable points Q i max , Q i min with corresponding coordinates ω max = max i ω i , ω min = min i ω i and saddles S , i max S i min . In this case the saddlenodes also have the coordinates /2 i max ϕ π = , ϕ π = /2 i min , respectively. It is possible that more than two simultaneous saddle-node bifurcations take place if more than one natural frequency reaches the same maximal (minimal) value max i ω i (min i ω i ) (this bifurcation transition is schematically represented in Figs 2(a-c) and 3(a-c)). 4) At least one pair of the natural frequencies satisfies the inequality This implies that the point Q i max Q (or ) i min disappears after the saddle-node bifurcation and a new stable limit cycle appears. The local saddle-node bifurcation is a global SNIC bifurcation by which a limit cycle appears from the 1-dimensional unstable invariant manifold of the saddle S i max S (or ) i min . The appearance of the limit cycle is possible because the mentioned 1-dimensional invariant manifold of S i max reaches the appropriate stable point Q i max along the coordinate i max ϕ that is closed on the torus. Stable fixed points Q 2 , …, Q n−1 and two stable limit cycles LC 1 , LC n co-exist after the first SNIC bifurcation. The system can have up to n − 1 limit cycles LC i . 5) The first SNIT bifurcation and the appearance of the first 2-dimensional limit torus LT 1 2 (or LT n 2 ) occur when |(ω 1 − ω n−1 )/b| = 1 (or when |(ω 2 − ω n )/b| = 1). In this case a stable limit torus LT 1 2 and stable limit cycle LC n−1 appear simultaneously and exist for |(ω 1 − ω n−1 )/b| > 1 (or the limit torus LT n 2 and stable limit cycle LC 2 appear simultaneously).  Fig. 6. The figure shows that in the WTA regime the CO can phase lock 1, 2, or 3 POs. In the first case the phase differences of non-synchronized POs can move in the same direction (positive or negative) or they can move in opposite directions. In the second case the phase difference of one non-synchronized PO is nearly stable while the phase difference of the other non-synchronized PO moves in positive or negative direction. In the third case all phase differences are stable. In any case, in this example the amplitude of the winning PO tends to 12 and the amplitudes of loser POs are near 2.
Example 2. Figure 7 shows different types of WTA dynamics that can appear in the system with 4 POs. The corresponding types of attractors for system (12)- (14) with the natural frequencies (24) are also indicated: (a) P 1,l for 1 4  Massive simulations. In this section we present different types of dynamics in terms of the frequencies with which they may appear in system (1)-(4) with functions (8)- (10). To check whether the WTA regime can be implemented by the system, we assign to a 1 an initial value that is greater than initial values of the amplitudes of other POs. It is expected that in this case the first PO will have a greater chance to win the competition for synchronization with the CO. In simulations we set 4 ≤ a 1 (0) ≤ 12, while the initial values of other amplitudes are a i (0) = 2, i ≠ 1. During a simulation, the connection strengths are adapted according to equation (4). The values of the parameters in (4) are c = 2, γ = 10. This means that in the stationary case the amplitude of the winner will tend to the value c + γ = 12, while the other amplitudes will tend to c = 2. In the non-stationary case the amplitude of the winner varies in a narrow range around 12 and the amplitudes of the other POs vary around 2. The difference in asymptotic dynamics of POs can be determined by checking whether the trajectory of a PO's amplitude approaches the lower or higher boundary. Before checking this, some time should pass for the   Two thresholds are introduced, H high and H low ., such that c < H low < H high < c + γ. In computations we set H high = 10 and H low = 3. If during the time interval (T 1 , T 2 ) the trajectory of a i (t) lies above H high , then the i-th PO is identified as the winner.
All parameter values for different simulation examples are summarized in Table 1. Stochasticity is included in the simulations through initial values. First, the natural frequencies of POs are randomly distributed in the interval (4.9, 5.1) in Examples 3, 4, 6, 7 and in the interval (4.25, 5.75) in Example 4. Second, initial phases of POs are randomly distributed in one of the ranges (−π/8, π/8), (−π/4, π/4), (−π/2, π/2), (−π, π) in Example 6. In other examples all initial phases of the POs are 0. The initial value of the natural frequency of the CO is equal to 5 in all examples, which is the mean of the distribution of the natural frequencies of POs. The initial value of the phase of the CO is always 0.
Massive computational experiments show that the following three types of dynamics can be observed: A. a 1 (t) > H high for T 1 < t < T 2 and the first PO is the only oscillator for which this inequality is valid. This means that the first PO is the winner of the competition for the synchronization with the CO. B. a i (t) > H high for T 1 < t < T 2 and a single value of the index i ≠ 1. This means that the first PO lost the competition despite its initial advantage in influencing on the CO. C. a i (t) < H low for T 1 < t < T 2 and all i = 1, …, n. This is the case when there is no winner of the competition.
Note that due to Proposition 2 there cannot be more than one winner of the competition in the stationary case. Simulations show that the same situation is valid for a non-stationary case too. Another important fact is that no chaotic behavior has been observed in the system for the parameter values presented in Table 1.
It is reasonable to think that a WTA system demonstrates better performance if in most cases its dynamics correspond to situation (A), while the other cases occur only rarely (if at all). In particular, the probability of the (C)-type dynamics should be minimized, since in this case the result of the competition is indefinite.
In each of the following examples we ran 1000 simulations and determined the number of cases when one of the types of dynamics (A)-(C) took place. The results for different parameter values are presented in the form of histograms in Figs 8,9,10,11 and 12. Example 3. In this example we vary the initial value of the amplitude of the first PO a 1 (0) consecutively assigning to it the values 4, 6, 8, 10, 12. It is expected that the higher the value of a 1 (0), the greater the chance that the first PO will win the competition. This expectation is confirmed by Fig. 8. The probability that the first PO wins the competition increases from 34.1% for a 1 (0) = 4 to 91.8% for a 1 (0) = 12. Respectively, the probability that another Parameter γ of equation (4) 10 10 10 10 10 Parameter β of equation (4)  PO wins the competition decreases from 60.3% to 7.7%. The probability that case (C) will be the outcome of the competition is always low and decreases when a 1 (0) increases. In the following examples we set a 1 (0) = 8 (to give the first PO a good chance to win the competition) and study how the variation of other parameters influences this probability.

Example 4.
In this example we vary the parameter ν of the function f. This parameter controls the steepness of the function f at 0. The results presented in Fig. 9 show that greater values of ν slightly increase the probability that the first PO will be the winner. More importantly, greater values of ν significantly decrease the probability that case (C) occurs.

Example 5.
This example is similar to Example 4 with the only difference that the range of the distribution of the natural frequencies of POs is 1.5, which prevents the possibility of stationary dynamics if the first PO is the winner. However, the system preserves the capability of WTA behavior with approximately the same efficiency as in Example 4. This is illustrated by Fig. 10. Example 6. This example shows that the width of the range of the distribution of initial phases of POs plays a crucial role for determining the WTA regime. Figure 11 shows that the best results for WTA can be obtained if this range is small. The probability that the first PO wins the competition decreases rapidly if the width of the distribution increases. Also the probability of case (C) rapidly increases if the range of this distribution. exceeds 0.5π.

Example 7.
This example shows that the speed of adaptation of the amplitudes of POs β is important for successful implementation of the WTA regime. The best results are obtained when β is small (Fig. 12). This represents the fact that for proper functioning of the system the process of phase synchronization should be much faster than that of amplitude adaptation.

Discussion
The objective of this paper was to suggest a WTA system whose functioning is based on oscillatory principles. As we have shown, the WTA can be realized in a system built from generalized phase oscillators with a central oscillator (CO). The WTA appears as a result of the competition between peripheral oscillators (POs) for synchronization with the CO. The amplification or inhibition of the activity of oscillators is a function of the phase relations between the CO and POs. In this respect our implementation of WTA differs from the traditional approach, where excitation or inhibition is directly conditioned by excitatory and inhibitory connections between the elements of the system. Our approach is consistent with the putative role of oscillatory processes in brain functions such as decision making and selective attention. It also opens the possibility of new hardware implementations of WTA.
Our results show that there are two different types of dynamics that can be associated with the WTA regime for the system in phase differences: 1. Stationary states, when there are constant relations between the phases of the CO and POs. In this case one of the POs (the winner) works inphase with the CO, while the difference between the phases of the CO and other POs is greater than π/2. A special case is represented by the system with identical natural frequencies of oscillators. In this case the "losers" work in anti-phase relative to the CO. 2. Non-stationary states. In this case one of the POs (the winner) works nearly coherently with the CO, while the set of POs is divided into two subsets. One subset contains the POs whose phases indefinitely move in the positive or negative direction relative to the phase of the CO. In the other subset the difference between the phase of a PO and the CO oscillates in a small range around a value greater than π/2.
According to Propositions 2 and 4, in the stationary case there cannot be more than one winning PO. This has been proved under two types of conditions. If the POs are identical, the statement is true for rather arbitrary functions of system (1)-(4). If POs have different natural frequencies, the statement is true at least for the functions specified by formulas (8)- (10).
Analytical analysis and numerical simulations show that the latter statement is also valid in the non-stationary case. The conditions for the existence of the stationary WTA regime are described by Lemmas 1-4 in the Methods section, Remark 3, and inequalities (35). The bifurcations that lead to the transformation of the WTA behavior from the stationary to non-stationary form are described by Propositions 5 and 6. The propositions describe a new type of bifurcation (saddle-node on invariant torus (SNIT) bifurcation) which has not been previously described in the literature.
Computer simulations under conditions (8)-(10) demonstrate that our WTA implementation is sensitive to the parameters of the system. The following conditions can increase the probability that the PO with the highest initial amplitude is the winner: 1. The highest initial amplitude is significantly greater than the initial amplitudes of the other POs. The second condition in this list is in line with the neurobiological role of phase reset in information processing in the brain 75 . The probability of an indefinite result of the competition can be reduced by the proper choice of parameter values.

Methods
Proofs of the main results. In this section we present the ideas of the proofs of the main results.

Stability of equilibria (Proof of Proposition 2).
For the functions f(x), g(x), h(x) satisfying conditions (5)-(7) (2n + 1)-dimensional system (16) has the Jacobian matrix bg a n f a n f n f n f a n f b g a n f n f n f a n f a n f where ϕ ϕ A a n f a n f a n f a n f where I n = diag{1, …, 1} is the n × n identity matrix, E n,m is the n × m-matrix of 1, 0 n,m is the n × m-matrix of 0. Note that the stability of a fixed point does not depend on the common frequency ω. Denote Then for the fixed point P k we have A a n p a n p a n p a n p J bq a n p a n p a n p b q a n p a n p a n p b q a n p a n p b q a n p a n p Other eigenvalues can be found as the eigenvalues of the (n + 1) × (n + 1)-block Using these eigenvalues we can formulate stability conditions for fixed points of different types. The fixed point P 0 corresponds to the antiphase regime of all POs relative to the CO (no-winner state). Linear neutrality of this point along one direction corresponds to the bifurcation line in the two-parametric bifurcation plane. The eigenvalues corresponding to P 0 can also describe the dynamics when a continuous set of periodic orbits exists close to this point. by formula (21). Each coordinate ϕ i , i ≠ l, in this formula is one of the two solutions of the equation bsinφ i = ω i − ω l , namely the one that is closer to π. The necessary condition for the existence of such fixed points is max i≠l |(ω l − ω i )/b| ≤ 1. One can check (in a similar way that has been used for the stable fixed point Q l ) that there exist n − 1 other fixed points S l,j , j ≠ l, which have the same coordinates as Q l except one coordinate . The point S l,j has the same Jacobian (34) as Q l but with  ϕ j instead of ϕ j . In the conditions of stability each fixed point S l,j has one positive eigenvalue λ ϕ = −b cos( ) j j  . This means that S l,j is a saddle point that belongs to a 1-dimensional unstable invariant manifold W u (S l,j ). The eigenvector corresponding to the eigenvalue λ j belongs to the coordinate line of the variable ϕ j . This means that both branches of W u (S l,j ) have the same coordinate as S l,j (and Q l ) except ϕ j and these branches spread along the variable ϕ j (at least locally) (Fig. 2(a)). The divergence in the transversal direction to W u (S l,j ) along ϕ j is negative (the Jacobian shows that W u (S l,j ) is stable in the transversal direction), which implies that both branches of W u (S l,j ) reach the stable point Q l from two opposite sides along ϕ j , forming a locked circle. The local saddle-node bifurcation of Q l and S l,j ( Fig. 2(b)) is simultaneously a global (SN) bifurcation on the invariant circle, (SNIC bifurcation), and it leads to the appearance of a stable limit cycle LC l (Fig. 2(c)). This bifurcation can happen if only |(ω l − ω j )/b| = 1 at the point ϕ j = ±π/2. Due to (24) SNIC bifurcations and disappearance of Q l can only happen if j = 1 or j = n. It is obvious that the first (in the sense of the continuous movement of ω i away from ω l ) such bifurcation occurs when ω 1 − ω n = |b| and l = 1 or l = n. This bifurcation simultaneously happens with both of the two POs, the first and the last in numeration (24).
Thus, we have two simultaneous SNIC bifurcations with the points Q 1 , S 1,n and Q n , S n,1 , where Q 1 is the winning PO with the largest frequency ω 1 and = S S : n 1 1 , is the corresponding saddle with the unstable manifold W u (S 1 ) along the phase variable ϕ n , (the same for Q n and W u (S n ), where = S S : n n ,1 ). As a result of the bifurcations, we have stable limit cycles LC 1 , LC n , where each of the cycles has one unbounded phase variable ϕ n and ϕ 1 , respectively. Summing up the above considerations we can formulate Proposition 5.
The cycle LC 1 is non-homologous to the zero periodic orbit on the torus  n in phase space, the variable ϕ n of the cycle is unbounded while all other phase variables are bounded (the same is true for LC n and ϕ 1 ). SNIC bifurcations with the fixed points Q 1 , Q n are not only possible ones. Other SNIC bifurcations at the point Q l occur when the distance ω l − ω n (or the distance ω 1 − ω l ) is equal to |b|, while the distance ω 1 − ω l (respectively, ω l − ω n ) is smaller than |b|. Such bifurcations lead to the appearance of new cycles LC l , l = 2, …, n − 1, and they usually happen simultaneously with more complex bifurcation of the existing cycles LC 1 , LC n which we will describe later.
Besides limit cycles with one running variable ϕ i (i ≠ l), the non-stationary WTA regime can be associated with more complex solutions when n ≥ 3 and there are several POs (more then one) whose variables ϕ i are boundlessly running in the positive or negative direction. Let us describe how such regimes may appear.
Consider a 2-dimensional hyperplane with free variables ϕ j , ϕ m and all other variables being the constant coordinates of the fixed point Q l (Fig. 2(a)). In this hyperplane the point Q l has the coordinates ( , ) j m ϕ ϕ and the saddle S l,j has the coordinates π ϕ ϕ − ( , ) j m . The SNIC bifurcation of Q l and S l,j leads to the emergence of a saddle-node point SN l,j (Fig. 2(b)) and then to a stable limit cycle LC l (Fig. 2(c)). According to the m-th equation (23), the system has also two additional fixed points S l,m belonging to the 1-dimensional manifold W u (S l,m ) and the saddle S l,j,m belonging to the 2-dimensional unstable manifold W uu (S l,j,m ) (S l,j,m is an unstable node on the 2-dimensional plane ϕ ϕ ∈  ( , ) j m l 2 ). The points S l,m and S l,j,m have a common 1-dimensional invariant manifold which is stable for the first point and unstable for the second one. In other worlds, these points lie on the circle formed by two branches of the manifold. The SNIC bifurcation of S l,m and S l,j,m leads to the appearance of a saddle-node point SN l,j,m ( Fig. 2(b)) and then to a saddle limit cycle SC l (Fig. 2(c)) with one unstable invariant manifold, spreading along the variable ϕ m . Two SNIC bifurcations and the appearance of LC l and SC l happen almost simultaneously ( Fig. 2(a-c)). In this case we have stable and unstable limit cycles on the torus . l 2  In the case of the non-simultaneous SNIC bifurcations we have two variants of Cherry flows 68,69,74,76 on the torus  l 2 : 1) a stable cycle LC l and two fixed points S l,m , S l,j,m or 2) a saddle (unstable on  l 2 ) limit cycle SC l and two fixed points Q l , S l,m . Note that the manifold l 2  of system (12)- (14) (that is represented here by the torus l 2  ) is a 2-dimensional unstable invariant manifold of the saddle cycle SC l which is stable in the transversal 2n − 1 directions in global (2n + 1)-dimensional phase space × × T R R n n (when SC l exists). The saddle-node bifurcation of LC l and SC l on the invariant torus  l 2 (SNIT) is similar to the SNIC bifurcations considered above (one can also call this bifurcation a saddle-node of periodic orbits with global reinjection). Locally, it is a saddle-node (fold) bifurcation of the stable cycle LC l and the saddle cycle SC l (which happens along the variable ϕ m ) that leads to the transitivity of trajectories through the torus. As a result of the SNIT bifurcation, we obtain first saddle-node limit cycle SNC l (Fig. 2(d)) and then a stable (in transversal direction) limit torus that is denoted as LT l 2 (Fig. 2(e)). The described SNIT bifurcation is a global bifurcation on the torus l 2  (as SNIC is a global bifurcation on a circle) and it leads to transitivity of the trajectories on the whole 2D torus. We note that the transitivity on the whole two-torus is true only for the fat Cantor set of parameter values with irrational winding ratio. The bifurcation occurs when the phase difference |ω l − ω m | reaches the value |b| but |ω l − ω j | ≥ |b| and |ω l − ω i | < |b| for i ≠ l, i ≠ j, i ≠ m. We already mentioned that if the range of the natural frequencies of POs increases the first two SNIC bifurcations occur simultaneously when the 1-st or the n-th PO becomes the winner. This leads to the appearance of the limit cycles LC 1 and LC n . One can check that enlarging further the frequency distribution we obtain a SNIT bifurcation simultaneously with a SNIC bifurcation. As a result of such simultaneous bifurcations, the pair of LT 1 2 and LC n−1 or the pair of LT n 2 and LC 2 appears. The emergence of other limit tori is associated with more complex bifurcations of a similar type.
In the same way it is possible to show that there exists a saddle torus ST l 2 that appears as a result of a SNIT bifurcation of two saddle limit cycles with 1-dimensional and 2-dimensional unstable manifolds. Two simultaneous bifurcations of the appearance of stable and saddle 2D invariant tori LT l 2 and ST l 2 are shown in Fig. 3(c-e).