A versatile class of prototype dynamical systems for complex bifurcation cascades of limit cycles

A general class of prototype dynamical systems is introduced, which allows to study the generation of complex bifurcation cascades of limit cycles, including bifurcations breaking spontaneously a symmetry of the system, period doubling and homoclinic bifurcations, and transitions to chaos induced by sequences of limit cycle bifurcations. The prototype systems are adaptive, with friction forces f(V(x)) being functionally dependent exclusively on the mechanical potential V(x), characterized in turn by a finite number of local minima. We discuss several low-dimensional systems, with friction forces f(V) which are linear, quadratic or cubic polynomials in the potential V. We point out that the zeros of f(V) regulate both the relative importance of energy uptake and dissipation respectively, serving at the same time as bifurcation parameters, hence allowing for an intuitive interpretation of the overall dynamical behavior. Starting from simple Hopf- and homoclinic bifurcations, complex sequences of limit cycle bifurcations are observed when the energy uptake gains progressively in importance.

The term 'prototype dynamical system' is employed for generic, but otherwise reduced systems, allowing to study and to understand a certain relevant phenomenon (like dynamical behavior and/or bifurcation scenario). For this, the dynamical behavior of the system should be dominated by the prime phenomenon of interest, with the system being otherwise simple enough to allow for straightforward numerical and (at least partial) analytic investigations [1][2][3][4] . Additionally, their dynamical behavior can often be understood in terms of general concepts, such as energy balance, symmetry breaking, etc.
Examples of prototype systems are the normal forms of standard bifurcation analysis 5,6 and classical systems, like the Van der Pol oscillator 5 , or the Lorenz model 7 , which have been of central importance for the development of dynamical systems (systems) theory. As an example we consider the Liénard equation, a generic adaptive mechanical system, which includes the Van der Pol oscillator and the Takens-Bogdanov system 8,9 . The periodically forced extended Liénard systems with a double-well potential have also been studied by many authors (see e.g. the double-well Duffing oscillator [10][11][12]. In this paper we propose a new class of autonomous Liénard-type systems, which allow to study cascades of limit cycle bifurcations, using a bifurcation parameter controlling directly the balance between energy dissipation and uptake, and hence the underlying physical driving mechanism. Though there are a range of alternative construction methods for dynamical systems in the literature (see e.g. [13][14][15] ), they generally involve abstract concepts, such as implicitly defined manifolds, or mathematical tools accessible only to researchers with an in-depth math training. In contrast to these methods, we provide here a mechanistic design procedure, based on the construction of attractors through the interaction of generalized friction forces with potential forces, an intuitive concept especially suitable for interdisciplinary investigations (e.g. in modeling cardiovascular systems 16 or for solving optimization problems 17 ), making it easily accessible and implementable for other scientific communities (such as neuroscience, biology etc.) as well.
As an introductory example for the role of balance between energy uptake and dissipation, in both local and global bifurcations, we reconsider the Bogdanov-Takens system, which is often used as a prototype system for homoclinic bifurcations 5 . Here, the mechanical potential is a third order polynomial, as illustrated in Fig. 1. The friction force is directly proportional to the velocity y, hence fixpoints of (2) correspond to the minima and the maxima of the potential V(x). The dynamics of the Bogdanov-Takens system is controlled by the parameter μ, defining, in terms of the mechanical energy E, the regions of dissipation and energy uptake in the potential valley, compare Fig. 1. The region x > μ of energy uptake increases when the bifurcation parameter μ is decreased, leading to two consecutive transitions. Initially the potential minimum becomes repelling, undergoing a supercritical Hopf bifurcation and a stable limit cycle emerges. Decreasing μ further, the extension of the limit cycle increases, merging at μ c with the stable and unstable manifolds of the saddle, resulting in a homoclinic bifurcation.

Results
The key mechanism leading to the bifurcations in the Bogdanov-Takens systems is the availability of a parameter, allowing to change the balance between energy uptake and energy dissipation along limit cycles. Our aim is to generalize this idea to the case of mechanical systems characterized by an arbitrary number of potential minima. For this purpose we consider x y x x y y 4 2 which describes a 2d-dimensional system, with d spatial coordinates, and friction forces f(V(x)) depending functionally only on the mechanical potential V(x), allowing for a fine-tuned control of the energy dissipation and uptake around the respective potential minima. A well known example of a system of type (4) is the Van der Pol oscillator: for which the regions of energy uptake and dissipation remain fixed, with ε regulating the overall influence of the velocity-dependent force.
The simplest generic class of friction functions f(V) entering (4) are polynomial: where α regulates the overall strength of the friction, and the individual μ 1 < μ 2 < μ 3 are the respective zeros, the points at which dissipation changes to anti-dissipation and vice versa, compare Fig. 2. When using f 1 (V) and the mechanical potential V(x) = x 3 /3 − x 2 /2, the resulting flow in phase space is equivalent to the one of the Bogdanov-Takens system (2), as shown in Fig. 1.

Generalized mechanical potentials.
We are interested in using (4) as prototype dynamical systems, especially for the case of non-trivial mechanical potentials V(x) having an arbitrary number M of local minima. One could in principle consider higher-order polynomials for this purpose, however these do not allow to control the overall height of the potential and the relative width of the local minima in as simple a fashion.
For this purpose we use throughout this study potential functions of the kind where the z n > 0 determine the half-width of the respective local minima, and the p n satisfy the self-consistent condition: For deep minima, with (z i + z j ) ≪ |x i − x j |, the positions and the heights of the local minima are close to x n and V n respectively. We found that a relative accuracy of 10 −2 for V n can already be achieved in general after three or four iterations.
Limit cycle bifurcation cascades. The system of type (4) allows to describe complex cascades of limit cycle bifurcations. In Fig. 2 we show some illustrative examples, using a symmetric double-well potential and linear/quadratic/cubic friction functions respectively, see Eq. (6). We used numerical methods (see the Methods section) to obtain the respective full bifurcation diagrams, with solid/dashed lines denoting stable/unstable fixpoints and limit cycles. The corresponding flow in phase space is illustrated in Fig. 3. For negative μ 1 the two fixpoints (± 1, 0) are stable, for the case of f 1 (V) and f 3 (V), and stable limit cycles evolve via two supercritical Hopf bifurcations (bifurcations). For f 2 (V), on the other hand, a subcritical Hopf bifurcation is observed at μ 1 = 0. The respective stable/unstable limit cycles merge for f 1 (V) and f 2 (V) in a homoclinic bifurcation, whereas a more complex bifurcation diagram emerges for f 3 (V). Saddle node bifurcations of limit cycles are present for both f 2 (V) and f 3 (V).

Chaos via period doubling of limit cycles.
We consider now a prototype system (4) with a two-dimensional symmetric potential function V(x), as defined in (7), having two minima x 1,2 = ± (1, − 1), and a linear friction term f 1 (V) = 0.5(μ − V). Both diagonals in the (x 1 , x 2 ) plane are symmetry axes of the system, as discussed in the Methods section. In Fig. 4 we present examples of stable limit cycles and of a chaotic trajectory, as projected to the (x 1 , x 2 ) plane. In Fig. 5 the corresponding bifurcation diagram is presented. The diagram shows Hopf bifurcations (H), homoclinic bifurcations (HO), branching of limit cycles via spontaneous symmetry breaking (SSB), period doubling of limit cycles (PD), and a transition to chaotic behavior:  the two potential minima become unstable, just as for the one-dimensional spatial system presented in Fig. 2, resulting in two equivalent supercritical Hopf bifurcations. We note that, as a result of the symmetric potential function (9), a second branch of limit cycles is created by the two Hopf bifurcations (see the discussion in the Methods section and the Supplementary Information). However, since in the parameter region of interest these limit cycles are mostly unstable, we have not investigated them in detail.
( ) 0 143 HO 1 the limit cycles merge, as in Fig. 4(a,b), in a homoclinic transition. The limit cycle stays, however, exactly on the diagonal x 1 + x 2 = 0.
SSB At the first branching point of limit cycles, μ ≈ .
( ) 0 171 SSB 1 the symmetry with respect to the diagonal (1, − 1) is spontaneously broken, as in Fig. 4(b,c), with the two limit cycles still being symmetric with respect to the (1, 1) diagonal. The latter symmetry is broken at the second branching point μ ≈ .
, as in Fig. 4(c,d), creating four symmetry related stable limit cycles. PD For larger values of the bifurcation parameter μ 1 a series of period-doubling of limit cycles is observed, with the first occurring at μ ≈ .
, as in Fig. 4(d,e). The next period-doubling transition occurs at μ ≈ . we observe seemingly chaotic trajectories, as illustrated in Fig. 4(f). Studying the transition to chaos is not the subject of the present investigation and we leave it to future work. We presume however, that the transition occurs via an accumulation of an infinite number of of period-doubling transitions of limit cycles, similar to the ones observed for the Lorenz system 18 and for the Rössler attractor 19,20 .
Our prototype system (4) is not generically dissipative. We have evaluated the average contraction rate σ, as defined by (22) in the Methods section, and presented the results in Fig. 5. Phase space contracts trivially along the attracting limit cycles, but also, on average, in the chaotic region, where the average Lyapunov exponent λ becomes positive. λ is negative for μ 1 < 0, when only stable fixpoints are present, vanishing for intermediate values of μ 1 , when stable limit cycles are present. The later is due to the fact, see Fig. 7 and the corresponding Methods section, that two initially close trajectories will generally flow to the same limit cycle with the relative distance becoming constant. V(x)). The bifurcation parameter μ 1 is 0.1, 0.15, 0.25, 0.265, 0.2698, 0.3 from (a) to (f). In (d) the four limit cycles can be mapped onto each other by using the symmetry operations σ 1,2 or σ 3,4 , as discussed in the Methods section. For (e) only a single of the four stable limit cycles is shown. This needs to circle the two potential minima twice in order to retrace itself. In (f) an example of a chaotic trajectory is given.

