Geometric analysis of oscillations in the Frzilator model

(1)

Geometric analysis of oscillations in the Frzilator model

Taghvafard, Hadi; Jardon Kojakhmetov, Hildeberto; Szmolyan, Peter; Cao, Ming

Published in:

Journal of Mathematical Analysis and Applications DOI:

10.1016/j.jmaa.2020.124725

IMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite from it. Please check the document version below.

Document Version

Publisher's PDF, also known as Version of record

Publication date: 2021

Link to publication in University of Groningen/UMCG research database

Citation for published version (APA):

Taghvafard, H., Jardon Kojakhmetov, H., Szmolyan, P., & Cao, M. (2021). Geometric analysis of oscillations in the Frzilator model. Journal of Mathematical Analysis and Applications, 495(1), [124725]. https://doi.org/10.1016/j.jmaa.2020.124725

Copyright

Other than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).

Take-down policy

If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately and investigate your claim.

Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons the number of authors shown on this cover page is limited to 10 maximum.

(2)

Contents lists available atScienceDirect

Journal

of

Mathematical

Analysis

and

Applications

www.elsevier.com/locate/jmaa

Geometric

analysis

of

oscillations

in

the

Frzilator

model

Hadi Taghvafarda,∗,1,2, Hildeberto Jardón-Kojakhmetovb,3, Peter Szmolyanc,4, Ming Caoa,2

a

EngineeringandTechnologyInstitute,FacultyofScienceandEngineering,UniversityofGroningen, 9747AGGroningen,theNetherlands

b

FacultyofScienceandEngineering,BernoulliInstitute,UniversityofGroningen,9747AG,Groningen, theNetherlands

c_Institute_for_Analysis_and_Scientiﬁc_Computing,_Technische_Universität_Wien,₁₀₄₀_Wien,_Austria

a r t i cl e i n f o a b s t r a c t

Articlehistory:

Received2December2019 Availableonline27October2020 SubmittedbyS.A.Fulling Keywords:

Slow-fastsystem Relaxationoscillations Blow-upmethod

Geometricsingularperturbation theory

Myxobacteria

Abiochemical oscillator model, describingdevelopmental stage of myxobacteria, isanalyzedmathematically.Observationsfromnumericalsimulationsshowthatin a certain range of parameters, the correspondingsystem of ordinary diﬀerential equations displays stable and robust oscillations. In this work, we use geomet-ric singularperturbation theoryandblow-upmethod to provetheexistenceof a stronglyattractinglimitcycle.Thiscyclecorrespondstoarelaxationoscillationof anauxiliarysystem,whosesingularperturbationnatureoriginatesfromthesmall Michaelis-Mentenconstantsofthebiochemicalmodel.Inaddition,wegiveadetailed descriptionofthestructureofthelimitcycle,andthetimescalesalongit.

1. Introduction

Oscillators are ubiquitous in diﬀerent ﬁelds of science such as biology [33], biochemistry [8,11], neuro-science [15], medicine [12,30,29],and engineering [31].Inparticular, biochemicaloscillationsoftenoccur in

several contexts including signaling, metabolism, development, and regulation of important physiological

cellfunctions [25]. Inthispaper, westudyabiochemicaloscillatormodelthatdescribesthedevelopmental

stageofmyxcobacteria.Myxcobacteriaismulticellular organismsthatarecommoninthetopsoil [3].During

vegetationgrowth,i.e. when foodis ample, myxobacteriaconstitute small swarms byamechanism called

* Correspondingauthor.

E-mailaddress:taghvafard@gmail.com(H. Taghvafard).

1

Current address: Division of SystemsBiomedicine and Pharmacology, LeidenAcademic Center forDrug Research, Leiden University,Leiden,theNetherlands.

2 _The_work_of_H._Taghvafard_and_M._Cao_was_supported_in_part_by_the_European_Research_Council_{(ERC-CoG-771687)}_and_the

NetherlandsOrganisationforScientiﬁcResearch (NWO-vidi-14134).

3

H.Jardón-KojakhmetovwassupportedbytheAlexandervonHumboldtFoundation.

4

TheresearchofP.SzmolyanwasfundedbyWWTFunderthegrantMA14-049.

https://doi.org/10.1016/j.jmaa.2020.124725

(3)

“gliding” [16].Incontrast,underastarvationcondition,theyaggregateandinitiateacomplex

developmen-tal cycle during which small swarms are transformed into a multicellular single body known as “fruiting

body”,whose roleisto producesporesfornextgeneration ofbacteria [3].Duringtheaforementioned

tran-sition, myxobacteria passthroughadevelopmentalstage called the“ripplephase” [16,3],characterized by

complexpatterns ofwavesthatpropagatewithinthewholecolony.

Two geneticallydistinct molecular motors are concentratedat the cell poles ofmyxobacteria, 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 cellsforward, 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 oppositecellpoles [16]. As aresult, byforward and

backward motionofmyxobacteria, complexspatial wavepatterns arecreated.Inparticular,wavepatterns

are producedbythecoordination ofmotionof individual cellsthroughadirectend-to-endcontactsignal,

the“C-signal”.Duringtheripplephaseofdevelopment,theC-signalinginduces reversals,whilesuppresses

them duringtheaggregationstage ofdevelopment.Observationsfrom experimentsresultedinproposinga

biochemicaloscillatorin [16],knownastheFrzilator,whichactsas a“clock”to controlreversals.

TheFrzilatorisdetailedinSection2.1.Fromournumericalsimulations,itappearsthatthisbiochemical

oscillatorisrobustundersmallvariationofparameters.Moreimportantly,itseemsthat(almost)allsolutions

converge to a“unique”limitcycle. Regardingthepreviousproperty,in [28] ithasbeen shownthatwithin

acertainrangeofparametervalues,(almost)alltrajectoriesareoscillatory,thesystemhasaﬁnite number

ofisolatedperiodicorbits,atleastoneofwhichisasymptoticallystable.Althoughsomebiologicalsystems

mayproduce morethanonestable periodicsolutionfor acertainrange ofparameters [4], thecoexistence

betweenmultiplestable solutionshasnotyetbeenobservedexperimentally [9].

The main contributionofthis paper isto provethat,withinacertain rangeofparameter values,there

existsastronglyattracting periodicorbitfortheFrzilator.Moreover,thedetaileddescriptionofthestructure

of such periodic orbit is given. The methodology used to prove the aforementioned result consists ﬁrst

on an appropriate rescaling of the original model, which leadsto aslow-fast (or two timescales) system;

next, we take advantageof the two timescales of the rescaled systemto develop ageometric analysis via

techniques ofmulti-timescaledynamicalsystems.Fromthemulti-timescalenatureoftheproblem,itturns

outthatthelimitcycleisinfactarelaxationoscillator,meaningthatthereareseveraltimescalesalongthe

orbit oftheoscillator. From ananalyticalpoint ofview,themain diﬃcultyofthis analysis isthedetailed

description of atransition along two non-hyperboliclines (seedetails inSection3). Ouranalysis is based

ontheapproachdevelopedin [18,19] wheresimilar mechanisms,leadingtoanattracting limitcycle inthe

Goldbeter minimal model [8], have been studied. The C-signal of the biological oscillator plays a crucial

role both in ouranalysis and inthe behavior of the limit cycle. Thus, besides proving the existenceof a

strongly attractinglimitcycleusing geometricsingularperturbationtheory andblow-upmethod,wehave

performedatwo-parameterbifurcation analysisto showwhichcombination oftheC-signal and

Michaelis-Mentenconstantsofthesystemleadstooscillatorymotions.Inaddition,wehavecomputedacertainrange

forγ overwhichourmainresultisvalid.

The rest of this paper is organized as follows. In Section 2 we introduce the model, perform some

preliminaryanalysis onthemodel,andbrieﬂy introducethetoolswhich wearegoingto useinthepaper.

In Section3 we give the slow-fast analysis of an auxiliary system, corresponding to the original system.

Moreprecisely, wediscussthebehavior ofthedynamicswhenε→ 0.InSection4,we presenttheblow-up

analysis ofthenon-hyperbolicparts.Weconcludethepaper withadiscussionand outlookinSection5.

2. Detailedmodelandpreliminaryanalysis

Inthissectionweprovide apreliminaryanalysisofthebiochemicaloscillator proposedin[16].Westart

(4)

Fig. 1. Essential components of the Frzilator.

thebehaviorofthetrajectoriesandtheroleofparameters,andproposeauniﬁcationofthem.Afterwards,in

Subsection2.2,wepresentatwo-parameterbifurcationanalysiswhereweclarifythenatureandtheroleof

