Quit[];

Néhány egyszerű sztochasztikus biológiai modell

A modelleket lásd például Érdi és Tóth Könyvének 5. fejezetében.

1. modell: Egyszerű beáramlás (Poisson-folyamat)

inflow[ini_:0,limit_,λ_]:=NestList[{#[[1]]+RandomReal[ExponentialDistribution[λ]],#[[2]]+1}&,{0,ini},limit]

Egymást λ paraméterű exponenciális eloszlású távolságban követő véletlen időpontokban növekszik a részecskék száma eggyel. Az időpontok és a darabszámok tipikus listája ilyen lehet:

inflow[0,15,3.6]

Mennyi ideig jutunk el a 15 sorsolással?

Mean[GammaDistribution[15,1/3.6]]

Show[ListPlot[inflow[0,15,3.6],Joined->True],Plot[3.6 t,{t,0,1.1Mean[GammaDistribution[15,1/3.6]]}],AxesLabel->{"idő","mennyiség"},PlotLabel->Style["Egyszerű beáramlás\n sztochasztikus és determinisztikus modelljének összevetése",16],ImageSize->800]

Show[ListPlot[inflow[0,150,3.6],Joined->True],Plot[3.6 t,{t,0,1.05Mean[GammaDistribution[150,1/3.6]]}],AxesLabel->{"idő","mennyiség"},ImageSize->800,PlotLabel->Style["Egyszerű beáramlás\n sztochasztikus és determinisztikus modelljének összevetése\nhosszabb időtartamon",16]]

Érdemes megfigyelni, hogy az egymás után kapott realizációk mennyire különböznek egymástól.

GraphicsGrid[Table[Show[ListPlot[inflow[0,150,3.6],Joined->True],Plot[3.6 t,{t,0,1.05Mean[GammaDistribution[150,1/3.6]]}]],{i,3},{j,3}],ImageSize->800]

Kiegészítő magyarázat

Determinisztikus megfelelője az x'=λ differenciálegyenlet.

A formális reakciókinetika jelöléseivel a 0(Overscript[->, λ])X folyamatról van szó.

2. modell: Születési-halálozási folyamat

A populáció  létszáma (vagy a részecskék száma) két ok miatt változhat: születés és halálozás.

Ha a folyamat az n állapotban van (n számú részecske van jelen), akkor (λ+μ)n paraméterű exponenciális eloszlású idő múlva növekszik vagy csökken a részecskék száma eggyel.  A kétféle változás valószínűsége (λ/(λ+μ)), illetve  (μ/(λ+μ)). Az időpontok és a darabszámok tipikus listáját előállító függvény és egy tipikus lista ilyen lehet:

birthanddeath[ini_:0,limit_,λ_,μ_]:=NestList[{#[[1]]+RandomReal[ExponentialDistribution[(λ+μ)#[[2]]]],#[[2]]+RandomChoice[{λ,μ}->{1,-1}]}&,{0,ini},limit]

birthanddeath[1,15,12,2]

ListPlot[birthanddeath[1,125,12,2],ImageSize->600,AxesStyle->Directive[Large],Joined->True,Epilog->First@Plot[(E^(10t)),{t,0,0.6},PlotStyle->Red],AxesLabel->{"idő","létszám"},PlotLabel->Style["Születési-halálozási folyamat\n sztochasztikus és determinisztikus modelljének összevetése",16]]

GraphicsGrid[Table[ListPlot[birthanddeath[5,125,12,4],ImageSize->600,AxesStyle->Directive[Large],Joined->True,AxesLabel->{"time","particle No"},Epilog->First@Plot[(E^(10t)),{t,0,0.6},PlotStyle->Red]],{i,3},{j,3}],ImageSize->1000]

Amikor λ=μ, akkor egyensúly körüli fluktuációt látunk:

ListPlot[birthanddeath[100,125,12,12],ImageSize->600,AxesStyle->Directive[Large],Joined->True,AxesOrigin->{0,0},Epilog->First@Plot[100,{t,0,0.6},PlotStyle->Red]]

GraphicsGrid[Table[ListPlot[birthanddeath[100,125,12,12],ImageSize->600,AxesStyle->Directive[Large],Joined->True,AxesLabel->{"time","particle No"},AxesOrigin->{0,50},Epilog->First@Plot[100,{t,0,0.6},PlotStyle->Red]],{i,3},{j,3}],ImageSize->1000]

Kiegészítő magyarázat

A formális reakciókinetika jelöléseivel az 0(Overscript[←, μ])X(Overscript[->, λ])2X folyamatról van szó.

Determinisztikus megfelelője az x'=λ x-μ x differenciálegyenlet.

Ez a folyamat tekinthető akár folytonos idejű, állapotfüggő bolyongásnak is.

3. modell: Volterra–Lotka-modell

Ha a folyamat az (m,n) állapotban van (m, illetve n számú részecske van jelen), akkor (Subscript[k, 1])m+(Subscript[k, 2])m n+(Subscript[k, 3])n paraméterű exponenciális eloszlású idő múlva változik a részecskék száma. A háromféle változás valószínűségét a táblázat mutatja:

Az időpontok és a darabszámok tipikus listáját előállító függvény és egy tipikus lista ilyen lehet:

volterra[ini_:0,limit_,k1_,k2_,k3_]:=NestList[{#[[1]]+RandomReal[ExponentialDistribution[k1#[[2,1]]+k2#[[2,1]]#[[2,2]]+k3#[[2,2]]]],#[[2]]+RandomChoice[{k1#[[2,1]],k2#[[2,1]]#[[2,2]],k3#[[2,2]]}->{{1,0},{-1,1},{0,-1}}]}&,{0,ini},limit]

v=volterra[{150,200},150,500,2,100];

Elkészítjük az idősorokat.

times=First/@v;

Transpose[Last/@v];

prey=Transpose[{First/@v,First@Transpose[Last/@v]}];

pred=Transpose[{First/@v,Last@Transpose[Last/@v]}];

ListPlot[Evaluate[{prey,pred}],Joined->True,PlotStyle->{Directive[Red],Directive[Blue]},AxesLabel->{"t","m,n"}]

Graphics[{Red,Point/@Transpose[Transpose[Last/@v]]}]

Összegezve és mellétéve a determinisztikust is:

nd=NDSolve[{x'[t]==x[t]-x[t] y[t],y'[t]==-y[t]+x[t] y[t],x[0]==0.7`,y[0]==1.2`},{x,y},{t,0,25}];

vlabra[ini_,limit_,k1_,k2_,k3_,opts___]:=Module[{v},v=volterra[ini,limit,k1,k2,k3];
prey=Transpose[{First/@v,First@Transpose[Last/@v]}];
pred=Transpose[{First/@v,Last@Transpose[Last/@v]}];
GraphicsRow[{ListPlot[Evaluate[{prey,pred}],Joined->True,PlotStyle->{Directive[Red],Directive[Blue]},AxesLabel->{"t","m,n"}],Graphics[{Red,Point/@Transpose[Transpose[Last/@v]]}],Plot[Evaluate[{x[t],y[t]}/.nd],{t,0,25}]},opts]
]

Kisebb részecskeszám esetén a realizációk kvalitatíve erősen különböznek egymástól

kicsi[sr_]:=(SeedRandom[sr];vlabra[{42,72},8000,1,1/60,1,ImageSize->800])

GraphicsColumn[Table[kicsi[i],{i,{30,90,20,70,80,100}}],Spacings->{0,-500}]

SeedRandom[100];
vlabra[{420,720},60000,1,1/600,1,ImageSize->800]

Az alábbi igen sokáig megy.

SeedRandom[2];
GraphicsGrid[Table[vlabra[{420,720},60000,1,1/600,1],{i,2},{j,2}],ImageSize->900]

Kiegészítő magyarázat

A formális reakciókinetika jelöléseivel az X(Overscript[->, (Subscript[k, 1])])2X, X+Y(Overscript[->, (Subscript[k, 2])])2Y, Y(Overscript[->, (Subscript[k, 3])])0, folyamatról van szó.

Determinisztikus megfelelője az x'=(Subscript[k, 1]) x-(Subscript[k, 2]) x y, y'=(Subscript[k, 2]) x y=(Subscript[k, 3])y differenciálegyenletrendszer.

Ez a folyamat tekinthető akár síkbeli, folytonos idejű, állapotfüggő bolyongásnak is.

Bebizonyítható, hogy a sztochasztikus modell stacionárius eloszlása az origóra koncentrált degenerált eloszlás,  tehát lényegesen viselkedik a hosszútávú viselkedés a determinsztikus modellétől.

Tetszőleges reakciókinetikai modell sztochasztikus modellje szimulálható a ReactionKinetics.m csomaggal.

Mathematica problémák

A szimulációt táblázatban végezve sokkal hosszabbak a futási idők.

Több szimulációt párhuzamosan illene végezni.

Hivatkozások

Érdi P.; Sipos T.; Tóth J.: Összetett kémiai reakciók sztochasztikus szimulálása számítógéppel, Magyar Kémiai Folyóirat 79 (3) (1973), 97–108.

Érdi, P.; Tóth, J.: Mathematical Models of Chemical Reactions. Theory and Applications of Deterministic and Stochastic Models, Manchester University Press, Manchester, Princeton University Press, Princeton, 1989.

Nagy, A. L.; Papp, D.; Tóth, J.: ReactionKinetics.m

Rényi A.: Kémiai reakciók tárgyalása a sztochasztikus folyamatok elméletével, Magyar Tud. Akad. Alkalm. Mat. Int. Kozl. 2 (1954), 83–101.