Geometric analysis of Oscillations in the Frzilator model

A biochemical oscillator model, describing developmental stage of myxobacteria, is analyzed mathematically. Observations from numerical simulations show that in a certain range of parameters, the corresponding system of ordinary differential equations displays stable and robust oscillations. In this work, we use geometric singular perturbation theory and blow-up method to prove the existence of a strongly attracting limit cycle. This cycle corresponds to a relaxation oscillation of an auxiliary system, whose singular perturbation nature originates from the small Michaelis-Menten constants of the biochemical model. In addition, we give a detailed description of the structure of the limit cycle, and the timescales along it.


Introduction
Oscillators are ubiquitous in different fields of science such as biology [33], biochemistry [8,11], neuroscience [15], medicine [12], and engineering [32]. In particular, biochemical oscillations often occur in several contexts including signaling, metabolism, development, and regulation of important physiological cell functions [26]. In this paper, we study a biochemical oscillator model that describes the developmental stage of myxcobacteria. Myxcobacteria are multicellular organisms that are common in the topsoil [18]. During vegetation growth, i.e. when food is ample, myxobacteria constitute small swarms by a mechanism called "gliding" [16]. In contrast, under a starvation condition, they aggregate and initiate a complex developmental cycle during which small swarms are transformed into a multicellular single body known as "fruiting body", whose role is to produce spores for next generation of bacteria [18]. During the aforementioned transition, myxobactria pass through a developmental stage called the "ripple phase" [16,18], characterized by complex patterns of waves that propagate within the whole colony.
Two genetically distinct molecular motors are concentrated at the cell poles of myxobacteria, allowing them to glide on surfaces; these two motors are called Adventurous (A-motility) and Social (S-motility) motors, respectively [16]. The role of the former is to push the cells forward, while the role of the latter is to pull them together. So, in order for a cell to reverse its direction, it has to alternatively activate its A-motility (push) and S-motility (pull) motors at opposite cell poles [16]. As a result, by forward and backward motion of myxobacteria, complex spatial wave patterns are created. In particular, wave patterns are produced by the coordination of motion of individual cells through a direct end-to-end contact signal, the "C-signal". During the ripple phase of development, the C-signaling induces reversals, while suppresses them during the aggregation stage of development. Observations from experiments resulted in proposing a biochemical oscillator in [16], known as the Frzilator, which acts as a "clock" to control reversals.
The Frzilator is detailed in Section 2.1. From our numerical simulations, it appears that this biochemical oscillator is robust under small variation of parameters. More importantly, it seems that (almost) all solutions converge to a "unique" limit cycle. Regarding the previous property, in [30] it has been shown that within a certain range of parameter values, (almost) all trajectories are oscillatory, the system has a finite number of isolated periodic orbits, at least one of which is asymptotically stable. Although some biological systems may produce more than one stable periodic solution for a certain range of parameters [4], the coexistence between multiple stable solutions has not yet been observed experimentally [9].
The main contribution of this paper is to prove that, within a certain range of parameter values, there exists a strongly attracting periodic orbit for the Frzilator. Moreover, the detailed description of the structure of such periodic orbit is given. The methodology used to prove the aforementioned result consists first on an appropriate rescaling of the original model, which leads to a slow-fast (or two timescales) system; next, we take advantage of the two timescales of the rescaled system to develop a geometric analysis via techniques of multi-timescale dynamical systems. From the multitimescale nature of the problem, it turns out that the limit cycle is in fact a relaxation oscillator, meaning that there are several timescales along the orbit of the oscillator. From an analytical point of view, the main difficulty of this analysis is the detailed description of a transition along two non-hyperbolic lines (see details in Section 3). Our analysis is based on the approach developed in [19,20] where similar mechanisms, leading to an attracting limit cycle in the Goldbeter minimal model [8], have been studied.
The rest of this paper is organized as follows. In Section 2 we introduce the model, perform some preliminary analysis on the model, and briefly introduce the tools which we are going to use in the paper. In Section 3 we give the slow-fast analysis of an auxiliary system, corresponding to the original system. More precisely, we discuss the behavior of the dynamics when ε → 0. In Section 4 we present the blow-up analysis of the non-hyperbolic parts. We conclude the paper with a discussion and outlook in Section 5.

Detailed model and preliminary analysis
In this section we provide a preliminary analysis of the biochemical oscillator proposed in [16]. We start in Subsection 2.1 by presenting a detailed description of the model under study. Furthermore, we describe the behavior of the trajectories and the role of parameters, and propose a unification of them. Afterwards, in Subsection 2.2, we present a two-parameter bifurcation analysis where we clarify the nature and the role of two distinct parameters of the system. Finally, in Subsection 2.3 we provide a brief introduction to slow-fast systems and the main techniques for their analysis.