twodistinctparametersofthesystem.Finally,inSubsection2.3weprovideabriefintroductiontoslow-fast

systemsandthemaintechniques fortheiranalysis.

2.1. Modeldescription

Westudyabiochemicaloscillatormodelwhichdescribesthesocial-behaviortransitionphaseof

myxobac-teria [16]. This model, which is known as the Frzilator (or simply “Frz”) model, is based on a negative

feedbackloop.IntheFrzmodel,therearethreeproteins,namely,amethyltransferase(FrzF),the

cytoplas-mic methyl-accepting protein (FrzCD), and aprotein kinase (FrzE). A direct and end-to-end collision of

two myxobacteriaresultsinproducingasignal,so-called “C-signal”,underwhich aproteincalledFruA is

phosphorylated.Thesignalfrom phosphorylatedFruA(FruA-P)activates theFrzproteinsas follows[16]:

(i) the methyltransferaseFrzF (FrzF∗) is activatedby 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. Fig. 1 shows a

schematicrepresentation of interactionsbetween proteinsof theFrz system.Foramoredetailed

explana-tionof the modelandits biologicalbackground, see [16].Denote f,c and e respectively as thefraction of

activatedFrzF,methylatedFrzCD,andphosphorylatedFrzE. Thesefractions aregivenby[16]

f = [FrzF ∗_] [FrzF∗] + [FrzF], c = [FrzCD-M] [FrzCD] + [FrzCD-M], e = [FrzE-P] [FrzE] + [FrzE-P].

Theinteraction between theFrz proteinsis modeledbyMichaelis-Menten kineticsand henceleadsto the

dynamicalsystem df dτ = ka(1− f) − kdf e, dc dτ = km(1− c)f − kdmc, de dτ = kp(1− e)c − kdpe, (1) where

(5)

ka= kmaxa Ka+ (1− f) , kd= k_dmax Kd+ f , km= kmax m Km+ (1− c) , kdm = kmax dm Kdm+ c , kp= kmax p Kp+ (1− e) , kdp = k_dpmax Kdp+ e . (2)

Remark1.Dueto thefactthatf,c ande representfractions ofactiveproteinconcentrations,theirvalues

are restrictedto [0,1].Sothe fraction ofinactive proteinconcentrationsare given by(1− f), (1− c) and

(1− e).Therefore,hereafter, ouranalysis isrestrictedtotheunitcube

Q =(f, c, e)∈ R3| f ∈ [0, 1], c ∈ [0, 1], e ∈ [0, 1]. (3)

As mentionedin[16],theFrzsystemhasthewell-knownproperty of“zero-orderultrasensitivity”which

requires that the Michaelis-Menten constants Ka,Kd,Km,Kdm,Kp and Kdp have to be small [10]. It is

observed numerically in [16] that forthe parametervaluesKa = 10−2,Kd = Km = Kdm = Kp = Kdp =

5× 10−3, kmax

d = 1 min−1, kmmax = kpmax = 4 min−1, kmaxdm = kdpmax = 2 min−1, and kmaxa = 0.08 min−1,

system(1) hasanattractingperiodicsolution.Owingtothefactthattheunitofkmax

d ,kmaxm ,kpmax,kdmmaxand

kmax

dp ismin−1,tomakethemodelnon-dimensional,wedivideallequationsofsystem (1) bykdmax.Notethat

this doesnotchangethequalitativebehavior ofthe system,whilemakes allitsparametersdimensionless.

For simplicity,weunify allthedimensionlessMichaelis-Mentenconstants byKa= 2Kd= 2Km= 2Kdm=

2Kp = 2Kdp = ε 1. After unifying all Michaelis-Menten constants by ε, denoting γ := kamax, and

substituting(2) in (1),weobtainthefollowingdynamicalsystem

df dτ = γ(1− f) ε + (1− f)− 2f e ε + 2f, dc dτ = 8(1− c)f ε + 2(1− c)− 4c ε + 2c, (4) de dτ = 8(1− e)c ε + 2(1− e)− 4e ε + 2e.

Figs.2and3shownumericallycomputed attractinglimitcycleas wellas timeevolution ofsystem (4) for

ε= 10−3 andγ = 0.08.

Remark2.Forouranalysis inthis paper,we ﬁxγ = 0.08,while laterweshow thatthis parametercanbe

relaxed tosomeextent,seeRemark3andAppendixA.

Thedynamics alongthelimitcycle,showninFig.2,canbesummarized asfollows. Initially,allprotein

ratiosf,c ande areclosetozero,underthedynamics(4),thevariablef increases(duetotheactionofthe

C-signal), while c and e stay closeto zero. Once thevariable f passesthe activationthreshold f∗ := 0.5,

thevariable c increasesveryfast.Next,oncethevariable c passesthethresholdc∗:= 0.5,thevariable e is

activated andalso increasesveryfast until itreaches itsmaximumvalue, i.e., e= 1.Due to thefactthat

there isanegative feedbackfrom e tof , theincreaseine resultsinthedegradationof variablef .Once f

reaches thethresholdf∗,variablec decreases,andoncec reachesthethresholdc∗,thevariablee decreases

varyfast.Asaresult,thevariablesf andc reachtheirlowestvalues(i.e.veryclosetozero),butthevariable

e reachesthethresholde∗:= γ.Oncethevariablee dropsbelowthethresholde∗,thevariablef isactivated

and increases.Thisbehaviorisrepeated inaperiodicmannerandalimitcycleisformed(seeFig.2).

For system (4), a parameter-robustness analysis with respect to ε and γ = 0.08 is presented in [28].

(6)

Fig. 2. Numericallycomputed attractinglimitcycle ofsystem(1) forε= 10−3 andγ = 0.08.Thearrowsindicatethedirection andspeedoftheﬂowalongthelimitcycle.Asinglearrowcorrespondstotheslowesttimescale,whilethreearrowsindicatethe fastestone.

Fig. 3. Numerically computed periodic solution of system (1) for ε = 10−3 and γ = 0.08.

ε∈ (0,ε∗) withε∗:= 0.05517665.Moreover,itisproventhatforε∈ (0,ε∗),almostalltrajectoriesconverge to aﬁnite numberof periodicsolutions, oneofwhichis orbitally asymptoticallystable.Inthis article,we

provethe existenceof astrongly attracting limitcycle whichexplains the numericallycomputed periodic

orbit,forsuﬃcientlysmallε> 0. 2.2. Two-parameter bifurcationanalysis

Thissectionis devotedto thetwo-parameterbifurcationanalysis of (4).Inparticular,we areinterested

in understandingthe behavior of system (4) under the variation of parameters (ε,γ). To this end, letus

represent (4) by

˙x = G(x; ε, γ), (5)

where x = [ f c e ], and G(x;ε,γ) denotes the right-hand side of (4). We have used the numerical

continuationsoftware Matcont [5] tocomputethetwo-parameterbifurcationdiagramof (5) withrespect

to(ε,γ),presented inFig.4,where theverticalandthehorizontalaxesshow, respectively,thebehavior of

G(x;ε,γ) with respect to ε and γ.The bluecurve indicatesthatfor any 0< γ < 1 and any ε below the

curve, the system has unstable equilibria. Owing to the fact thatsystem (4) is cyclic [28, Theorem 5.7],

according to [28, Theorems 5.5 & 5.7] almost all trajectories of system (4) converge to a limit cycle for

any 0< γ < 1 andany ε below the bifurcation curve,depicted inFig.4. Forthose valuesof ε whichare

above the blue curve, the system is not oscillatory anymore, i.e. the equilibriumpoint is stable. In fact,

thebluecurveisacurveofHopfbifurcationswheretheequilibriaofthesystemswitchesfrombeingstable

to unstable: with ﬁxed 0 < γ < 1, as ε passes through the curve from above to below, a limit cycle is

generated.

As shown in Fig. 4, there are two points, denoted by “GH”, which are generalized Hopf (or Bautin)

(7)

Fig. 4. Two-parameterbifurcationanalysisof (4) withrespectto(ε,γ).(Forinterpretationofthecolorsintheﬁgure(s),thereader isreferredtothewebversionofthisarticle.)

theﬁrstLyapunovexponentcoeﬃcientoftheHopfbifurcation vanishes [23]. ComputedbyMatcont,the

valuesof(ε,γ) at“GH” pointsareasfollows:

(ε1, γ1) = (0.060907128, 0.086423772), (ε2, γ2) = (0.043172692, 0.949470320). (6)

