GLOBAL BIFURCATION ANALYSIS OF TOPP SYSTEM

In this paper, we study the 3-dimensional Topp model for the dynamics of diabetes. First, we reduce the model to a planar quartic system. In particular, studying global bifurcations, we prove that such a system can have at most two limit cycles. Next, we study the dynamics of the full 3-dimensional model. We show that for suitable parameter values an equilibrium bifurcates through a Hopf-saddle-node bifurcation. Numerical analysis suggests that near this point Shilnikov homoclinic orbits ex-ist. In addition, chaotic attractors arise through period doubling cascades of limit cycles.


Introduction
In this paper, we carry out a global qualitative analysis, first of all, of a reduced planar quartic Topp system which models the dynamics of diabetes [Goel, 2015;Topp et al., 2000].
Diabetes mellitus is a disease of the glucose regulatory system characterized by fasting or postprandial hyperglycemia. There are two major classifications of diabetes based on the etiology of the hyperglycemia. Type 1 diabetes (also referred to as juvenile onset or insulindependent diabetes) is due to an autoimmune attack on the insulin secreting β cells. Type 2 diabetes (also referred to as adult onset or non-insulin-dependent diabetes) is associated with a deficit in the mass of β cells, reduced insulin secretion, and resistance to the action of insulin; see [Topp et al., 2000] and the references therein.
Blood glucose levels are regulated by two negative feedback loops. In the short term, hyperglycemia stimulates a rapid increase in insulin release from the pancreatic β cells. The associated increase in blood insulin levels causes increased glucose uptake and decreased glucose production leading to a reduction in blood glucose. On the long term, high glucose levels lead to increase in the number of β-cells. An increased β-cell mass represents an increased capacity for insulin secretion which, in turn, leads to a decrease in blood glucose. Type 2 diabetes has been associated with defects in components of both the short-term and long-term negative feedback loops [Topp et al., 2000].
Mathematical modeling in diabetes research has focused predominately on the dynamics of a single variable, usually blood glucose or insulin level, on a timescale measured in minutes [Topp et al., 2000]. Generally, these models are used as tools for measuring either rates (such as glucose production and uptake rates or insulin secretion and clearance rates) or sensitivities (such as insulin sensitivity, glucose effectiveness, or the sensitivity of insulin secretion rates to glucose). Two model-based studies have examined coupled glucose and insulin dynamics [Topp et al., 2000]. In each of these studies, multiple parameter changes, representing multiple physiological defects, were required to simulate glucose and insulin dynamics observed in humans with diabetes. In doing so, three distinct pathways were found to the diabetic state: regulated hyperglycemia, bifurcation and dynamical hyperglycemia [Topp et al., 2000].
In our study, we reduce the 3D Topp diabetes dynamics model [Goel, 2015;Topp et al., 2000] to a planar quartic dynamical system and study global bifurcations of limit cycles that could occur in this system, applying the new bifurcation methods and geometric approaches developed in [Broer and Gaiko, 2010;Gaiko, 2003;Gaiko, 2012a;Gaiko, 2012b;Gaiko, 2012c;Gaiko, 2015;Gaiko, 2016;Gaiko, 2018;Gaiko and Vuik, 2018]. In Section 2, we consider the Topp model of diabetes dynamics. In Section 3, we carry out the global qualitative analysis of the reduced Topp system. Finally, in Section 4, we perform a numerical study of the full 3-dimensional Topp model getting new qualitative phenomena for this model.

Topp Model of Diabetes Dynamics
In [Topp et al., 2000], a novel model of coupled βcell mass, insulin, and glucose dynamics was presented, which is used to investigate the normal behavior of the glucose regulatory system and pathways into diabetes. The behavior of the model is consistent with the observed behavior of the glucose regulatory system in response to changes in blood glucose levels, insulin sensitivity, and β-cell insulin secretion rates.
In the post-absorptive state, glucose is released into the blood by the liver and kidneys, removed from the interstitial fluid by all the cells of the body, and distributed into many physiological compartments, e. g., arterial blood, venous blood, cerebral spinal fluid, interstitial fluid [Topp et al., 2000].
Since we are primarily concerned with the evolution of fasting blood glucose levels over a time-scale of days to years, glucose dynamics are modeled with a singlecompartment mass balance equatioṅ G = a − (b + cI)G.
(2.1) Insulin is secreted by pancreatic β-cells, cleared by the liver, kidneys, and insulin receptors, and distributed into several compartments, e. g., portal vein, peripheral blood, and interstitial fluid. The main concern is the long-time evolution of fasting insulin levels in peripheral blood. Since the dynamics of fasting insulin levels on this time-scale are slow, we use a single-compartment equation given byİ Despite a complex distribution of pancreatic β cells throughout the pancreas, β-cell mass dynamics have been successfully quantified with a single-compartment modelβ Finally, the Topp model isĠ = a − (b + cI)G, with parameters as in [Topp et al., 2000].

