Vakcinálási modellek számítógépes vizsgálata

Diplomamunka

Vizi Zsolt
Alkalmazott matematikus MSc

Témavezetők:

Dr. Röst Gergely
egyetemi adjunktus
Alkalmazott és Numerikus Matematika Tanszék

Dr. Karsai János
egyetemi docens
Analízis Tanszék

Bevezetés

1. Előzmények

1.1. Járványterjedési modellek

1.2. Backward bifurkáció

1.3. Impulzív vakcinálási modellek

2. SIVS modell folytonos vakcinálással

2.1. A modell

2.2. A betegségmentes egyensúlyi helyzet

2.3. Bifurkáció-analízis

3. SIVS modell impulzív vakcinálással

3.1. A modell

3.2. Periódus-leképezés és dinamika a betegségmentes altéren

3.3. A betegségmentes periodikus megoldás stabilitása

3.4. A betegségmentes periodikus megoldás globális stabilitása

3.5. Perzisztencia-tulajdonság (Subscript[R, c])>1 esetben

3.6. Bifurkáció- analízis

4. Numerikus vizsgálatok

4.1. Fázisképek, megoldásgörbék, attraktivitási tartományok

4.2. A bifurkációs diagram

5. Járványügyi kérdések

5.1. Gyakorlati példák backward bifurkációra

5.2. Impulzív modellhez tartozó reprodukciós szám vizsgálata

5.3. Backward bifurkáció vizsgálata a (φ,T) síkon

6. Összegzés

Hivatkozások

A. Melléklet

A .1. Az Impulse programcsomag

Ezúton szeretnék köszönetet mondani tanáraimnak, dr. Röst Gergelynek, akinek iránymutatásai és tanácsai révén formálódott és fejlődött a dolgozat, és dr. Karsai Jánosnak, a számítógépes kísérletek megoldása során nyújtott ötleteiért és kreativitásáért. Köszönöm fáradhatatlan türelmük és segítőkészségük. 

A kutatás az Európai Unió és Magyarország támogatásával a TÁMOP 4.2.4.A/2-11-1-2012-0001 azonosító számú „Nemzeti Kiválóság Program – Hazai hallgatói, illetve kutatói személyi támogatást biztosító rendszer kidolgozása és működtetése konvergencia program” című kiemelt projekt keretei között valósult meg.

A dolgozat célja az SIVS típusú járványterjedési modellek elméleti és számítógéppel segített vizsgálata volt. A vizsgálatok elsősorban a modellekben fellépő bifurkációkra fókuszáltak, de járványügyi kérdéseket is megfogalmaztunk, amelyek elméleti és numerikus módon is megválaszolásra kerültek. A következőkben a dolgozat fejezeteinek tartalmát foglaljuk össze.
Az első fejezetben a célkitűzésekhez szükséges előzményeket gyűjtöttük össze: először a járványtani modellek kompartment rendszerekkel történő felírását tisztáztuk, bevezettük a legfontosabb epidemiológiai paramétereket. Ezután a modellekben megjelenő bifurkációk típusait tekintettük át, kihangsúlyozva a különbségeket és az ebből adódó gyakorlati problémákat. Végül az impulzív vakcinálási stratégiához szükséges háttér került kidolgozásra: bevezettük a rögzített időpillanatban ható impulzust tartalmazó impulzív differenciálegyenletet.
A második fejezetben a folytonos SIVS-modellt vizsgáltuk. Az eredmények jó része Brauer [2] cikkéből valók, amelyben először bizonyított erre a modellre backward bifurkációt. Felírtuk az egyensúlyi helyzeteket, a betegségmentes egyensúly lokális aszimptotikus stabilitására vonatkozó feltételből pedig definiáltuk a modellhez tartozó kontroll reprodukciós számot. Az endemikus egyensúlyi helyzetek vizsgálatából pedig szükséges és elégséges feltételt tudtunk adni a backward bifurkáció megjelenésére: ehhez egy másodfokú egyenlet gyökvizsgálatát végeztük el.
A harmadik fejezet önálló eredményeket tartalmaz az impulzív SIVS-modellel kapcsolatban - ezek többsége a XXXI. Országos Tudományos Diákköri Konferenciára írt dolgozatból [1] valók. A legfontosabb különbség, hogy egyensúlyi pontok helyett periodikus megoldásokat kellett vizsgálnunk. Először a betegségmentes altéren adtunk meg az egyetlen periodikus megoldást (röviden BMPM, vagyis betegségmentes periodikus megoldás), ennek globális stabilitását bizonyítottuk a betegségmentes altéren. Ezután a lokális aszimptotikus stabilitásra adtunk feltételt, amiből definiáltuk az impulzív modellhez tartozó reprodukciós számot: kiderült, hogy eredményünk epidemiológiai szempontból is interpretálható. A globális stabilitáshoz (Subscript[R, c])<1 esetén újabb feltételre volt szükség, majd a betegség egyenletesen erős perzisztenciáját bizonyítottuk (Subscript[R, c])>1 esetben. Végül az egész téren adtunk szükséges és elégséges feltételt a backward bifurkáció megjelenésére, amelyhez a legfontosabb eszköz a Ljapunov-Schmidt redukciós eljárás volt. Tudomásunk szerint ez az első eredmény, amely impulzív modellben backward bifurkációt igazol.
A negyedik fejezetben a fentiekhez implementáltunk számítógépes kísérleteket: mindkét vakcinálási stratégia esetén felrajzoltuk a fázisképet, a megoldásgörbéket és egy konkrét paraméteregyüttes esetén az attraktivitási tartományokat. Az impulzív rendszer kezeléséhez az egyik témavezetőm, dr. Karsai János által fejlesztett Impulse package-et használtuk. Végül a bifurkációs diagramokat rajzoltuk fel: folytonos esetben egyszerűen csak gyökvizsgálat eredményeként kapott egyszerű feltételt kellett implementálnunk, így egy interaktív ábrát kaptunk. Az impulzív modell esetén a feltétel sokkal bonyolultabb volt, emellett az endemikus egyensúlyokra sincs explicit képlet, azonban numerikus úton, a peremérték problémáknál használt shooting method-hoz hasonló eljárással sikerült az instabil ágat is kirajzolnunk.
Az ötödik fejezetben járványügyi kérdésekkel foglalkoztunk, főleg az impulzív modellre koncentrálva. Először a backward bifurkáció különböző fertőzésekkel, betegségekkel kapcsolatos modellekben való felbukkanására mutattunk példát. Egy járvány esetén impulzív esetben általában a vakcinálási tényezőt (φ) és az impulzusok közötti időtartamot (T) tudjuk szabályozni, ezért megvizsgáltuk az (Subscript[R, c])=(Subscript[R, c])(φ,T) szintvonalait és interaktív ábrát készítettünk a “biztonságos” tartományról (vagyis azon (φ,T) párok tartományáról, amelyre az (Subscript[R, c])<1). Végül megvizsgáltuk, hogy egy backward bifurkációt okozó paraméter együttes esetén, ha a (φ,T)-t egy kritikus (Subscript[R, c]) értékhez tartozó szintvonal mentén változtatom, megmarad-e a bistabilitás: kiderült, hogy a szinvonalak mentén ugyanúgy viselkedik a rendszer, vagyis kihal a betegség vagy nagy kezdeti betegszám esetén (Subscript[R, c])<1 esetén is beáll egy endemikus egyensúlyba.
A dolgozat egy a Mathematica-ban készült interaktív dokumentum, amely dinamikus kísérletezéseket is lehetetővé tesz az “olvasó” számára. A használathoz szükséges Wolfram CDF-Player telepítő programja a mellékelt CD-n megtalálható, illetve a http://www.wolfram.com/cdf-player/ site-ról letölthető.

Előzmények