InFig.4,theredcurvesarethecurvesof“limitpoints”(orsaddle-nodebifurcation)ofcycles.Forparameter

values (ε,γ) between theblueand redcurves inFig.4, at leasttwo limit cyclesexist simultaneously, i.e.,

for γ closeto0or γ closeto1,with asuitable0< ε 1, atleast onestableand oneunstable limitcycle 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,forsuﬃciently smallε,system (4) hassimilarqualitativebehaviors whenγ belongsto certain

boundswhicharecloseto0and1.Inthisregard,weemphasizethatalthoughthepositionofthelimitcycle

changeswhen γ iscloseto 1(see,forinstance,Fig.5),thegeometric analysisof thedynamicsisthesame

as thecasethatγ iscloseto0,forsuﬃcientlysmallε.

Remark4.InSection2.1, wehave uniﬁedallthe Michaelis-Mentenconstants of system (1) byε, resulted

in system (4). Although γ has a similar size as the Michaelis-Menten constants, we have not uniﬁed it

with them.Thereasonisthatγ is theC-signalof thebiologicaloscillator, i.e.,the inputunderwhichthe

dynamics are triggered (see Fig.1). Since suchasignalcoordinates themovementsof the individual cells

and inﬂuencestheshapeofthefruitingbodies [16],itiscrucialtohaveγ inacertainrangebetween0and

1.Ouranalysisthroughout thepaperclearlyshowstheroleofγ.

2.3. Preliminariesonslow-fast systems

Ourgoalistounderstandthedynamicsof(4) forsmallε inthelimitε→ 0.However,asitisseenin(4),

when the variables f,c and e arevery closeto the boundaryof Q, thelimiting behavior is diﬀerentfrom

the casethattheyare awayfrom theboundary.To resolvethe aforementionedproblem,one possibility is

to consideranauxiliarysystemwhichissmoothlyequivalent to(4).Tothis end,letusdeﬁne

Hε(f, c, e) := H₁ε(f )H₂ε(c)H₃ε(e),

(8)

Fig. 5. Numericallycomputedattractinglimitcycleofsystem (1) forε= 10−3andγ = 0.9.Thearrowsindicatethedirectionand speedoftheﬂowalongthelimitcycle.Asinglearrowcorrespondstotheslowesttimescale,whilethreearrowsindicatethefastest one.

H₁ε(f ) := (ε + 1− f)(ε + 2f),

H₂ε(c) := (ε + 2− 2c)(ε + 2c),

H₃ε(e) := (ε + 2− 2e)(ε + 2e).

(7)

Note that Hε_(f,_c,_e) _{> 0 for} _any _ε _{> 0 and} _any _(f,_c,_e) _{∈ Q.} _Therefore, _we _can _{reparametrize} _time _of

system(4) bymultiplyingbothsidesof(4) inHε(f,c,e),whichleadsto thefollowing dynamicalsystem

df dτ = γ(1− f) ε + (1− f)− 2f e ε + 2f Hε(f, c, e), dc dτ = 8(1− c)f ε + 2(1− c)− 4c ε + 2c Hε(f, c, e), (8) de dτ = 8(1− e)c ε + 2(1− e)− 4e ε + 2e Hε(f, c, e),

where,forsimplicity,werecycleτ todenote thereparametrizedtime.Onecanrewrite(8) asfollows

Xε: ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ df dτ = [γ(1− f)(ε + 2f) − 2fe(ε + 1 − f)] H2ε(c)H3ε(e), dc dτ = [8(1− c)f(ε + 2c) − 4c(ε + 2 − 2c)] H ε 1(f )H3ε(e), de

dτ = [8(1− e)c(ε + 2e) − 4e(ε + 2 − 2e)] H

ε

1(f )H2ε(c).

(9)

Thevector ﬁeld(9) is smoothly equivalent to (4) for ε> 0 [1], which from now onis theobject ofstudy. Themainreasontorewritesystem(4) intotheformofsystem(9) isthatthelatterisasingularlyperturbed

ODE whichallows us toanalyze thesystem usinggeometric methods. Moreover,note thatin contrastto

(4),system(9) ispolynomial,whichisanotherofitsadvantages.

2.3.1. Slow-fastsystems

System (9) isaslow-fastsysteminnon-standardform.Thus,inthissection,wepresentabriefsummary

of thebasicdeﬁnitions andresults regardingslow-fastsystems ofsuch aform.Our aimis to provide just

necessary terminologiesneededforthis paper.Foradetailedexposition,thereaderisreferredto [32].

Aslow-fastsystem(SFS)in standardform isasingularlyperturbedordinarydiﬀerentialequationwith

twotimescalesoftenpresented as

˙x = F (x, y, ε),

(9)

where the dot ˙ denotes derivativewith respect to the slow time t,F and G are assumed to be smooth,

x∈ Rns_, _y ∈ Rnf_,_and ₀_{< ε} 1 is_a_small _parameter_that_describes_the_timescale _separation_between_x and y.Fordetailssee[21].

Incontrast,aslow-fastsystemin non-standardform isavectorﬁeldas

z= Z(z, ε), (11)

where z∈ Rn_,_n_{≥ 2,}_the_prime _denotes_the_derivative_with_respect _to_the_fast_{time-parameter}_{τ ,}_and_the

smoothvectorﬁeldZ : Rn× R→ Rn isassumedtosatisfyZ(z,ε)= Z0(z)+ εW (z,ε).Oneusuallyassumes

thatε∈ (0,ε0),with ε 1.Themain diﬀerencebetween (10) and (11) is thatinsystem (11) thereisno

evidenttimescaleseparationbetweenthecomponentsofz.As(11) isaperturbationproblem,oneusually

startsitsanalysis byconsideringitsunperturbedversion.

Deﬁnition 1.Thelimitε→ 0 of(11),thatis

z = Z(z, 0) = Z0(z), (12)

is calledthelayerequation.

The set of singularities of Z0 has a crucial role in the analysis of (11), and allows one to distinguish

whether(11) deﬁnesaregularor asingularperturbationproblem:

Deﬁnition 2.Let

C = {x ∈ Rn_{| Z(z, 0) = Z}

0(z) = 0∈ Rn} . (13)

Then:

• System (11) is a regular perturbation problem if the set C is either empty or it consists entirely of

isolatedsingularities.

• System (11) is a singular perturbation problem if there exists a subset C0 ⊆ C which forms a

k-dimensional,1≤ k < n, diﬀerentiablemanifold.ThesetC0 iscalledthecritical manifold of(11).

Remark5.Ingeneral,andasithappensinthispaper,thesetC0isnotnecessarilyamanifold.Forexample,

itcouldbeformedbytheunion,alonglowerdimensionalsubmanifolds,ofdisjointk-dimensionalmanifolds.

Nevertheless, itiscustomary to keepreferringto C0 as thecritical manifold.Onecould referto itas “the

critical set”toavoidsuchambiguity.

We note that generically,given thatthe critical manifold C0 is k-dimensional, themap Z0 : Rn → Rn

hasconstantrank(n− k).This alsoimpliesthattheJacobianDZ0,evaluatedateachpointz∈ C0,hasat

least k zeroeigenvalues(calledtrivialeigenvalues)thatcorrespondto thetangentspaceTzC0,and(n− k)

nontrivial eigenvalues.

Deﬁnition 3.LetC0,n⊆ C0 denotethesubsetwhere allnontrivialeigenvaluesof DZ0|C0,n arenonzero, and

letC0,h⊆ C0,n denotethesubsetwhereallnontrivialeigenvaluesofDZ0|C0,h havenonzerorealpart.Inthe latter case,C0,his callednormally hyperbolic. Onthe otherhand,pointsinC0\C0,n or inC0,n\C0,h,ifthey

exist, are called non-hyperbolic, and one says that“the critical manifold C0 loses normal hyperbolicity at

(10)

Ifin(11) werescale timebydeﬁningtheslowtimeparametert= ετ ,onecanrewrite(11) as ˙z = 1

εZ0(z) + W (z, ε). (14)

The limitof (14) as ε→ 0, is called the reduced problem and is well deﬁnedif z ∈ C0 and the projection

ofW (z,ε) toC0 iswell deﬁned(see[32,section 3.2]).TheoverallideaofGeometricSingularPerturbation

Theory (GSPT) is to analyze the layer equation and the reduced problem, and then use perturbation

argumentstodescribe thedynamics of(11).

OneofthemainresultsofGeometricSingularPerturbationTheory(GSPT)isconcernedwithslow-fast

systemswithinasmallneighborhoodofanormallyhyperboliccriticalmanifold.