Bifurcation Analysis of the Reduced System
Consider system (2.6). Its finite singularities are determined by the algebraic system which can give us at most three singular points in the first quadrant: a saddle S and two antisaddles (nonsaddles) -A 1 and A 2 -according to the second Poincaré index theorem [Bautin and Leontovich, 1990;Gaiko, 2003]. Suppose that with respect to the x-axis they have the following sequence: A 1 , S, A 2 . System (2.6) can also have one singular point (an antisaddle) or two singular points (an antisaddle and a saddle-node) in the first quadrant.
To study singular points of (2.6) at infinity, consider the corresponding differential equation Dividing the numerator and denominator of the righthand side of (3.2) by x 4 (x ̸ = 0) and denoting y/x by u (as well as dy/dx), we will get the equation 3) for all infinite singularities of (3.2) except when x = 0 (the "ends" of the y-axis); see [Bautin and Leontovich, 1990;Gaiko, 2003]. For this special case we can divide the numerator and denominator of the right-hand side of (3.2) by y 4 (y ̸ = 0) denoting x/y by v (as well as dx/dy) and consider the equation See Table 1 for the color coding. The Hopf-saddle-node bifurcation is located at the point (a, l) = (0.2, 1).
Applying the definition of a field rotation parameter [Bautin and Leontovich, 1990;Gaiko, 2003], to system (2.6), let us calculate the corresponding determinants for the parameters a, b, c, α, and β, respectively: . It follows that in the first quadrant the signs of ∆ a , ∆ b , ∆ c depend on the sign of βx 2 − α y(1 + x 2 ) and that the signs of ∆ α and ∆ β depend on the sign of a − (b + c y)x on increasing (or decreasing) the parameters a, b, c, α, and β, respectively. Therefore, to study limit cycle bifurcations of system (2.6), it makes sense together with (2.6) to consider also the auxiliary system (2.7) with field-rotation parameter γ : ∆ γ = P 2 + Q 2 ≥ 0.
Theorem 3.1. The reduced Topp system (2.6) can have at most two limit cycles.
Proof. In [Broer et el., 2007;Broer and Gaiko, 2010;Li and Xiao, 2007;Zhu et al., 2002], where a similar quartic system was studied, it was proved that the cyclicity of singular points in such a system is equal to two and that the system can have at least two limit cycles; see also [Gaiko, 2016;Gaiko and Vuik, 2018;Gonzalez-Olivares et al., 2011;Lamontagne, 2008] with similar results. Consider systems (2.6)-(2.7) supposing that the cyclicity of singular points in these systems is equal to two and that the systems can have at least two limit cycles. Let us prove now that these systems have at most two limit cycles. The proof is carried out by contradiction applying Catastrophe Theory; see [Gaiko, 2003;Perko, 2002].
We will study more general system (2.7) with three parameters: α, β, and γ (the parameters a, b, and c can be fixed, since they do not generate limit cycles). Suppose that (2.7) has three limit cycles surrounding the singular point A 1 , in the first quadrant. Then we get into some domain of the parameters α, β, and γ being restricted by definite conditions on three other parameters, a, b, and c. This domain is bounded by two fold bifurcation surfaces forming a cusp bifurcation surface of multiplicity-three limit cycles in the space of the parameters α, β, and γ.
The corresponding maximal one-parameter family of multiplicity-three limit cycles cannot be cyclic, otherwise there will be at least one point corresponding to the limit cycle of multiplicity four (or even higher) in the parameter space.
Extending the bifurcation curve of multiplicity-four limit cycles through this point and parameterizing the corresponding maximal one-parameter family of multiplicity-four limit cycles by the field rotation parameter, γ, according to the Perko monotonicity theorem [Gaiko, 2003;Perko, 2002], we will obtain two monotonic curves of multiplicity-three and one, respectively, which, by the Wintner-Perko termination principle [Gaiko, 2003;Perko, 2002], terminate either at the point A 1 or on a separatrix cycle surrounding this point. Since on our assumption the cyclicity of the singular point is equal to two, we have obtained a contradiction with the termination principle stating that the multiplicity of limit cycles cannot be higher than the multiplicity (cyclicity) of the singular point in which they terminate.
If the maximal one-parameter family of multiplicitythree limit cycles is not cyclic, using the same principle, this again contradicts the cyclicity of A 1 not admitting the multiplicity of limit cycles to be higher than two. This contradiction completes the proof in the case of one singular point in the first quadrant.
Suppose that system (2.7) with three finite singularities, A 1 , S, and A 2 , has two small limit cycles around, for example, the point A 1 (the case when limit cycles surround the point A 2 is considered in a similar way).  Then we get into some domain in the space of the parameters α, β, and γ which is bounded by a fold bifurcation surface of multiplicity-two limit cycles. The corresponding maximal one-parameter family of multiplicity-two limit cycles cannot be cyclic, otherwise there will be at least one point corresponding to the limit cycle of multiplicity three (or even higher) in the parameter space. Extending the bifurcation curve of multiplicity-three limit cycles through this point and parameterizing the corresponding maximal one-parameter family of multiplicity-three limit cycles by the field rotation parameter, γ, according to the Perko monotonicity theorem [Gaiko, 2003;Perko, 2002], we will obtain a monotonic curve which, by the Wintner-Perko termination principle [Gaiko, 2003;Perko, 2002], terminates either at the point A 1 or on some separatrix cycle surrounding this point. Since we know at least the cyclicity of the singular point which on our assumption is equal to one in this case, we have obtained a contradiction with the termination principle.
If the maximal one-parameter family of multiplicitytwo limit cycles is not cyclic, using the same principle, this again contradicts the cyclicity of A 1 not admitting the multiplicity of limit cycles higher than one. Moreover, it also follows from the termination principle that either an ordinary (small) separatrix loop or a big loop, or an eight-loop cannot have the multiplicity (cyclicity) higher than one in this case. Therefore, according to the same principle, there are no more than one limit cycle in the exterior domain surrounding all three finite singularities, A 1 , S, and A 2 .
Thus, taking into account all other possibilities for limit cycle bifurcations (see [Broer et el., 2007;Broer and Gaiko, 2010;Li and Xiao, 2007;Zhu et al., 2002]), we conclude that system (2.7) (and (2.6) as well) cannot have either a multiplicity-three limit cycle or more than two limit cycles in any configuration. The theorem is proved.