Model description
We study a biochemical oscillator model which describes the social-behavior transition phase of myxobacteria [16]. This model, which is known as the Frzilator (or simply "Frz") model, is based on a negative feedback loop. In the Frz model, there are three proteins, namely, a methyltransferase (FrzF), the cytoplasmic methyl-accepting protein (FrzCD), and a protein kinase (FrzE). A direct and end-to-end collision of two myxobacteria results in producing a signal, so-called "C-signal", under which a protein called FruA is phosphorylated. The signal from phosphorylated FruA (FruA-P) activates the Frz proteins as follows [16]: (i) the methyltransferase FrzF (FrzF * ) is activated by the protein FruA-P; (ii) in response to FrzF * , the protein FrzCD is methylated (FrzCD-M); (iii) the phosphorylation of FrzE (FrzE-P) is activated by the methylated form of FrzCD; (iv) FrzF * is inhibited by the phosphorylated form of FrzE. Figure 1 shows a schematic representation of interactions between proteins of the Frz system. For a more detailed explanation of the model and its biological background, see [16]. Denote f, c and e respectively as the fraction of activated FrzF, methylated FrzCD, and phosphorylated FrzE. These fractions are given by [16] The interaction between the Frz proteins is modeled by Michaelis-Menten kinetics and hence leads to the dynamical system where , k dp = k max dp K dp + e .
(2) Remark 1. Due to the fact that f, c and e represent fractions of active protein concentrations, their values are restricted to [0, 1]. So the fraction of inactive protein concentrations are given by . Therefore, hereafter, our analysis is restricted to the unit cube As mentioned in [16], the Frz system has the well-known property of "zero-order ultrasensitivity" which requires that the Michaelis-Menten constants K a , K d , K m , K dm , K p and K dp have to be small [10]. It is observed numerically in [16] that for the parameter values K a = 10 −2 , dm = k max dp = 2 min −1 , and k max a = 0.08 min −1 , system (1) has an attracting periodic solution. For simplicity, we unify all the dimensionless Michaelis-Menten constants by K a = 2K d = 2K m = 2K dm = 2K p = 2K dp = ε 1. After unifying all Michaelis-Menten constants by ε, denoting γ := k max a , and substituting (2) in (1), we obtain the following dynamical system Figures 2 and 3 show numerically computed attracting limit cycle as well as time evolution of system (4) for ε = 10 −3 and γ = 0.08.
Remark 2. For our analysis in this paper, we fix γ = 0.08, while later we show that this parameter can be relaxed to some extent, see Remark 3 and Appendix Appendix A.
The dynamics along the limit cycle, shown in Fig. 2, can be summarized as follows. Initially, all protein ratios f, c and e are close to zero, under the dynamics (4), the variable f increases (due to the action of the C-signal), while c and e stay close to zero. Once the variable f passes the activation threshold f * := 0.5, the variable c increases very fast. Next, once the variable c passes the threshold c * := 0.5, the variable e is activated and also increases very fast until it reaches its maximum value, i.e., e = 1. Due to the fact that there is a negative feedback from e to f , the increase in e results in the degradation of variable f . Once f reaches the threshold f * , variable c decreases, and once c reaches the threshold c * , the variable e decreases vary fast. As a result, the variables f and c reach their lowest values (i.e. very close to zero), but the variable e reaches the threshold e * := γ. Once the variable e drops below the threshold e * , the variable f is activated and increases. This behavior is repeated in a periodic manner and a limit cycle is formed (see Figure  2).
For system (4), a parameter-robustness analysis with respect to ε and γ = 0.08 is presented in [30]. More precisely, using bifurcation analysis, it is shown that system (4) is robust under the variation of ε for ε ∈ (0, ε * ) with ε * := 0.05517665. Moreover, it is proven that for ε ∈ (0, ε * ), almost all trajectories converge to a finite number of periodic solutions, one of which is orbitally asymptotically stable. In this article, we prove the existence of a strongly attracting limit cycle which explains the numerically computed periodic orbit, for sufficiently small ε > 0.

Two-parameter bifurcation analysis
This section is devoted to the two-parameter bifurcation analysis of (4). In particular, we are interested in understanding the behavior of system (4) under the variation of parameters (ε, γ). To this end, let us represent (4) byẋ = G(x; ε, γ), where x = f c e , and G(x; ε, γ) denotes the right-hand side of (4). We have used the numerical continuation software Matcont [5] to compute the two-parameter bifurcation diagram of (5) with respect to (ε, γ), presented in Fig. 4, where the vertical and the horizontal axes show, respectively, the behavior of G(x; ε, γ) with respect to ε and γ. The blue curve indicates that for any 0 < γ < 1 and any ε below the curve, the system has unstable equilibria and hence exhibits oscillatory behavior. For those values of ε which are above the blue curve, the system is not oscillatory anymore, i.e. the equilibrium point is stable. In fact, the blue curve is a curve of Hopf bifurcations where the equilibria of the system switches from being stable to unstable: with fixed 0 < γ < 1, as ε passes through the curve from above to below, a limit cycle is generated.
In Fig. 4, the red curves are the curves of "limit points" (or saddle-node bifurcation) of cycles. For parameter values (ε, γ) between the blue and red curves in Fig. 4, at least two limit cycles exist simultaneously, i.e., for γ close to 0 or γ close to 1, with a suitable 0 < ε 1, at least one stable and one unstable limit cycle coexist.
Remark 3. As we mentioned in Section 2.1, due to the property of "zero-order ultrasensitivity", the Michaelis-Menten constants and hence ε have to be small. Our observation from numerical simulations shows that, for sufficiently small ε, system (4) has similar qualitative behaviors when γ belongs to certain bounds which are close to 0 and 1. In this regard, we emphasize that although the position of the limit cycle changes when γ is close to 1 (see, for instance, Fig. 5), the geometric analysis of the dynamics is the same as the case that γ is close to 0, for sufficiently small ε.
Remark 4. In Section 2.1, we have unified all the Michaelis-Menten constants of system (1) by ε, resulted in system (4). Although γ has similar size as the Michaelis-Menten constants, we have not unified it with them. One reason is that the unit of γ is "min −1 ", while the Michaelis-Menten constants are unitless. The other reason is that the simultaneous limit (ε, γ) → (0, 0) is very singular because a certain point (0, 0, γ), playing crucial role in our analysis, approaches (0, 0, 0) which is the intersection of three critical manifolds f = 0, c = 0, and e = 0. It would be interesting to study this limit further, which could explain the coalescence of the Hopf curve and the saddle-node curve at (0, 0), see Fig. 4. Similar remark holds as (ε, γ) → (0, 1).

Preliminaries on slow-fast systems
Our goal is to understand the dynamics of (4) for small ε in the limit ε → 0. However, as it is seen in (4), when the variables f, c and e are very close to the boundary of Q, the limiting behavior is different from the case that they are away from the boundary. To resolve the aforementioned problem, one possibility is to consider an auxiliary system which is smoothly equivalent to (4). To this end, let us define Note that H ε (f, c, e) > 0 for any ε > 0 and any (f, c, e) ∈ Q. Therefore, we can reparametrize time of system (4) by multiplying both sides of (4) in H ε (f, c, e), which leads to the following dynamical system H ε (f, c, e), where, for simplicity, we recycle τ to denote the reparametrized time. One can rewrite (8) as follows The vector field (9) is smoothly equivalent to (4) for ε > 0 [2], which from now on is the object of study. The main reason to rewrite system (4) into the form of system (9) is that the latter is a singularly perturbed ODE which allows us to analyze the system using geometric methods. Moreover, note that in contrast to (4), system (9) is polynomial, which is another of its advantages.