Theorem 1(Fenichel[7]).Let S0⊆ C0,h be a compactand normally hyperbolic critical manifoldof a SFS

(11). Then,forε> 0 suﬃciently small,thefollowingholds:

• There exists a locally invariant manifold Sε which is diﬀeomorphic to S0 and lies within distance of

orderO(ε) fromS0.

• Theﬂow on Sεconvergestothereduced ﬂow onS0as ε→ 0.

• Sε hasthesame stabilitypropertiesasS0.

Inanutshell,Fenichel’stheorem impliesthatthedynamicsoftheslow-fastsystem (11) nearacompact

and normallyhyperbolic critical manifold isaregular perturbation ofthe reducedproblem.On theother

hand,thedynamicsnearnon-hyperbolicpointscanbequiteintricate.Inthenextsectionwebrieﬂypresent

apowerfultechniquetostudy thedynamics ofslow-fastsystemsnear aclassof non-hyperbolicpoints.

2.3.2. Theblow-up method

Theblow-upmethodwasintroduced todescribethedynamics ofSFSs nearnon-hyperbolicpoints, and

isthemain mathematicaltechniqueusedinforthcomingsectionofthisarticle.Herewejustprovideabrief

descriptionofthemethod,formoredetailstheinterestedreaderisreferredto [6,17,20,21].

Deﬁnition4.Considerageneralizedpolarcoordinatetransformation

Φ : Sn× I → Rn+1

Φ(¯z, ¯ε, ¯r)→ (¯rα¯z, ¯rγε) = (z, ε),¯ (15) wheren_i=1z¯2

i+¯ε2= 1 andr¯∈ I whereI isa(possiblyinﬁnite)intervalcontaining0∈ R.Thecorresponding

(quasi-homogeneous)5 _{blow-up is}_deﬁned_by_(¯_z,_ε,_¯_r)_¯ _{= Φ}−1_(z,_ε)._The_map_{Φ is}_called _blow_down.6

Forthepurposes ofthisarticle,itissuﬃcienttoletr¯∈ [0,ρ),withρ> 0.Themainideaoftheblow-up

methodistoconstructanew,butequivalent,vectorﬁeldtoZ (11),whichisdeﬁnedinahigherdimensional

manifold,butwhosesingularitiesaresimplercompared tothoseofZ.

Deﬁnition5.Theblownupvector ﬁeldZ is¯ inducedbytheblow-upmap asZ = DΦ¯ −1◦ Z ◦ Φ, whereDΦ

denotesthederivativeofΦ.IfZ vanishes¯ onSn× {0} withorderm∈ N,wedeﬁnethedesingularizedvector

ﬁeld ˜Z = _r_¯1mZ.¯

5

Ahomogeneousblow-up(orsimplyblow-up)referstoalltheexponentsα,β,γ setto1.

6

Notethattheblow-upmapstheorigin0∈ Rn+1tothesphereSn_{× {0} while}

theblowdown doestheopposite,hencethe names.

(11)

Note thatthevectorﬁelds Z andZ are˜ equivalentonSn_{× {¯r > 0}.}_Moreover,_if_the_weights _(α,_{γ) are}

well chosen,thesingularitiesofZ˜|¯r=0arepartiallyhyperbolicorevenhyperbolic,makingtheanalysisofZ˜

simpler thanthatofZ. Due to theequivalence betweenZ andZ,˜ oneobtainsallthelocal information of

Z around 0∈ Rn+1from theanalysisofZ around˜ Sn× {¯r ≥ 0}.

While doing computations, it is more convenient to study the vector ﬁeld Z in˜ charts. A chart is a

parametrization of a hemisphere of Sn_{× I and} _is _obtained _by _setting _one _of _the _coordinates_(¯_z,_ε)_¯ _{∈ S}n

to ±1 inthe deﬁnitionofΦ.For example,oneofthemostimportant chartsinthe blow-upmethodisthe

centralchartdeﬁnedbyKε¯={¯ε = 1}.Afterwestudythedynamicsintherelevantcharts,weconnectthe

ﬂow together via transition maps, allowing us acomplete descriptionof the ﬂow of Z near˜ Sn _{× {0}.} _In

turn,and asmentioned above,theflowof Z is˜ equivalent to theflowof Z forε> 0 sufficientlysmall. For

moredetailsseeSection4and[17,21].

Remark6.Itisalsopossibletoblow-uponlysome ofthevariablesinthesystem(11),andkeeptheothers

unchanged.Inthispaper, weblow-upanon-hyperboliclineofequilibria toacylinder, seeSection4.

3. Geometricsingularperturbationanalysis

Thegoalofthissectionisto givethedetailedanalysis oftheslow-faststructure oftheauxiliarysystem

(9).

3.1. Layerproblemandthecritical manifold

Settingε= 0 in(9) resultsinthelayer problem

df dτ = (γ− e) H 0_{(f, c, e),} dc dτ = 2 (2f− 1) H 0_{(f, c, e),} de dτ = 2(2c− 1)H 0_{(f, c, e),} (16) with H0(f, c, e) = 32f ce(1− f)(1 − c)(1 − e).

Apart form the isolated equilibriumpoint P := (0.5,0.5,γ), whichis inside thecube Q, the boundaryof

Q, which consists of six planes, is the equilibria set of the layer problem (16). We denote each plane of

equilibria byS0,i(i= 1,2,...,6) asfollows:

S0,1:= (f, c, e)∈ R3| f = 0, c ∈ [0, 1], e ∈ [0, 1], S0,2:= (f, c, e)∈ R3| f ∈ [0, 1], c = 0, e ∈ [0, 1], S0,3:= (f, c, e)∈ R3| f ∈ [0, 1], c ∈ [0, 1], e = 0, S0,4:= (f, c, e)∈ R3| f = 1, c ∈ [0, 1], e ∈ [0, 1], S0,5:= (f, c, e)∈ R3| f ∈ [0, 1], c = 1, e ∈ [0, 1], S0,6:= (f, c, e)∈ R3| f ∈ [0, 1], c ∈ [0, 1], e = 1. (17) ThereforeS0:= 6

i=1S0,iisthecriticalmanifold.Thestabilityofsystem(9) changesatlinesf ∈ S0,2,f ∈

(12)

Fig. 6. ThecriticalmanifoldS0=6i=1S0,i,non-hyperboliclinesf,c,e,f,c,einred,all12non-hyperbolicedgesinblue,and

inparticular,thetwonon-hyperbolicedges1and2shallplayanimportantroleinouranalysis.

Moreover,the12edgesoftheunitcube,wherethe6planesS0,iintersect,arenon-hyperboliclines aswell.

However,forour analysis,only thelines 1=S0,1∩ S0,2 and 2=S0,2∩ S0,3 arecrucial (see Fig.6). The

stabilityof pointsinS0 issummarizedinthefollowinglemma.

Lemma1.The criticalmanifold S0 ofthelayer problem(16) hasthefollowingproperties:

• S0,1isattracting fore> e∗ andrepellingfore< e∗.

• S0,2isattracting forf < f∗ andrepellingforf > f∗.

• S0,3isattracting forc< c∗ andrepellingforc> c∗.

• S0,4isattracting fore< e∗ andrepellingfore> e∗.

• S0,5isattracting forf > f∗ andrepellingforf < f∗.

• S0,6isattracting forc> c∗ andrepellingforc< c∗.

• Theequilibrium P := (0.5,0.5,γ) is asaddle-focuspoint.

• The linesf ∈ S0,1,c ∈ S0,2,e∈ S0,3,f ∈ S0,4,c ∈ S0,5,e∈ S0,6, all12 edgesof theunit cube,and

inparticular,theedges1=S0,1∩ S0,2 and2=S0,2∩ S0,3 arenon-hyperbolic.

Proof. Theeigenvaluesofthelinearizationofsystem(16) atpoints, e.g.,intheplaneS0,1aregivenby

λ1= λ2= 0, λ3=−32ce(c − 1)(e − 1)(e − γ).

Itisclearthatλ3iszeroattheboundaryofS0,1,andalsoalongthelinelegivenbye= e∗.Therefore,S0,1is

attractingfore> e∗anditisrepellingfore< e∗.Theproofoftheothercasesisperformedanalogously.

We denote the interior of the cube Q by Q.˚ Note that when (f,c,e)∈ ˚Q, the layer problem (16) can

be dividedby thepositive termH0_(f,_c,_e)_{= 32f ce(1}_{− f)(1}_{− c)(1}_{− e).} _Therefore_away_from _the _critical

manifold S, all the variables evolve on thefast time scale τ and the orbits of thelayer problem (16) are

