Quit[];

Járvány terjedése térben

Farkas Miklós: Dynamical models in biology (Academic Press 2001), 73. old.

A modell

eq={(∂_t)S[t,x]==-r S[t,x]I[t,x]+d ΔS[t,x],(∂_t)I[t,x]==r S[t,x]I[t,x]-a I[t,x]+d ΔI[t,x]};

(alkalmas kezdeti és peremfeltételekkel),

Jelölések

ahol:

S(t,x) a veszélyeztetett (susceptible) vagy fertőzésnek kitett populáció térbeli sűrűsége a t időpontban

I(t,x) a fertőzőtt (infective) populáció térbeli sűrűsége a t időpontban  a t időpontban

r>0 a fertőzési ráta (egy beteg egységnyi idő alatt ennyi százalékkal csökkenti az egészségesek számát)

a>0 a fertőzöttek mortalitási rátája (pontosabban a fertőzők száma egységnyi idő alatt ennyi százalékkal csökken, ahol a csökkenés oka csak elhalálozás lehet, pl. pestis)

(1/a) a fertőzöttek átlagos élettartama

(r/a) kontakt ráta

d>0 az (azonosnak feltételezett) diffúziós együttható

Δ:=((∂^2)/(∂(x^2)))+((∂^2)/(∂(y^2))) (pl. síkon)

Egyszerűsítő feltevések

A fenti modell a következő szituációt írja le:

a tér egydimenziós, a számegyenes pozitív fele (vagyis végtelen)

kiindulunk egy pontból, és nézzük valamilyen irányban a terjedést

S(0,x)=:(Subscript[S, 0])>0,I(0,x)=0

Az egyszerűsített modell dimenziómentes alakja

eqdim={D[h[τ,y],τ]==- h[τ,y]v[τ,y]+ D[h[τ,y],{y,2}],D[v[τ,y],τ]==h[τ,y]v[τ,y]-b v[τ,y]+ D[v[τ,y],{y,2}],
h[0,y]==1,v[0,y]==0};

ahol

τ:=r (Subscript[S, 0])t

y:=(Sqrt[((r (Subscript[S, 0]))/d)])x

h[τ,y]:=((S[t,x])/(Subscript[S, 0]))

v[τ,y]:=((I[t,x])/(Subscript[S, 0]))

b:=(a/(r (Subscript[S, 0])))

Térben homogén egyensúlyi helyzetek

(Subscript[E, 0]):=(0,0) (érdektelen)

E((Subscript[h, 0])):=((Subscript[h, 0]),0) ahol (Subscript[h, 0])>0. (Subscript[h, 0])=1 felel meg az ((Subscript[S, 0]),0)  kezdeti feltételnek

Haladóhullám-megoldások

Ilyen alakú megoldásokat keresünk:

h(τ,y)=H(y-cτ), v(τ,y)=V(y-cτ)

z:=y-cτ

Jelentősége, hogy a következő ábrán szereplő egyenesek mentén ugyanaz h és v értéke. Ezt úgy kell érteni, hogy például a (Subscript[τ, 1]) időpontban az (Subscript[y, 1]) helyen ugyanaz, mint a (Subscript[τ, 2]) időpontban az (Subscript[y, 2]) helyen. Sebessége, amivel terjed, számítható:

(Subscript[y, 2])-c (Subscript[τ, 2])=(Subscript[y, 1])-c (Subscript[τ, 1])=:(Subscript[z, 1])

c =(((Subscript[y, 2])-(Subscript[y, 1]))/((Subscript[τ, 2]) -(Subscript[τ, 1])))

ahol c a járvány terjedési sebessége.

Plot[{y-3,y-1,y+2},{y,0,5},AspectRatio->Automatic,PlotRange->{-1,4.5},AxesLabel->{τ,y},PlotStyle->Directive[Thick],Epilog->{Dashed,Line[{{1.5,0},{1.5,0.5},{0,0.5}}],Text["(Subscript[τ, 1])",{1.5,-0.1}],Text["(Subscript[y, 1])",{0.2,0.6}],Line[{{3.5,0},{3.5,2.5},{0,2.5}}],Text["(Subscript[τ, 2])",{3.5,-0.1}],Text["(Subscript[y, 2])",{0.2,2.6}],Text["(Subscript[z, 1])",{4.2,3.5}]}]