Slow-fast systems
A Slow-Fast System (SFS) is a singularly perturbed ordinary differential equation with two timescales often presented asẋ = F (x, y, ε), where the "dot" denotes derivative with respect to the slow time t, F and G are assumed to be smooth, x ∈ R ns , y ∈ R n f , and 0 < ε 1 is a small parameter that describes the timescale separation between x and y. The SFS presented in (10), where the timescale separation is explicitly given, is said to be in the standard form. To study standard SFSs we usually define a new fast time τ = t ε with which system (10) can be rewritten as where now the "prime" denotes d dτ . Since ε is a small parameter, we would like to draw conclusions on the overall behavior of the trajectories of a SFS from limiting systems obtained by taking the limit ε → 0. In such a limit (10) becomes a Differential Algebraic Equation 5 (DAE) of the forṁ 5 Also known as Constrained Differential Equation [31].
which is called the reduced problem. On the other hand (11) becomes the layer problem Remark 5. The SFS (11) is a regular ε-perturbation problem. Therefore, its solutions can be expressed as O(ε) perturbations of solutions of (13). However, the previous is valid only for time τ of order O(1), or equivalently for time t of order O(ε) in (10). To describe trajectories for longer time, techniques outside the scope of regular perturbation theory are needed.
Systems (12) and (13) are not equivalent. However, the critical manifold provides a relationship between the two.
Definition 1 (Critical manifold). The critical manifold is defined as Note that the critical manifold C 0 serves as the phase space of the DAE and as the set of equilibrium points of the layer problem. An important property that a critical manifold may posses is normal hyperbolicity.
Definition 2 (Normal Hyperbolicity). Consider a SFS (10) and its associated critical manifold C 0 . A point p ∈ C 0 is said to be hyperbolic if the matrix D y G(p), where D y denotes the total derivative with respect to y, has all its eigenvalues with non-zero real part. The critical manifold C 0 is said to be normally hyperbolic (NH) if every point p ∈ C 0 is hyperbolic.
Fenichel theory [22] describes the dynamics of a SFS with a normally hyperbolic critical manifold.
Theorem 1 (Fenichel). Let S 0 ⊆ C 0 be a compact and normally hyperbolic critical manifold of an SFS. Then, for ε > 0 sufficiently small, the followings hold: • There exists a locally invariant manifold S ε which is diffeomorphic to S 0 and lies within distance of order O(ε) from S 0 .
• The vector field X ε restricted to S ε is a smooth perturbation of the reduced problem.
• S ε has the same stability properties as S 0 .
In words, Fenichel theory says that if a SFS has a compact and normally hyperbolic critical manifold S 0 , the dynamics of the slow-fast systems can be inferred from the reduced flow along S 0 and the flow of the layer equation, which provide the stability properties of S 0 .
Often slow-fast systems have critical manifolds which lose normal hyperbolicity at certain points. In fact, like the system studied in this article, many interesting phenomena in several timescales such as relaxation oscillations and canards, are associated to the loss of normal hyperbolicity [23,20,1,27].
Remark 6. In a more general context, SFSs do not have to be given in standard form as in (11). That is, SFSs can be defined by and ODE of the form z = H(z, ε). In such a case the corresponding critical manifold S 0 is defined by S 0 = {z ∈ R ns+n f | H(z, 0) = 0}, while the layer equation reads as z = H(z, 0). Under normal hyperbolicity of the critical manifold, all the Fenichel theory results hold for the aforementioned general case [22]. When the critical manifold has nonhyperbolic points, a careful combination of Fenichel theory and the blow-up method can be employed for a detailed analysis of the dynamics of the SFS. In the following subsection we briefly describe the blow-up method. We later show that (9) is indeed a general SFS, and provide a detailed geometric analysis of (9) by means of the blow-up method and Fenichel theory.

The blow-up method
The blow-up method was introduced to describe the dynamics of SFSs near non-hyperbolic points, and is the main mathematical technique used in forthcoming section of this article. Here we just provide a brief description of the method, for more details the interested reader is referred to [6,21,22,17].
First of all, note that a SFS written in the fast-time scale is an ε-parameter family of vector fields. Thus, it is convenient to lift such family up and instead consider a single vector field of the form X : Definition 3. Consider a generalized polar coordinate transformation where ns i=1x For the purposes of this article, it is sufficient to letr ∈ [0, ρ), with ρ > 0. The main idea of the blow-up method is to construct a new, but equivalent, vector field to X, which is defined in a higher dimensional manifold, but whose singularities are simpler compared to those of X.

Definition 4. The blown up vector fieldX is induced by the blow-up map asX
Note that the vector fieldsX andX are equivalent on S ns+n f ×{r > 0}. Moreover, if the weights (α, β, γ) are well chosen, the singularities ofX|r =0 are partially hyperbolic or even hyperbolic, making the analysis ofX simpler than that of X. Due to the equivalence between X andX, one obtains all the local information of X around 0 ∈ R ns+n f +1 from the analysis ofX around While doing computations, it is more convenient to study the vector fieldX in charts. A chart is a parametrization of a hemisphere of S ns+n f × I and is obtained by setting one of the coordinates (x,ȳ,ε) ∈ S ns+n f to ±1 in the definition of Φ. For example, one of the most important charts in the blow-up method is the central chart defined by Kε = {ε = 1}. After we study the dynamics in the relevant charts, we connect the flow together via transition maps, allowing us a complete description of the flow ofX near S ns+n f × {0}. In turn, and as mentioned above, the flow ofX is equivalent to the flow of X for ε > 0 sufficiently small. For more details see Section 4 and [22].
Remark 7. It is also possible to blow-up only some of the variables in the system (15), and keep the others unchanged. In this paper, we blow-up a non-hyperbolic line of equilibria to a cylinder, see Section 4.

Geometric singular perturbation analysis
The goal of this section is to give the detailed analysis of the slow-fast structure of the auxiliary system (9).