identicalto theorbitsofthelinearsystem

df dτ = γ− e, dc dτ = 2(2f− 1), de dτ = 2(2c− 1). (18)

(13)

Table 1

Foreachrowi weshowtheintervalofdeﬁnitionofSa

ε,iandtherelationbywhichitisdeﬁned,

allanalogoustoLemma2.

i Ia i Sε,ia 2 (f, e)∈δ, 1 2− δ × [δ, 1 − δ] c =₁_−2ff ε + O(ε2) 3 (f, c)∈ [δ, 1 − δ] ×δ, 1 2− δ e = c 1−2cε + O(ε2) 6 (f, c)∈ [δ, 1 − δ] ×1 2+ δ, 1− δ e = 1 + 1 2(1−2c)ε + O(ε2)

Remark7.System(18) isthelimitof(4) whenε→ 0 and(f,c,e)∈ ˚Q. 3.2. Reducedproblem, slowmanifolds, andslowdynamics

From Subsection 3.1,we knowthattheboundaryof Q isthecriticalmanifold S0.Any compactsubset

ofS0thatdoesnotcontainanynon-hyperbolicpointisnormallyhyperbolic,andhenceFenicheltheory[7]

is applicable.Inother words,this theoryimpliesthatthenormally hyperbolicparts ofS0perturbto slow

manifolds,whichliewithinadistanceoforderO(ε) ofthecriticalmanifoldS0.Inthefollowing,wecompute

theslowmanifoldsand analyzethereducedﬂowsintheplanesS0,1,S0,2,S0,3 andS0,6whichareessential

forouranalysis.

Lemma 2.For suﬃciently small δ > 0, there exist ε0 > 0 and a smooth function hε,1(c,e) deﬁned on

Ia

1 = [δ, 1− δ] × [γ + δ, 1 − δ] suchthat themanifold

Sa

ε,1={(f, c, e) ∈ Q | f = hε,1(c, e), (c, e)∈ I1a} , (19)

is alocallyinvariantattractingmanifoldof(9) forε∈ (0,ε0].The functionhε,1(c,e) hastheexpansion

hε,1(c, e) =

γ

2(e− γ)ε + O(ε

2_). ₍₂₀₎

Proof. Since theset Ia

1 is hyperbolic,Fenichel theory impliesthat there exists asuﬃciently small ε0 > 0

such that the functionhε,1(c,e) hasthe expansion hε,1(c,e) = η(c,e)ε+ O(ε2) for all ε∈ (0,ε0]. Due to

invariance,wecansubstitutehε,1(c,e) intotheequationof df_dτ in(9) andidentifycoeﬃcientsofε.Bydoing

so,weobtain

η(c, e) = γ

2(e− γ). (21)

Note that(21) reﬂects thefactthatthemanifold S_ε,1a isnotwell-deﬁned when e= γ. Thus, theinvariant

manifold Sa

ε,1 isgivenasstatedinthelemma,whichcompletestheproof.

For the sake of brevity, we summarize theanalysis inthe planes S0,2, S0,3 and S0,6 in Table 1, which

is shown byfollowing the sameline of reasoning as the oneof Lemma2. For moredetails, theinterested

readeris referredto [27].

Remark 8.Similarresults canbe obtainedforthe “repelling” parts Sr

ε,i,i= 1,2,...,6.However,these are

not neededinouranalysis. Nonetheless,we point outthattheslow manifoldsS_ε,ir wouldbe expressed by

thesamefunctionshε,i and appropriateintervalsIir.

Remark9.Theexpansionsofthefunctionshε,i(·,·),i= 1,2,3,6,alsoexplainwhyitisnecessarytorestrict

thedomain ofdeﬁnitionof theslowmanifoldsto Ia

(14)

Fig. 7. Flow of the slow vector ﬁeld inS0,1, non-hyperbolic line ein red, and sections c1, c2in cyan.

WenowturntotheanalysisofthereducedﬂowsintheplanesS0,1,S0,2, S0,3andS0,6which,respectively,

means the planesf = 0,c = 0,e = 0 ande = 1. We know thatsystem (9) has the fast time scale τ . By

substituting the functionshε,i, i = 1,2,3,6 into (9), transforming the fast time variable to the slow one

byt = ετ , and settingε = 0, theequations governingthe slowdynamics on thecritical manifoldS0,iare

computed.Inthefollowing,wegivetheanalysis intheplaneS0,1.

After substituting hε,1 into system (9), the dynamics of thereduced system in S0,1, i.e., onthe plane

f = 0,isgovernedby c =−32ce 2_(c_{− 1)(e − 1)} e− γ ε + O(ε 2_), e =32ce 2_(c_{− 1)(e − 1)(2c − 1)} e− γ ε + O(ε 2_), (22)

where denotesthediﬀerentiationwithrespectto τ .Nowbydividingoutafactorofε, whichcorresponds

toswitching fromthefasttimevariableto theslowone,wehave

˙c = −32ce 2_(c_{− 1)(e − 1)} e− γ + O(ε), ˙e =32ce 2_(c_{− 1)(e − 1)(2c − 1)} e− γ + O(ε), (23)

where the overdot represents diﬀerentiation with respect to t = ετ . Now, by setting ε = 0 in (23), the

reducedﬂowonS0,1isgiven by

˙c = −32ce 2_(c_{− 1)(e − 1)} e− γ , ˙e = 32ce 2_(c_{− 1)(e − 1)(2c − 1)} e− γ . (24)

Asitisclear,thevectorﬁeld(24) issingularatthelinee,givenbye= e∗.Inotherwords,theﬂow(24) is

notdeﬁnedonthelinee.Thelinesc= 0,e= 0,c= 1,ande= 1,showninFig.7,arelines ofequilibria.

Theline c = 0 isattracting for e> e∗ and itis repelling fore< e∗, whilethe line c= 1 is attractingfor

e< e∗ andrepellingfore> e∗.

Bydividingoutthefactor 32ce2(c_e−γ−1)(e−1) in(24),theorbitsofthereducedﬂowcanbederivedfrom the

(15)

˙c =−1,

˙e = 2c− 1, (25)

whichcanbeintegrated explicitly.

Remark 10. For e > e∗, systems (24) and (25) have qualitatively the same dynamics when c,e ∈ (0,1).

In particular, the vector ﬁeld (25) is C∞-equivalent but not C∞-conjugate to the vector ﬁeld (24). For

the case that e < e∗, the direction of the vector ﬁeld (24) is not preserved in the vector ﬁeld (25).

However, for our analysis, it suﬃces to study the ﬂow of system (24) when e > e∗, or equivalently

onSa 0,1.

Fig. 8. The reducedﬂowsalongSa

0,2,S0,3a andS0,6a .Thevectorﬁeldsare singularatred lines.Thethickbluelinesarelines of

(16)

Lemma 3.For e > e∗, the reduced ﬂow (24) on S0,1 and hence the slow ﬂow (23) on Sε,1a maps section

{c = c1} to {c = c2},where0< c2< c1< 1₂;thismapiswell-deﬁnedanditsﬁrstderivative withrespectto

e isequaltoone.

Proof. Itsuﬃcestoconsider(25).LetΠ(e) denotethemap from{c = c1} to{c = c2} induced bytheﬂow

of(25).Then,itisstraightforwardtogetΠ(e)= e+ c2− c22− c1+ c21,fromwhichthestatementfollows.

In order to obtain the equations governing the slow ﬂow along Sa

ε,2,Sε,3a and Sε,6a , a similar analysis

can be done by inserting the functions hε,2,hε,3 and hε,6 into (9) and dividing out afactor of ε, which

correspondstoswitching totheslowtimescalet= ετ .Next,bysettingε= 0 oneobtainsthereducedﬂow

on the critical manifolds S0,2, S0,3 and S0,6. For the sakeof brevity, we have summarized the slow ﬂows

alongS0,2, S0,3andS0,6inFig.8.Formoredetails,theinterestedreaderisreferredto [27].

3.3. Singular cycle

Inthissection,wepresenttheoverallbehaviorofthesingularcycle,whichisaclosedcurveconsistingof

alternatingparts of thelayer problem,and thecritical manifoldS0. However,bythe informationthatwe

havesofarfromthecriticalmanifoldandthelayerproblem,wecannotfullydescribethesingularcycleclose

tothenon-hyperboliclines1 and2.Afulldescriptionofthesingularcycleforthosepartsthatcannotbe

derivedfromthecriticalmanifoldandthelayerproblem ispresentedinSection4bytheblow-upmethod.