Analysis of 3-Dimensional Topp Model
In this section, we study numerically the dynamics of the 3-dimensional Topp model (2.4) getting new qualitative phenomena for this model. Our particular interest is to identify the bifurcations leading to chaotic dynamics. We fix the following parameter values: b = 1, c = 1, m = 2, n = 1.
The remaining parameters α, a, and l will be used for bifurcation analysis.
We start by studying equilibrium solutions and their stability. The Topp system (2.4) has at most three equilibria which are given by , where ξ ± = 1 ± √ 1 − l. Note that E 2,− and E 2,+ coalesce in a saddle-node bifurcation which occurs for l = 1. Now assume that l = 1. In this case it follows that E 2,+ = E 2,− = (1, a − 1, 2α(a − 1)).
A straightforward calculation shows that the characteristic polynomial of the Jacobian matrix of (2.4) evaluated at E 2,± is given by p(λ) = −λ(λ 2 − T λ + D), where T = α+a and D = α(2a−1). Note that λ = 0 is a zero of p(λ); indeed this is the eigenvalue associated with the saddle-node bifurcation. For 0 < a < 1 2 and α = −a it follows that T = 0 and D > 0, which implies that p(λ) also has two imaginary zeros λ = ±i √ −a(2a − 1). In conclusion, in the three-dimensional (α, a, l)-parameter space there is a plane of saddle-node bifurcations given by l = 1 and a line segment of Hopf-saddle-node bifurcations given by (−α, α, 1) where − 1 2 < α < 0. The possible unfoldings of the Hopf-saddle-node (HSN) bifurcation are presented in Kuznetsov (2004). The HSN bifurcation is a codimension-two bifurcation which forms an organizing centre in the twodimensional (a, l)-parameter plane. From the HSN  point typically other bifurcation curves emanate, such as Hopf-Neȋmark-Sacker bifurcations which lead to quasiperiodic attractors. In addition, Shilnikov homoclinic bifurcations can occur subordinate to a HSN bifurcation [Broer & Vegter (1984)]. In certain cases, Shilnikov homoclinic are associated with the existence of chaotic dynamics and strange attractors. The HSN bifurcation and related Shilnikov bifurcations occur in many atmospheric models [Broer et al. (2002), Broer & Vitolo (2008), Crommelin et al. (2004), Sterk et al. (2010), Van Veen (2003]. We take cross sections in the parameter space by fixing α and study bifurcations and routes to chaos in the (a, l)-plane. The Lyapunov diagram in Figure 1 shows a classification of the dynamical behaviour of the Topp model in different regions of the (a, l)-parameter plane where α = −0.2 is kept fixed. The diagram suggests that periodic attractors and chaotic attractors with a positive Lyapunov exponent occur for regions in the parameter plane with positive Lebesgue measure. For other values of − 1 2 < α < 0 the Lyapunov diagrams look qualitatively similar, see Figure 2 for the case α = −0.1. Now we fix the parameters α = −0.2 and a = 0.33 and perform a more detailed bifurcation analysis by varying the parameter l. For l = 0.9999 the equilibrium E 2,− is stable. Continuation with decreasing l shows that E 2,− becomes unstable through a supercritical Andronov-Hopf bifurcation which occurs for l ≈ 0.99852. Next, we continue the periodic orbit born at the Andronov-Hopf bifurcation. For l ≈ 0.995641 the periodic orbit becomes unstable through a period doubling bifurcation. Presumably this is the first period doubling of an infinite cascade.
Continuation of the periodic orbit beyond the first period doubling bifurcation reveals the following phenomenon. The unstable periodic orbit bifurcates further through a rapid succession of saddle-node bifurcations. Presumably, infinitely many saddle-node bifurcations occur. The newly born periodic orbits themselves may bifurcate through period doubling bifurcations. Figure 3 shows a bifurcation diagram in which the period of the orbit is plotted as a function of the continuation parameter l. Clearly, the diagram suggests that the periods of the periodic orbits born through the saddle-node bifurcations tend to infinity.
The phenomenon depicted in Figure 3 can be explained as follows. During the continuation the periodic orbits born through the saddle-node bifurcations become arbitrarily close to an equilibrium. Hence, this bifurcation sequence leads to a homoclinic orbit. Figure 4 shows a periodic orbit which has a striking resemblance to a Shilnikov homoclinic orbit which is formed by an intersection of the one-dimensional unstable manifold and the two-dimensional stable manifold of the equilibrium E 2,+ . Indeed, it is expected that these Shilnikov homoclinic orbits occur along a curve in the (a, l)plane which emanates from the HSN bifurcation point [Kuznetsov (2004)]. Likewise, there may also be curve emanating from the HNS point along which there are Shilnikov homoclinic orbits which are formed by the one-dimensional stable manifold and two-dimensional unstable manifold of the equilibrium E 2,− . The numerical computation of these curves and performing a more detailed bifurcation analysis will be pursued in our forth- coming work. Finally, we explore the chaotic regime for 0.994775 < l < 0.993466. From the flow of the Topp model we numerically compute a Poincaré map by computing the intersections of the integral curves with the plane G = 0.9. Figure 5 shows a bifurcation diagram of the Poincaré map. The period doubling bifurcations of periodic attractors are clearly visible. After what is presumably an infinite cascade of period doublings we find chaotic attractors. Figure 6 shows a chaotic attractor for the parameter values (α, a, l) = (−0.2, 0.35, 0.994). In Fig-ure 7 the corresponding attractor of the Topp flow is shown.
The attractor in Figure 6 seems to have the geometric structure of a "fattened curve". In fact, we conjecture that the attractor is Hénon-like, which means that the attractor is the closure of the 1-dimensional unstable manifold of a fixed point. For the classical Hénon map the existence of such attractors has been proven by Benedicks and Carleson (1991). In turn this would imply that the attractor in Figure 7 is formed by the closure of the unstable manifold of a periodic orbit of saddle type.
Hénon-like attractors appear in many applications which range from climate models [Broer et al. (2002), Broer et al. (2011)] to control systems [Ghane et al. (2019)]. Their occurrence in the Topp model will be investigated in more detail by the authors in forthcoming work.