Layer problem and the critical manifold
Setting ε = 0 in (9) results in the layer problem with Apart form the isolated equilibrium point P := (0.5, 0.5, γ), which is inside the cube Q, the boundary of Q, which consists of six planes, is the equilibria set of the layer problem (17). We denote each plane of equilibria by S 0,i (i = 1, 2, ..., 6) as follows: S 0,i is the critical manifold. The stability of system (9) changes at lines f ∈ S 0,2 , f ∈ S 0,5 (given by f = f * ); c ∈ S 0,3 , c ∈ S 0,6 (given by c = c * ); and e ∈ S 0,1 , e ∈ S 0,4 (given by e = e * ). Moreover, the 12 edges of the unit cube, where the 6 planes S 0,i intersect, are nonhyperbolic lines as well. However, for our analysis, only the lines 1 = S 0,1 ∩S 0,2 and 2 = S 0,2 ∩S 0,3 are crucial (see Figure 6). The stability of points in S 0 is summarized in the following lemma.
, e, f , c , e in red, all 12 non-hyperbolic edges in blue, and in particular, the two non-hyperbolic edges 1 and 2 shall play an important role in our analysis.
• S 0,1 is attracting for e > e * and repelling for e < e * .
• S 0,2 is attracting for f < f * and repelling for f > f * .
• S 0,3 is attracting for c < c * and repelling for c > c * .
• S 0,4 is attracting for e < e * and repelling for e > e * .
• S 0,5 is attracting for f > f * and repelling for f < f * .
• S 0,6 is attracting for c > c * and repelling for c < c * .
Proof. The eigenvalues of the linearization of system (17) at points, e.g., in the plane S 0,1 are given by It is clear that λ 3 is zero at the boundary of S 0,1 , and also along the line l e given by e = e * . Therefore, S 0,1 is attracting for e > e * and it is repelling for e < e * . The proof of the other cases is performed analogously.
We denote the interior of the cube Q byQ. Note that when (f, c, e) ∈Q, the layer problem (17) can be divided by the positive term . Therefore away from the critical manifold S, all the variables evolve on the fast time scale τ and the orbits of the layer problem (17) are identical to the orbits of the linear system Remark 8. System (19) is the limit of (4) when ε → 0 and (f, c, e) ∈Q.
3.2. Reduced problem, slow manifolds, and slow dynamics From Subsection 3.1, we know that the boundary of Q is the critical manifold S 0 . Any compact subset of S 0 that does not contain any non-hyperbolic point is normally hyperbolic, and hence Fenichel theory [7] is applicable. In other words, this theory implies that the normally hyperbolic parts of S 0 perturb to slow manifolds, which lie within a distance of order O(ε) of the critical manifold S 0 . In the following, we compute the slow manifolds and analyze the reduced flows in the planes S 0,1 , S 0,2 , S 0,3 and S 0,6 which are essential for our analysis.
Note that (22) reflects the fact that the manifold S a ε,1 is not well-defined when e = γ. Thus, the invariant manifold S a ε,1 is given as stated in the lemma, which completes the proof.
For the sake of brevity, we summarize the analysis in the planes S 0,2 , S 0,3 and S 0,6 in Table 1, which is shown by following the same line of reasoning as the one of Lemma 2. For more details, the interested reader is referred to [29].
Remark 9. Similar results can be obtained for the "repelling" parts S r ε,i , i = 1, 2, ..., 6. However, these are not needed in our analysis. Nonetheless, we point out that the slow manifolds S r ε,i would be expressed by the same functions h ε,i and appropriate intervals I r i . Remark 10. The expansions of the functions h ε,i (·, ·), i = 1, 2, 3, 6, also explain why it is necessary to restrict the domain of definition of the slow manifolds to I a i to exclude their singularities. We now turn to the analysis of the reduced flows in the planes S 0,1 , S 0,2 , S 0,3 and S 0,6 which, respectively, means the planes f = 0, c = 0, e = 0 and e = 1. We know that system (9) has the fast time scale τ . By substituting the functions h ε,i , i = 1, 2, 3, 6 into (9), transforming the fast time variable to the slow one by t = ετ , and setting ε = 0, the equations governing the slow dynamics on the critical manifold S 0,i are computed. In the following, we give the analysis in the plane S 0,1 .
After substituting h ε,1 into system (9), the dynamics of the reduced system in S 0,1 , i.e., on the plane f = 0, is governed by where denotes the differentiation with respect to τ . Now by dividing out a factor of ε, which corresponds to switching from the fast time variable to the slow one, we havė where the overdot represents differentiation with respect to t = ετ . Now, by setting ε = 0 in (24), the reduced flow on S 0,1 is given bẏ As it is clear, the vector field (25)   which can be integrated explicitly.
Remark 11. For e > e * , systems (25) and (26) have qualitatively the same dynamics when c, e ∈ (0, 1). In particular, the vector field (26) is C ∞ -equivalent but not C ∞ -conjugate to the vector field (25). For the case that e < e * , the direction of the vector field (25) is not preserved in the vector field (26). However, for our analysis, it suffices to study the flow of system (25) when e > e * , or equivalently on S a 0,1 .

Lemma 3.
For e > e * , the reduced flow (25) on S 0,1 and hence the slow flow (24) on S a ε,1 maps section {c = c 1 } to {c = c 2 }, where 0 < c 2 < c 1 < 1 2 ; this map is well-defined and its first derivative with respect to e is equal to one.
Proof. It suffices to consider (26). Let Π(e) denote the map from {c = c 1 } to {c = c 2 } induced by the flow of (26). Then, it is straightforward to get Π(e) = e + c 2 − c 2 2 − c 1 + c 2 1 , from which the statement follows.
In order to obtain the equations governing the slow flow along S a ε,2 , S a ε,3 and S a ε,6 , a similar analysis can be done by inserting the functions h ε,2 , h ε,3 and h ε,6 into (9) and dividing out a factor of ε, which corresponds to switching to the slow time scale t = ετ . Next, by setting ε = 0 one obtains the reduced flow on the critical manifolds S 0,2 , S 0,3 and S 0,6 . For the sake of brevity, we have summarized the slow flows along S 0,2 , S 0,3 and S 0,6 in Fig. 8. For more details, the interested reader is referred to [29].

Singular cycle
In this section, we present the overall behavior of the singular cycle, which is a closed curve consisting of alternating parts of the layer problem, and the critical manifold S 0 . However, by For f > f * , the reduced flow (28) contracts the variable f between sections c = c 1 and c = c 2 with 0 < c 1 < c < c 2 < c * .
For f < f * , the reduced flow (29) contracts the variable f between sections c = c 1 and c = c 2 with c * < c 1 < c < c 2 < 1. the information that we have so far from the critical manifold and the layer problem, we cannot fully describe the singular cycle close to the non-hyperbolic lines 1 and 2 . A full description of the singular cycle for those parts that cannot be derived from the critical manifold and the layer problem is presented in Section 4 by the blow-up method.

Remark 13.
At the singular level, there is no visible flow on the segments ω 6 and ω 8 . The blow-up analysis, carried out in Section 4, will reveal a hidden flow for such segments.