The construction of the singular cycle Γ0 starts at the point pf := (0.5,0,0). This point is connected

to the point p1 := (1+

√_γ

2 ,0.5,0) ∈ c through the orbit ω1 of the reduced ﬂow (27). Starting at p1, the

layer problem (18) intersectsthe attractingpartofthe planeSa

0,6 inapoint, denotedbyp2.This pointis

connected to apoint, denotedby qe _∈c_, _through_the _orbit _ω

3 of the reduced ﬂow (28). Starting at qe,

throughthelayer problem(18),theorbitω4 intersects theplaneS0,1a atapoint,denotedbyqe. Theorbit

ω5 ofthereducedﬂow(22) connectsqetoapoint, denotedbype∈ 1,whichistheintersectionofS0,1a and

Sa

0,2;peisconnectedtothepointpe:= (0,0,γ) byasegmentontheline1,denotedbyω6.Theorbitω7of

thereducedﬂow(26)connectspeto thepointpf := (γ

2

4,0,0);Finally,pf isconnectedto pf byasegment

ontheline2,denotedbyω8.Hence,thesingularcycleΓ0∈ R3ofsystem(9) forε= 0 isdeﬁnedasfollows

(seeFig.9):

Γ0:= ω1∪ ω2∪ ω3∪ ω4∪ ω5∪ ω6∪ ω7∪ ω8. (29)

Remark11.Alltheorbitsωj (j = 1,2,...,8)areknownanalytically.

Owing to the fact that the layer problem is linear, all the points that connect ωj to ωj+1 are

ex-plicitly known. For the particular quantity γ = 0.08, we have pf _{= (0.5,}_0,_0), _p

1 ≈ (0.6414,0.5,0),

p2 ≈ (0.3638,0.8485,1), qe ≈ (0.0771,0.5,1),qe ≈ (0,0.3438,0.9743), pe ≈ (0,0,0.7487),pe= (0,0,0.08),

andpf = (0.0016,0,0).

Remark12.Atthesingularlevel,thereisnovisibleﬂowonthesegmentsω6andω8.Theblow-upanalysis,

carriedoutinSection4,will revealahiddenﬂowforsuchsegments.

3.4. Main result

InviewofthesingularcycleΓ0,introduced intheprevioussubsection,wearenowreadyto presentthe

(17)

Fig. 9. Schematic diagram of the singular cycle Γ0.

Theorem2.AssumethatΓ0isthesingularcycledescribedinSection3.3.Thenforsuﬃcientlysmallε> 0,

there exists a unique attracting periodic orbit Γε of the auxiliary system (9), which tends to the singular

cycle Γ0 asε→ 0.

Inorderto proveTheorem 2,weneedto introducethefollowing sections

Σ1:={(f, c, e) ∈ R3 | (f, e) ∈ R1, c = δ1},

Σ2:={(f, c, e) ∈ R3 | (c, e) ∈ R2, f = δ2},

Σ3:={(f, c, e) ∈ R3 | (f, e) ∈ R3, c = δ3},

(30)

where Rj (j = 1,2,3) are suitable small rectangles, and δj are chosen suﬃciently small. Note thatΣ1 is

transversal toω4,Σ2 istransversaltoω6,and Σ3 istransversalto ω8,seeFig.9.

AccordingtothedeﬁnitionofthesectionsΣi,introducedin(30),wedeﬁne thefollowingPoincarémaps

fortheﬂow ofthesystem(9)

π1: Σ1→ Σ2,

π2: Σ2→ Σ3,

π3: Σ3→ Σ1,

(31)

wherethemapπ1describesthepassagefromΣ1toΣ2alongthenon-hyperbolicline1,themapπ2describes

thepassagefrom Σ2 toΣ3alongthenon-hyperbolicline2,andthemapπ3 describesthepassagefromΣ3

to Σ1.Themap π3 consists ofslowﬂow alongSε,3a , followedbythefastdynamics from aneighborhood of

p1 toaneighborhood ofp2,followedbytheslowﬂowalongSε,6a toaneighborhoodofqe. Throughthefast

dynamics,thisneighborhoodismappedtoaneighborhoodofqe,followedbytheslowﬂowalongSε,1a toΣ1.

Wesummarizethepropertiesoftheabovemaps inthefollowinglemmas.

Lemma 4.If thesection Σ1 ischosensuﬃciently small,then thereexistsε0> 0 such that themap

(18)

iswell-deﬁnedforε∈ [0,ε0] andsmoothforε∈ (0,ε0].The mapπ1isastrongcontractionwithcontraction

rateexp(−K/ε) forsomeK > 0.TheimageofΣ1isatwo-dimensionaldomainofexponentiallysmallsize,

whichconvergestothepointq2:= Σ2∩ ω7 asε→ 0.

Lemma5.If thesectionΣ2 ischosensuﬃciently small,thenthere existsε0> 0 suchthat themap

π2: Σ2→ Σ3, (c, e)→ (πf2(c, e, ε), πe2(c, e, ε)), (33)

iswell-deﬁnedforε∈ [0,ε0] andsmoothforε∈ (0,ε0].The mapπ2isastrongcontractionwithcontraction

rateexp(−K/ε) forsomeK > 0.TheimageofΣ2isatwo-dimensionaldomainofexponentiallysmallsize,

whichconvergestothepointq3:= Σ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,

whichwillbe presentedinSubsections4.1 and4.2.

Remark13. Thepoints on the line c when 0.5 < f < 1,and on theline c when 0< f < 0.5 are jump

foldpoints,i.e., thetrajectoryswitchesfrom theslowdynamics tothefastdynamics.Oneway ofshowing

the aforementioned is by following [24, Lemma 6]. For the reader’s convenience, we have included such

computations inAppendixB.This alsoexplainswhythebehaviorof thetrajectorynearc isverysimilar

to thebehavior ofstandardslow-fastsystemswith twoslow variablesand onefast variablenear ageneric

“fold” line, studiedin [26] based on the blow-up method. The critical manifolds S0,3 and S0,6 of system

(9) can be viewed as a standard folded critical manifold, which has been straightened out by a suitable

diﬀeomorphism.Thisleadstothecurved ﬁbersofthelayerproblem(16).Therefore,wecanusetheresults

of[26] tounderstandthebehaviorof (9) closetothenon-hyperboliclines c andc.

Thefollowinglemmadescribesthemap fromthesectionΣ3 tothesectionΣ1,deﬁnedin(31).

Lemma6.If thesectionΣ3 ischosensuﬃciently small,thenthere existsε0> 0 suchthat themap

π3: Σ3→ Σ1, (f, e)→ (π3f(f, e, ε), π

e

3(f, e, ε)), (34)

is well-deﬁned for ε ∈ [0,ε0] and smooth for ε ∈ (0,ε0]. The image of Σ3 is an exponentially thin strip

lying exponentiallyclose toS1

a,ε∩ Σ1, i.e., itswidth in thef -direction is O(exp(−K/ε)) for some K > 0.

Moreover, π3(Σ3) convergestoasegment ofS0,1a ∩ Σ1 asε→ 0.

Proof. Thebasicideaoftheproofisbasedonthemapthathasbeen alreadydescribedinFig.9forε= 0,

denotedbyπ0

3,and thentreatπ3 asan ε-perturbationofπ03.If thesectionΣ3 ischosen suﬃcientlysmall,

then the trajectories starting in Σ3 can be described by the slow ﬂow along the manifold Sε,3a combined

with the exponential contractiontowards the slowmanifold until they reacha neighborhood of thejump

pointsonthelinec.Applying[26,Theorem1] closetothejumppints,thetrajectoriesswitchfromtheslow

dynamics to the fast dynamics, and hence pass the non-hyperbolic line c; this transition is well-deﬁned

for ε∈ [0,ε1], and smooth forε∈ (0,ε1] for someε1 > 0.Note that[26, Theorem 1] guarantees thatthe

contractionof the solutionsinthee-direction persistsduringthe passage throughthefold-line c, as itis

at most algebraically expanding. After that, the solutionsfollow the fastdynamics ω2 until theyreach a

neighborhood of thepoint p2, see Fig.9. Next,the solutionsfollow the slow ﬂow alongthe manifold Sε,6a

combined with theexponential contraction towardsthe slowmanifold until theyreach aneighborhood of

