Quit[];

Kolmogorov-rendszerek 2.

Egy zsákmányfajból és két ragadozófajból álló modell (m=1,n=2)

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

A modell

ClearAll[eq,i,j,m,n,ini];eqpray[i_,m_,n_]:=(Subscript[x, i])'[t]== (Subscript[x, i])[t](Subscript[F, i])@@Join[Array[(Subscript[x, #])[t]&,m],Array[(Subscript[y, #])[t]&,n]];
eqpred[j_,m_,n_]:=(Subscript[y, j])'[t]== (Subscript[y, j])[t](Subscript[G, j])@@Join[Array[(Subscript[x, #])[t]&,m],Array[(Subscript[y, #])[t]&,n]];
inipray[i_]:=(Subscript[x, i])[0]==(Subscript[x0, i]);inipred[j_]:=(Subscript[y, j])[0]==(Subscript[y0, j]);

Például 2 ragadozó és 3 zsákmányfaj esetén az egyenletek:

Join[Table[eqpray[i,2,3],{i,2}],Table[eqpred[j,2,3],{j,3}]]//TableForm

és a kezdeti feltételek:

Join[Table[inipray[i],{i,2}],Table[inipred[j],{j,3}]]//TableForm

{eqpray[#,2,3],eqpred[#,2,3],inipray[#]}&/@Range[2]

Tömören:

kep[m_,n_]:=Flatten[Reverse[
Partition[Sort[Flatten[Join[{eqpray[#,m,n],inipray[#]}&/@Range[m],{eqpred[#,m,n],inipred[#]}&/@Range[n]]]],m+n]]]

kep[2,3]

Még tömörebben:

ivp[m_,n_]:=Join[Array[eqpray[#,m,n]&,m],Array[eqpred[#,m,n]&,n],
Array[inipray[#]&,m],Array[inipred[#]&,n]]

ivp[2,3]

Ennek speciális esetével fogunk foglalkozni:

ivp[1,2]

Jelölések

ahol:

(Subscript[x, 1])(t) a zsákmányfaj létszáma a t időpontban

(Subscript[y, j])(t) a j-edik ragadozófaj létszáma a t időpontban (j=1,2)

F a zsákmányfaj létszámának változási rátája (F∈(𝒞^1)((ℝ^3),ℝ))

(Subscript[G, j]) a j-edik  ragadozófaj létszámának változási rátája ((Subscript[G, j])∈(𝒞^1)((ℝ^3)(Subscript[,, i])ℝ);   j=1,2)

Az F és a (Subscript[G, j]) függvényekre vonatkozó természetes feltételek  (j=1,2)

F(0,0)>0, azaz a zsákmányfaj önállóan fennmarad

∃ K>0∀x>0&x≠K: (x-K)F(x,0)<0, azaz a zsákmányfaj nem növekedhet korlátlanul

(Subscript[D, 1+j])F(x,y)<0, (ahol (Subscript[D, 1+j]) a j-edik  ragadozófaj létszáma (vagyis (Subscript[y, j])) szerinti deriválást jelöli), azaz a j-edik  ragadozófaj létszámának növekedése csökkenti a zsákmányfaj változási rátáját

(Subscript[G, j])(0,y)<0, azaz  a j-edik  ragadozófaj létszáma táplálék hiányában csökken

(Subscript[D, 1])(Subscript[G, j])(x,y)>0, (ahol (Subscript[D, 1]) a zsákmányfaj létszáma (azaz x) szerinti deriválást jelöli), azaz a zsákmányfaj létszámának növekedése növeli a j-edik  ragadozófaj változási rátáját

(Subscript[D, 1+j])(Subscript[G, j])(x,y)≤0, azaz a j-edik  ragadozófajnál lehet belső (intraspecifikus) verseny

ϕ[x_,y1_,y2_,K_,a1_,a2_,μ1_,μ2_,r_]:=r(1-(x/K))-((μ1 y1)/(a1+x))-((μ2 y2)/(a2+x));
(Subscript[F, 1])[x_,y1_,y2_]:=ϕ[x,y1,y2,K,a1,a2,μ1,μ2,r];
γ1[x_,y1_,a1_,μ1_,d1_]:=((μ1 x)/(a1+x))-d1;
γ2[x_,y2_,a2_,μ2_,d2_]:=((μ2 x)/(a2+x))-d2;
(Subscript[G, 1])[x_,y1_,y2_]:=γ1[x,y1,a1,μ1,d1];(Subscript[G, 2])[x_,y1_,y2_]:=γ2[x,y2,a2,μ2,d2];

{(Subscript[F, 1])[x,y1,y2],(Subscript[G, 1])[x,y1,y2],(Subscript[G, 2])[x,y1,y2]}

akkor az egyenletek és a kezdeti feltételek:

ivp[1,2]

Vegyük észre, hogy (Subscript[γ, 1]) és (Subscript[γ, 2]) (és így (Subscript[G, 1]) és (Subscript[G, 2]) is) csak x-től függ.

Jelölések

ahol:

r(1-((Subscript[x, 1])/K)) a zsákmányfaj logisztikus növekedési rátája, és

r>0 a belső növekedési ráta (azaz születési és a halálozási ráta különbsége)

K>0 a környezet fenntartóképessége a zsákmányfaj szempontjából

(((Subscript[μ, j]) (Subscript[x, 1])(Subscript[y, j]))/((Subscript[a, j])+(Subscript[x, 1]))) az úgynevezett Holling-féle funkcionális válasz, és

(Subscript[μ, j])>0 a j-edik ragadozófaj maximális növekedési rátája

(Subscript[a, j])>0 a j-edik  ragadozófaj félszaturációs állandója, ugyanis x = (Subscript[a, j]) esetén lesz a j-edik ragadozó növekedési rátája a maximális érték fele

a ragadozó létszámával arányosan csökkenti a zsákmányfaj növekedési rátáját

(Subscript[d, j])>0 a j-edik ragadozófaj mortalitása

Az Allee-hatás zónája

A zsákmány nullklína-felülete

Solve[ϕ[x,y1,y2,K,a1,a2,μ1,μ2,r]==0,y1][[1]]

%/.{y2->0}

Az alább használandó tipikus paraméterek:

{r->3, μ1->7/10 ,μ2->5/10,a1->1/6,a2->3/2,d1->6/10,d2->2/10,K->3,β1->1/10,β2->3/10,λ->1};

ContourPlot[Evaluate[ϕ[x,y1,0,K,a1,a2,μ1,μ2,r]/.{K->3,a1->1/6,a2->3/2, μ1->7/10 ,μ2->5/10,r->3}]==0,{x,0,4},{y1,0,4},FrameLabel->{Style[x,Italic,16],Style["(Subscript[y, 1])",Italic,16]},PlotLabel->"Ez a nullavonal-felület metszete,\n amikor (Subscript[y, 2])=0.\n"]

ContourPlot3D[Evaluate[ϕ[x,y1,y2,K,a1,a2,μ1,μ2,r]/.{K->3,a1->1/6,a2->3/2, μ1->7/10 ,μ2->5/10,r->3}]==0,{y1,0,4},{y2,0,4},{x,0,4},AxesLabel->{Style["(Subscript[y, 1])",Italic,16],Style["(Subscript[y, 2])",Italic,16],Style[x,Italic,16]},PlotLabel->"Ez a hagymaalakú nullavonal-felület.\n"]

Solve[γ1[x,y1,a1,μ1,d1]==0,x][[1]]

Solve[γ2[x,y2,a2,μ2,d2]==0,x][[1]]

Pozitív stacionárius pontot csak akkor fogunk kapni, amikor a két ragadozófaj ugyanannál a zsákmánymennyiségnél kezd el növekedni, ha  μi>di (i,1,2), ami egyben a ragadozók túlélésének feltétele, valamint ha a fenti két érték megegyezik.

Milyen μ2 értéknél fordul ez elő?

ClearAll[legyen];
legyen[a1_,a2_,μ1_,d1_,d2_]:=First@ Solve[((a1 d1)/(d1-μ1))==((a2 d2)/(d2-μ2)),μ2];
legyen[a1,a2,μ1,d1,d2]

És valóban:

{((a1 d1)/(d1-μ1)),((a2 d2)/(d2-μ2))}/.legyen[a1,a2,μ1,d1,d2]//Simplify

Konkrét paramétereinkkel:

%/.{a1->1/6,a2->3/2, μ1->7/10 ,μ2->5/10,d1->6/10,d2->2/10}

Hogyan függ μ2 értéke d2-étől, ha a többi paraméter értékét rögzítjük?

legyen[a1,a2,μ1,d1,d2]/.{a1->1/6,a2->3/2, μ1->7/10 ,d1->6/10}

Most számoljuk ki a stacionárius pontokat:

ClearAll[stac];
stac[y1_,K_,a1_,a2_,μ1_,r_,d1_,d2_]:=FullSimplify@First@Solve[{ϕ[x,y1,y2,K,a1,a2,μ1,μ2,r]==0,γ1[x,y1,a1,μ1,d1]==0,
γ2[x,y2,a2,μ2,d2]==0}/.legyen[a1,a2,μ1,d1,d2],{x,y2}]

stac[y1,K,a1,a2,μ1,r,d1,d2]

Ezek a pontok tehát egy egyenesen vannak, aminek egyenletrendszere kedvenc paramétereink mellett:

stac[y1,K,a1,a2,μ1,r,d1,d2]/.{a1->1/6,a2->3/2, μ1->7/10 ,d1->6/10,d2->2/10,r->3}//Expand

Ezeknél a konkrét paraméterértékeknél az egyenes minden K mellett benne halad az x=1 síkban.

ClearAll[pp];pp[d2_]:=ParametricPlot3D[Evaluate[{1,y1,y2/.stac[y1,K,a1,a2,μ1,r,d1,d2]/.{K->3,a1->1/6,a2->3/2, μ1->7/10 ,d1->6/10,r->3}}],{y1,0,3},PlotStyle->Directive[{Thick,Yellow}],AxesLabel->{Style[x,Italic,16],Style["(Subscript[y, 1])",Italic,16],Style["(Subscript[y, 2])",Italic,16]},PlotLabel->"A stacionárius pontok halmaza\naz x = 1 síkban haladó egyenes,\n K = 3 és (Subscript[d, 2]) = "<>ToString[TraditionalForm[d2]]]

GraphicsColumn[{pp[1/6],pp[1/5],pp[1/2]},ImageSize->400]

A fenti ábrát úgy kell érteni, hogy minden d2-höz meghatározzuk azt a μ2-t, amelyre a két ragadozó nullklína-sík egybeesik.

Akár folytonosan is változtathatjuk d2 értékét.

Manipulate[pp[d2],{d2,1/10,2,1/10},SaveDefinitions->True]

Mivel

μ2/.legyen[a1,a2,μ1,d1,d2]/.{a1->1/6,a2->3/2,μ1->7/10,r->3,d1->3/5}

ϕ[x,y1,y2,K,a1,a2,μ1,5d2/2,r]/.{a1->1/6,a2->3/2,μ1->7/10,r->3}//Simplify

ClearAll[sh];
sh[d2_,K_,xig_,y1ig_,y2ig_,plp_:2,opts___]:=ContourPlot3D[Evaluate[ϕ[x,y1,y2,K,a1,a2,μ1,5d2/2,r]/.{a1->1/6,a2->3/2,μ1->7/10,r->3}]==0,{x,0,xig},{y1,0,y1ig},{y2,0,y2ig},AxesLabel->{x,(Subscript[y, 1]),(Subscript[y, 2])},PlotPoints->plp,PlotLabel-> "Hagymán elhelyezkedő stacionárius pontok",opts]

Az alábbi manipulációnál lépésenként érdemes változtatni a paramétereket, mert lassan készül a rajz.

Manipulate[Show[sh[d2,3,10,5,20,4],pp[d2]],{d2,1/10,2,1/10},SaveDefinitions->True]

Zipzár-bifurkáció

Ha a paraméterek teljesítik az egyensúlyi helyzet létezésének fent már említett

λ:=((a1 d1)/(μ1-d1))=((a2 d2)/(μ2-d2))>0

szükséges feltételét, akkor az egyenletrendszer így alakítható át:

ujivp=FullSimplify[ivp[1,2][[{2,3}]]/.{μ1->β1+d1,μ2->β2+d2}]

ered2=Simplify[ujivp/.{d1->λ β1/a1, d2->λ β2/a2}]

equ=Prepend[ered2,ivp[1,2][[1]]]

Az (1) feltétel azt jelenti, hogy mindkét ragadozó ugyanannál a zsákmánymennyiségnél kezd el növekedni. Például a következő paraméterértékeknél a feltétel úgy teljesül, hogy λ értéke 1:

{((a1 d1)/(μ1-d1)),((a2 d2)/(μ2-d2))}/.{ μ1->7/10,μ2->5/10, a1->1/6, a2->3/2,d1->3/5,d2->1/5}

A { μ1->7/10,μ2->5/10, a1->1/6, a2->3/2,d1->3/5,d2->1/5} paraméterkombináció olyan, hogy  (Subscript[y, 1]) az úgynevezett K-stratéga, míg  (Subscript[y, 2])  az úgynevezett r-stratéga:

{a1<a2,μ1/d1<μ2/d2}/.{μ1->7/10,μ2->5/10, a1->1/6, a2->3/2,d1->3/5,d2->1/5}

Az r-stratéga a túlélését bőséges táplálékkal és nagy növekedési rátával biztosítja, míg a K stratéga növekedési rátája alacsony, és szűkös táplálékviszonyok mellett is túlél.

ini={(Subscript[x, 1])[0]==1,(Subscript[y, 1])[0]==2,(Subscript[y, 2])[0]==3};

nds=NDSolve[Join[equ/.{K->3, μ1->7/10,μ2->5/10, a1->1/6, a2->3/2,r->3, β1->1/10,β2->3/10,λ->1},ini],{(Subscript[x, 1])[t],(Subscript[y, 1])[t],(Subscript[y, 2])[t]},{t,0,100}]

Needs["PlotLegends`"]

Plot[Evaluate[{(Subscript[x, 1])[t],(Subscript[y, 1])[t],(Subscript[y, 2])[t]}/.nds],{t,0,100},PlotStyle->Directive[Thick],PlotRange->{0,14},PlotLegend->{"zsákmány", "K-stratéga","r-stratéga"},LegendPosition->{0.7,0},AxesLabel->{t,Row[{"(Subscript[x, 1])(t), ","(Subscript[y, 1])(t), ","(Subscript[y, 2])(t)"}]}]

ParametricPlot[Evaluate[{(Subscript[x, 1])[t],(Subscript[y, 1])[t]}/.nds],{t,30,100},PlotStyle->Directive[Thick],AxesLabel->{"(Subscript[x, 1])(t)","(Subscript[y, 1])(t)"},AspectRatio->1]

Ismét szükségünk lesz a stacionárius pontok egyenesére:

stac[y1,K,a1,a2,μ1,r,d1,d2]/.{a1->1/6,a2->3/2, μ1->7/10 ,d1->6/10,d2->2/10,r->3}//Expand

Manipulate[DynamicModule[{nds,S,x,y},nds=NDSolve[{S'[t]==r S[t] (1-((S[t])/K))-((μ1 x[t]S[t])/(a1+S[t]))-((μ2 y[t]S[t])/(a2+S[t])),
x'[t]==x[t]((β1(S[t]-λ))/(a1+S[t])),y'[t]==y[t]((β2(S[t]-λ))/(a2+S[t])),S[0]==1,x[0]==2,y[0]==3}/.{ μ1->7/10,μ2->5/10, a1->1/6, a2->3/2,r->3, β1->1/10,β2->3/10,λ->1},{S,x,y},{t,0,300}][[1]];
Dynamic@ParametricPlot3D[{Evaluate[{S[t],x[t],y[t]}/.nds],Evaluate[{λ,t,15-(15/K)-3 t}/.{λ->1}]},{t,0,300},BoxRatios->{1,1,1},PlotRange->{{0,2},{0,7},{2,6}},AxesLabel->{"S[t]","x[t]","y[t]"},  PlotStyle->{Directive[Thick,Red],Directive[Thick,Blue]},PlotLabel->Style[K,"Section"]]],{{K,3},2,4,0.05},LocalizeVariables->False,SaveDefinitions->True]

Az ábra a zipzár bifurkációt mutatja. Az egyensúlyi helyzetek egy egyenesen helyezkednek el, kis, de λ-nál nagyobb K értékekre ezek Ljapunov-értelemben stabilisak, és az egyenes attraktora a rendszernek. K növekedésével a  K=a2+2λ értéknél az ℓ egyenes (S,y) síkon lévő pontja elveszíti stabilitását. Ha a2+2λ<K<a1+2λ, akkor van az egyenesnek egy pontja, amelytől balra az egyenes pontjai instabilisak, míg jobbra stabilisak. K növekedésével ez a pont az egyenesen balról jobbra mozog, mint ahogyan egy zipzár kinyílik, úgy veszítik el a hátrahagyott pontok stabilitásukat. Ez a jelenség a bőség paradoxona, mert a táplálékbőség destabilizálja a rendszert. Az r- és K-stratéga versenyében az r-stratéga előnyre tesz szert a K-stratégával szemben. A K=a1+2λ értéknél az (x,y) síkon Hopf-bifurkáció történik.

{a2+2λ,a1+2λ}/.{a1->1/6, a2->3/2,λ->1}