Main result
In view of the singular cycle Γ 0 , introduced in the previous subsection, we are now ready to present the main result.
Theorem 2. Assume that Γ 0 is the singular cycle described in Section 3.3. Then for sufficiently small ε > 0, there exists a unique attracting periodic orbit Γ ε of the auxiliary system (9), which tends to the singular cycle Γ 0 as ε → 0.
In order to prove Theorem 2, we need to introduce the following sections where R j (j = 1, 2, 3) are suitable small rectangles, and δ j are chosen sufficiently small. Note that Σ 1 is transversal to ω 4 , Σ 2 is transversal to ω 6 , and Σ 3 is transversal to ω 8 , see Fig. 9.
According to the definition of the sections Σ i , introduced in (31), we define the following Poincaré maps for the flow of the system (9) where the map π 1 describes the passage from Σ 1 to Σ 2 along the non-hyperbolic line 1 , the map π 2 describes the passage from Σ 2 to Σ 3 along the non-hyperbolic line 2 , and the map π 3 describes the passage from Σ 3 to Σ 1 . The map π 3 consists of slow flow along S a ε,3 , followed by the fast dynamics from a neighborhood of p 1 to a neighborhood of p 2 , followed by the slow flow along S a ε,6 to a neighborhood of q e . Through the fast dynamics, this neighborhood is mapped to a neighborhood of q e , followed by the slow flow along S a ε,1 to Σ 1 . We summarize the properties of the above maps in the following lemmas.
Lemma 5. If the section Σ 2 is chosen sufficiently small, then there exists ε 0 > 0 such that the map is well-defined for ε ∈ [0, ε 0 ] and smooth for ε ∈ (0, ε 0 ]. The map π 2 is a strong contraction with contraction rate exp(−K/ε) for some K > 0. The image of Σ 2 is a two-dimensional domain of exponentially small size, which converges to the point q 3 := Σ 3 ∩ ω 1 as ε → 0.
The proofs of Lemmas 4 and 5 are based on the blow-up analysis of the lines 1 and 2 , respectively, which will be presented in Subsections 4.1 and 4.2.
Remark 14. The points on the line c when 0.5 < f < 1, and on the line c when 0 < f < 0.5 are jump points, i.e., the trajectory switches from the slow dynamics to the fast dynamics. Further, it is shown that this behavior is very similar to the behavior of standard slow-fast systems with two slow variables and one fast variable near a generic "fold" line, studied in [28] based on the blow-up method. The critical manifolds S 0,3 and S 0,6 of system (9) can be viewed as a standard folded critical manifold, which has been straightened out by a suitable diffeomorphism. This leads to the curved fibers of the layer problem (17). Therefore, we can use the results of [28] to understand the behavior of (9) close to the non-hyperbolic lines c and c .
Proof. The basic idea of the proof is based on the map that has been already described in Fig. 9 for ε = 0, denoted by π 0 3 , and then treat π 3 as an ε-perturbation of π 0 3 . If the section Σ 3 is chosen sufficiently small, then the trajectories starting in Σ 3 can be described by the slow flow along the manifold S a ε,3 combined with the exponential contraction towards the slow manifold until they reach a neighborhood of the jump points on the line c . Applying [28,Theorem 1] close to the jump pints, the trajectories switch from the slow dynamics to the fast dynamics, and hence pass the non-hyperbolic line c ; this transition is well-defined for ε ∈ [0, ε 1 ], and smooth for ε ∈ (0, ε 1 ] for some ε 1 > 0. Note that [28, Theorem 1] guarantees that the contraction of the solutions in the e-direction persists during the passage through the fold-line c , as it is at most algebraically expanding. After that, the solutions follow the fast dynamics ω 2 until they reach a neighborhood of the point p 2 , see Fig. 9. Next, the solutions follow the slow flow along the manifold S a ε,6 combined with the exponential contraction towards the slow manifold until they reach a neighborhood of the point q e . Again applying [28, Theorem 1] close to the jump points, the solutions which are very close to the non-hyperbolic line c switch from the slow dynamics to the fast dynamics, and hence pass the non-hyperbolic line c , where the corresponding transitions are well-defined for ε ∈ [0, ε 2 ], and smooth for ε ∈ (0, ε 2 ] for some ε 2 > 0, and then follow the fast dynamics (ω 4 ) until they reach a neighborhood of the point q e . Finally, the solutions follow the slow flow along the manifold S a ε,1 combined with the exponential contraction towards the slow manifold until they reach the section Σ 1 . Theorem 1 of [28] implies that the map π 3 is at most algebraically expanding in the direction of e when Σ 3 is chosen sufficiently small. On the other hand, the slow manifold S a ε,1 is exponentially contracting in the direction of f (Fenichel theory). Therefore, the image of Σ 3 is a thin strip lying exponentially close to S a ε,1 ∩ Σ 1 . Hence, the statements of the lemma follow. Now we are ready to give the proof of the main result. Proof of Theorem 2. Let us define the map π : Σ 3 → Σ 3 as a combination of the maps π j (j = 1, 2, 3), described in Lemmas 4, 5 and 6. More precisely, we define π = π 2 • π 1 • π 3 : Σ 3 → Σ 3 .

Blow-up analysis
The slow-fast analysis that we have done in Section 3 does not explain the dynamics of system (9) close to the non-hyperbolic lines 1 and 2 . As the segments ω 5 and ω 7 lie on these lines (see Fig.  9), we need a detailed analysis close to the lines 1 and 2 , which is carried out in this section via the blow-up method [22,14,21]. To apply this, we extend system (9) by adding ε as a trivial dynamic variable and obtain df dτ where H ε 1 (f ), H ε 2 (c) and H ε 3 (e) are defined in (7). Note that for the extended system (36), the lines 1 × {0} and 2 × {0} are sets of equilibria. Due to the fact that the linearization of (36) around these lines has quadruple zero eigenvalues, system (36) is very degenerate close to 1 × {0} and 2 × {0}. To resolve these degeneracies, we use the blow-up method, given in next subsections.
For the analysis of system (36) near the line 1 × {0}, we define three charts K 1 , K 2 and K 3 by settingc = 1,ε = 1, andf = 1 in (37), respectively: The changes of coordinates for the charts K 1 to K 2 , and K 2 to K 3 in the blown-up space are given in the following lemma.
Lemma 7. The changes of coordinates K 1 to K 2 , and K 2 to K 3 are given by The goal of this subsection is to construct the transition map π 1 : Σ 1 → Σ 2 , defined in (32), and prove Lemma 4. Before going into the details, let us briefly describe our approach. We describe the transition map π 1 : Σ 1 → Σ 2 via an equivalent one in the blown-up space. More specifically we define is the cylindrical blow-up defined by (37), the maps Π i are local transitions induced by the blown-up vector fields which are detailed below, and κ 12 and κ 23 denote the changes of coordinates, given in Lemma 7.π 1 is the transition map in the blown-up space and due to the fact that Φ is a diffeomorphism, it is equivalent to π 1 . A schematic of the problem at hand is shown in Fig. 10.
The left picture in Fig. 10 illustrates the critically manifolds S a 0,1 and S a 0,2 , and the corresponding flows in blue. The non-hyperbolic line 1 is shown in orange. For e > γ, the reduced flows on both critically manifolds approach the line 1 . At the point on the line 1 with e = γ, a transition from S a 0,1 to S a 0,2 is possible as indicated in the figure. The right picture in Fig. 10 schematically shows the configuration in the blown-up space. The cylinder corresponding to r = 0 is show in orange. The part of the phase space corresponding toε = 0 and r > 0 are shown outside of the the cylinder. Here we recover the layer problem, the critically manifolds, and the reduced flows inS a 0,1 andS a 0,2 . In the blown-up space, the manifolds S a 0,1 and S a 0,2 are separated and hence gained hyperbolicity, in particular they are attractive, as indicated below in Fig. 11a. All these assertions will be proven in this section.
Roughly speaking, in chart K 1 we continue the attracting slow manifoldS a 0,1 onto the cylinder. Chart
2. ε 1 = 0: in this case, the dynamics (44) is represented by (47) From (47), one concludes that the plane f 1 = 0 is the plane of equilibria which is denoted byS a 0,1 , see Fig. 11a. The non-zero eigenvalue alongS a 0,1 is given by λ = 32e 1 (1 − e 1 )(1 − r 1 )(γ − e 1 ). For 0 ≤ r 1 < 1 and e 1 > γ, the planeS a 0,1 is attracting. As the e 1 -axis is a part ofS a 0,1 , we denote that part of the e 1 -axis that γ ≤ e 1 ≤ 1 by e1 . We also have another curve of equilibria which is defined by r 1 = 0, and f 1 = e1−γ 2 , denoted by M r 1 , see Fig. 11a. This curve of equilibria is of saddle-type with the eigenvalues λ = ±32e 1 (e 1 − 1)(e 1 − γ). Note that we have recovered the information of the previous case here. 3. r 1 = 0: in this case, the dynamics (44) is represented by By setting ε 1 = 0, we again have the line e1 and the curve M r 1 . The Jacobian matrix at a point in e1 has two eigenvalues: one zero and the other one is λ = 32e 1 (1 − e 1 )(γ − e 1 ). So the line e1 is attracting when e > γ. As in this case we have two zero eigenvalues, it implies that there exists a two-dimensional center manifold, namely, C a,1 .
Remark 15. In chart K 1 , the most important role is played by the two-dimensional center manifold C a,1 , see Lemma 9. In fact, this is the continuation of the critical manifoldS a 0,1 . We summarize the analysis performed in this subsection in the following lemmas. 1. The linearization of (44) alongS a 0,1 has three zero eigenvalues, and the nonzero eigenvalue λ = 32e 1 (1 − e 1 )(1 − r 1 )(γ − e 1 ), which for r 1 = 0 corresponds to the flow in the invariant plane (f 1 , e 1 ).