thepoint qe_._Again _applying _[₂₆_,_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-hyperbolicline c_, _where_the _{corresponding}_transitions _are _{well-deﬁned}_for_ε_{∈ [0,}_ε

(19)

ε∈ (0,ε2] for someε2> 0, andthenfollow thefastdynamics (ω4)until theyreachaneighborhood ofthe

point qe.Finally,thesolutionsfollowtheslowﬂow alongthemanifoldSε,1a combinedwith theexponential

contractiontowardstheslowmanifolduntil theyreachthesectionΣ1.

Theorem 1of[26] impliesthatthemapπ3 isatmostalgebraicallyexpandinginthedirectionofe when

Σ3 is chosen suﬃciently small. On the other hand,the slow manifold Sε,1a is exponentially contractingin

the directionof f (Fenicheltheory). Therefore,theimage ofΣ3 isathinstriplying exponentiallyclose to

Sa

ε,1∩ Σ1.Hence,thestatementsofthelemmafollow.

Now wearereadytogivetheproofofthemain result.

Proof of Theorem 2. Letus deﬁne the map π : Σ3 → Σ3 as a combination of the maps πj (j = 1,2,3),

described inLemmas 4,5and6.Moreprecisely,wedeﬁne

π = π2◦ π1◦ π3: Σ3→ Σ3.

IfthesectionΣ3ischosensuﬃcientlysmall,Lemma6impliesthatthereexistsε3> 0 suchthatthemapπ3

iswell-deﬁnedforε∈ [0,ε3] andsmoothforε∈ (0,ε3],andtheimageofΣ3isathinstriplyingexponentially

close toSa

ε,1∩ Σ1,i.e.,π3(Σ3) isexponentiallycontractingwithrateexp(−K3/ε),forsomeK3> 0,inthe

f -directionwhileitisbounded inthee-direction.

Next,iftheentrysectionΣ1 ischosensuchthatΣ1⊃ π3(Σ3),Lemma4impliesthatthereexists ε1> 0

suchthatthemapπ1iswell-deﬁned foranyε∈ [0,ε1] and smoothforε∈ (0,ε1],andπ1is anexponential

contraction with rate exp(−K1/ε) for someK1 > 0. Finally, if the entry section Σ2 is chosen such that

Σ2⊃ π1(Σ1),Lemma5impliesthatthereexistsε2> 0 suchthatthemapπ2iswell-deﬁnedforanyε∈ [0,ε2]

and smooth for any ε ∈ (0,ε2], and further, π2 is an exponential contraction with rate exp(−K2/ε), for

someK2> 0,suchthatΣ3⊃ π2(Σ2).

Denotingε0:= min{ε1,ε2,ε3} andK := min{K1,K2,K3},themap π : Σ3→ Σ3 iswell-deﬁnedfor any

ε∈ [0,ε2],andsmoothforε∈ (0,ε0].Further,basedonthecontractingpropertiesofthemapsπi,i= 1,2,3,

weconcludethatπ(Σ3)⊂ Σ3iscontractionwithrateexp(−K/ε).TheBanachﬁxed-pointtheoremimplies

the existence ofa uniqueﬁxed point for themap π, corresponding to theattracting periodicorbit ofthe

system (9). Moreover, due to the last assertion of Lemmas 4, 5 and 6, the periodicorbit Γε tends to the

singularcycleΓ0as ε→ 0.Thiscompletestheproof.

4. Blow-upanalysis

The slow-fast analysis that we have done in Section 3 does not explain the dynamics of system (9)

close to the non-hyperboliclines 1 and 2. As thesegments ω5 and ω7 lieon these lines (see Fig. 9), we

need adetailed analysis close to the lines 1 and 2, which is carried outin this section via the blow-up

method [21,14,20].Toapplythis,weextendsystem(9) byaddingε asatrivialdynamicvariableandobtain

df dτ = [γ(1− f)(ε + 2f) − 2fe(ε + 1 − f)] H ε 2(c)H3ε(e), dc dτ = [8(1− c)f(ε + 2c) − 4c(ε + 2 − 2c)] H ε 1(f )H3ε(e), de

dτ = [8(1− e)c(ε + 2e) − 4e(ε + 2 − 2e)] H

ε 1(f )H2ε(c), dε dτ = 0, (35) whereHε

1(f ),H2ε(c) andH3ε(e) aredeﬁnedin(7).Notethatfortheextendedsystem(35),thelines1× {0}

(20)

quadruplezeroeigenvalues,system(35) isverydegeneratecloseto1× {0} and2× {0}.Toresolvethese

degeneracies,weusetheblow-upmethod,giveninnextsubsections.

4.1. Blow-up ofthenon-hyperbolic line1× {0}

Theblow-upofthenon-hyperbolicline1×{0} ispresentedinthissubsection.Tothisend,wetransform

thenon-hyperbolicline ofsteadystates1× {0} by

f = r ¯f , c = r¯c, ε = r ¯ε, e = ¯e, (36)

where f¯2_{+ ¯}_c2_{+ ¯}_ε2 _{= 1 and}_r_{≥ 0.}_Note _that_since_(f,_c,_e)_{∈ Q,} _we_may_further_assume_that_{f ,}¯_¯_c_{≥ 0 and} ¯

e∈ [0,1]. Sinceall weights are equalto 1in (36), this is ahomogeneous blow-up.For ﬁxed e,¯ eachpoint (0,0,e) is¯ blown-up toasphereS2, andtheline1× {0} isblown-upto acylinderS2× [0,1],seeFig.10.

Fortheanalysis ofsystem(35) near theline 1× {0},we deﬁnethree chartsK1,K2 andK3 bysetting

¯

c = 1,ε = 1,¯ andf = 1 in¯ (36),respectively:

K1: f = r1f1, c = r1, ε = r1ε1, e = e1, (37)

K2: f = r2f2, c = r2c2, ε = r2, e = e2, (38)

K3: f = r3, c = r3c3, ε = r3ε3, e = e3. (39)

ThechangesofcoordinatesforthechartsK1to K2, andK2 toK3 intheblown-upspace aregiveninthe

followinglemma.

Lemma7.The changesof coordinatesK1 to K2,andK2 toK3 aregivenby

κ12: f2= f1 ε1 , c2= 1 ε1 , ε2= r1ε1, e2= e1, ε1> 0, (40) κ23: r3= r2f2, c3= c2 f2 , ε3= 1 f2 , e3= e2, f2> 0. (41)

Thegoalof thissubsectionis toconstruct thetransition mapπ1: Σ1 → Σ2,deﬁnedin(31),and prove

Lemma4.Beforegoingintothedetails,letusbrieﬂydescribeourapproach.Wedescribethetransitionmap

π1: Σ1→ Σ2viaanequivalentoneintheblown-upspace.Morespeciﬁcallywedeﬁne

π1:= Φ◦ ¯π1◦ Φ−1, (42)

where

¯

π1:= Π3◦ κ23◦ Π2◦ κ12◦ Π1,

andΦ: S2_×[0,_1]_×[0,_r

0)→ R4isthecylindricalblow-updeﬁnedby(36),themapsΠiarelocaltransitions

induced by the blown-up vector ﬁelds which are detailed below, and κ12 and κ23 denote the changes of

coordinates,giveninLemma7.π¯1 isthetransitionmap intheblown-up spaceanddue tothefactthatΦ

isadiﬀeomorphism,itisequivalent toπ1.A schematicoftheproblem athandisshowninFig.10.

TheleftpictureinFig.10illustratesthecritically manifoldsSa

0,1and S0,2a ,and thecorrespondingﬂows

in blue. The non-hyperbolic line 1 is shown in orange. For e > γ, the reduced ﬂows on both critically

manifolds approach the line 1. Atthe point on the line 1 with e = γ, a transition from S0,1a to S0,2a is

possible as indicated inthe ﬁgure. The right picture in Fig.10 schematicallyshows the conﬁguration in

(21)

Fig. 10. Theleftﬁgureshowsthedynamicsclosetothenon-hyperbolicline1.Therightﬁgureshowsthecorrespondingdynamics

intheblown-upspace.

corresponding toε = 0 and¯ r > 0 areshownoutsideofthe cylinder.Herewerecoverthelayerproblem,the

critically manifolds,and thereducedﬂowsinS¯a

0,1andS¯0,2a .Intheblown-up space,themanifoldsS0,1a and

Sa

0,2are separatedand hencegainedhyperbolicity,inparticular theyare attractive,as indicated belowin

Fig.11a.Allthese assertionswill beproveninthis section.

Roughly speaking, inchart K1 we continue theattracting slow manifold S¯0,1a onto thecylinder. Chart

K2 isusedtotracktheﬂowacrossthecylinder.Theexitoftheﬂowfromthecylinderanditstransition to

¯

Sa

0,2isstudiedinchartK3,see Figs.10and15.ThedetailedanalysisofthemapsΠi introducedin(42),is

given intheforthcomingsubsections.

4.1.1. Analysis inchartK1

Aftersubstituting(37) into(35) anddividingoutalltheequationsbythecommonfactorr1,theequations

governing thedynamics inchartK1aregiven by

f₁ =−4f1Γ1G11+ [γ(1− r1f1)(ε1+ 2f1)− 2f1e1(r1ε1+ 1− r1f1)]G12, r₁= 4r1Γ1G11, e₁= 4r1[2r1(1− e1)(r1ε1+ 2e1)− e1(r1ε1+ 2− 2e1)]G13, ε₁=−4ε1[2r1f1(1− r1)(ε1+ 2)− (r1ε1+ 2− 2r1)]G11, (43) where wedenote Γ1:= [2r1f1(1− r1)(ε1+ 2)− (r1ε1+ 2− 2r1)], G11:= (r1ε1+ 1− r1f1)(r1ε1+ 2− 2e1)(ε + 2f1)(r1ε1+ 2e1), G12:= (r1ε1+ 2− 2r1)(r1ε1+ 2− 2e1)(ε1+ 2)(r1ε1+ 2e1), G13:= (r1ε1+ 1− r1f1)(r1ε1+ 2− 2r1)(ε1+ 2f1)(ε1+ 2). (44)

From(43) itisclearthattheplanesr1= 0 andε1= 0 areinvariant.Hence,weconsiderthefollowingcases:

1. r1= ε1= 0:inthis case,thedynamics (43) issimpliﬁedto

e₁= 0,

f1 = 32f1e1(1− e1)[2f1+ γ− e1].

(22)

Fig. 11. Dynamics of (43) restricted to invariant subspaces.

For ﬁxede, theequilibriaof thesystem(45) are theattractingpointpa

1 = (f1,r1,e1,ε1)= (0,0,e1,0),

andtherepellingpointpr

1= (f1,r1,e1,ε1)= (e1₂−γ,0,e1,0).Notethatthetwohyperbolicpointspa1and

pr

1 intersectat thenon-hyperbolicpoint(f1,r1,e1,ε1)= (0,0,γ,0),seeFig.11a.

2. ε1= 0:inthis case,thedynamics (43) isrepresentedby

f₁ = 32f1e1(1− e1)(1− r1)(1− r1f1)[(γ− e1)− 2f1(2r1f1− 1)],

r₁ = 64r1f1e1(1− e1)(1− r1)(1− r1f1)[2r1f1− 1],

e₁= 64r1f1e1(1− e1)(1− r1)(1− r1f1)[2r1− 1].

(46)

From (46),oneconcludesthattheplanef1= 0 is theplaneofequilibria whichisdenotedbyS¯0,1a ,see

Fig.11a.Thenon-zeroeigenvaluealongS¯a

0,1isgivenbyλ= 32e1(1− e1)(1− r1)(γ− e1).For0≤ r1< 1

and e1 > γ, the plane S¯0,1a is attracting. As the e1-axis is apart of S¯0,1a , we denote that part of the

e1-axisthatγ≤ e1≤ 1 bye1.Wealsohaveanothercurveofequilibriawhichisdeﬁnedbyr1= 0,and

f1=e1₂−γ,denotedbyMr1,see Fig.11a.Thiscurveofequilibriaisofsaddle-typewiththeeigenvalues

λ=±32e1(e1− 1)(e1− γ).Note thatwehaverecoveredtheinformationofthepreviouscasehere.

3. r1= 0:inthiscase,thedynamics(43) isrepresentedby

e₁= 0,

f₁ = 8e1(1− e1)[(ε1+ 2)(γ(ε1+ 2f1)− 2f1e1) + 4f1(ε1+ 2f1)],

ε₁= 32e1ε1(1− e1)(ε1+ 2f1).

(47)

Bysettingε1= 0,weagainhavethelinee1 andthecurveM

r

1.TheJacobianmatrixatapoint ine1 hastwoeigenvalues:onezeroandtheotheroneisλ= 32e1(1− e1)(γ− e1).Sothelinee1 isattracting whene> γ.Asinthiscasewehavetwozeroeigenvalues,itimpliesthatthereexists atwo-dimensional centermanifold,namely,Ca,1.

Remark14.InchartK1,themostimportantroleisplayedbythetwo-dimensionalcentermanifoldCa,1,

see Lemma9.Infact,this isthecontinuationofthecriticalmanifoldS¯_0,1a .

Wesummarizetheanalysisperformedinthis subsectioninthefollowinglemmas.

Lemma8.System (43) hasthefollowingmanifoldsofequilibria:

1. Theplane S¯a

(23)

2. Mr

1={(f1,r1,e1,ε1)| f1= e1−γ₂ , r1= 0, e1∈ [γ,1], ε1= 0}.

Lemma 9.The followingpropertiesholdforsystem(43):

1. Thelinearizationof(43) alongS¯a

0,1hasthreezeroeigenvalues,andthenonzeroeigenvalueλ= 32e1(1−

e1)(1− r1)(γ− e1), whichforr1= 0 correspondstotheﬂow in theinvariantplane (f1,e1).

2. Thereexists athree-dimensionalcentermanifoldWc

a,1 of thelinee1 which containstheplane of

equi-libria S¯a

0,1 andthe two-dimensional center manifold Ca,1.The manifoldWa,1c is attracting, and in the

setD1,deﬁnedby

D1:={(f1, r1, e1, ε1)| 0 ≤ r1≤ δ1, e1∈ I1, 0≤ ε1≤ α1},

isgivenby thegraph

f1= ha,1(r1, e1, ε1),

where I1 isasuitable interval,and α1,δ1 > 0 are suﬃciently small.Fortheparticular pointpa,1∈ e1

wheree0∈ I1,thefunctionha,1(r1,e0,ε1) hastheexpansion

ha,1(r1, e0, ε1) =

γ

2(e0− γ)ε1+ O(ε 2

1). (48)

3. ThereexistsK > 0 suchthat theorbitsthatare nearthecentermanifoldWc

a,1are attractedtoWa,1c by

an exponentialrate oforder O(exp(−Kt1)).

Proof. A straightforward calculation shows the ﬁrst claim. Due to the fact thatthe linearization of (43)

along S¯a

0,1 has three zero eigenvalues, there exists [2,13] an attracting three-dimensional center manifold

Wc

a,1 at thepointpa,1.To deriveequation(48),weﬁrstexpandf1 totheﬁrstorderofvariables r1,e1 and

ε1, andthen plug into (43). Bycomparingthe coeﬃcientsofr1,e1 andε1, equation(48) is obtained.The

last claimisprovenbythecentermanifoldtheoryappliedatthepointpa,1.

Remark 15.The attracting center manifold Wc

a,1 recovers parts of the slow manifold Sε,1a away form the

line 1, and extendsit into an O(ε) neighborhood of 1. The slow manifold Sε,1a is obtained as a section

ε= constant ofWc

a,1.InchartK1,thiscentermanifoldisgivenbythegraph(48).

Note thatin chartK1, our goalisto understand thedynamics (43) close to thecentermanifold Wa,1c ,

whichcorrespondstoasuﬃcientlysmallneighborhoodoftheslowmanifoldS¯a

0,1.Assumethatδ1,α1,β1> 0

are smallconstants.Letus deﬁnethesections

Δin₁ :={(f1, r1, e1, ε1)| (f1, r1, e1, ε1)∈ D1, r1= δ1},

Δout₁ :={(f1, r1, e1, ε1)| (f1, r1, e1, ε1)∈ D1, ε1= α1},

R₁in:={(f1, r1, e1, ε1)| (f1, r1, e1, ε1)∈ D1, r1= δ1,|f1| ≤ β1}.

(49)

NotethatbythewaywehavedeﬁnedΔin

1 ,weinfacthaveΔin1 = ¯Σ1:= Φ−1(Σ1×{[0,ρ1]}) forsomeρ1> 0,

see Fig.10. Furthermore, theconstants δ1,α1,β1 arechosen suchthatRin1 ⊂ Δin1 , andthe intersectionof

thecentermanifoldWc

a,1 withΔin1 liesinRin1 ,i.e.,Wa,1c ∩ Δin1 ⊂ Rin1 .

Let us denote Π1 as the transition map from Δin1 to Δout1 , induced by the ﬂow of (43). In order to

construct map Π1, we reduce system (43) to the center manifold Wa,1c and analyzethe system based on

the dynamicsonWc

a,1.Tothisend,bysubstituting(48) into(43) andrescalingtime,theﬂowofthecenter