Quit[]

Térbeli mintázatok: Turing-szerkezetek

Az alábbiakban a térben elhelyezkedő populációkra és térben végbemenő kémiai reakciókra vonatkozó számítások szemléletes kiegészítése következik. (Szóhasználatunk gyakrabban utal az első esetre.) Itt elsősorban Murray alább megadott könyvére támaszkodtunk.

A modell

Evolúciós egyenletek

((∂u(t,x))/(∂ t))=(Subscript[d, u]) (((∂^2)u(t,x))/(∂ (x^2)))+f(u(t,x),v(t,x))

((∂v(t,x))/(∂ t))=(Subscript[d, v]) (((∂^2)v(t,x))/(∂ (x^2)))+g(u(t,x) ,v(t,x))

x∈Ω⊂(ℝ^1) vagy   (ℝ^2),     t∈(Subscript[ℝ, 0]^(+))

Dirichlet-féle peremfeltételek

u(t,x)=v(t,x)=0     (x∈∂Ω,  t∈(Subscript[ℝ, 0]^(+)))

Neumann-féle peremfeltételek

((∂u(t,x))/(∂ n))=((∂v(t,x))/(∂ n))=0     (x∈∂Ω,  t∈(Subscript[ℝ, 0]^(+)))

A megfelelő kinetikai differenciálegyenletrendszer

((d U(t))/(d t))=f(U(t),V(t))

((d V(t))/(d t))=g(U(t) ,V(t))

Jelölések

u(t,x) és v(t,x) valamilyen faj sűrűsége a t∈(Subscript[ℝ, 0]^(+)) időpontban az x∈Ω helyen

(Subscript[d, u]),(Subscript[d, v])>0 a diffúziós együttható

Ω⊂(ℝ^1) vagy   (ℝ^2) a vizsgált térrész, sima határú tartomány; határa: ∂Ω

n az Ω tartomány határáról kifelé mutató normális vektor

U(t) és V(t) valamilyen faj mennyisége a t∈(Subscript[ℝ, 0]^(+)) időpontban az egész Ω tartományban

A jelenség

A fenti (1)–(2) modell a következő szituációt írja le:

a tér egy- vagy kétdimenziós;

két faj, két kémiai anyagfajta, emberi populáció két részpopulációja egymásra reagál és a magasabb sűrűségű helyről az alacsonyabb sűrűségű helyre vándorolnak;

a (3) Dirichlet-féle peremfeltételek jelentése: a tartomány határán nincs egyed;

a (4) Neumann-féle peremfeltételek jelentése: a tartomány határán nincs ki- és beáramlás.

A fenti (5)–(6) modell eltekint a térbeli változásoktól.

Turing-bifurkáció

A diffúziótól általában azt várjuk, hogy tompítja az egyenlőtlenségeket. Előfordulhat azonban, hogy egy tében homogén stabilis dinamika a  térbeli változások hatására instabilissá válik, és egy inhomogén térbeli mintázat válik stabilissá. Ezt a forgatókönyvet Turing írta le először annak modellezésére, ahogyan a kezdetben térben homogén szerkezetű megtermékenyített petesejtből létrejön a bonyolult szerkezetű embrió. A jelenséget olyan kísérleti (kémiai) rendszerben, ahol valóban csak ezek a hatások érvényesülnek, Castel és mtsai mutatták ki majd 40 évvel később.

Egy speciális eset egy térdimenzióban

Előkészületek

Legyen

f1[u_,v_]:=u(10+(u/2)-((u^2)/200)-v)
g1[u_,v_]:=v(-(1/10)+u-(3/2)v)

és határozzuk meg az (5)–(6) kinetikai egyenletrendszer pozitív egyensúlyi helyzetét:

pozsol=First@Select[Solve[{f1[u,v]==0,g1[u,v]==0}],(u>0&&v>0/.#)&]

valamint nullavonalait:

ContourPlot[{f1[u,v]==0,g1[u,v]==0},{u,10,40},{v,10,40},PlotLabel->"Nullavonalak és \na pozitív egyensúlyi helyzet",FrameLabel->{Style[u,Italic,16],Style[v,Italic,16]},Epilog->First@Graphics[{Red,PointSize->Large,Point[{u,v}/.pozsol]}]]

A Turing-bifurkáció kialakulásának feltételei és alkalmazásuk

az (5)–(6) kinetikai differenciálegyenletrendszernek egy pozitív egyensúlyi helyzete van, E:=((OverBar[u]),(OverBar[v])).  (Ekkor az (u(t,x):=(OverBar[u]),v(t,x):=(OverBar[v])) konstans megoldás (!) az eredeti (1)–(2) rendszernek is megoldása.)

(∂_1)f((OverBar[u]),(OverBar[v]))(∂_2)g((OverBar[u]),(OverBar[v]))<0

ha (Subscript[d, u])=1, akkor kell, hogy d:=(Subscript[d, v])>(Subscript[d, crit]) legyen

A feltételeket ellenőrző program:

Tcheck[f_,g_]:=Module[{det,psol=First@Select[Solve[{f[u,v]==0,g[u,v]==0},{u,v}],(u>0&&v>0/.#)&],d,dq},
det=Det[D[{f[u,v],g[u,v]},{{u,v}}]/.psol];dq=Reduce[((((D[g[u,v],v]+d D[f[u,v],u]))^2)/.psol)>4 d det&&(D[g[u,v],v]+d D[f[u,v],u]/.psol)>0&&d>1,{d}];
If[(D[g[u,v],v]+D[f[u,v],u]/.psol)<0&&det>0,dq,False]]

Tcheck[f1,g1]//N

Ezekre a függvényekre igazak a feltételek; akkor most d értékét rögzítjük.

d=70;

A térben homogén egyensúlyi helyzet akkor instabilis, amikor az alábbi polinom negatív.

poly[λ_]:=Det[D[{f1[u,v],g1[u,v]},{{u,v}}]]-λ(D[g1[u,v],v]+d D[f1[u,v],u])+(λ^2)d/.pozsol

Plot[Evaluate[poly[λ]],{λ,0,6}]

Egy dimenzióban, azaz amikor Ω=[0,ℓ], a Laplace-egyenlet megoldásai az x↦(Subscript[R, k])(x):=cos(k π (x/ℓ)) függvények, a hozzájuk tartozó sajátérték: (Subscript[λ, k])=((((k/ℓ)π))^2)     (k∈ℤ). ℓ növelésével a sajátértékek sűrűbben helyezkednek el, elérhető a Turing-bifurkáció, viszont rögzített k mellett (Subscript[λ, k])végul kiesik az intervallumból, egyre több es egyre nagyobb frekvenciájú hullám jelenik meg.

Reduce[poly[((((k/1)π))^2)]<0&&k>0,{k},Integers]

Reduce[poly[((((k/5)π))^2)]<0&&k>0,{k},Integers]

Reduce[poly[((((k/32)π))^2)]<0&&k>0,{k},Integers]

Az egyenleteket egységesen a [0,1] szakaszra írjuk fel, az ℓ hossz az egyenletben jelenik meg, ennek értékét változtatjuk, az eredmények összhangban vannak a linearizált modell alapján tett következtetésekkel.

ClearAll[pde1];
pde1[ℓ_]:={
D[u[t,x],t]==(1/(ℓ^2)) D[u[t,x],{x,2}]+f1[u[t,x],v[t,x]],
D[v[t,x],t]== (d/(ℓ^2)) D[v[t,x],{x,2}]+g1[u[t,x],v[t,x]]};

cosinusok[x_,tagok_:7]:=(1/3)Sum[Cos[i π x],{i,1,tagok}]

neumann1={
Derivative[0,1][u][t,0]==Derivative[0,1][u][t,1]==0,
u[0,x]==24 (x^2)(((1-x))^2)+3cosinusok[x]+31,
Derivative[0,1][v][t,0]==Derivative[0,1][v][t,1]==0,
v[0,x]==25 (x^2)(((1-x))^2)-2cosinusok[x]+20};

pdsol1[ℓ_,maxtime_]:=NDSolve[{pde1[ℓ],neumann1},{u,v},{t,0,maxtime},{x,0,1}]

Plot[Evaluate[{u[30,x],v[30,x]}/.pdsol1[1,30]],{x,0,1},PlotRange->{0,50},AxesLabel->{x/ℓ,{u,v}},PlotLabel->"u és v mint a hely függvénye,\n hosszú idő után; ℓ=1"]

Plot[Evaluate[{u[50,x],v[50,x]}/.pdsol1[5,50]],{x,0,1},PlotRange->{0,50},AxesLabel->{x/ℓ,{u,v}},PlotLabel->"u és v mint a hely függvénye,\n hosszú idő után; ℓ=5"]//Quiet

Plot[Evaluate[{u[50,x],v[50,x]}/.pdsol1[32,50]],{x,0,1},PlotRange->{0,50},AxesLabel->{x/ℓ,{u,v}},PlotLabel->"u és v mint a hely függvénye,\n hosszú idő után; ℓ=32"]//Quiet

Feladat: Kísérletezzünk az alábbi kinetikával:

f2[x_,y_]:=a-x-((ρ x y)/(1+x+k (x^2)))/.{k->0.1,ρ->18.5,a->92}
g2[x_,y_]:=α(b-y)-((ρ x y)/(1+x+k (x^2)))/.{k->0.1,ρ->18.5,α->1.5,b->56}

Síkbeli jelenségek: a térdimenziók száma kettő

Síkban bonyolultabbak a sajátfüggvények, a mintázatok és a számolások. A példákban Dirichlet-peremfeltétel szerepel a rövidebb futásidő érdekében.

Előkészületek

pde2[ℓ_,w_]:={
(∂_t)u[t,x,y]==(1/(ℓ^2))(∂_(x,x))u[t,x,y]+(1/(w^2))(∂_(y,y))u[t,x,y]+f1[u[t,x,y],v[t,x,y]],
(∂_t)v[t,x,y]==d((1/(ℓ^2))(∂_(x,x))v[t,x,y]+(1/(w^2))(∂_(y,y))v[t,x,y])+g1[u[t,x,y],v[t,x,y]]};

dirichlet={
u[t,x,1]==u[t,x,0]==u[t,0,y]==u[t,1,y]==31,u[0,x,y]==31+3Sin[π x] Sin[2π y],
v[t,x,1]==v[t,x,0]==v[t,0,y]==v[t,1,y]==20,v[0,x,y]==20+3Sin[2π x] Sin[3π y]};

pdsol2[ℓ_,w_,maxtime_,ini_:dirichlet,opts___]:=NDSolve[{pde2[ℓ,w],ini},{u,v},{t,0,maxtime},{x,0,1},{y,0,1},opts]

Példa: ℓ = 6, w = 3

Plot3D[Evaluate[{u[50,x,y],v[50,x,y]}/.pdsol2[6,3,50][[1]]],{x,0,1},{y,0,1},PerformanceGoal->"Quality",AxesLabel->{x,y,{u,v}},PlotLabel->"ℓ = 6, w = 3"]

És itt látszik a mintázat.

RegionPlot[Evaluate[u[50,x,y]>39/.pdsol2[6,3,50][[1]]],{x,0,1},{y,0,1},FrameLabel->{x,y},PlotLabel->"ℓ = 6, w = 3"]

Példa: ℓ = 15, w = 9

elore=pdsol2[15,9,50,dirichlet,MaxStepSize->{0.1,0.01,0.01}][[1]];

GraphicsRow[{Plot3D[Evaluate[{u[50,x,y],v[50,x,y]}/.elore],{x,0,1},{y,0,1},PerformanceGoal->"Quality",AxesLabel->{"\n"x,"\n"y,{u,v"  "}},PlotLabel->"ℓ = 15, w = 9"],RegionPlot[Evaluate[u[50,x,y]>39/.elore],{x,0,1},{y,0,1},FrameLabel->{x,y},PlotLabel->"ℓ = 15, w = 9"]}]

Síkbeli mintázatok: konkrétabb biológiai alkalmazás

Jelölések

u(t,x) és v(t,x) valamilyen kémiai anyag (aktivátor, illetve inhibítor) sűrűsége a t időpontban az x helyen

(∂_1)f((OverBar[u]),(OverBar[v]))>0 mutatja, hogy u katalizálja saját szaporodását (autokatalízis)

(∂_2)g((OverBar[u]),(OverBar[v]))<0 mutatja, hogy v gátolja saját szaporodását (autoinhibíció)

(Subscript[d, u]),(Subscript[d, v])>0 a diffúziós együttható

ℓ>0 a vizsgált egydimenziós térrész mérete

A jelenség

A fenti (1)–(2) modell a következő szituációt írja le:

a vizsgált térrész kétdimenziós: Ω:={(x,y);  -ℓ<x<ℓ, -w b(x)<y<w b(x)}, ahol az alábbiakban b lineáris függvény

rp[ℓ_,w_,b_,opts___]:=RegionPlot[{-ℓ<x<ℓ,-w b[x]<y<w b[x]},{x,-5,5},{y,-20,20},opts]

domain1[x_]:=((3-x)/2)

rp[3,6,domain1,FrameLabel->{x,y},PlotLabel->"A sötét térrészben\n zajlik a folyamat",Epilog->First@Graphics[{Text[Style["Ω",16],{-1,0}]}]]

Egy másfajta ábrázolási lehetőség

bp[ℓ_,w_,b_,opts___]:=Plot3D[Boole[-ℓ<x<ℓ∧-w b[x]<y<w b[x]],{x,-5,5},{y,-20,20},opts]

bp[3,6,domain1,AxesLabel->{x,y,1},PlotPoints->33]

Láttuk korábban, hogy a tartomány mérete sokat számít. Ha túl kicsi, akkor nem alakul ki minta, ha túl nagy, akkor az eredmény sok hullám szuperpozíciója, ami nagyon összetett, a különbségek nem térnek el nagyon egymástól. Ennek következménye, hogy kis állatok nem mintásak. Kicsit nagyobbakon megjelenhetnek keresztirányu csíkok, sávok. Ha elég széles az állat, hosszirányban is mintás lesz, ami együttesen foltokat, pöttyöket eredményez. Majd a hossziranyú csíkok maradnak, végül ismét nem alakul ki minta. A két véglet egyezik a tapasztalattal, az egér, az elefánt egyszínű. De egyes zebrafajok között is van különbség aszerint, hogy mikor megy végbe a feltételezett folyamat. Ha később, akkor nagyobb az embrió, több csík jelenik meg, amelyek sűrűbben helyezkednek el.

Az Ω  tartományt a [-1,1]×[-1,1] négyzetbe transzformáló függvények:

p[x_,y_]:=(x/ℓ)
q[x_,y_]:=(y/(w b[(x/ℓ)]))

Transzformálódnak a koncentrációk is: u(t,x,y)=:uu(t,p(x,y),q(x,y))     v(t,x,y)=:vv(t,p(x,y),q(x,y)), így az egyenletek új alakjához:

Δuu=Simplify[(D[uu[t,p[x,y],q[x,y]],{x,2}]+D[uu[t,p[x,y],q[x,y]],{y,2}])/.{x->ℓ p,y->w q b[ p]}];

Δvv=Δuu/.{uu->vv};

pde3={(∂_t)uu[t,x,y]==Δuu+f[uu[t,x,y],vv[t,x,y]],
(∂_t)vv[t,x,y]==d Δvv+g[uu[t,x,y],vv[t,x,y]]}/.{p->x,q->y};

dirichlet3={
uu[t,x,1]==uu[t,x,-1]==uu[t,-1,y]==uu[t,1,y]==31,uu[0,x,y]==31+3Sin[π x] Sin[2π y],
vv[t,x,1]==vv[t,x,-1]==vv[t,-1,y]==vv[t,1,y]==20,vv[0,x,y]==20-2Sin[2π x] Sin[3π y]};

dirichlet3periodic={
uu[t,x,1]==uu[t,x,-1],uu[t,-1,y]==uu[t,1,y]==31,uu[0,x,y]==31+3Sin[π x] Sin[2π y],
vv[t,x,1]==vv[t,x,-1],vv[t,-1,y]==vv[t,1,y]==20,vv[0,x,y]==20-2Sin[2π x] Sin[3π y]};

zaj=Sum[Sin[RandomInteger[10]π x]Sin[RandomInteger[10]π y],{i,1,5}];

neumann3periodic={
uu[t,x,1]==uu[t,x,-1],
Derivative[0,1,0][uu][t,-1,y]==Derivative[0,1,0][uu][t,1,y]==0,
uu[0,x,y]==31+3Cos[π x] Cos[2π y],
vv[t,x,1]==vv[t,x,-1],
Derivative[0,1,0][vv][t,-1,y]==Derivative[0,1,0][vv][t,1,y]==0,
vv[0,x,y]==20-2Cos[2π x] Cos[3π y]};

Egy állat farka leírható egy kúppal, ami egy trapéz periodikus peremfeltétel mellett. Az eredmény sok pötty.

adat1={ℓ->15,w->9,b->domain1,f->f1,g->g1,d->70};

Timing[pdsol3=NDSolve[{pde3/.adat1,dirichlet3},{uu,vv},{t,0,50},{x,-1,1},{y,-1,1},MaxStepSize->{1,0.1,0.1}]]//Quiet

Manipulate[Plot3D[Evaluate[{uu[t,x,y],vv[t,x,y]}/.pdsol3],{x,-1,1},{y,-1,1},AxesLabel->{"\n"x," "y,"u,v  "},PerformanceGoal->"Quality",Mesh->None],{t,0,50,1}]

Plot3D[Evaluate[{uu[50,x,y],vv[50,x,y]}/.pdsol3],{x,-1,1},{y,-1,1},AxesLabel->{"\n"x," "y,"u,v  "},PerformanceGoal->"Quality",Mesh->None]

Timing[pdsol3=NDSolve[{pde3/.adat1,dirichlet3},{uu,vv},{t,0,50},{x,-1,1},{y,-1,1},MaxStepSize->{1,0.02,0.02}]]//Quiet

Plot3D[Evaluate[{uu[50,x,y],vv[50,x,y]}/.pdsol3],{x,-1,1},{y,-1,1},AxesLabel->{"\n"x," "y,"u,v  "},PerformanceGoal->"Quality",Mesh->None]

Ezért kár volt sokkal hosszabb ideig számolni.

RegionPlot[Evaluate[(uu[50,x,y]>35)/.pdsol3],{x,-1,1},{y,-1,1}]

ParametricPlot3D[{p ℓ,w q b[p],uu[50,p,q]}/.adat1/.pdsol3,{p,-1,1},{q,-1,1},PerformanceGoal->"Quality",Mesh->None,ColorFunction->Function[{x,y,z},Hue[z]],ViewPoint->{0,0,+∞}]

Feladatok

Kísérletezzünk (türelmesen!) az alábi adatokkal:

domain2[x_]:=1.05-x

adat2={ℓ->25,w->20,b->domain2,f->f1,g->g1,d->80};

dirichlet3periodiczaj={uu[t,x,1]==uu[t,x,-1],uu[t,-1,y]==uu[t,1,y]==31,uu[0,x,y]==31+3zaj,
vv[t,x,1]==vv[t,x,-1],vv[t,-1,y]==vv[t,1,y]==20,vv[0,x,y]==20-2zaj};

Az alábbi eseetben a hosszirányú csíkok megmaradnak. Ujjlenyomatra jobban hasonlít.

adat3={ℓ->100,w->10,b->domain2,f->f1,g->g1,d->80};

ParametricPlot3D[{p ℓ,w q b[p],uu[50,p,q]}/.adat2/.p3sol2,{p,-1,1},{q,-1,1},PerformanceGoal->"Quality",Mesh->None,ColorFunction->Function[{x,y,z},Hue[z]],ViewPoint->{0,0,∞}]

Egy állat elnagyolt testfelépítése.

domain=Interpolation[{{-0.4,0.5},{-0.5,1},{-0.6,1},{-0.7,0.3},{-1,0.2},
{0,0.5},{0.4,0.5},{0.5,1},{0.6,1},{0.7,0.2},{1,0.1}},
InterpolationOrder->1];

adat4={ℓ->20,w->6,b->domain,f->f1,g->g1,d->80};

Problémák a Mathematicával

Miért ne lennének kompatibilisek a kezdeti és a peremfeltételek az első feladatnál

Hivatkozások

Castets, V.; Dulos, E., Boissonade, J.; De Kepper, P.: Experimental evidence of a sustained standing Turing-type nonequilibrium chemical pattern, Phys. Rev. Lett. 64 (1990), 2953–2957.

Farkas Miklós: Dynamical Models in Biology

Murray, J. D.: Mathematical Biology, Biomathematics, Volume 19, Springer Verlag

Köszönetnyilvánítás

Ez a notebook jelentősen támaszkodik Kovács Dániel matematikus hallgató munkájára.