2.
There exists a three-dimensional center manifold W c a,1 of the line e1 which contains the plane of equilibriaS a 0,1 and the two-dimensional center manifold C a,1 . The manifold W c a,1 is attracting, and in the set D 1 , defined by is given by the graph where I 1 is a suitable interval, and α 1 , δ 1 > 0 are sufficiently small. For the particular point p a,1 ∈ e1 where e 0 ∈ I 1 , the function h a,1 (r 1 , e 0 , ε 1 ) has the expansion h a,1 (r 1 , e 0 , ε 1 ) = γ 2(e 0 − γ)   3. There exists K > 0 such that the orbits that are near the center manifold W c a,1 are attracted to W c a,1 by an exponential rate of order O(exp(−Kt 1 )).
Proof. A straightforward calculation shows the first claim. Due to the fact that the linearization of (44) alongS a 0,1 has three zero eigenvalues, there exists [3,13] an attracting three-dimensional center manifold W c a,1 at the point p a,1 . To derive equation (49), we first expand f 1 to the first order of variables r 1 , e 1 and ε 1 , and then plug into (44). By comparing the coefficients of r 1 , e 1 and ε 1 , equation (49) is obtained. The last claim is proven by the center manifold theory applied at the point p a,1 .
Remark 16. The attracting center manifold W c a,1 recovers parts of the slow manifold S a ε,1 away form the line 1 , and extends it into an O(ε) neighborhood of 1 . The slow manifold S a ε,1 is obtained as a section ε = constant of W c a,1 . In chart K 1 , this center manifold is given by the graph (49).
Note that in chart K 1 , our goal is to understand the dynamics (44) close to the center manifold W c a,1 , which corresponds to a sufficiently small neighborhood of the slow manifoldS a 0,1 . Assume that δ 1 , α 1 , β 1 > 0 are small constants. Let us define the sections Note that by the way we have defined ∆ in 1 , we in fact have ∆ in Fig. 10. Furthermore, the constants δ 1 , α 1 , β 1 are chosen such that R in 1 ⊂ ∆ in 1 , and the intersection of the center manifold W c a,1 with ∆ in 1 lies in R in 1 , i.e., W c a,1 ∩ ∆ in 1 ⊂ R in 1 .
Let us denote Π 1 as the transition map from ∆ in 1 to ∆ out 1 , induced by the flow of (44). In order to construct map Π 1 , we reduce system (44) to the center manifold W c a,1 and analyze the system based on the the dynamics on W c a,1 . To this end, by substituting (49) into (44) and rescaling time, the flow of the center manifold is given by where the derivative is with respect to the new timescale, namely, t 1 . Now let us consider a solution of (51), namely, (r 1 (t 1 ), e 1 (t 1 ), ε 1 (t 1 )) which satisfies the following conditions: From equation ε 1 = ε 1 with the conditions ε 1 (0) = ε in 1 and ε 1 (T out ) = α 1 , we can calculate the time that (r 1 (t 1 ), e 1 (t 1 ), ε 1 (t 1 )) needs to travel from ∆ in 1 to ∆ out 1 , which is given by Since we can estimate the time evolution of e 1 (t 1 ), which is given by Hence, in view of (53), one has e 1 (T out ) = e out 1 := r in We summarize the analysis performed for chart K 1 in the following theorem.
Theorem 3. For system (44) with sufficiently small δ 1 , α 1 , β 1 and R in 1 ⊂ ∆ in 1 , the transition map is well-defined and has the following properties: is a three-dimensional wedge-like region in ∆ out 1 .

The transition map Π 1 is given by
where e out 1 is given in (55), Ψ(·) is an exponentially small function, and h a,1 (·) is of order O(ε 1 ), due to (49).