A haladó hullámok differenciálegyenlete

H''+c H'-H V=0

V''+c V'+(H-b) V=0

Az ((Subscript[S, 0]),0) kezdeti feltételnek megfelelő egyensúlyi helyzet (1,0).

(A h tengely pozitív felén lévő összes pont egyensúlyi helyzet.)

Cauchy-féle normál alak

Alkalmazva az átviteli elvet:

(Subscript[(OverDot[x]), 1])=(Subscript[x, 2])
(Subscript[(OverDot[x]), 2])=-c (Subscript[x, 2])+(Subscript[x, 1])(Subscript[x, 3])
(Subscript[(OverDot[x]), 3])=(Subscript[x, 4])
(Subscript[(OverDot[x]), 4])=-c (Subscript[x, 4])+(b-(Subscript[x, 1]))(Subscript[x, 3])

(Az érdekes) egyensúlyi helyzet: (1,0,0,0). Itt a Jacobi-mátrix:

(({{0, 1, 0, 0}, {0, -c, 1, 0}, {0, 0, 0, 1}, {0, 0, b-1, -c}}))

CharacteristicPolynomial[(({{0, 1, 0, 0}, {0, -c, 1, 0}, {0, 0, 0, 1}, {0, 0, b-1, -c}})),λ]//Simplify

A 0 sajátérték arra utal, hogy a v=0 egyenes minden pontja egyensúlyi helyzet, így a h tengely centrumsokaság.

Stabilitás

A stabilitás szükséges feltétele: c>2(Sqrt[1-b])>0.

Pontosan ekkor teljesül ugyanis, hogy a másodfokú tényező gyökei negatív valós számok. (A komplex sajátértékek esetén a változók negatív értékeket vehetnek föl, ezért - is - ezt az esetet kizárjuk.)

A stabilitás feltételének jelentése az eredeti változókban: c>2(Sqrt[1-(a/(r (Subscript[S, 0])))])>0⟺(Subscript[S, 0])>(a/r).

Ez a járvány – pontszerű modellből jól ismert – feltétele.

A haladó hullám minimális sebessége

(Subscript[c, min])=2(Sqrt[1-(a/(r (Subscript[S, 0])))])

Az eredeti koordinátákkal:

(Subscript[c, emin])=(Sqrt[d r (Subscript[S, 0])])2(Sqrt[1-(a/(r (Subscript[S, 0])))])=2(Sqrt[d (r (Subscript[S, 0])-a)])

Ha a,r,(Subscript[S, 0]),d ismert, akkor ez számítható.

Alkalmazás a pestisre (Farkas (2001), 76. old.)

A pestis esetére:

(1/a)∼0.033 év (kb. ennyi idő alatt belehalnak) ⟹ a=30 ((év)^(-1))

r∼0.4 ((mérföld)^2)/év

(Subscript[S, 0])∼100/((mérföld)^2)

b∼0.75

d∼4000 ((mérföld)^2)/év

param={a->30,r->0.4,S0->100,b->0.75,d->4000};

{S0,(a/r)}/.param

Tlimit=35;

nds=NDSolve[{H''[z]+c H'[z]-H[z]V[z]==0,V''[z]+c V'[z]+(H[z]-b)V[z]==0,H[0]==0.75,H'[0]==.09,V[0]==-0.04,V'[0]==.08}/.{c->1,b->0.75},{H[z],V[z]},{z,0,Tlimit}];
ParametricPlot[{H[z],V[z]}/.nds,{z,0.69,Tlimit},AxesLabel-> {h,v},AspectRatio->1,AxesOrigin->{0.7,0},Epilog->Arrow[{{0.838,0.0325},{0.828,0.03}}]]
GraphicsRow[{Plot[H[z]/.nds,{z,0.69,Tlimit},AxesLabel->{"z","h(z)"},AxesOrigin->{0,0}],Plot[V[z]/.nds,{z,0.69,Tlimit},AxesLabel->{"z","v(z)"},PlotRange->All,AxesOrigin->{0,0}]}]

A haladó hullám képe

Plot3D[Evaluate[V[z]/.nds/.z->y-τ],{τ,0,1},{y,τ,6},AxesLabel->{τ,y,v},PlotPoints->50,BoxRatios->{1,4,1}]

Hivatkozások

Farkas, M.: Dynamical models in biology, Academic Press, New York, 2001.