Figure 4. Stable limit cycles and chaotic orbits of the Liénard prototype system (4) in a two-dimensional symmetric double well potential V(x) (color coded, as defined by Eq. (9)) and with a linear friction term f 1 (V(x)) = 0.5 (μ 1 −
For larger values of μ 1 > 0.322 the chaotic region transforms into a phase of intermittent chaos as illustrated in Fig. 6, in which an extended quasi-regular flow along the (−1, 1) diagonal is interseeded by a roughly perpendicular bursting flow. This behavior is, to a certain extend, reminiscent to a scenario of intermittent chaos 21 , in which a strange attractor is embedded in a higher-dimensional space with partly unstable directions. We have, however, not investigated the observed intermittent dynamics in detail.

Discussion
We have proposed and discussed a prototype dynamical system (4) in which the friction forces ∝ f(V) depend functionally only on the mechanical potential V(x). We have shown that complex cascades of limit cycle bifurcations can be obtained even for two dimensional phase spaces, when the friction function f(V) alternates between regions of energy uptake and dissipation.
We have also introduced a generic class of potential functions (7), which allows to define, in a relative straightforward manner, mechanical potentials with an arbitrary number of local minima and varying depth. Considering a simple double-well prototype system with two spatial dimensions (and with a four-dimensional phase space), we have shown that symmetry induced bifurcations of limit cycles and period-doubling of limit cycle transition to chaotic behavior can occur.
As discussed in the Methods section, the presence of stable and unstable fixpoints, the birth of limit cycles through transition from dissipation to energy uptake in the neighborhood of the local minima, or the symmetry properties of the prototype system do not depend on the particular method used to  Fig. 4. The second branches of limit cycles emerging from the two destabilized minima (see Supplementary Fig. S1) are not shown here. One observes Hopf and homoclinic bifurcations (H and HO), branching of limit cycles via spontaneous symmetry breaking and period doubling (SSB and PD), as well as a transition to chaos for μ 1 > 0.2705. The red/green lines indicate the maximal/minimal x 1 -values of the respective limit cycles. The blue curve is the second x 1 -maxima after period doubling. The right diagram represents a zoom-in of the transition to a chaotic region, indicated by the shaded green and red areas. Only the first two period doubling bifurcations are shown. Bottom row: The average Lyapunov exponent λ and the contraction rate σ, calculated as described in the Methods section, for the corresponding μ 1 parameter intervals. For the left figure Δ μ 1 = 0.001 parameter stepsize was used. Increasing the resolution more and more periodic windows (with = λ 0) become visible, as shown on the right plot, where Δ μ 1 is decreased by a factor of ten. construct the potential function. The only requirement it has to fulfill is the existence of a certain number of local minima. Hence, other potential functions could also be considered. For example, one could study the biquadratic version of the potential (9), used in our study of chaotic behavior with the prototype system (4). We did not study in detail the bifurcation diagram for the potential function (10), however, we have checked that one would get similar results to the ones presented in Figs 4 and 5, having the same underlying driving mechanism in terms of a linear friction function f(V) ∝ (μ 1 − V). For increasing values of the μ 1 control parameter, first stable fixpoints, then spatially separated-, merging-and symmetry breaking limit cycles can also be observed. Furthermore, using the potential function (10) chaotic behavior has also been found.
As a future perspective, we note that by changing the depth of the minima, one could control the order in which the fixpoints are going to be destabilized, which might lead to other interesting phenomena. Adding an extra (maybe slow) dynamics to the positions or depths of the minima, the metadynamics  of the attractors 22 may also be considered. In this case the p n parameters should be recalculated in each time-step using the self-consistent equations (8), which proved to be fast enough for practical purposes.
Models, for which the equations of motion are derived from higher order principles, provide promising results for the understanding of many different phenomena, such as the optimization hardness of boolean satisfiability problems 17 or the complex dynamics of biological neural networks 23,24 . Generally, these methods involve the construction of a generating functional, such as the cost function or energy functional [25][26][27][28] , with the dynamics of the system being defined by a gradient decent rule. When all equations are derived from the same generating functional, the system corresponds mathematically to a gradient system for which the asymptotic behavior is determined by stable fixpoints (nodes). They can thus not produce limit cycles or oscillatory behaviors. To by-pass this problem, additional equations of motions are usually defined, derived either from a second generating functional, to induce objective function stress 29,30 , or from other considerations. In our model, the system has an inherent inertia, which, in the presence of dissipation, leads to damped oscillations around the equilibria (minima of the potential). By creating regions of antidissipation, stable oscillatory dynamics and chaotic behavior is stabilized. Considering nonsymmetric and/or higher dimensional potential functions, we expect to find an even richer set of dynamical behaviors (see the Methods section), a scenario worth to be investigated in the future. As a possible application one could use the system for modeling the dynamics of various, complex and adaptive dynamical systems, for which the generalized potential function (energy landscape, cost functional, etc.) is approximated by (8) or is found from some other considerations. An alternative type of prototype dynamical system has been shown to be useful for understanding the coexistence of spiking and bursting neural activity observed in electrophysiological experiments 2 .
Finally, we note that the concepts of dynamical systems theory, such as attractors, slow points and bifurcations have been used recently to understand phenomena of surprisingly diverse fields. We mention here the modeling of birdsongs 31 and migraine dynamics 32 , and the control mechanisms for the movements of humanoid robots 33 . The common approach, considered in these works, is the aim to construct simple and, to a certain extend, idealized dynamical systems, which allow for an in-depth understanding of certain dynamical behaviors. We hence believe that prototype systems allowing, in an intuitive manner, for the construction of models with a predefined set of attractors, as presented here, could offer a useful tool for understanding the behavior of a range of interesting interdisciplinary problems.

Methods
The bifurcation diagrams shown in Figs 2 and 5 have been constructed by using the PyDSTool 34 software package. In this section we provide the analytic calculations for the study of fixpoint stability for 2-, 3and 2d-dimensional prototype systems respectively. We note that, these properties are valid irrespective of the particular shape considered for the potential function. This is followed by a discussion of symmetry properties, and presentation of numerical methods used to estimate the average Lyapunov exponent and the contraction rate.
Hopf bifurcations in the prototype system. 2 fixpoints of the 4-dimensional prototype systems (4) correspond to critical points of the V(x 1 , x 2 ) potential function. Classification of the local minima and saddle critical points with respect to their stability can be achieved by evaluating the eigenvalues of the Jacobian of the system in terms of the Hessian of the potential function: are also solutions. As one could expect from the definition of the class of prototype systems introduced here (see Eq. (4)), the symmetry properties of the system are closely related to the particular symmetries of the potential function considered for modeling a certain behavior. Thus, finding the corresponding σ i symmetry operations, could reveal new limit cycle solutions related by symmetry.
Lyapunov exponent and contraction rate. The local Lyapunov exponent λ is determined from the growth rate of the distance Δ r(t) = Δ r 0 e λt , between point pairs with an initial displacement, which we have taken to be Δ r 0 = 10 −8 . The measurement of the Lyapunov exponent was started after a transient of t tr = 1.5 ⋅ 10 4 . Considering 100 random initial conditions the average Lyapunov exponent λ is then given by the slope of the initial linear part of the 〈 ln(Δ r)〉 curve (as given by the brown lines in Fig. 7).
The contraction rate σ, is defined as the average of local contraction rates along a set of trajectories Γ for different initial conditions: where L = ∫ Γ ds is the length of the trajectory, and f is the flow, viz the right-hand side of the evolution equations (4). σ is negative for dissipative systems, in which the phase space contracts 5,6 .