Analysis in chart K 2
After substituting (39) into (36) and dividing out all the equations by the common factor r 2 , the equations governing the dynamics in chart K 2 are given by Due to the fact that r 2 = ε in chart K 2 , we have presented (56) in terms of ε. Note that since r 2 = ε = 0, system (56) is a family of three-dimensional vector fields which are parametrized by ε. Moreover, system (56) is a slow-fast system in the standard form, i.e., e 2 is the slow variable, and f 2 and c 2 are the fast variables. Since the differentiation in (56) is with respect to the fast time variable, namely τ 2 , by transforming it to the slow time variable we have t 2 = ετ 2 , and hence where the derivative is with respect to t 2 . Now by setting ε = 0 in (56) we obtain the corresponding layer problem which has the associated critical manifold c 2 = 0 and f 2 = γ 2(e2−γ) , denoted by N 0 2 (see Fig. 12). The Jacobian matrix corresponding to (58) along this critical manifold has the eigenvalues As it is clear form (59), the critical manifold restricted to e 2 ∈ (γ, 1) is normally hyperbolic, and specially, is fully attracting since both of the eigenvalues are negative. As e 2 approaches γ from above, f 2 develops a singularity along N 0 2 . Thus, the behavior of N 0 2 as e → γ has to be studied in chart K 3 . Using Fenichel theory and the dynamics in chart K 2 for ε = 0, one is able to describe the dynamics for 0 < ε 1 in this chart, i.e., there exists a slow manifold N ε 2 which is the ε-perturbation of N 0 2 . We summarize the properties of the critical manifold of chart K 2 in the following lemma.
We summarize the analysis performed in chart K 2 in the following lemma.
Proof. The transition map Π 2 : ∆ in 2 → ∆ out 2 is described by Fenichel theory, i.e., all orbits starting from ∆ in 2 are attracted by the slow manifold N ε 2 , with a contraction rate exp(−K/ε) for some K > 0, and after some time they reach the section ∆ out 2 .
Remark 18. The slow manifold N ε 2 corresponds to the perturbation of N 0 2 when ε = constant. The family of all such manifolds is denoted by N 2 .

Analysis in chart K 3
Solutions in chart K 2 which reach the section ∆ out 2 must be continued in chart K 3 . For this reason, we continue our analysis in chart K 3 . After substituting (40) into (36) and dividing out all the equations by the common factor r 3 , the equations governing the dynamics in chart K 3 are given by where we denote System (63) has three invariant subspaces, namely, r 3 = 0, ε 3 = 0 and their intersection. Recall that by definition e = e 3 and thus 0 < e 3 < 1. 2. ε 3 = 0 and r 3 ≥ 0: In the invariant plane ε 3 = 0, the dynamics is governed by . Recall that c = r 3 c 3 and therefore V (r 3 , c 3 , e 3 ) > 0. The equilibria of the system are the plane c 3 = 0, denoted byS 0,2 , and the curve of equilibria given by c 3 = 2 e3−γ , denoted by M r 3 . The change of stability of the points inS 0,2 occurs at r 3 = 0.5, i.e., for r 3 < 0.5 the points are attracting, while for r 3 > 0.5 they are repelling. We denote the attracting part ofS 0,2 byS a 0,2 . The e 3 -axis, which we denote by e3 , is a boundary ofS a 0,2 , which is a line of equilibria.
The equilibria of the system are the planes c 3 = 0, and the line ε 3 = 2(e3−γ) γ , denoted by N 0 3 . The Jacobian of (66) along the curve N 0 3 has the eigenvalues and hence N 0 3 is fully attracting. In fact, N 0 3 is exactly the critical manifold N 0 2 that we found in chart K 2 . In other words, N 0 3 is the image of N 0 2 under the transformation κ 23 , defined in (42).
Remark 20. Note that the attracting manifold N 0 2 that is unbounded in chart K 2 , is now bounded in chart K 3 . So the behavior of the critical manifold that is not visible in chart K 2 when e → γ, is now visible in chart K 3 . For e 3 = γ, the critical manifold N 0 3 intersects the line e3 at the non-hyperbolic point q e3 = (e 3 , c 3 , ε 3 ) = (γ, 0, 0).
We summarize the analysis of the invariant planes, performed in this subsection, in the following Lemma.
Remark 21. Note that the dynamics in the invariant plane ε = 0 corresponds to the reduced flow on S 2 a in the original system. Summarizing the analysis, we have the following lemma.
Lemma 13. The following properties hold for system (70): 1. The curve N 0 3 has a one dimensional stable manifold, and a two dimensional center manifold away from the point q e3 . 2. The linearization of (70) at the points in e3 is given by 3. The point q e3 is nilpotent.
As we already mentioned, our goal in chart K 3 is to describe the dynamics (63) close to the line e3 , and especially at the point q e3 . To this end, we defined the map is transversal to the slow manifold in the plane ε 3 = 0 for e < γ. From Lemma 13 we know that the point q e3 is nilpotent. Thus, in order to describe the transition map Π 3 we need to blow-up the point q e3 . For such a point, a similar analysis has been carried out in [19,Theorem 5.8], in view of which we have the following theorem.
3 is a small rectangle centered at the intersection point N 0 3 ∩∆ in 3 . For sufficiently small α 3 , the transition map Π 3 : R 3 → ∆ out 3 induced by the flow of (70) is welldefined and satisfies the following properties:

2.
Restricted to the lines r 3 = constant in R 3 , the map is contracting with the rate exp(−K/r 3 ) for some K > 0. 3. The image Π 3 (R 3 ) is an exponentially thin wedge-like containing the curve σ out 3 .
Finally, if we set α 3 = δ 2 (recall the definition of Σ 2 ) we actually have that ∆ out In the above subsections, we have presented the detailed analysis of the blow-up of the nonhyperbolic line 1 × {0} in charts K 1 , K 2 and K 3 , which has been summarized in Fig. 15. A summary of the analysis, carried out in such charts, is as follows. First of all, the critical manifolds S 0,1 (i.e., f = 0) and S 0,2 (i.e., c = 0) intersect in the non-hyperbolic line 1 , which is replaced by the orange cylinder, see Figs. 10 and 15. Note that in Fig. 15, the orbitsω 5 andω 7 in the blown-up space correspond to the orbits ω 5 and ω 7 , respectively. The point at whichω 5 reaches the cylinder is denoted byq e , and the point at whichω 7 starts is denoted byq e . Starting from the sectionΣ 1 , the trajectory follows the orbitω 5 onf = 0 until it reaches the pointq e . Our analysis in chart K 1 (Lemma 9) shows that there exists a three-dimensional attracting center manifold which is the continuation of the family of orbits (indexed by ε) of the attracting slow manifold S a ε,1 . This allows us to connect the family S a ε,1 into the chart K 2 which is inside the cylinder (see the thick orange manifold fromq e toN 0 in Fig. 15). Our analysis in chart K 2 (Lemma 10) shows that the slow manifold N 0 2 is normally hyperbolic and stable. Therefore, the family S a ε,1 is exponentially attracted by the slow manifold N ε 2 . Next, our analysis in chart K 3 shows that the unbounded critical manifold N 0 2 (see Figs. 12,14) limits in the point q e3 , which is exactly the pointq e in Fig.  15. Next, our analysis in chart K 3 (see Lemma 13 and Fig. 14) demonstrates that the unbounded critical manifold N ε 2 (see Figs. 12 and 14) limits at the point q e3 , which is exactly the pointq e in Fig. 15. In addition, we have proven that the point q e3 is degenerate, i.e., the linearization of the dynamics at q e3 has a nonzero (stable) eigenvalue and a triple zero eigenvalue (see Lemma 13), which allows us to construct a three-dimensional center manifold at the point q e3 . Now, by following the family N 2 along such a center manifold, we conclude (Lemma 4) that the continuation of N ε 2 for a sufficiently small ε > 0 intersects the sectionΣ 2 in a point, namely, (α 3 , c 3 (ε 3 ), e 3 (ε 3 ), ε 3 ) ∈Σ 2 , for some ε 3 ∈ [0, β 3 ], which is exponentially close to the slow manifold S a ε,2 . Note that the point (α 3 , c 3 (ε 3 ), e 3 (ε 3 ), ε 3 ) converges to the point q 2 := Σ 2 ∩ ω 7 as ε 3 → 0. All these analyses in charts K 1 , K 2 , and K 3 show that the transition mapπ 1 :Σ 1 →Σ 2 is well-defined for ε ∈ [0, ε 0 ] and is smooth for ε ∈ (0, ε 0 ], for some ε 0 > 0.
We are now ready to prove Lemma 4.
Proof of Lemma 4. The proof is carried out by constructing the map π 1 : Σ 1 → Σ 2 for ε > 0 as where Φ is given by (37), Φ −1 is the corresponding blown-up transformation, andπ 1 :Σ 1 →Σ 2 is a transition map which can equivalently be regarded as The proof is based on the corresponding transition mapπ 1 :Σ 1 →Σ 2 in the blown-up space and interpreting the result for fixed ε ∈ [0, ε 0 ] with ε 0 > 0. Recall that the transitionπ 1 :Σ 1 →Σ 2 is equivalent to the transition map π 1 : Σ 1 → Σ 2 in the sense that it has the same properties. Furthermore, via the matching maps κ ij defined in Lemma 7, we have appropriately identified the relevant sections in each of the charts, allowing us to follow the flow of the blown-up vector field along the three charts.
As summarized above, the transition mapπ 1 :Σ 1 →Σ 2 is well-defined for ε ∈ [0, ε 0 ] and smooth for ε ∈ (0, ε 0 ] for some ε 0 > 0. It remains to prove thatπ 1 is a contraction. From Lemma 3 we know that the solutions started inΣ 1 are contracting, see (Fig. 7). This family of orbits is continued to chart K 2 by spending an O(1)-time on the time scale t 2 of system (57). This continuation persists (Theorem 4) during the passage near the point q e3 in chart K 3 until it reaches the sectionΣ 2 . As the contraction persists fromΣ 1 toΣ 2 , one concludes thatπ 1 is a contraction. This completes the proof.
Recall that the goal of Lemma 5 is to describe the map π 2 : Σ 2 → Σ 3 in the original space. In this subsection, we present a sketch of the proof of Lemma 5 by constructing the corresponding map π 2 :Σ 2 →Σ 3 in the blown-up space, and interpreting the results for fixed ε ∈ [0, ε 0 ] for some ε 0 > 0. For the sake of brevity, we have summarized the analysis of the blow-up of the non-hyperbolic line 2 × {0} in Fig. 16.
First of all, note that the non-hyperbolic line 2 , which is the intersection of the critical manifolds c = 0 and e = 0, has been blown-up to the orange cylinder (see Fig. 16). We have illustrated the slow flows in the planes c = 0 and e = 0 in blue. Note that the orbitsω 7 andω 1 which are in the blown-up space correspond, respectively, to the orbits ω 7 and ω 1 in the original space (see Figs. 9 and 16). As it is shown in Fig. 16, the intersection of the cylinder withω 7 andω 1 is denoted byp f andp f , respectively.
Our analysis in chartK 1 proves that there exists a three-dimensional attracting center manifold at the pointp f , which is the continuation of the family indexed by ε of the attracting slow manifold S a 0,2 . In view of such a center manifold, the family of the slow manifoldS a 0,2 enters the chartK 2 . Our analysis in chartK 2 proves that there exists a hyperbolic attracting one-dimensional slow manifoldÑ 0 , which attracts the interior of the cylinder. Our analysis in chartK 3 shows that the critical manifoldÑ 0 limits at the pointp f (see Fig. 16). Note thatp f is a degenerate point, i.e., the linearization of the blown-up dynamics in chartK 3 atp f has a stable eigenvalue, and a triple zero eigenvalues which allows us to construct a three-dimensional attracting center manifold. Therefore the family of flows follows such a center manifold and then intersects the section Σ 3 in a point (f (ε), δ 3 , e(ε)), for some δ 3 > 0, which is exponentially close to the slow manifold S a ε,3 and converges to the point q 3 := Σ 3 ∩ ω 1 as ε → 0. This proves that the transition mapπ 2 :Σ 2 →Σ 3 and hence π 2 : Σ 2 → Σ 3 are well-defined for ε ∈ [0, ε 0 ] and also are smooth for ε ∈ (0, ε 0 ], for some ε 0 > 0. The proof of contraction of the transition map π 2 follows the same line of reasoning as that of the map π 1 , and hence is omitted for brevity.

Conclusions
In this work we have studied a model that describes several important properties of myxobacteria during development [16]. This model, which is in line with observation from experiments [16], acts as an internal clock to control the gliding motions in myxobacteria. When two cells collide with each other, the speed of the clock in both cells is affected, some spatial wave patterns are created, and hence leads to synchronization of cells, i.e., fruiting body. The model presented in [16] can reproduce observed spatial patterns in experiments, and furthermore, it can explain both the cellular oscillations and the developmental stage of myxobacteria from vegetative swarming to the rippling phase and hence to the formation of the fruiting body.
The model, described by a system of three ordinary differential equations, has oscillatory behavior for certain parameter values, and sufficiently small Michaelis-Menten constants which we have unified them by a parameter ε. We have analyzed the dynamics of this oscillator in the limits of ε, and proven that for sufficiently small ε, there exists a strongly attracting limit cycle. The geometric method could be pushed to analyze the global uniqueness of the limit cycle which is clearly of great interest from both the mathematical and biological point of view. This requires a more global analysis of the singular flows, and in particular, connecting orbits between the critical manifolds S 0,i by orbits of the layer problem. As the layer problem is linear, this is possible. Our approach has been based on the geometric perturbation analysis and blow-up method. The geometric perturbation theory and geometric desingularization by several blow-ups allow us to fully understand the structure of the limit cycle. We emphasize that the approach and tools presented in this paper, i.e. geometric singular perturbation theory and the blow-up method, are not limited