({{{, ({{x'(t)=f(x(t)), ha t≠(Subscript[t, i]),}, {x((SuperPlus[(Subscript[t, i])]))=(Subscript[J, i])(x((SuperMinus[(Subscript[t, i])]))), i=1,2,3,...}})}})

rendszer egy impulzív rendszer rögzített pillanatokban ható impulzusokkal. 
Valamely

({{{, ({{x'(t)=f(x(t)), ha t≠(Subscript[t, i]),}, {x((SuperPlus[(Subscript[t, i])]))=(Subscript[J, i])(x((SuperMinus[(Subscript[t, i])]))), i=1,2,3,...,}, {x((SuperPlus[(Subscript[t, 0])]))=(Subscript[x, 0]),  }})}})

kezdetiérték-probléma x(t,(SuperPlus[(Subscript[t, 0])]),(Subscript[x, 0])) megoldását az alábbi módon kapjuk. Legyen k:=inf{i∈ℕ: (Subscript[t, i])>(Subscript[t, 0])}. Vegyük az

({{{, ({{x'(t)=f(x(t)),  }, {x((SuperPlus[(Subscript[t, 0])]))=(Subscript[x, 0.]),  }})}})

kezdetiérték-probléma nem folytatható x(t) megoldását, amely értelmezve van valamilyen [(Subscript[t, 0]),τ) intervallumon. Ha τ≤(Subscript[t, k]), akkor a feladat megoldott, ha τ>(Subscript[t, k]), akkor tekintsük a megoldást a [(Subscript[t, 0]),(Subscript[t, k])] intervallumon, és folytassuk az

({{{, ({{x'(t)=f(x(t)),  }, {x((SuperPlus[(Subscript[t, k])]))=(Subscript[J, k])(x((SuperMinus[(Subscript[t, k])]))),  }})}})

kezdetiérték-probléma megoldásával. Fontos, hogy a megoldást mindig az x((SuperPlus[(Subscript[t, 0])]))=(Subscript[x, 0]) kezdeti feltétellel indítjuk és mindig a differenciálegyenlet megoldásával kezdjük. A konstrukció miatt a rendszer megoldásai szakaszonként folytonos függvények, elsőfajú szakadásuk van a (Subscript[t, i]) pillanatokban. Mindenütt balról folytonosak, és folytonosan differenciálhatók értelmezési tartományukon, kivéve a  t=(Subscript[t, i]) helyeket. Ha  t≠(Subscript[t, i]), akkor megoldásai az x'=f(t,x) egyenletnek, ha t=(Subscript[t, i]), akkor eleget tesznek az x((SuperPlus[(Subscript[t, i])]))=(Subscript[I, i])(x((SuperMinus[(Subscript[t, i])]))) ugrásfeltételnek. További részletek a [] könyvben.

(Subscript[x, n+1])=F((Subscript[x, n])),

ahol F((Subscript[x, n]))=(Subscript[I, n+1])(x((SuperMinus[(Subscript[t, n+1])])))=(Subscript[I, n+1])(Φ((Subscript[x, n]),(Subscript[t, n+1])-(Subscript[t, n]))), amelyben Φ((Subscript[x, 0]),t) a x'=f(x) differenciálegyenlethez tartozó, (Subscript[x, 0]) kezdeti értékből induló folyam. Ennek a leképezésnek vizsgálhatók az (SuperStar[x]) fixpontjai (ha léteznek), meghatározható a stabilitásuk. Az impulzív rendszer periodikus megoldásai, az F leképezés fixpontjai, ugyanolyan stabilitási tulajdonsággal.

SIVS modell folytonos vakcinálással

Az így megalkotott modellt az . ábra szemlélteti, és a következő egyenletrendszer írja le:

({{{, ({{S' (t)= μ-β S(t) I(t)-μ S(t)+γ I(t),,  }, {I'(t) = β S (t)I(t)-(μ+γ) I(t).,  }})}})

A modellre röviden SIS modellként hivatkozunk.

({{{, ({{S'(t)  = μ-β S(t)  I (t)-μ S(t) +γ I(t)-φ S(t)+θ V(t),,  }, {I' (t) = β S(t)  I (t)-(μ + γ) I(t)+σ β V(t)  I(t),,  }, {V'(t) =φ S(t)-σ β V(t)  I(t)-(μ+θ) V(t). ,  }})}})

Ezt a modellt SIVS modellnek nevezzük (. ábra).

Mivel fennáll, hogy  1 = S(t) + I(t)+V(t), ezért a rendszert kétdimenziósra redukálhatjuk:

({{{, ({{S'(t)  = μ-β S(t)  I (t)-μ S(t) +γ I(t)-φ S(t)+θ(1-S(t)-I(t)),,  }, {I' (t) = β S(t)  I (t)-(μ + γ) I(t)+σ β (1-S(t)-I(t)) I(t).,  }})}})

Az () modell egyensúlyi helyzeteit a következő egyenletrendszer megoldásai adják:

(Piecewise[{{ 0=μ-β S I -μ S +γ I-φ S+θ(1-S-I),}, {0=β S  I -(μ + γ) I+σ β (1-S-I) I.}}])

A betegségmentes egyensúlyi helyzetben I=0, ekkor a S-re kapjuk, hogy

(SuperStar[S])=((μ+θ)/(μ+θ+φ)).

Ha a () egyenletet linearizáljuk egy (S,I) egyensúlyi pontban, akkor a következő mátrixot kapjuk:

(({{-β I-μ-φ-θ, -β S+γ-θ}, {β I-σ β, β S-(μ+γ)- σ β (1-S-2 I)}})).

A betegségmentes egyensúlyban kiértékelve:

(({{-μ-φ-θ, -β (SuperStar[S])+γ-θ}, {-σ β, β (SuperStar[S])-(μ+γ)- σ β (1-(SuperStar[S]))}})).

Az (SuperStar[S])=((μ+θ)/(μ+θ+φ)) feltételt felhasználva megmutatható, hogy a sajátértékek negatívak, ha

(β/(μ+γ))((μ+θ+σ φ )/(μ+θ+φ))<1,

így definiálhatjuk a modellhez tartozó kontroll reprodukciós számot:

(Subscript[R, c])(φ):=(β/(μ+γ))((μ+θ+σ φ )/(μ+θ+φ)).

Ezzel a következő tételhez jutottunk:

A () egyenlet betegségmentes egyensúlyi helyzete lokálisan aszimptotikusan stabil, ha (Subscript[R, c])(φ)<1, instabil, ha (Subscript[R, c])(φ)>1.

Felhasználva, hogy I=0 esetén S+V=1 (Subscript[R, c])(φ) a (β/(μ+γ))((SuperStar[S])+σ (SuperStar[V])) alakot ölti. Látjuk, hogy a kontroll reprodukciós szám azt fejezi ki, hogy egy fertőzött hány másodlagos fertőzést generál, mivel I=0 esetén a fogékonyak száma (SuperStar[S]), fertőzhető vakcináltak száma σ (SuperStar[V]), a ténylegesen fertőzhető populáció (SuperStar[S])+σ (SuperStar[V]) méretű, ezt az átviteli rátával (β) és az átlagos betegségben töltött idővel ((1/(μ+γ))) kell szoroznunk.

Tekintsük ismét a () egyenletrendszert. Ha 0≤σ<1, akkor az egyenletrendszer további megoldásai az endemikus egyensúlyi helyzeteket szolgáltatja. Fejezzük ki S-et az első egyenletből, majd helyettesítsük a második egyenletbe, így az

A (I^2)+B I+C=0

egyenletet kapjuk, ahol

A=σ β,
B=(μ+θ+σ φ)+σ(μ+γ)-σ β,
C=(((μ+γ) (μ+θ+φ))/(β (μ+θ)))-(μ+θ+σ φ).

Speciálisan, ha σ=0, akkor a () egyenlet elsőfokú és egyetlen megoldása (Subscript[I, e])=1-(((μ+γ) (μ+θ+φ) )/(β (μ+θ+σ φ )))=1-(1/((Subscript[R, c])(φ))), ami pozitív akkor és csakis akkor, ha (Subscript[R, c])(φ)>1. Tehát σ=0 esetén egyetlen endemikus egyensúly létezik, ha (Subscript[R, c])(φ)>1 és (Subscript[I, e])->0, amint  (Subscript[R, c])(φ)->(SuperPlus[1]). Ebben az esetben tehát nem lehet endemikus egyensúly  (Subscript[R, c])(φ)<1 esetén, vagyis  (Subscript[R, c])(φ)-t véve bifurkációs paraméternek, ebben az esetben nem jelenhet meg backward bifurkáció  (Subscript[R, c])(φ)=1-nél.

Látható, hogy a C előjelét az  (Subscript[R, c])(φ) egyértelműen meghatározza: ha  (Subscript[R, c])(φ)>1, akkor C<0,  (Subscript[R, c])(φ)=1 esetén C=0 és ha  (Subscript[R, c])(φ)<1, akkor C>0. Ha σ>0, akkor () egy másodfokú egyenlet és (Subscript[R, c])(φ)>1 esetén egyetlen pozitív endemikus egyensúlyi helyzet létezik. Ha (Subscript[R, c])(φ)=1, akkor C=0 és az egyetlen nemnulla megoldás (Subscript[I, e])=(B/A), ami éppen akkor pozitív, ha B<0. Összegezve, (Subscript[R, c])(φ)=1 esetén van egy pozitív endemikus egyensúly, ha B<0 és C=0. Mivel az egyensúlyi helyzetek folytonosan függenek φ-től, ezért szükségképpen létezik egy intervallum (Subscript[R, c])(φ)=1-től balra, ahol két pozitív egyensúlyi helyzet van, amelyek

(Subscript[I, e])=((-B±(Sqrt[(B^2)–  4A C]))/(2 A))

Így a következő tételt fogalmazhatjuk meg:

A () egyenlet esetén (Subscript[R, c])(φ)=1-ben backward bifurkáció lép fel akkor és csak akkor, ha B<0 és C=0.

Az előbbiekhez hasonló számolásokkal az eredeti SIS-modellhez tartozó alap reprodukciós szám is kiszámolható: (Subscript[R, 0])=(β/(μ+γ)). Ha az SIVS-modell esetén feltesszük, hogy képesek vagyunk a fogékony populáció egészét azonnal vakcinálni (formálisan, φ->∞), illetve, hogy a vakcinált nem veszítheti el védettségét (vagyis θ=0), akkor az () rendszer a következő egyenlettel ekvivalens:

I'=σ β I (1-I)-(μ+γ) I

Igazolható, hogy ehhez az egyenlethez definiálható reprodukciós szám (Subscript[R, 0]^( * ))=((σ β)/(μ+γ))=σ (Subscript[R, 0]) alakú; mivel 0≤σ≤1, ezért  (Subscript[R, 0]^( * ))≤(Subscript[R, 0]). Ezek alapján, ha (Subscript[R, 0])<1, akkor a vakcinálás felgyorsítja a fertőzés eltűnését, de az a vakcina nélkül is kihalna. Ha (Subscript[R, 0]^( * ))>1, akkor a vakcina valamelyest kontrollálja a fertőzők számát, azonban önmagában nem képes a betegséget eltüntetni. Közegészségügyi szempontból egy érdekes eset, ha (Subscript[R, 0]^( * ))<1<(Subscript[R, 0]), vagyis a fertőzés vakcinálás nélküli esetben perzisztens, de megfelelő vakcinával dolgozva a betegség kihal. Ebben az esetben szükségünk van arra a (Subscript[φ, 0]) vakcinálási tényezőre, ami képes az (Subscript[R, 0]) értékét 1-re csökkenteni.

(Subscript[φ, 0]):=(μ+θ) ((β-(μ+γ))/(μ+γ-σ β))=(μ+θ) (((Subscript[R, 0])-1)/(1-(Subscript[R, 0]^( * ))))

mértékű vakcinálási tényező szükséges ahhoz, hogy a fertőzést kontrollálhassuk. Matematikai szempontból is fontos kérdés, hogy ha (Subscript[R, c])(φ)=1-ben backward bifurkáció lép fel, akkor mekkora vakcinálási rátára van szükségünk, hiszen már korábban is láttuk, hogy ebben az esetben nem elegendő 1-re csökkenteni a reprodukciós számot.

SIVS modell impulzív vakcinálással

Az 1.3. részben bevezetett jelölésekkel a következőképpen írhatjuk fel az impulzív SIVS modellt (. ábra):

({{{, ({{S'(t)=μ-β S(t)  I (t)-μ S(t)+γ I(t)+θ V(t),  }, {I' (t) = β S(t)  I (t)-(μ+γ) I(t)+σ β V(t)  I(t), t≠n Tesetén, }, {V'(t) =-σ β V(t)  I(t)-(μ+θ) V(t) ,  }})}})
	({{{, ({{S(n (SuperPlus[T]))=(1-φ) S(n (SuperMinus[T])),  }, {I(n (SuperPlus[T]))=I(n (SuperMinus[T])), t=n T esetén. }, {V(n (SuperPlus[T]))=V(n (SuperMinus[T]))+φ S(n (SuperMinus[T])),  }})}})

Értelemszerűen a kezdeti értékeket is úgy választjuk, hogy S(0)+I(0)+V(0)=1 és S(0)≥0, I(0)≥0, V(0)≥0. Az egyenlet alakjából látjuk, hogy a megoldások végig nemnegatívak maradnak, vagyis az

Ω={(S,I,V)∈(ℝ^3):S+I+V=1, S≥0, I≥0, V≥0}

tartomány invariáns. Mivel S(t) + I(t)+V(t)=1, ezért a háromdimenziós modellt kétdimenziósra redukálhatjuk. A vakcináltakra vonatkozó egyenletet elhagyva, a többi egyenletbe a V(t)=1 -S(t) - I(t) összefüggést helyettesítve

({{{, ({{S'(t)=μ+θ-β S(t) I (t)-(μ+θ) S(t)+(γ-θ) I(t) ,,  }, {I'=β S (t) I (t)-(μ+γ) I(t)+σ β (1-S(t)-I(t) ) I(t) , t≠n Tesetén, }})}})
	({{{, ({{S(n (SuperPlus[T]))=(1-φ) S(n (SuperMinus[T])),,  }, {I(n (SuperPlus[T]))=I(n (SuperMinus[T])),  t=n T esetén. }})}})

A továbbiakban ezt a rendszert fogjuk vizsgálni.

S'(t)=μ+θ-(μ +θ)S(t)

S(n (SuperPlus[T]))=(1-φ)S(n (SuperMinus[T]))

({{{, ({{x'(t)+P x(t)=Q,,  }, {x((Subscript[t, n]))=(Subscript[x, n]).,  }})}})

Ekkor az általános megoldás a következő alakú:

x(t)=(Q/P)+(e^(-P(t-(Subscript[t, n]))))((Subscript[x, n])-(Q/P)).

S(t)=1+(e^(-(μ+θ)(t-n T)))(S(n T)-1).

Jelölje (Subscript[S, n+1]) az (n+1). beavatkozás utáni fertőzhető részpopuláció nagyságát, vagyis:

(Subscript[S, n+1])=S((n+1) (SuperPlus[T]))=(1-φ) S((n+1)(SuperMinus[T]))=(1-φ)(1+(e^(-(μ+θ)T))(S(n (SuperPlus[T]))-1))=(1-φ)(1+(e^(-(μ+θ)T))((Subscript[S, n])-1)).

Legyen

ψ(S)=(1-φ)(1+(e^(-(μ+θ)T))(S-1)).

Így (Subscript[S, n+1])=ψ((Subscript[S, n])). A periodikus megoldáshoz a ψ fixpontját kell megtalálnunk:

ψ((SuperStar[S]))=(SuperStar[S]).

Kapjuk, hogy

(SuperStar[S])=(((1-(e^(-(μ+θ)T)))(1-φ))/(1-(1-φ)(e^(-(μ+θ)T)))).

Általában egy

({{{, ({{x'(t)+P x(t)=Q,,  }, {x(n (SuperPlus[T]))=(1-φ)x(n (SuperMinus[T])),  }})}})

impulzív rendszer esetén:

(SuperStar[x])=(((1-φ)(Q/P)(1-(e^(-P T))))/(1-(1-φ)(e^(-P T)))).

A () egyenletrendszer betegségmentes periodikus megoldása a következő alakú:

(OverBar[S])(t)=1+(e^(-(μ+θ)(t-n T)))((((1-(e^(-(μ+θ)T)))(1-φ))/(1-(1-φ)(e^(-(μ+θ)T))))-1), ha n T<t<(n+1)T,
(OverBar[S])(n T)=(((1-(e^(-(μ+θ)T)))(1-φ))/(1-(1-φ)(e^(-(μ+θ)T)))),
(OverBar[I])(t)=0, t≥0.

Ez a megoldás a betegségmentes altéren globálisan aszimptotikusan stabil.

() T-periodikus, ugyanis  ha t≠n T, akkor

(OverBar[S])(t)=1+(e^(-(μ+θ)(t-n T)))((((1-(e^(-(μ+θ)T)))(1-φ))/(1-(1-φ)(e^(-(μ+θ)T))))-1)=1+(e^(-(μ+θ)(t+T-(n+1) T)))((((1-(e^(-(μ+θ)T)))(1-φ))/(1-(1-φ)(e^(-(μ+θ)T))))-1)=(OverBar[S])(t+T).

Ha t=n T, akkor (OverBar[S])(t)=(SuperStar[S])=(OverBar[S])((n+1)T).

Megmutatjuk, hogy () megoldása ()-nak. Mivel a függvények periodikusak, ezért a megoldást elég a [0,T] intervallumon vizsgálni, ahol

(OverBar[S])(t)=1+(e^(-(μ+θ)t))((((1-(e^(-(μ+θ)T)))(1-φ))/(1-(1-φ)(e^(-(μ+θ)T))))-1).

Ezért

((OverBar[S])(t))'=-(μ+θ)(e^(-(μ+θ)t))((((1-(e^(-(μ+θ)T)))(1-φ))/(1-(1-φ)(e^(-(μ+θ)T))))-1)=(μ+θ)(1-(OverBar[S])(t))= μ +θ- β (OverBar[S])(t) 0 - (μ +θ)(OverBar[S])(t)  + (γ-θ)0,

((OverBar[I])(t))'=0= β S (t)0 - (μ + γ)0 +σ β (1-S(t) -0 )0.

Tehát az ((OverBar[S])(t),(OverBar[I])(t)) megoldása a () rendszernek. Az (SuperStar[S]) fixpont konstrukciója miatt t=n T pontokban az impulzív feltétel is teljesül.

A stabilitáshoz a leképezés fixpontjának stabilitását kell ellenőriznünk. A ψ fixpontja pontosan akkor aszimptotikusan stabil, ha  |ψ'((SuperStar[S]))|<1. Elemi számolás adja, hogy

0<(ψ'((SuperStar[S])))'=(1-φ)((e^(-T(μ+θ))))<1.

Tehát az (SuperStar[S]) aszimptotikusan stabil fixpontja ψ-nek. Mivel a ψ(S) lineáris leképezés, ezért (SuperStar[S]) globálisan aszimptotikusan stabil fixpontja ψ-nek.

Egy interaktív ábra a ()-() rendszer megoldásaira:

A () egyenletrendszer betegségmentes periodikus megoldása lokálisan aszimptotikusan stabil, ha (1/T)(∫_0^T)(OverBar[S])(u)+σ(1-(OverBar[S])(u))du<((μ+γ)/β).

Linearizáljuk az egyenletrendszert a betegségmentes egyensúly körül. Legyen

({{{, ({{S(t)=s(t)+(OverBar[S])(t),  }, {I(t)=i(t)+(OverBar[I])(t).,  }})}})

Így a következő differenciálegyenletrendszert kapjuk:

({{{, ({{s'(t) = - β (OverBar[S])(t)  i (t) - μ s(t)  + γ i(t) +θ (1-s(t) -i(t) ),  }, {i' (t)= β (OverBar[S]) (t) i (t) - (μ + γ) i(t) +σ β (1-(OverBar[S])(t)) i(t).,  }})}})

A () egyenletrendszernek a (0,0) egyensúlyi helyzete. Vizsgáljuk meg ennek a stabi-litását! Ehhez definiáljuk az egyenletrendszer A(t) alapmátrixát:

A(t)=(({{(Subscript[s, 1])(t), (Subscript[s, 2])(t)}, {(Subscript[i, 1])(t), (Subscript[i, 2])(t)}})),

ahol ((Subscript[s, 1])(t),(Subscript[i, 1])(t)) és ((Subscript[s, 2])(t),(Subscript[i, 2])(t)) az () egyenletrendszer ((Subscript[s, 1])(0),(Subscript[i, 1])(0)) és ((Subscript[s, 2])(0),(Subscript[i, 2])(0)) kezdeti értékhez tartozó megoldásai.
Legyen

({{(Subscript[s, 1])(0)=1,, (Subscript[s, 2])(0)=0,}, {(Subscript[i, 1])(0)=0,, (Subscript[i, 2])(0)=1.}})

Floquet-elméletet fogunk alkalmazni, ezért először megoldjuk a () közönséges differenciálegyenlet-rendszert a 0≤t≤T intervallumon az adott kezdőfeltételekkel.

(Subscript[s, 1])(t)=(e^(-(μ+θ)t)), (Subscript[i, 1])(0)=0, (Subscript[i, 2])(t)=(e^((∫_0^t) β (OverBar[S]) (u)- (μ + γ)+σ β (1-(OverBar[S])(u)) )du)).

Ezeket kiértékelve a t=T-ben és alkalmazva az impulzust, megkapjuk az

M:=(({{(1-φ)(e^(-(μ+θ)T)), (Subscript[s, 2])(T)}, {0, (e^((∫_0^T) β (OverBar[S]) (u)- (μ + γ)+σ β (1-(OverBar[S])(u)) )du))}}))

monodrómia mátrixot. Az (Subscript[s, 2])(T) meghatározása már nehezebb, de nincs is rá szükségünk, mert a stabilitásvizsgálathoz csak az M mátrix sajátértékei kellenek és a karakterisztikus egyenletben az (Subscript[s, 2])(T) együtthatója 0. M karakterisztikus egyenlete:

((1-φ)(e^(-(μ+θ)T))-λ)((e^((∫_0^T) β (OverBar[S]) (u)- (μ + γ)+σ β (1-(OverBar[S])(u)) )du))-λ)=0,

amelynek gyökei: (Subscript[λ, 1])=(1-φ)(e^(-(μ+θ)T)) és (Subscript[λ, 2])=(e^((∫_0^T) β (OverBar[S]) (u)- (μ + γ)+σ β (1-(OverBar[S])(u)) )du)).

Az egyenletrendszer (0,0) megoldása lokálisan stabil, ha az előbb kiszámolt sajátértékek abszolútértékben 1-nél kisebbek. A (Subscript[λ, 1]) mindig 1 alatt marad, mivel az exponenciális kifejezés kitevője negatív. A (Subscript[λ, 2])<1 egyenlőtlenség pontosan akkor teljesül, ha a kitevő negatív, tehát

(∫_0^T) β (OverBar[S]) (u)- (μ + γ)+σ β (1-(OverBar[S])(u)) )du=
=(∫_0^T) β (OverBar[S]) (u)+σ β (1-(OverBar[S])(u)) )du- (μ + γ)T<0.

Átrendezve,

(1/T)(∫_0^T) (OverBar[S]) (u)+σ (1-(OverBar[S])(u)) )du< ((μ + γ)/β).

Tehát a betegségmentes megoldás lokálisan aszimptotikusan stabil, ha (1/T)(∫_0^T) (OverBar[S]) (u)+σ (1-(OverBar[S])(u)) )du< ((μ + γ)/β) .

Vegyük észre, hogy ha a tételben szereplő stabilitási feltétel jobb oldalát 1-re rendezzük, akkor a következőt ekvivalens egyenlőtlenséget kapjuk:

(Subscript[R, c]):=(β/(,μ+γ))((1/T)(∫_0^T) (OverBar[S]) (u)du+σ(1/T)(∫_0^T) (1-(OverBar[S])(u)) du)<1.

Mivel a betegségmentes altéren I(t)=0, ezért (OverBar[V])(t)=1-(OverBar[S])(t). Így az (Subscript[R, c]) a következő alakot ölti:

(Subscript[R, c])=(β/(,μ+γ))((1/T)(∫_0^T) (OverBar[S]) (u)du+σ(1/T)(∫_0^T) (OverBar[V])(u) du).

Látjuk, hogy a kontroll reprodukciós szám azt fejezi, hogy a járvány kezdetén egy fertőzött hány másodlagos fertőzést generál (a populációban átlagosan (1/T)(∫_0^T) (OverBar[S]) (u)du fogékony és σ(1/T)(∫_0^T) (OverBar[V])(u) du olyan vakcinált tartózkodik, amely elvesztheti védettségét; ezek β rátával találkoznak a fertőző egyedekkel) a betegség teljes ideje ((1/(,μ+γ))) alatt. Mivel pontosan ismerjük  (OverBar[S]) (t) alakját, ezért explicit meghatározható az (Subscript[R, c]) is:

(Subscript[R, c])=((β )/(μ+γ))(σ+(1/T)(1-σ)(T-((1 )/(μ+γ))(1-(e^(-T(μ+θ))))(((1-(e^(-T(μ+θ))))(1-φ))/(1-(1-φ)(e^(-T(μ+θ))))))).

Tehát a stabilitási feltétel azzal ekvivalens, hogy az (Subscript[R, c]) kontroll reprodukciós szám kisebb, mint 1.

Ha γ-θ≤0 , akkor (Subscript[R, c])<1 esetén a betegségmentes periodikus megoldás globálisan aszimptotikusan stabil.

Mivel S'(t)=μ+θ-β S(t) I (t)-(μ+θ) S(t)+(γ-θ) I(t) ≤μ+θ-(μ+θ) S(t) a tételben és a paraméterekre megfogalmazott feltételek miatt, ezért tekintsük a következő összehasonlító rendszert:

({{{, ({{x'(t)=μ-μ x(t)+θ-θ x(t),,  }, {x(n (SuperPlus[T]))=(1-φ)x(n (SuperMinus[T])).,  }})}})

Ennek a rendszernek a megoldása a betegségmentes periodikus megoldás, amiről már láttuk, hogy globálisan aszimptotikusan stabil a betegségmentes altéren. Az impulzív diffe-renciálegyenletekre vonatkozó összehasonlítási elv értelmében bármely (Subscript[ϵ, 1])>0 esetén létezik (Subscript[m, 1])≥0  egész szám, melyre

S(t)<(OverBar[S])(t)+(Subscript[ϵ, 1]), 	n T<t<(n+1)T, n T>(Subscript[m, 1])T.

A () másik egyenletéből:

I'(t)≤I(t) (β ((OverBar[S]) (t)+(Subscript[ϵ, 1]))-(μ+γ)+σ β-σ β ((OverBar[S])(t)+(Subscript[ϵ, 1]))))
n T<t<(n+1)T, n T>(Subscript[m, 1])T.

Ez alapján definiáljuk a következő összehasonlító rendszert:

({{{, ({{y'(t)=y(t) (β(1-σ) ((OverBar[S]) (t)+(Subscript[ϵ, 1]))-(μ+γ-σ β )),,  }, {y(n (SuperPlus[T]))=y(n (SuperMinus[T])).,  }})}})

így y(t)≥I(t), amennyiben y(0)≥0. Ha az egyenletet (n T,(n+1) T]-n integráljuk, a következőt kapjuk:

y((n+1) T) = y(n T) (e^((∫_(n T)^((n+1) T))β(1-σ) ((OverBar[S]) (t)+(Subscript[ϵ, 1])) - (μ + γ-σ β )dt)).

Ebből iterációval adódik, hogy

y(n T) = (Subscript[y, 0]) (e^(n(∫_0^T)β(1-σ) ((OverBar[S]) (t)+(Subscript[ϵ, 1])) - (μ + γ-σ β )dt)).

ahol (Subscript[y, 0])=y((SuperPlus[0]))>0. Ha (Subscript[R, c])<1, vagyis (1/T)(∫_0^T) (OverBar[S]) (u)+σ (1-(OverBar[S])(u)) )du< ((μ + γ)/β), akkor (Subscript[ϵ, 1])-et választhatjuk olyan kicsire, hogy (1/T)(∫_0^T) (OverBar[S]) (u)+σ (1-(OverBar[S])(u)) +(Subscript[ϵ, 1]))du< ((μ + γ)/β) is teljesüljön. Ekkor (lim_(n->∞))y(n T)=0, mivel a kitevő negatív. Emellett felírhatjuk a () rendszer megoldását is:

y(t)=y(n T) (e^((∫_(n T)^t)β(1-σ) ((OverBar[S]) (s)+(Subscript[ϵ, 1])) - (μ + γ-σ β )ds)) , 	n T<t<(n+1)T.

amiből (lim_(t->∞))y(t)=0. Így megkaptuk, hogy (lim_(t->∞))I(t)=0. 
A célunk megmutatni, hogy (Subscript[R, c])<1 esetén S(t)-> (OverBar[S]) (t), amint t->∞. Az előzőből következik, hogy minden (Subscript[ϵ, 2])>0 esetén létezik (Subscript[m, 2])≥0, (Subscript[m, 2])∈ℤ, (Subscript[m, 1])<(Subscript[m, 2]), melyre I(t)<(Subscript[ϵ, 2]), amikor  t>(Subscript[m, 2])T>(Subscript[m, 1])T.
Újra () első egyenletét felhasználásával kapjuk, hogy

S'(t)>μ-β S(t) (Subscript[ϵ, 2])-μ S(t)+γ (Subscript[ϵ, 2])+θ (1-S(t)-(Subscript[ϵ, 2]) ).

Ezzel definiáljuk a következő összehasonlító rendszert:

({{{, ({{z'(t)=μ-β z(t) (Subscript[ϵ, 2])-μ z(t)+γ (Subscript[ϵ, 2])+θ (1-z(t)-(Subscript[ϵ, 2]) ),,  }, {z(n (SuperPlus[T]))=(1-φ)z(n (SuperMinus[T])).,  }})}})

Felhasználva a () és () képletet, kiszámolható a (OverBar[z])(t) T-periodikus megoldás és (SuperStar[z]) fixpont, amelyről belátható, hogy globálisan aszimptotikusan stabil. (OverBar[z])(t) a következő alakú:

(OverBar[z])(t)=((μ+γ (Subscript[ϵ, 2])+θ(1-(Subscript[ϵ, 2])))/(μ+θ+β (Subscript[ϵ, 2])))+(e^(-(μ+θ+β (Subscript[ϵ, 2]))(t-n T)))((((1-φ)((μ+γ (Subscript[ϵ, 2])+θ(1-(Subscript[ϵ, 2])))/(μ+θ+β (Subscript[ϵ, 2])))(1-(e^(- (μ+θ+β (Subscript[ϵ, 2]))T))))/(1-(1-φ)(e^(-(μ+θ+β (Subscript[ϵ, 2])) T))))+((μ+γ (Subscript[ϵ, 2])+θ(1-(Subscript[ϵ, 2])))/(μ+θ+β (Subscript[ϵ, 2])))),n T<t<(n+1)T.

Legyen z(t) az összehasonlító rendszer tetszőleges (Subscript[z, 0])=z((SuperPlus[0]))>0 kezdeti értékből induló megoldása. Az impulzív differenciálegyenletekre vonatkozó összehasonlítási elvből tudjuk, hogy a () rendszer bármelyik megoldására, amelyik az  (Subscript[S, 0])=S((SuperPlus[0]))=(Subscript[z, 0])>0 és (Subscript[I, 0])=I((SuperPlus[0]))>0 kezdeti értékekből indul, létezik (Subscript[m, 3])≥0, (Subscript[m, 3])∈ℤ, (Subscript[m, 1])<(Subscript[m, 2])<(Subscript[m, 3]), hogy

S(t)>(OverBar[z])(t)-(Subscript[ϵ, 2]), n T<t<(n+1)T, n T>(Subscript[m, 3])T.

Így (OverBar[z])(t)-(Subscript[ϵ, 2])<S(t)<(OverBar[S])(t)+(Subscript[ϵ, 1]), ahol (OverBar[z])(t)->(OverBar[S])(t) egyenletesen, amint (Subscript[ϵ, 2])->0. Mivel (Subscript[ϵ, 1]) és (Subscript[ϵ, 2]) tetszőlegesen kicsik, S(t)->(OverBar[S])(t), amint t->∞. Tehát az (OverBar[S])(t) globálisan attraktív,  vagyis a rendszer betegségmentes megoldása globálisan aszimptotikusan stabil.

Ha (Subscript[R, c])>1, akkor létezik egy pozitív (Subscript[m, I]) konstans, hogy a () rendszer minden megoldása esetén (Underscript[liminf, t->∞]) I(t)≥(Subscript[m, I]) , vagyis a betegségre erős egyenletes perzisztencia teljesül.

Először tekintsünk egy tetszőleges (SuperStar[I])∈(0,1) számot. Tegyük fel, hogy létezik (Subscript[t, 0])>0, hogy I(t)<(SuperStar[I]) minden t≥(Subscript[t, 0]) esetén. Ekkor () első egyenletéből

S'(t)>μ-β (SuperStar[I])-μ S(t)+θ (1-S(t)-(SuperStar[I]) )

t≥(Subscript[t, 0]) esetén. Tekintsük a következő összehasonlító rendszert

({{{, ({{w'(t)=μ-β (SuperStar[I])-μ w(t)+θ (1-w(t)-(SuperStar[I]) ),,  }, {w(n (SuperPlus[T]))=(1-φ)w(n (SuperMinus[T])),,  }, {w((SuperPlus[(Subscript[(OverTilde[t]), 0])]))=0,  }})}})

t≥(Subscript[(OverTilde[t]), 0]) esetén, ahol (Subscript[(OverTilde[t]), 0]):=⌈((Subscript[t, 0])/T)⌉T. Ennek a rendszernek is kiszámolható a (OverBar[w]) T-periodikus, globálisan aszimptotikus megoldása

(OverBar[w])(t)=((μ+θ-(β+θ) (SuperStar[I]))/(μ+θ))+(e^(-(μ+θ)(t-n T)))((((1-φ) ((μ+θ-(β+θ) (SuperStar[I]))/(μ+θ)) (1-(e^(- (μ+θ)T))) )/(1-(1-φ)(e^(-(μ+θ) T))))+((μ+θ-(β+θ) (SuperStar[I]))/(μ+θ)))
n T<t<(n+1)T,

így minden ϵ>0 esetén létezik olyan  (Subscript[t, 1])>(Subscript[(OverTilde[t]), 0]), hogy

w(t)>(OverBar[w])(t)-ϵ

t>(Subscript[t, 1]) esetén. Az impulzív differenciálegyenletekre vonatkozó összehasonlítási elvet felhasználva

S(t)≥w(t)>(OverBar[w])(t)-ϵ

t>(Subscript[t, 1]) esetén. () második egyenletét felhasználva:

I'(t)≥I(t)(β ((OverBar[w])(t)-ϵ)-(μ + γ)+σ β (1-((OverBar[w])(t)-ϵ))),

t>(Subscript[t, 1]) esetén. Legyen (SuperStar[K])∈(SuperPlus[ℤ]) és (SuperStar[K])T>(Subscript[t, 1]) és integráljuk (n T,(n+1) T]-on az előző egyenlőtlenséget

I((n+1)T)≥I(n T)(e^((∫_(n T)^((n+1) T))β ((OverBar[w])(t)-ϵ)-(μ + γ)+σ β (1-((OverBar[w])(t)-ϵ))dt)),

amiből

I(((SuperStar[K])+n)T)≥I((SuperStar[K]) T) (e^(n(∫_0^T)β ((OverBar[w])(t)-ϵ)-(μ + γ)+σ β(1- ((OverBar[w])(t)-ϵ))dt)).

Legyen α:=(∫_0^T)β ((OverBar[w])(t)-ϵ)-(μ + γ)+σ β(1- ((OverBar[w])(t)-ϵ))dt. Látjuk, hogy ϵ->0,  (SuperStar[I])->0 esetén (OverBar[w])(t)->(OverBar[S])(t). Ha ϵ és (SuperStar[I]) elég kicsi, akkor az (Subscript[R, c])>1 feltétel átrendezéséből kapjuk, hogy α>0. Így

I(((SuperStar[K])+n)T)≥I((SuperStar[K]) T) (e^(n α))->∞, amint n->∞.

Ez ellentmond I(t) korlátosságának, tehát ha (SuperStar[I]) elegendően kicsi, akkor nincs olyan (Subscript[t, 0]) , hogy I(t)≤(SuperStar[I]) t≥(Subscript[t, 0]) esetén. Emiatt két eset lehetséges: létezik (Subscript[t, 2])>(Subscript[t, 1]), melyre I(t)>(SuperStar[I]) minden t≥(Subscript[t, 2])-re vagy I(t) oszcillál (SuperStar[I]) körül. Az első eset épp a perzisztenciát igazolja, (Subscript[m, I]):=(SuperStar[I]), ezért a továbbiakban az oszcilláló esettel foglalkozunk.

Az előző levezetésből látszik, hogy van olyan (Subscript[t, 3]), hogy I((Subscript[t, 3]))>(SuperStar[I]). Legyen (Subscript[t, 4])>(Subscript[t, 3])  tetszőleges úgy, hogy I((Subscript[t, 4]))<(SuperStar[I]).  Legyen (SubStar[t]):=(sup_(t<(Subscript[t, 4]))){I(t)≥(SuperStar[I])} és (SuperStar[t]):=(inf_(t>(Subscript[t, 4]))){I(t)≥(SuperStar[I])}, ekkor I(t)<(SuperStar[I]), ha t∈((SubStar[t]),(SuperStar[t])), és I(t) folytonossága miatt  I((SubStar[t]))=I((SuperStar[t]))=(SuperStar[I]). Felső becslést fogunk adni a (SuperStar[t])-(SubStar[t]) különbségre. Legyen (Subscript[τ, 1]):=⌈((SubStar[t])/T)⌉T-(SubStar[t]), (Subscript[τ, 2]):=(Subscript[t, 1])-(Subscript[(OverTilde[t]), 0]), (Subscript[τ, 3]):=⌈(((Subscript[t, 1])+(Subscript[τ, 1])+(Subscript[τ, 2]))/T)⌉T-((Subscript[t, 1])+(Subscript[τ, 1])+(Subscript[τ, 2])). Legyen (Subscript[n, 0])=min{n∈N: (e^(-(2+(Subscript[τ, 2]))(μ + γ+σ β)+n α))>1}.   Megmutatjuk, hogy (SuperStar[t])-(SubStar[t])<2+(Subscript[τ, 2])+(Subscript[n, 0])T. (Vegyük észre, hogy (Subscript[τ, 1])<1, (Subscript[τ, 3])<1, (Subscript[τ, 2])-t pedig az () rendszer a megoldástól függetlenül definiálja, így (Subscript[n, 0]) is független a megoldástól.) Mivel α>0, ezért ilyen (Subscript[n, 0]) létezik. Tegyük fel, hogy (SuperStar[t])-(SubStar[t])>2+(Subscript[τ, 2])+(Subscript[n, 0])T. Ekkor

I'(t)≥-(μ + γ+σ β) I(t)

miatt I((SubStar[t])+(Subscript[τ, 1])+(Subscript[τ, 2])+(Subscript[τ, 3]))≥(SuperStar[I])(e^(-((Subscript[τ, 1])+(Subscript[τ, 2])+(Subscript[τ, 3]))(μ + γ+σ β)))>(SuperStar[I])(e^(-(2+(Subscript[τ, 2]))(μ + γ+σ β))). 
Ha t>(SubStar[t])+(Subscript[τ, 1])+(Subscript[τ, 2])+(Subscript[τ, 3]), akkor a ()-ban leírt egyenlőtlenség teljesül, így

I((SubStar[t])+(Subscript[τ, 1])+(Subscript[τ, 2])+(Subscript[τ, 3])+(Subscript[n, 0])T)≥(SuperStar[I])(e^(-((Subscript[τ, 1])+(Subscript[τ, 2])+(Subscript[τ, 3]))(μ + γ+σ β)))(e^((Subscript[n, 0])α))>(SuperStar[I])(e^(-(2+(Subscript[τ, 2]))(μ + γ+σ β)))(e^((Subscript[n, 0])α))>(SuperStar[I]), ami ellentmondás. Tehát (SuperStar[t])-(SubStar[t])<2+(Subscript[τ, 2])+(Subscript[n, 0])T és I(t)≥(SuperStar[I])(e^(-(μ + γ+σ β)(2+(Subscript[τ, 2])+(Subscript[n, 0])T)))=:(Subscript[m, I]) minden t∈((SubStar[t]),(SuperStar[t])) esetén, ezért speciálisan (Subscript[t, 4]) esetén is, ami tetszőleges volt. Ezzel a tételt bebizonyítottuk.

A dolgozat elsődleges célja már az elején is az volt, hogy választ kapjunk arra, vajon a SIVS modellben impulzív vakcinálás esetén is fellép-e a backward bifurkáció jelensége. A numerikus szimulációk azt mutatták, hogy a két vakcinálási modell lényegesen nem különbözik egymástól, így találtunk a paramétereknek egy olyan együttesét, amely esetén bistabi-litás lépett fel: ha a kezdeti értékek nem voltak közel egymáshoz, akkor a megoldások különböző periodikus megoldáshoz konvergáltak. Tehát a sejtésünk az volt, hogy az impulzív SIVS modellben is felléphet a backward bifurkáció. Ennek matematikai bizonyítása azonban már közel sem triviális. A következőkben egy A. Lakmeche és O. Arino által 2000-ben publikált cikk [] technikáit módosítjuk és alkalmazzuk, amelynek alapgondolata a Ljapunov-Schmidt módszer és az implicitfüggvény-tétel használata. A [] cikkel ellentétben nem a impulzusok közti időtartamot, hanem a modell interpretációjához jobban illeszkedve, a φ vakcinálási tényezőt választottuk bifurkációs paraméternek.

A következő jelöléseket vezetjük be: 
X=X(t)=(S(t),I(t)) a megoldások vektora, 
F(X)=((Subscript[F, 1])(S,I),(Subscript[F, 2])(S,I))=
=(μ-β S I-μ S+γ I+θ (1-S-I),β S I-(μ+ γ) I+σ β (1-S-I) I)
a redukált SIVS modell jobboldala,
Θ(φ,X)=((Subscript[Θ, 1])(φ,(S,I)),(Subscript[Θ, 2])(φ,(S,I)))=((1-φ)S,I) a vakcinálási impulzusok.
(Subscript[t, i])=i T az impulzusok időpontjai, i=0,1,2,...
Ekkor a modell a következő alakot ölti:

S'=(Subscript[F, 1])(S,I)

I'=(Subscript[F, 2])(S,I)

S((SuperPlus[(Subscript[t, i])]))=(Subscript[Θ, 1])(φ,(S((SuperMinus[(Subscript[t, i])])),I((SuperMinus[(Subscript[t, i])]))))

I((SuperPlus[(Subscript[t, i])]))=(Subscript[Θ, 2])(φ,(S((SuperMinus[(Subscript[t, i])])),I((SuperMinus[(Subscript[t, i])])))).

Látjuk, hogy (Subscript[F, 2])(S,0)≡(Subscript[Θ, 2])(φ,(S,0))≡(Subscript[Θ, 1])(φ,(0,I))≡0 és (Subscript[Θ, 1])(φ,X)≠0, ha S≠0, illetve (Subscript[Θ, 2])(φ,X)≠0, ha I≠0.
Legyen (Subscript[Φ, t]) az ()-()-hez tartozó folyam

X(t)=(Subscript[Φ, t])((Subscript[X, 0])), 0<t≤T,

ahol (Subscript[X, 0])=X(0). Ekkor X(T)=(Subscript[Φ, T])((Subscript[X, 0]))=:Φ((Subscript[X, 0])). A vakcina azonnal kifejti hatását, legyen X((SuperPlus[T])) a vakcina beadása utáni állapota a populációnak. Így X((SuperPlus[T]))=Θ(φ,Φ((Subscript[X, 0]))). 
Legyen Ψ a következőképpen definiált operátor:

Ψ(φ,X)=((Subscript[Ψ, 1])(φ,X),(Subscript[Ψ, 2])(φ,X))=Θ(φ,Φ(X))

és jelöljük (Subscript[D, X])Ψ-vel a Ψ deriváltját X szerint. 
Ekkor (Subscript[X, 0]) egy T-periodikus megoldása ()-()-nek pontosan akkor, ha

Ψ(φ,(Subscript[X, 0]))=(Subscript[X, 0]),

illetve (Subscript[X, 0]) aszimptotikusan stabil, ha (Subscript[D, X])Ψ spektrálsugara kisebb, mint 1.
Már láttuk, hogy a (), () problémának (Subscript[x, 2])=0 feltétellel van egy aszimptotikusan stabil, T-periodikus megoldása, ezt (OverBar[S])(t)-vel jelöltük. Rögzítsük a vakcinálási tényező kivételével a paramétereket, ekkor az (Subscript[R, c])=1 egyenlőség által meghatározott φ értéket jelöljük (Subscript[φ, 0])-val. A továbbiakban a (Subscript[φ, 0])-hoz tartozó T-periodikus megoldással foglalkozunk, jelöljük ezt (OverBar[(Subscript[S, 0])])(t)-vel.
Végül legyen ζ(t)=((OverBar[(Subscript[S, 0])])(t),0) a T-periodikus megoldása ()-()-nek a kétdimenziós térben.

Mivel

(Subscript[D, X])Ψ(φ,X)=(Subscript[D, X])Θ(φ,Φ(X)) (Subscript[D, X])Φ(X)=(({{((∂ (Subscript[Θ, 1]))/(∂S)), ((∂ (Subscript[Θ, 1]))/(∂I))}, {((∂ (Subscript[Θ, 2]))/(∂S)), ((∂ (Subscript[Θ, 2]))/(∂I))}}))(({{((∂ (Subscript[Φ, 1]))/(∂S)), ((∂ (Subscript[Φ, 1]))/(∂I))}, {((∂ (Subscript[Φ, 2]))/(∂S)), ((∂ (Subscript[Φ, 2]))/(∂I))}}))(X),

emiatt (Subscript[X, 0])=((Subscript[x, 0]),0)-ban (ahol (Subscript[x, 0])=(OverBar[(Subscript[S, 0])])(0)=(SuperStar[S]))

(Subscript[D, X])Ψ((Subscript[φ, 0]),(Subscript[X, 0]))=(Subscript[D, X])Θ((Subscript[φ, 0]),Φ((Subscript[X, 0]))) (Subscript[D, X])Φ((Subscript[X, 0]))=(({{1-(Subscript[φ, 0]), 0}, {0, 1}}))(({{((∂ (Subscript[Φ, 1]))/(∂S)), ((∂ (Subscript[Φ, 1]))/(∂I))}, {0, ((∂ (Subscript[Φ, 2]))/(∂I))}}))((Subscript[X, 0])).

A betegségmentes periodikus megoldás stabilitási feltételeit most a következő alakban tudjuk felírni:

|(1-(Subscript[φ, 0]))((∂ (Subscript[Φ, 1])((Subscript[X, 0])))/(∂S)) |<1 és |((∂ (Subscript[Φ, 2])((Subscript[X, 0])))/(∂I)) |<1.

A ()-()-hez tartozó

(d/(dt))((Subscript[D, X])(Subscript[Φ, t])((Subscript[X, 0])))=(Subscript[D, X])F((Subscript[Φ, t])(X)) (Subscript[D, X])(Subscript[Φ, t])(X)

variációs egyenletből a következő feltétellel

(Subscript[D, X])(Subscript[Φ, t])((Subscript[X, 0]))=(Subscript[Id, (ℝ^2)])

kiszámítható ((∂ (Subscript[Φ, 1])((Subscript[X, 0])))/(∂S)), ((∂ (Subscript[Φ, 1])((Subscript[X, 0])))/(∂I)) és ((∂ (Subscript[Φ, 2])((Subscript[X, 0])))/(∂I)). Ezekre a

((∂ (Subscript[Φ, 1])((Subscript[X, 0])))/(∂S))=(e^((∫_0^T)((∂ (Subscript[F, 1])(ζ(r)))/(∂S))dr)),
((∂ (Subscript[Φ, 1])((Subscript[X, 0])))/(∂I))=(∫_0^T)(e^((∫_u^T)((∂ (Subscript[F, 1])(ζ(r)))/(∂S))dr))((∂ (Subscript[F, 1])(ζ(u)))/(∂I))(e^((∫_0^u)((∂ (Subscript[F, 2])(ζ(r)))/(∂I))dr))du ,
((∂ (Subscript[Φ, 2])((Subscript[X, 0])))/(∂S))=0,
((∂ (Subscript[Φ, 2])((Subscript[X, 0])))/(∂I))=(e^((∫_0^T)((∂ (Subscript[F, 2])(ζ(r)))/(∂I))dr))

képletek adódnak. Ezért () második egyenlőtlensége a következő formában is felírható:

(e^((∫_0^T)((∂ (Subscript[F, 2])(ζ(r)))/(∂I))dr)) <1.

A következőkben azt fogjuk vizsgálni, milyen feltételek mellett léphet fel backward bifurkáció a ζ közelében. A számoláshoz a φ és az X változók helyett (OverBar[φ]) és (OverBar[X]) változókra térünk át, melyekre φ=(Subscript[φ, 0])+(OverBar[φ]) és X=(Subscript[X, 0])+(OverBar[X]). 
Az új változókkal a fixpont-probléma a következőképpen is olvasható

N((OverBar[φ]),(OverBar[X]))=0,

ahol

N((OverBar[φ]),(OverBar[X]))=((Subscript[N, 1])((OverBar[φ]),(OverBar[X])),(Subscript[N, 2])((OverBar[φ]),(OverBar[X])))=(Subscript[X, 0])+(OverBar[X])-Ψ((Subscript[φ, 0])+(OverBar[φ]),(Subscript[X, 0])+(OverBar[X])).

Ha ((OverBar[φ]),(OverBar[X])) az N zérushelye, akkor (Subscript[X, 0])+(OverBar[X]) a Ψ((Subscript[φ, 0])+(OverBar[φ]), .) fixpontja. 
Mivel ζ egy T-periodikus, (Subscript[φ, 0])-hoz tartozó betegségmentes megoldása ()-()-nek, ezért ehhez tartozik egy (Subscript[X, 0]) fixpontja Ψ( (Subscript[φ, 0]),.)-nak. Az (OverBar[(Subscript[S, 0])])(t) stabil megoldás, ezért kapjuk, hogy

1-|(1-(Subscript[φ, 0]))((∂ (Subscript[Φ, 1])((Subscript[X, 0])))/(∂S)) |≠0

Jelöljük N deriválját a következő mátrixszal:

(Subscript[D, X])N((OverBar[φ]),(OverBar[X]))=(({{a', b'}, {c', d'}}))

Legyen a'=(Subscript[a, 0]), b'=(Subscript[b, 0]), c'=(Subscript[c, 0]), d'=(Subscript[d, 0]) az ((OverBar[X]);(OverBar[φ]))=((0,0),0) helyen. A mátrix elemeit kiszámolva

a'=1-(((∂ (Subscript[Θ, 1]))/(∂S))((∂ (Subscript[Φ, 1]))/(∂S))+((∂ (Subscript[Θ, 1]))/(∂I))((∂ (Subscript[Φ, 2]))/(∂S)))((Subscript[X, 0])+(OverBar[X]);(Subscript[φ, 0])+(OverBar[φ])),
b'=-(((∂ (Subscript[Θ, 1]))/(∂S))((∂ (Subscript[Φ, 1]))/(∂I))+((∂ (Subscript[Θ, 1]))/(∂I))((∂ (Subscript[Φ, 2]))/(∂I)))((Subscript[X, 0])+(OverBar[X]);(Subscript[φ, 0])+(OverBar[φ])),
c'=-(((∂ (Subscript[Θ, 2]))/(∂S))((∂ (Subscript[Φ, 1]))/(∂S))+((∂ (Subscript[Θ, 2]))/(∂I))((∂ (Subscript[Φ, 2]))/(∂S)))((Subscript[X, 0])+(OverBar[X]);(Subscript[φ, 0])+(OverBar[φ])),
d'=1-(((∂ (Subscript[Θ, 2]))/(∂S))((∂ (Subscript[Φ, 1]))/(∂I))+((∂ (Subscript[Θ, 2]))/(∂I))((∂ (Subscript[Φ, 2]))/(∂I)))((Subscript[X, 0])+(OverBar[X]);(Subscript[φ, 0])+(OverBar[φ])),

és ((OverBar[φ]),(OverBar[X]))=(0,(0,0))-ra, felhasználva Θ(X)-et, illetve, hogy ((∂ (Subscript[Φ, 2]))/(∂S))((Subscript[X, 0]))=0, kapjuk hogy

(Subscript[a, 0])=1-(1-(Subscript[φ, 0]))(((∂ (Subscript[Φ, 1]))/(∂S)))((Subscript[X, 0])),
(Subscript[b, 0])=-(1-(Subscript[φ, 0]))(((∂ (Subscript[Φ, 1]))/(∂I)))((Subscript[X, 0])),
(Subscript[c, 0])=0,
(Subscript[d, 0])=1-(((∂ (Subscript[Φ, 2]))/(∂I)))((Subscript[X, 0])).

Így (Subscript[c, 0])=0 és (Subscript[a, 0])>0, mivel (1-(Subscript[φ, 0]))(((∂ (Subscript[Φ, 1]))/(∂S)))((Subscript[X, 0]))<1 (OverBar[S])(t) stabilitása miatt.

Ahhoz, hogy a betegségmentes periodikus megoldástól különböző periodikus megoldás lépjen fel, a (Subscript[D, X])N(0,(0,0)) Jacobi-mátrix determinánsának szükségszerűen nullának kell lennie. Ebből következik, hogy (Subscript[d, 0])=0. Ez esetünkben pontosan azt jelenti, hogy bifurkáció ott léphet fel, ahol a () éppen 1-gyel egyenlő, vagyis az (Subscript[R, c])=1.

A továbbiakban tegyük fel, hogy (Subscript[d, 0])=0, vagyis (Subscript[R, c])=1. A betegségmentes periodikus megoldástól különböző periodikus megoldás az (Subscript[R, c])=1-ből ágazik el, ami előreágazásnál egy endemikus egyensúly megjelenését jelentheti, hátrafelé ágazásnál pedig éppen a backward bifurkációt okozná. Csak a pozitív megoldásokat vizsgálva akarunk elégséges feltételt adni a bifurkáció létrejöttére.

A következőkben a fixpontegyenlettel ekvivalens () egyenlet analízisét végezzük el. Mivel feltesszük, hogy det((Subscript[D, X])N(0,(0,0)))=0, így nem tudjuk alkalmazni az implicitfügg-vény-tételt az X változó φ függvényeként való megadására. Ahhoz, hogy a problémát kezelni tudjuk, a Ljapunov-Schmidt redukciós eljárást használjuk. Az eljárás lényege, hogy megadjunk egy olyan ekvivalens egyenletrendszert, amelynél már alkalmazható az implicitfüggvény-tétel. Eddig azt tudjuk, hogy N(0,(0,0))=0, legyen (Subscript[D, X])N(0,(0,0))=E egy lineáris leképezés mátrixa.

1. lépés: Felbontjuk a fázisteret az E leképezés mag- és képtere segítségével.
Látjuk, hogy dim Ker(E)=1=co-dim Im(E). Legyen P és Q a Ker(E)-re és Im(E)-re vonatkozó projekció, melyekre:

P((ℝ^2))=Ker(E)=Span{(Subscript[Y, 0])} , ahol (Subscript[Y, 0])=(-((Subscript[b, 0])/(Subscript[a, 0])),1) ,  
Q((ℝ^2))=Im(E)=Span{(1,0)}, (I-Q)((ℝ^2))=Span{(0,1)}.

2. lépés: A felbontást alkalmazzuk az egyenletre - így két egyenletet kapunk, amelyek együttesen ekvivalensek az eredetivel.

Nyilvánvaló, hogy u∈(ℝ^2) esetén u=0⇔Q(u)=0 és (I-Q)(u)=0.  Ebből kapjuk, hogy

N((OverBar[φ]),(OverBar[X]))=0
⇔Q(N((OverBar[φ]),(OverBar[X])))=0 és (I-Q)(N((OverBar[φ]),(OverBar[X])))=0⇔(Subscript[N, 1])((OverBar[φ]),(OverBar[X]))=0 és (Subscript[N, 2])((OverBar[φ]),(OverBar[X]))=0.

Ezen kívül alkalmazzuk az (ℝ^2)=Ker(E)⊕Im(E) felbontást (OverBar[X])-ra, vagyis (OverBar[X])=α (Subscript[Y, 0])+Z, ahol  α∈ℝ, Z=(Subscript[z, 1]) (1,0)=((Subscript[z, 1]),0), (Subscript[z, 1])∈ℝ egyértelmű előállítás.

Így () ekvivalens a következővel:

(Subscript[N, 1])((OverBar[φ]),α (Subscript[Y, 0])+Z)=0,
(Subscript[N, 2])((OverBar[φ]),α (Subscript[Y, 0])+Z)=0.

3. lépés: Megmutatjuk, hogy az egyik egyenlet esetén már alkalmazható az implicitfügg-vény-tétel. 
() első egyenletéből

((∂ (Subscript[N, 1])(0,(0,0)))/(∂(Subscript[z, 1])))=((∂ (Subscript[N, 1])(0,(0,0)))/(∂S))((∂ S)/(∂(Subscript[z, 1])))=(Subscript[a, 0])≠0.

Tehát az implicitfüggvény-tétel szerint létezik egy elegendően kicsi δ>0 és egy egyértelmű, folytonos (SuperStar[Z]) függvény, hogy

(SuperStar[Z])((OverBar[φ]),α)=((SuperStar[(Subscript[z, 1])])((OverBar[φ]),α),0),
				(SuperStar[Z])(0,0)=(0,0),
				(Subscript[N, 1])((OverBar[φ]),α (Subscript[Y, 0])+(SuperStar[Z])((OverBar[φ]),α))=0

minden ((OverBar[φ]),α) -ra, melyre |α|<δ és |(OverBar[φ])|<δ és

((∂ (SuperStar[(Subscript[z, 1])])(0,0))/(∂α))=-(((∂ (Subscript[N, 1])(0,(0,0)))/(∂α))/((∂ (Subscript[N, 1])(0,(0,0)))/(∂S)))=-((((∂ (Subscript[N, 1]))/(∂S))((∂S)/(∂α))+((∂ (Subscript[N, 1]))/(∂I))((∂ I)/(∂α)))/((∂ (Subscript[N, 1]))/(∂S)))=-((((∂ (Subscript[N, 1]))/(∂S))(((-(Subscript[b, 0]))/(Subscript[a, 0])))+((∂ (Subscript[N, 1]))/(∂I))((∂ I)/(∂α)))/((∂ (Subscript[N, 1]))/(∂S)))=-(((((∂ (Subscript[N, 1]))/(∂S))))^(-1))((∂ (Subscript[N, 1]))/(∂I))+((Subscript[b, 0])/(Subscript[a, 0]))=0,

vagyis ((∂ (SuperStar[Z])(0,0))/(∂α))=0.

4. lépés: A tétel eredményeképp kapott megoldást a másik egyenletbe helyettesítjük.

Tehát N((OverBar[φ]),(OverBar[X]))=0 akkor és csakis akkor, ha

f((OverBar[φ]),α)=(Subscript[N, 2])((OverBar[φ]),α (-((Subscript[b, 0])/(Subscript[a, 0])),1)+((SuperStar[(Subscript[z, 1])])((OverBar[φ]),α),0)))=(Subscript[N, 2])((OverBar[φ]),(-((Subscript[b, 0])/(Subscript[a, 0]))α+(SuperStar[(Subscript[z, 1])])((OverBar[φ]),α),α))=0.

Eddig azt tudjuk, hogy f((OverBar[φ]),α)-nak zérushelye a (0,0), ezért szükség van f((OverBar[φ]),α) magasabb rendű deriváltjaira egy olyan i rendig, melyre (D^i)f(0,0)≠0. Ha bevezetjük a η((OverBar[φ]))=(Subscript[φ, 0])+(OverBar[φ]), (Subscript[η, 1])((OverBar[φ]),α)=(Subscript[x, 0])-((Subscript[b, 0])/(Subscript[a, 0]))α+(SuperStar[(Subscript[z, 1])])((OverBar[φ]),α), (Subscript[η, 2])((OverBar[φ]),α)=α új változókat, hosszas számolás után (lásd []) a következőt kapjuk:

(D^2)f(0,0)((((OverBar[φ]),α))^2)=2B α(OverBar[φ])+C (α^2),

ahol

B=(((∂^2) f((OverBar[φ]),α) )/(∂α ∂(OverBar[φ])))=(((∂^2) (Subscript[Φ, 2])((Subscript[η, 1]),(Subscript[η, 2])) )/(∂I ∂S))(((Subscript[Φ, 1])((Subscript[x, 0]),0) )/(Subscript[a, 0])),

C=(((∂^2) f((OverBar[φ]),α) )/(∂(α^2)))=2(((∂^2) (Subscript[Φ, 2])((Subscript[η, 1]),(Subscript[η, 2])) )/(∂S∂I))((Subscript[b, 0])/(Subscript[a, 0]))-(((∂^2) (Subscript[Φ, 2])((Subscript[η, 1]),(Subscript[η, 2])) )/(∂(I^2))).

Ebből következik, hogy

f((OverBar[φ]),α)=(α/2)(OverTilde[f])((OverBar[φ]),α),

ahol

(OverTilde[f])((OverBar[φ]),α)=2 B(OverBar[φ])+C α+(1/α)(Subscript[ℴ, ((OverBar[φ]),α)])((((|(OverBar[φ])|+|α|))^2)),

Továbbá ((∂(OverTilde[f])(0,0))/(∂(OverBar[φ])))=2B (hasonlóan ((∂(OverTilde[f])(0,0))/(∂α))=C). Így, ha B≠0 (hasonlóan, ha C≠0), akkor az implicitfüggvény-tétel értelmében található egy  (OverBar[φ])=ξ(α) (hasonlóan  α=λ((OverBar[φ]))), melyre minden α esetén (hasonlóan (OverBar[φ]) esetén) a 0 környezetében (OverTilde[f])(ξ(α),α)=0 (hasonlóan (OverTilde[f])((OverBar[φ]),λ((OverBar[φ])))=0 ). Ekkor, ha BC≠0, akkor (α/(OverBar[φ]))≃-(C/(2B)). Tehát, ha f((OverBar[φ]),α)=0, akkor (α/(OverBar[φ]))≃-(C/(2B)).

BC=0 nem határozza meg a viselkedést; ekkor az f harmadik deriváltjait is ki kell számolni a sorfejtésben. Ezek alapján a következő tételt fogalmazhatjuk meg:

Tekintsük a (21)-ben definiált  Ψ(φ,X) operátorcsaládot. Amint a φ paraméter áthalad (Subscript[φ, 0])-n, az (Subscript[X, 0]) fixpont egy környezetében egy nemtriviális fixpont jelenik meg, ami a () rendszer endemikus periodikus megoldásának felel meg. Ha BC<0, akkor a bifurkáció szuperkritikus, míg a BC>0 esetben szubkritikus.

Ha φ=(Subscript[φ, 0]), akkor éppen (Subscript[d, 0])=0, azaz (Subscript[R, c])=1. Tehát a bifurkáció pontosan akkor történik, amikor a reprodukciós szám áthalad a kritikus értéken. A bifurkáció irányát az  (α/(OverBar[φ]))≃-(C/(2B)) reláció magyarázza, az α>0 ág tartozik a (3) endemikus periodikus megoldásához. A BC=0 degenerált esetben magasabb rendű deriváltak kiszámítására van szükség.

Mivel a φ növekedése esetén az (Subscript[R, c]) csökken, így a tétel állítása a következőképpen fogalmazható meg: a ((OverBar[φ]),α) síkon lévő szuperkritikus bifurkáció a modellben egy backward bifurkációt jelent. A tételbeli szubkritikus bifurkáció a modellben megjelenő transzkritikus bifurkációt magyarázza. 
A tétel illusztrálásaképp két konkrét esetben megnéztük a BC szorzat előjelét. Mivel a C-ben három, egymásban lévő integrállal van dolgunk, amit a Mathematica is nehezen kezel, így egy egyszerű közelítéssel dolgoztunk. A σ=0.1, γ=0.89, μ=0.01, θ=0.01, T=1, φ=0.174, β=5.32 paraméterek esetén az (Subscript[R, c])=1-ben backward bifurkáció alakul ki; ekkor a BC szorzat értéke -30.4828. Az előbbi paraméterek közül a μ értékét 0.31-re, a β értékét 2.23-ra változtatva már forward bifurkáció lép fel; a BC szorzat ekkor 2684.44, tehát pozitív.

Numerikus vizsgálatok

Folytonos vakcinálási modell

A folytonos modellel kapcsolatban először a fázisképet (. ábra) érdemes megvizsgálni. Az alábbi ábrán a három egyensúlyi helyzetet látunk, amiből az instabil egy nyeregpont. Ez azért is fontos észrevétel, mivel az impulzív rendszer esetén fellépő instabil egyensúlyi helyzetről is úgy sejtjük, hogy nyeregpont. Ez azt a problémát veti fel, hogy a differenciál egyenletben az idő megfordításával nem fogjuk megtalálni az eredeti rendszer instabil egyensúlyi helyzetét, így a bifurkációs diagram kirajzolásához más eszközöket kell használnunk.

A fáziskép (. ábra) kapcsán implementáltunk egy rutint, amely a 0≤S+I≤1 tartományt négyzetekre osztja fel (az átfogó mentén háromszögekre); a négyzetek középpontjaiból indítottunk megoldásokat, és egy nagy t-re rögzítettük az értékpárt (hasonlóan járunk el a bifurkációs diagramnál is). Egy négyzetet kékre színeztünk, ha a betegségmentes egyensúlyhoz konvergált, pirosra, ha az endemikus egyensúlyhoz (az egyensúlyi helyzeteket szintén feltüntettük). Így kirajzoltuk a két stabil periodikus megoldás vonzási tartományát (. ábra). A paraméterek a . ábra paraméterei.

Impulzív vakcinálási modell

Az impulzív rendszer esetén is rajzolható fáziskép. A  . ábrán a differenciálegyenlet jobb oldala által definiált vektormezőt, az impulzus mezőt (részletekért lásd []) és a trajektóriákat ábrázoltuk.

A megoldásgörbék ábrájában a függőleges vonalak a Mathematica függvényábrázoló rutinja miatt jelentek meg: a program a függvény képletéből pontokat generál le, amelyeket ezután szakaszokkal köt össze. A szakadásokat nem tudja kezelni, ezért az ugrásban egy függőleges vonallal folytatja.  A φ-t ismét rögzítettnek tekintjük, ezért a β értéke igazodik a rögzített paraméterekhez. A β mellett az impulzív modellhez tartozó (Subscript[R, c]) értékét is feltüntetjük. A megoldásokat az (((SuperStar[S]),ϵ), ((SuperStar[S]),1-(SuperStar[S])-ϵ)) pontokból indítjuk.

11. ábra.

Folytonos vakcinálási modell

Folytonos vakcinálás esetén láttuk, hogy a bifurkáció vizsgálatához a () másodfokú egyenlet gyökeit kell vizsgálnunk. A következő interaktív ábrában a másodfokú egyenlet megoldásait ábrázoljuk az (Subscript[R, c])(φ) függvényében. Korábban láttuk, hogy (Subscript[R, c])(φ) a β paraméter konstansszorosa, ezért (Subscript[R, c])(φ) és a többi paraméter már egyértelműen meghatározza β-t. Egy megoldórutint írtunk a másodfokú egyenlet megoldására, aminek bemenő paraméterei a modell paraméterei; így minden (Subscript[R, c])(φ)-re már a β is adott volt, így az eljárás kiszámolhatta a gyököket. Az impulzív vizsgálatainknál nagy hasznát vettük ennek az interaktív ábrának, mivel azt reméltük, hogy a folytonos modellel “analóg” az impulzív rendszer, így a folytonos vakcinálás esetén backward bifurkációt okozó paraméteregyüttes az impulzív esetén is ugyanezt a jelenséget idézi elő.

Impulzív vakcinálási modell

Már az első interaktív ábrák azt mutatták, hogy a folytonos modellhez hasonlóan, az impulzív modellben is beállíthatók a paraméterek úgy, hogy a kezdeti értékektől függően két stabil periodikus megoldás alakuljon ki. A folytonos esetben a bifurkációs diagramon ez egy visszahajló görbét eredményezett - természetes volt a célkitűzés, hogy az impulzív esetben is rajzoljunk egy bifurkációs diagramot, amelyen ugyanez a backward jelenség figyelhető meg. Mivel az egyenlő időközönként bekövetkező impulzusokkal ható rend-szereknél a periodikus megoldások stabilitását vizsgáltuk, így első ránézésre nem volt egyértelmű, mit ábrázoljunk a bifurkációs diagramon. Több alternatív lehetőség kínálkozott: az impulzus előtti vagy utáni fertőzők száma, a T-periódusban a fertőzők átlagos száma. Látható, hogy ezek között a szemléltetés szempontjából nem lehet igazán különbséget tenni, ezért a könnyebb implementálhatóság érdekében az impulzus előtti fertőzők száma mellett döntöttük. 
A célunk az volt, hogy egy konkrét esetben felrajzoljunk egy bifurkációs diagramot, amely-nek vízszintes tengelyén az (Subscript[R, c]) változik. A terv az volt, hogy egy adott (Subscript[R, c]) esetén különböző kezdeti értékekből indítunk megoldásokat, és egy nagy (OverBar[t]) időpontban megnézzük, mennyi fertőző van a populációban. Mivel a backward jelenséget akartuk megfigyelni, ezért közvetlenül a betegségmentes periodikus megoldás közeléből indultunk (ez az ((SuperStar[S]),0) pont volt), és a fertőzők számának növelésével egyszercsak elhagytuk a betegségmentes periodikus megoldás attraktivitási tartományát - ha ez bekövetkezett (Subscript[R, c])<1 esetén, akkor a diagramon két pont jelent meg. Tudjuk, hogy (Subscript[R, c])<1 esetén a betegségmentes megoldás aszimptotikusan stabil, így a másik pontnak egy endemikus periodikus megoldásnak kellett lennie. Már az első szimulációk eredményei is azt mutatták, amit szerettünk volna: (Subscript[R, c])<1 esetén is megjelent egy másik stabil ág az I=0 mellett. 
Ezek az eredmények azonban még nem a teljes képet mutatták, ugyanis a bifurkációs diagramon az instabil periodikus megoldásnak is rajta kell lennie. Az instabil rész kirajzolása általában nehéz feladat, erre egy lehetőség, hogy az időben visszafelé haladva tekint-jük a megoldásokat, hogy konvergálnak-e az instabil állapothoz. A folytonos modell fázisképéről azonban megállapítottuk, hogy egy nyeregpontról van szó - emiatt ezt az ötletet el kellett vetni, mivel numerikusan nagyon nehéz lett volna “eltalálni” a stabil sokaságot. Emiatt egy másik megközelítést választottunk, ami a peremérték problémáknál használt “shooting method” elvén alapul. A numerikus vizsgálatokat Wolfram Mathema-tica-val végeztük és a program széleskörben alkalmazható programozási eszközöket kínál listák (vektorok) és mátrixok kezelésére, ezért egy szemléletes módszert alkalmaztunk. Már a korai eredményekből is ki lehetett olvasni, körülbelül hol jelenik meg az első, (Subscript[R, c])=1 előtti endemikus periodikus megoldás (ez egy ((SuperStar[(Subscript[R, c])]),(SuperStar[I])) pont volt). Az [(SuperStar[(Subscript[R, c])]),1] intervallumot  kicsi, egyenlő hosszúságú intervallumokra osztottuk be, és minden (Subscript[R, c]) értékhez ((SuperStar[S]),k (Subscript[δ, k])) alakú pontokból indítottunk megoldásokat (k=1,2,..., a lépésköz a k-val együtt változhatott). Így adott (Subscript[R, c])-hez tartozó párok listájában volt egy (SuperStar[k]) index, amely előtt még (OverBar[t])-ra az ((SuperStar[S]),k (Subscript[δ, k]))-ból indított megoldás ((SuperStar[S]),ϵ)  alakúvá vált (ϵ<(10^(-10))), a (SuperStar[k]) -tól kezdve pedig a második komponens már legalább (SuperStar[I]) nagyságú volt. Ekkor a (SuperStar[k])(Subscript[δ, (SuperStar[k])]) körülbelül azt a pontot adja meg, ahonnantól távolodva már az endemikus periodikus megoldáshoz konvergálunk, előtte pedig az ((SuperStar[S]),0)-hoz. Az első ilyen ábrákon látszott, hogy a megfelelő beosztással ki tudjuk rajzolni az instabil részt. A következőkben az implementálásra vonatkozó legfontosabb információkat írjuk le.
A paraméterekre a következő értékeket rögzítettük: σ=0.1, γ=0.89, μ=0.01, θ=0.01, T=1, φ=0.174 (a rögzített φ-re úgy gondolunk, mintha a kampányban ekkora a ma-ximális rész, amelyet be tudunk oltani). Az (Subscript[R, c])-tengelyen a [0,1.5] intervallumot 0.05 hosszú egységekre osztottuk fel, minden osztópont esetén kezdeti értékek egy listáját írtuk fel és a lista elemeiből indítottuk a megoldásokat; a kapott megoldásokat a [0.9,1] intervallumon a (OverBar[t])=400 időpontban, a többi intervallumon a (OverBar[t])=300 időpontban értékeltük ki (a [0.9,1] intervallumon azért volt szükség későbbi időpontra, mert a rendszer már (10^(-4)) nagyságrendű I esetén sem állt be a betegségmentes periodikus megoldásba). A kezdeti értékek listájában a [0.55,1]  intervallum kivételével δ=0.05 lépésközt használtunk, a kezdeti értékek 0 és 0.9 közöttiek voltak (a fenti paraméterek esetén (SuperStar[S])=0.0859, nagyobb számot az I-re nem volt érdemes venni). Mivel az instabil ág a [0.55,1]  intervallumon volt megfigyelhető, ezért itt a kezdeti értékek között kisebb lépésközzel dolgoztunk; a program felgyorsítása érdekében folyamatosan figyeltük az eddigi eredményeinket, és kicsi lépésközt csak az instabil ágig alkalmaztunk, ezen érték felett már jóval nagyobb volt a kezdeti értékek közötti különbség. Minden intervallumon feljegyeztük a feltételezett instabil periodikus megoldásig alkalmazott lépésközt és ezzel súlyoztuk azt az indexet, amelytől kezdve a megoldáspárok már egy endemikus periodikus megoldást mutattak. A kellően finom beosztással egy folytonos görbe alakult ki. A piros színű pontok a stabil periodikus megoldásokat, a kék színű pontok az instabil periodikus megoldásokat mutatják.

Járványügyi kérdések

A következőkben néhány járványügyi szempontból fontos gyakorlati példát és numerikus vizsgálatot mutatunk be.

Már az 1.2. részben is láttuk, hogy fontos tudnunk, milyen bifurkáció lép fel a modellben. Ennek gyakorlati jelentőségét érzékelteti, hogy számos létező betegségre felállított modellben felbukkant backward bifurkáció. A következő cikkek a MathSciNet adatbázisából valók, a “backward bifurkáció” keresésre adott találatok közül (a teljesség igénye nélkül):

- Ngonghala, Calistus N. et al, Periodic oscillations and backward bifurcation in a model for the dynamics of malaria transmission
- Garba, S.M. et al, Cross-immunity induced backward bifurcation for a model of transmission dynamics of two strains of influenza
- Podder, C.N. et al, Risk-induced backward bifurcation in HSV-2 transmission dynamics
- Qesmi, R. et al, A hepatitis B and C virus model with age since infection that exhibits backward bifurcation
- Sharomi, O. et al, Re-infection-induced backward bifurcation in the transmission dynamics of Chlamydia trachomatis
- Garba, S.M. et al, Backward bifurcations in dengue transmission dynamics
- Wang, X., Backward bifurcation in a mathematical model for tuberculosis with loss of immunity

Egy konkrét cikk bevezetőjéből ragadnék ki egy részletet:

“Hepatitis-B virus infection is globally ubiquitous, but its distribution is very heterogeneous, with prevalence of serological markers in various nations ranging from less than 1% to more than 90%. We propose an explanation for this diversity using a mathematical model of hepatitis B virus transmission dynamics that shows, for the first time, ‘catastrophic’ behavior using realistic epidemiological processes and parameters.”
Medley, F.M. et al, Hepatitis-B virus endemicity: heterogeneity, catastrophic dynamics and control, Nature Medicine, 2001.

A részlet arra világít rá, hogy backward bifurkációval magyarázták a hepatitis-B vírus nagy varianciájú endemicitását; vagyis egy populáción belül két, egymáshoz közel lévő és azonos egészségügyi protokollt alkalmazó területen teljesen más lehet a fertőzöttek aránya. Ennek lehet az oka, hogy az egyik területen a betegek száma még nem érte el az endemikus egyensúlyt, így a gyors reakció hatására (Subscript[R, c]) kicsivel egy alá csökkentése megoldotta a problémát, míg a másik területen már beállt az endemikus egyensúly, így jóval hatékonyabb protokollra lenne szükség, ami az (Subscript[R, c]) értékét még jobban mérsékli.

A valóságban a paramétereink többsége rögzített, mivel a populációra, a vakcinára vagy magára a betegséget okozó vírusra, baktériumra jellemző. Járványügyi szempontból impulzív vakcinálás esetén a φ-t és a T-t tudjuk befolyásolni, így érdemes az (Subscript[R, c])-t úgy is vizsgálnunk, mint a φ és a T függvénye. A . ábrán az (Subscript[R, c])(φ,T) függvény néhány szintvonalát láthatjuk (vastag vonallal az (Subscript[R, c])=1-hez tartozó szintvonal):

Az ábra nem mutatja, de az (Subscript[R, c])(φ,T)-nek szakadása van az origóban: (lim_(T->0)) (Subscript[R, c])(c T,T)=((β (θ+μ+c σ))/((γ+μ) (c+θ+μ))) (egyébként kiszámolhatók a következők: (lim_(φ->0)) (Subscript[R, c])(φ,T)=((β )/(γ+μ)) és (lim_(T->0)) (Subscript[R, c])(φ,T)=((β σ)/(γ+μ))).  
Mivel a φ csökkentése és/vagy a T növelése rontja a kampány sikerességét, ezért az (Subscript[R, c])=1 szintvonal felett olyan (φ,T) párok vannak, amelyek (Subscript[R, c])>1 értékekhez tartoznak. Ezt mutatja a . ábra, ahol piros a “veszélyes tartomány”, kék a “biztonságos tartomány”. Forward bifurkáció esetén ennyit elég tudnunk, hogy a járványt elfojtsuk, hiszen akármennyi fertőző is van a populációban az (Subscript[R, c]) 1-re csökkentése (azaz egy olyan (φ,T) kiválasztása, amelyre (Subscript[R, c])(φ,T)=1) betegség kihalását okozza.

Célszerű megvizsgálni, hogy egy backward bifurkációt okozó paraméteregyüttesnél a (φ,T) milyen értékeire kapunk még backward bifurkációt, vagy másképp: ha nagyon sok fertőző van a populációban és backward bifurkáció lépett fel, hogy fog kinézni a “biztonságos tartomány”? Impulzív esetben a bifurkációs tételben szereplő szükséges és elegendő feltétel nagyon bonyolult, így nem világos a φ-vel és T-vel való kapcsolata. A vizsgálatainkhoz a korábban is tekintett (σ=0.1, γ=0.89, μ=0.01, θ=0.01, β=5.32) paraméteregyüttest tekintettük és a következő interaktív ábrát implementáltuk (itt τ-val paramétereztük a szintvonalakat).

A tapasztalataink meglepőek voltak: bármilyen (Subscript[R, c])=1-et teljesítő (φ,T) párra nagy fertőző szám esetén a rendszer egy endemikus egyensúlyba állt be. Sőt, ahogy (Subscript[R, c])-t csökkentettük, a szintvonalakon ugyanaz történt: ha (Subscript[R, c]) nem volt elég kicsi, a rendszer egy endemikus egyensúlyba állt be, viszont elegendően alacsony (Subscript[R, c]) esetén már a betegség kihalt a szintvonal bármelyik (φ,T) párját is tekintettük.

Összegzés

1. A betegségmentes periodikus megoldást explicit meghatároztuk és megmutattuk, hogy globálisan aszimptotikusan stabil a betegségmentes altéren.

2. A betegségmentes periodikus megoldásra explicit stabilitási feltételt fogalmaztunk meg a teljes fázistéren; ebből definiáltunk egy (Subscript[R, c]) reprodukciós számot.

3. (Subscript[R, c])<1 esetén egy újabb feltétel mellett a betegségmentes periodikus megoldás globális aszimptotikus stabilitása teljesül.

4. (Subscript[R, c])>1 esetén a betegség perzisztens lesz.

5. Az (Subscript[R, c])=1 kritikus pontban fellépő szubkritikus bifurkációra elégséges feltételt adtunk, így endemikus periodikus megoldások létezését bizonyítottuk (Subscript[R, c])<1 esetekben is.

A következő, backward bifurkációt mutató ábrán szemléltetjük eddigi eredményeinket.

Zs. Vizi, Szubkritikus bifurkáció és globális dinamika impulzív imperfekt vakcinálási modellben, OTDK dolgozat, 2013.

F. Brauer, Backward bifurcations in simple vaccination models, J. Math. Anal. Appl. 298 (2004) 418-431

J. Karsai, Impulzív jelenségek modelljei, Typotex Kiadó, 2002

A. Lakmeche, O. Arino, Bifurcation of non trivial periodic solutions of impulsive differential equations arising chemotherpeutic treatment, Dyn. Con., Disc. Imp. Sys. 7 (2000) 265-287

A. Lakmeche, O. Arino, Nonlinear mathematical model of pulsed-therapy of heterogeneous tumors, Nonlin. Anal: Real World Appl. 2 (2001) 455-465

Y. Li, J. Cui, The effect of constant and pulse vaccination on SIS epidemic models incorporating media coverage, Commun Nonlinear Sci. Num. Sim. 14 (2009) 2353-2365

P. Georgescu, H. Zhang, L. Chen, Bifurcation of nontrivial periodic solutions for an impulsively controlled pest management model, Appl. Math. Comp. 202 (2008) 675-687

B. Shulgin, L. Stone, Z. Agur, Theoretical Examination of the Pulse Vaccination Policy in the SIR Epidemic model, Math. Comp. Mod. 31. (2000), 207-215

M. Golubitsky, D. G. Schaeffer, Singularities and Groups in Bifurcation Theory Vol. I., Springer-Verlag, 1985

V. Lakshmikantham, D.D. Bainov, P.S. Simeonov, Theory of Impulsive Differential Equations, World Scientific Publishing, 1989.

Y. A. Kuznetsov, O. Diekmann, W.J. Beyn, Lecture Notes on Dynamical Systems (egyetemi jegyzet), 2012

M, Szűcs., Járványterjedési modellek számítógépes vizsgálata, Diplomamunka, 2009

F. Brauer, C. Castillo-Chavez, Mathematical Models in Population Biology and Epidemiology, 2nd Edition Springer Science, 2012

http://www.wolfram.com/

B. Buonomo, D. Lacitignola, On the backward bifurcation of a vaccination model woth nonlinear incidence, Nonlin. Anal. Mod. Contr., Vol. 16. (2011) 30-46

B.H. Singer, D.E. Kirschner, Influence of backward bifurcation on interpretation of (Subscript[R, 0]) in a model of epidemic tuberculosis with reinfection, Math. Biosc. Eng., Vol. 1., (2004) 81-93

A numerikus vizsgálatok során a Karsai János által fejlesztett Impulse programcsomagot használtuk; ezt a [] könyv számára készítette, amelyben általánosabb impulzív rendszerekkel és interaktív példákkal is foglalkozik. A következőkben a könyv alapján leírjuk az IDESolve-val kapcsolatos legfontosabb tudnivalókat.

Az IDESolve program rögzített időpillanatokkal rendelkező impulzív rendszerek közelítő megoldását végzi több kezdeti értékből kiindulva. A program számára a rendszert az

x'=f(t,x),ha t≠(Subscript[t, i]) ,
x((SuperPlus[(Subscript[t, i])]))=I((Subscript[t, i]),x((SuperMinus[(Subscript[t, i])]))); i=1,2,...

alakban adjuk meg (ez eltér az 1.3 fejezetben használt jelöléstől). A rendszer megoldásához meg kell adni az egyenletet, az impuzusok (Subscript[t, k]) időpillanatait, az impulzusfüggvényt, a kezdeti értékeket. Két impulzus között a differenciálegyenletet a beépített NDSolve program oldja meg (lásd a differenciálegyenletek számítógépes vizsgálatának módszereivel foglalkozó fejezetet). Az impulzusok (10^(-10)) idő alatt hatnak. A program eredménye InterpolatingFunction objektumokból álló lista lesz.

A program indítása az alábbi módon történik:

IDESolve[xdotlist, varlist, timelist, Impulse, x0list,{t, t0, t1}, option]

Az egyes paraméterek jelentése:

({{varlist, {x1,x2,..., xn},függő változók listája}, {xdotlist, {f1[t,x1,...,xn],...,fn[t,x1,...,xn]},a differenciálegyenlet jobb oldala}, {timelist
Impulse, az impulzusok időpontjai:{t1,t2,...}
impulzus–függvény}, {x0lst, {{x10,...,xn0},...,{x1M,...,xnM}}kezdőértékek listája,M–db n–dimenziós vektor}, {t, független változó}, {[t0,t1], a megoldás intervalluma}, {option, az NDSolve opciói}})

Az Impulse függvény megadása az I((Subscript[t, k]),x) jelöléshez illeszkedik:

Impulse[i_,timelist_,xlist_]:= ...

ahol

({{i, az adott impulzus pillanat indexe timelist–ben}, {timelist, ugyanaz,mint fent}, {xlist, a megoldás értéke a timelist[[i]]helyen}})

Alulírott Vizi Zsolt kijelentem, hogy a dolgozatban foglaltak saját munkám eredményei, és csak a hivatkozott forrásokat használtam fel. Tudomásul veszem, hogy dolgozatomat a Szegedi Tudományegyetem könyvtárában a kölcsönözhető könyvek között helyezik el, és az interneten is nyilvánosságra hozhatják.


2013. május 15.

....................................................

Vizi Zsolt