Gillespieho algoritmus - Gillespie algorithm

v teorie pravděpodobnosti, Gillespieho algoritmus (nebo příležitostně Algoritmus Doob-Gillespie) generuje statisticky správnou trajektorii (možné řešení) a stochastický systém rovnic, pro který reakční rychlosti jsou známy. Vytvořil Joseph L. Doob a další (kolem roku 1945), předložil Dan Gillespie v roce 1976 a popularizován v roce 1977 v příspěvku, kde jej využívá k efektivní a přesné simulaci chemických nebo biochemických systémů reakcí s využitím omezené výpočetní síly (viz stochastická simulace )^{[Citace je zapotřebí ]}. Jak se počítače zrychlily, algoritmus se používal k simulaci stále složitějších systémů. Algoritmus je zvláště užitečný pro simulaci reakcí v buňkách, kde je počet činidla je nízká a sledování polohy a chování jednotlivých molekul je výpočetně proveditelné. Matematicky jde o variantu a dynamická metoda Monte Carlo a podobně jako kinetický Monte Carlo metody. Používá se těžce v biologie výpočetních systémů.^{[Citace je zapotřebí ]}

Dějiny

Proces, který vedl k algoritmu, uznává několik důležitých kroků. V roce 1931 Andrei Kolmogorov představil diferenciální rovnice odpovídající časovému vývoji stochastických procesů, které probíhají skoky, dnes známé jako Kolmogorovovy rovnice (Markovův skok) (zjednodušená verze je známá jako hlavní rovnice v přírodních vědách). to bylo William Feller, v roce 1940, který našel podmínky, za nichž Kolmogorovovy rovnice připouštěly (správné) pravděpodobnosti jako řešení. Ve své teorémě I (dílo z roku 1940) uvádí, že čas do dalšího skoku byl exponenciálně rozdělen a pravděpodobnost další události je úměrná rychlosti. Jako takový navázal vztah Kolmogorovových rovnic s stochastické procesy Později Doob (1942, 1945) rozšířil Fellerova řešení nad rámec procesů čistého skoku. Metodu v počítačích implementoval David George Kendall (1950) pomocí Manchester Mark 1 počítač a později používán Maurice S. Bartlett (1953) ve studiích propuknutí epidemií. Gillespie (1977) získává algoritmus jiným způsobem pomocí fyzického argumentu.

Nápad za algoritmem

Tradiční kontinuální a deterministická biochemie rychlostní rovnice nepředpovídají přesně buněčné reakce, protože se spoléhají na hromadné reakce, které vyžadují interakce milionů molekul. Obvykle se modelují jako sada spojených obyčejných diferenciálních rovnic. Naproti tomu Gillespieho algoritmus umožňuje diskrétní a náhodnou simulaci systému s několika reaktanty, protože každá reakce je výslovně simulována. Trajektorie odpovídající jedné Gillespie simulaci představuje přesný vzorek z funkce pravděpodobnostní hmotnosti, který je řešením hlavní rovnice.

Fyzickým základem algoritmu je srážka molekul v reakční nádobě. Předpokládá se, že kolize jsou časté, ale kolize se správnou orientací a energií se vyskytují jen zřídka. Všechny reakce v rámci Gillespieho systému proto musí zahrnovat nejvýše dvě molekuly. Reakce zahrnující tři molekuly se považují za extrémně vzácné a jsou modelovány jako posloupnost binárních reakcí. Rovněž se předpokládá, že reakční prostředí je dobře promícháno.

Algoritmus

Nedávný přehled (Gillespie, 2007) nastiňuje tři různé, ale ekvivalentní formulace; přímé, první reakce a metody první rodiny, přičemž první dva jsou zvláštními případy druhého. Formulace metod přímé a první reakce je zaměřena na provedení obvyklých kroků inverze Monte-Carlo na takzvaném „základním předpokladu stochastické chemické kinetiky“, což je matematicky funkce

${displaystyle p (au, j | {oldsymbol {x}}, t) = a_ {j} ({oldsymbol {x}}) exp (au sum _ {j} a_ {j} ({oldsymbol {x}}) )}$,

kde každý z ${displaystyle a}$ termíny jsou sklonové funkce elementární reakce, jejichž argumentem je ${displaystyle {oldsymbol {x}}}$se počítá vektor druhů. The ${displaystyle au}$ parametr je čas do další reakce (nebo doba pobytu) a ${displaystyle t}$ je samozřejmě aktuální čas. Abych parafrázoval Gillespie, tento výraz se čte jako „uvedená pravděpodobnost ${displaystyle {oldsymbol {X}} (t) = {oldsymbol {x}}}$, že k další reakci systému dojde v nekonečně časovém intervalu ${displaystyle [t + au, t + au + d au]}$a bude stechiometrie odpovídající ${displaystyle j}$Tato reakce ". Tato formulace poskytuje okno přímým metodám a metodám první reakce naznačením ${displaystyle au}$ je exponenciálně distribuovaná náhodná proměnná a ${displaystyle j}$ je „statisticky nezávislá celočíselná náhodná proměnná s pravděpodobností bodu ${displaystyle a_ {j} ({oldsymbol {x}}) / součet _ {j} a_ {j} ({oldsymbol {x}})}}$".

Metodou generování Monte Carla je tedy jednoduše nakreslit dvě pseudonáhodná čísla, ${displaystyle r_ {1}}$ a ${displaystyle r_ {2}}$ na ${displaystyle [0,1]}$a vypočítat

${displaystyle au = {frac {1} {sum _ {j} a_ {j} ({oldsymbol {x}})}} log ({frac {1} {r_ {1}}})}$,

a

${displaystyle j =}$ nejmenší celé číslo vyhovující ${displaystyle sum _ {j '= 1} ^ {j} a_ {j'} ({oldsymbol {x}})> r_ {2} sum _ {j} a_ {j} ({oldsymbol {x}})}$.

Použitím této generující metody pro dobu pobytu a další reakci je algoritmus přímé metody uveden Gillespie jako

1. Inicializujte čas ${displaystyle t = t_ {0}}$ a stav systému ${displaystyle {oldsymbol {x}} = {oldsymbol {x}} _ {0}}$2. Se systémem ve stavu ${displaystyle {oldsymbol {x}}}$ v čase ${displaystyle t}$, vyhodnotit všechny ${displaystyle a_ {j} ({oldsymbol {x}})}$ a jejich součet ${displaystyle sum _ {j} a_ {j} ({oldsymbol {x}})}$3. Vyměňte další reakci za novou ${displaystyle tleftarrow t + au}$ a ${displaystyle {oldsymbol {x}} leftarrow {oldsymbol {x}} + u _ {j}}$4. Záznam ${displaystyle ({oldsymbol {x}}, t)}$ podle přání. Vraťte se ke kroku 1, jinak simulaci ukončete.

Tato skupina algoritmů je výpočetně nákladná, a proto existuje mnoho modifikací a adaptací, včetně další metody reakce (Gibson & Bruck), tau-skákat, stejně jako hybridní techniky, kde jsou hojné reaktanty modelovány s deterministickým chováním. Přizpůsobené techniky obecně kompromitují přesnost teorie za algoritmem, jak se připojuje k hlavní rovnici, ale nabízejí rozumné realizace pro výrazně vylepšené časové rámce. Výpočtová cena přesných verzí algoritmu je určena třídou vazby reakční sítě. Ve slabě vázaných sítích je počet reakcí ovlivněných jakoukoli jinou reakcí omezen malou konstantou. V silně vázaných sítích může jedno odpálení reakce v zásadě ovlivnit všechny ostatní reakce. Byla vyvinuta přesná verze algoritmu s konstantním časovým měřítkem pro slabě vázané sítě, která umožňuje efektivní simulaci systémů s velmi velkým počtem reakčních kanálů (Slepoy Thompson Plimpton 2008). Zobecněný Gillespieho algoritmus, který zohledňuje nemarkovovské vlastnosti náhodných biochemických událostí se zpožděním, vyvinuli Bratsun et al. 2005 a nezávisle Barrio et al. 2006, jakož i (Cai 2007). Podrobnosti viz níže uvedené články.

Formulace s částečným sklonem, jak byly vyvinuty nezávisle Ramaswamy et al. (2009, 2010) a Indurkhya a Beal (2010) jsou k dispozici ke konstrukci rodiny přesných verzí algoritmu, jejichž výpočetní náklady jsou úměrné počtu chemických druhů v síti, spíše než (většímu) počtu reakcí. Tyto formulace mohou snížit výpočetní náklady na škálování v konstantním čase pro slabě vázané sítě a škálovat maximálně lineárně s počtem druhů pro silně vázané sítě. Byla také navržena varianta generalizovaného Gillespieho algoritmu pro reakce se zpožděním s částečným sklonem (Ramaswamy Sbalzarini 2011). Použití metod částečné náchylnosti je omezeno na elementární chemické reakce, tj. Reakce s nejvýše dvěma různými reaktanty. Každá neelementární chemická reakce může být ekvivalentně rozložena na sadu elementárních, na úkor lineárního (v pořadí reakce) zvětšení velikosti sítě.

Zvažte následující objektově orientovanou implementaci metod přímé a první reakce, které jsou obsaženy ve třídě Pythonu 3:

z matematika import logtřída SSA:    "" "Kontejner pro SSA" ""    def __init__(já, Modelka, semínko=1234):        "" "Inicializovat kontejner pomocí modelu a generátoru pseudonáhodných čísel" ""        já.Modelka = Modelka        já.náhodný = Mersenne(semínko=semínko)    def Přímo(já):        "" "Neurčitý generátor přímých trajektorií" ""        zatímco Skutečný:            zatímco ne já.Modelka.výstup():                # vyhodnotit váhy a rozdělení                závaží = [                    (rxn, sto, pro(já.Modelka))                    pro (rxn, sto, pro) v já.Modelka.reakce                ]                rozdělit = součet(w[-1] pro w v závaží)                # vyhodnotit čas pobytu (MC krok 1)                pobývat = log(1.0 / já.náhodný.plovoucí()) / rozdělit                já.Modelka["čas"].připojit(já.Modelka["čas"][-1] + pobývat)                # vyhodnotit reakci (MC krok 2)                rozdělit = rozdělit * já.náhodný.plovoucí()                zatímco rozdělit >= 0.0:                    rxn, sto, pro = závaží.pop(0)                    rozdělit -= pro                pro druh, delta v sto.položky():                    já.Modelka[druh].připojit(já.Modelka[druh][-1] + delta)                já.Modelka.kaplan()            výtěžek já.Modelka            já.Modelka.resetovat()    def first_reaction(já):        "" "Neurčitý generátor trajektorií první reakce" ""        zatímco Skutečný:            zatímco ne já.Modelka.výstup():                # vyhodnotit další reakční časy                krát = [                    (                        log(                            1.0 / já.náhodný.plovoucí()                        ) / pro(já.Modelka),                        sto                    )                    pro (rxn, sto, pro) v já.Modelka.reakce                ]                krát.třídit()                # vyhodnotit reakční čas                já.Modelka["čas"].připojit(                    já.Modelka["čas"][-1] + krát[0][0]                )                # vyhodnotit reakci                pro druh, delta v krát[0][1].položky():                    já.Modelka[druh].připojit(                        já.Modelka[druh][-1] + delta                    )                já.Modelka.kaplan()            výtěžek já.Modelka            já.Modelka.resetovat()

Jak je vidět, stejně jako všechny metody Monte-Carlo, SSA vyžaduje semeno pro reprodukovatelnost, které se předává jménu Mersenne, generátor pseudonáhodných čísel. Implementace Mersenne najdete v Mersenne Twister článek, ale lze použít i jiný, např náhodný. náhoda. The Přímo a first_reaction členové jsou neurčité generátory, což znamená, že budou nepřetržitě vytvářet trajektorie, přičemž každá trajektorie je úplnou simulací systému, ve smyčce, dokud tuto smyčku nerozbije signál. Chcete-li skutečně implementovat takovou smyčku a získat určitý počet trajektorií pro analýzu, vyžaduje tato třída model, který jí byl předán v instanci. Účelem modelové třídy je vyhnout se směšování kinetických vlastností konkrétního procesu simulovaného logikou Gillespieho algoritmu. Model by měl být podtřídou slovníku s veřejnými členy kaplan, výstup, a resetovat že příslušně

• určit, které reakce jsou a nejsou platné na konci každé iterace dané trajektorie;
• vrátit se Skutečný pokud již nejsou možné reakce;
• na konci dané trajektorie vrátí hodnoty hash modelu na jejich původní hodnoty (tj. jejich počáteční podmínky).

Zvažte následující implementaci modelu, která se používá v níže uvedeném příkladu SIR:

třída SSAModel(diktát):    "" "Kontejner pro model SSA" ""    def __init__(        já, počáteční podmínky, sklony, stechiometrie    ):        """        Inicializovat model pomocí slovníku počátečních podmínek (každý        """        super().__init__(**počáteční podmínky)        já.reakce = seznam()        já.vyloučené_reakce = seznam()        pro reakce,sklon v sklony.položky():            -li sklon(já) == 0.0:                já.vyloučené_reakce.připojit(                    (                        reakce,                        stechiometrie[reakce],                        sklon                    )                )            jiný:                já.reakce.připojit(                    (                        reakce,                        stechiometrie[reakce],                        sklon                    )                )    def výstup(já):        "" "Návrat True k vymanění se z trajektorie" ""        # return True, pokud už žádné reakce        -li len(já.reakce) == 0: vrátit se Skutečný        # return False, pokud existuje více reakcí        jiný: vrátit se Nepravdivé    def kaplan(já):        "" "Ověření a zneplatnění modelových reakcí" ""                # Vyhodnoťte možné reakce        reakce = []        zatímco len(já.reakce) > 0:            reakce = já.reakce.pop()            -li reakce[2](já) == 0:                já.vyloučené_reakce.připojit(reakce)            jiný:                reakce.připojit(reakce)        reakce.třídit()        já.reakce = reakce        # vyhodnotit nemožné reakce        vyloučené_reakce = []        zatímco len(já.vyloučené_reakce) > 0:            reakce = já.vyloučené_reakce.pop()            -li reakce[2](já) > 0:                já.reakce.připojit(reakce)            jiný:                vyloučené_reakce.připojit(reakce)        vyloučené_reakce.třídit()        já.vyloučené_reakce = vyloučené_reakce    def resetovat(já):        "" "Vymazat trajektorii" ""        # resetovat druhy na počáteční podmínky        pro klíč v já: del já[klíč][1:]        # reset reakcí za počátečních podmínek        já.kaplan()

Příklady

Reverzibilní vazba A a B za vzniku AB dimerů

Jednoduchý příklad může pomoci vysvětlit, jak funguje Gillespieho algoritmus. Zvažte systém molekul dvou typů, A a B. V tomto systému A a B reverzibilně se spojí a vytvoří AB dimery takové, že jsou možné dvě reakce: buď A a B reagují reverzibilně za vzniku an AB dimer, nebo AB dimer disociuje na A a B. Konstanta reakční rychlosti pro danou jedinou molekulu A reagující s danou jedinou molekulou B molekula je ${displaystyle k_ {mathrm {D}}}$a reakční rychlost pro AB rozpad dimeru je ${displaystyle k_ {mathrm {B}}}$.

Pokud v době t existuje jedna molekula každého typu, pak je rychlost tvorby dimeru ${displaystyle k_ {mathrm {D}}}$, zatímco pokud existují ${displaystyle n_ {mathrm {A}}}$ molekuly typu A a ${displaystyle n_ {mathrm {B}}}$ molekuly typu B, rychlost tvorby dimeru je ${displaystyle k_ {mathrm {D}} n_ {mathrm {A}} n_ {mathrm {B}}}$. Pokud existují ${displaystyle n_ {mathrm {AB}}}$ dimery pak je rychlost disociace dimeru ${displaystyle k_ {mathrm {B}} n_ {mathrm {AB}}}$.

Celková rychlost reakce, ${displaystyle R_ {mathrm {TOT}}}$, v čase t je pak dáno

${displaystyle R_ {mathrm {TOT}} = k_ {mathrm {D}} n_ {mathrm {A}} n_ {mathrm {B}} + k_ {mathrm {B}} n_ {mathrm {AB}}}$

Nyní jsme tedy popsali jednoduchý model se dvěma reakcemi. Tato definice je nezávislá na Gillespieho algoritmu. Nyní popíšeme, jak použít Gillespieho algoritmus na tento systém.

V algoritmu postupujeme vpřed v čase ve dvou krocích: výpočet času do další reakce a určení, která z možných reakcí je další reakcí. Reakce se považují za zcela náhodné, takže pokud jde o rychlost reakce najednou t je ${displaystyle R_ {mathrm {TOT}}}$, pak čas, δt, dokud nedojde k další reakci, je náhodné číslo čerpané z funkce exponenciálního rozdělení s průměrem ${displaystyle 1 / R_ {mathrm {TOT}}}$. Posunuli jsme tedy čas z t na t + δt.

Spiknutí čísla A molekuly (černá křivka) a AB dimery jako funkce času. Jak jsme začali s 10 A a B molekuly v čase t= 0, počet B molekul se vždy rovná počtu A molekuly a tak to není zobrazeno.

Pravděpodobnost, že tato reakce je A vazba molekuly na a B molekula je jednoduše zlomek celkové rychlosti způsobené tímto typem reakce, tj.

pravděpodobnost, že reakce bude ${displaystyle P ({ce {{A} + B -> AB}}) = k_ {D} n_ {A} n_ {B} / R _ {{ce {TOT}}}}$

Pravděpodobnost, že další reakce bude AB disociace dimeru je jen 1 minus. Takže s těmito dvěma pravděpodobnostmi buď vytvoříme dimer redukcí ${displaystyle n_ {mathrm {A}}}$ a ${displaystyle n_ {mathrm {B}}}$ o jeden, a zvýšit ${displaystyle n_ {mathrm {AB}}}$ o jednu, nebo disociujeme dimer a zvyšujeme ${displaystyle n_ {mathrm {A}}}$ a ${displaystyle n_ {mathrm {B}}}$ o jednu a snížit ${displaystyle n_ {mathrm {AB}}}$ jedním.

Nyní máme oba pokročilý čas t + δta provedl jedinou reakci. Algoritmus Gillespie právě opakuje tyto dva kroky tolikrát, kolikrát je potřeba k simulaci systému, jak dlouho chceme (tj. Pro tolik reakcí). Výsledek simulace Gillespie, která začíná ${displaystyle n_ {mathrm {A}} = n_ {mathrm {B}} = 10}$ a ${displaystyle n_ {mathrm {AB}} = 0}$ na t= 0 a kde ${displaystyle k_ {mathrm {D}} = 2}$ a ${displaystyle k_ {mathrm {B}} = 1}$, je zobrazen vpravo. Pro tyto hodnoty parametrů je v průměru 8 ${displaystyle n_ {mathrm {AB}}}$ dimery a 2 z A a B ale vzhledem k malému počtu molekul jsou fluktuace kolem těchto hodnot velké. Gillespieho algoritmus se často používá ke studiu systémů, kde jsou tyto výkyvy důležité.

Byl to jen jednoduchý příklad se dvěma reakcemi. Složitější systémy s více reakcemi jsou řešeny stejným způsobem. Všechny reakční rychlosti musí být vypočítány v každém časovém kroku a musí být vybrána jedna s pravděpodobností rovnou jejímu dílčímu příspěvku k rychlosti. Čas je poté pokročilý jako v tomto příkladu.

Epidemie SIR bez zásadní dynamiky

The Model SIR je klasický biologický popis toho, jak určité nemoci procházejí populací pevné velikosti. Ve své nejjednodušší podobě existují ${displaystyle N}$ příslušníci populace, přičemž každý člen může být v kterémkoli okamžiku v jednom ze tří států - vnímavý, infikovaný nebo uzdravený a každý takový člen nevratně prochází těmito státy podle níže uvedeného směrovaného grafu. Můžeme označit počet náchylných členů jako ${displaystyle n_ {S}}$, počet infikovaných členů jako ${displaystyle n_ {I}}$a počet obnovených členů jako ${displaystyle n_ {R}}$. Proto ${displaystyle N = n_ {S} + n_ {I} + n_ {R}}$ pro jakýkoli okamžik.

Daný náchylný člen dále přejde do infikovaného stavu kontaktem s kterýmkoli z těchto prvků ${displaystyle n_ {I}}$ infikovaných členů, takže infekce se vyskytuje s rychlostí ${displaystyle alpha n_ {I}}$. Daný člen infikovaného stavu se zotaví bez závislosti na kterémkoli ze tří stavů, což je určeno sazbou β. Vzhledem k tomuto základnímu schématu je možné sestrojit následující nelineární systém.

Trajektorie získané přímou metodou epidemie SIR (oranžové čáry) se superponovaným numerickým řešením (černé čáry).

${displaystyle {frac {dn_ {S}} {dt}} = - {frac {alpha n_ {S}} {N}} n_ {I}}$,

${displaystyle {frac {dn_ {I}} {dt}} = vlevo ({frac {alpha n_ {S}} {N}} - eta ight) n_ {I}}$,

${displaystyle {frac {dn_ {R}} {dt}} = eta n_ {I}}$.

Tento systém nemá žádné analytické řešení. Jedním z přístupů by mohlo být použít Gillespieho algoritmus, mnohokrát simulovat systém a použít regresní techniku, jako jsou nejmenší čtverce, aby se vešel na polynom přes všechny trajektorie. Jak se počet trajektorií zvyšuje, taková polynomiální regrese se bude asymptoticky chovat jako numerické řešení (černé čáry). Kromě odhadu řešení neřešitelného problému, jako je epidemie SIR, umožňuje stochastická povaha každé trajektorie vypočítat jiné statistiky než ${displaystyle mathrm {E} [n | t]}$Při spuštění níže uvedeného skriptu se někdy dosáhne realizace epidemie SIR, která se drasticky liší od numerického řešení. Například když jsou všechny osoby vyléčeny (náhodou) velmi brzy nebo velmi pozdě.

Trajektorie prezentované na výše uvedeném obrázku byly získány pomocí následující implementace přímé metody v Pythonu spolu s modelovou třídou, jejíž členové interagují s mechanizmem přímé metody a provádějí obecné simulace s teoriemi, které jsou základem Gillespieho algoritmu. Ty jsou uvedeny výše. Řešitel ODE navíc z SciPy je vyvolána k získání numerického řešení systému diferenciálních rovnic, tj. reprezentace ${displaystyle mathrm {E} [n | t]}$.

z matplotlib import pyplotz numpy import linspacez scipy.integrovat import odeintPočet počátečních druhů a doba pobytuinitital_conditions = {    „s“: [480],    "já": [20],    „r“: [0],    "čas": [0.0],}# funkce sklonusklony = {    0: lambda d: 2.0 * d[„s“][-1] * d["já"][-1] / 500,    1: lambda d: 1.0 * d["já"][-1],}# změna druhu pro každý sklonstechiometrie = {    0: {„s“: -1, "já": 1, „r“: 0},    1: {„s“: 0, "já": -1, „r“: 1},}# vytvořit instanci kontejneru modelu epidemie SSAepidemický = SSAModel(    initital_conditions,    sklony,    stechiometrie)# vytvořit instanci kontejneru SSA s modelemepidemic_generator = SSA(epidemický)# udělejte pěknou, velkou postavupyplot.postava(figsize=(10,10), dpi=500)# Vytvořte dílčí graf pro vnímavé, infikované a uzdravené jedinceaxes_s = pyplot.subplot(311)axes_s.set_ylabel(„vnímaví jedinci“)axes_i = pyplot.subplot(312)axes_i.set_ylabel(„infikovaní jedinci“)axes_r = pyplot.subplot(313)axes_r.set_ylabel("uzdravení jednotlivci")axes_r.set_xlabel("čas (libovolné jednotky)")# simulujte a zakreslete 30 trajektoriítrajektorie = 0pro trajektorie v epidemic_generator.Přímo():    axes_s.spiknutí(trajektorie["čas"], trajektorie[„s“], barva="oranžový")    axes_i.spiknutí(trajektorie["čas"], trajektorie["já"], barva="oranžový")    axes_r.spiknutí(trajektorie["čas"], trajektorie[„r“], barva="oranžový")    trajektorie += 1    -li trajektorie == 30:        přestávka# numerické řešení pomocí obyčejného solversiru diferenciální rovnicet = linspace(0, 14, počet=200)y0 = (480, 20, 0)alfa = 2.0beta = 1.0def diferenciální_SIR(n_SIR, t, alfa, beta):    dS_dt = -alfa * n_SIR[0] * n_SIR[1] / 500    dI_dt = ((alfa * n_SIR[0] / 500) - beta) * n_SIR[1]    dR_dt = beta * n_SIR[1]    vrátit se dS_dt, dI_dt, dR_dtřešení = odeint(diferenciální_SIR, y0, t, args=(alfa, beta))řešení = [[řádek[i] pro řádek v řešení] pro i v rozsah(3)]# vykreslit numerické řešeníaxes_s.spiknutí(t, řešení[0], barva="Černá")axes_i.spiknutí(t, řešení[1], barva="Černá")axes_r.spiknutí(t, řešení[2], barva="Černá")pyplot.ukázat()

Další čtení

• Gillespie, Daniel T. (1977). "Přesná stochastická simulace spojených chemických reakcí". The Journal of Physical Chemistry. 81 (25): 2340–2361. CiteSeerX  10.1.1.704.7634. doi:10.1021 / j100540a008.
• Gillespie, Daniel T. (1976). "Obecná metoda pro numerickou simulaci stochastického časového vývoje spojených chemických reakcí". Journal of Computational Physics. 22 (4): 403–434. Bibcode:1976JCoPh..22..403G. doi:10.1016/0021-9991(76)90041-3.
• Gibson, Michael A .; Bruck, Jehoshua (2000). „Efektivní přesná stochastická simulace chemických systémů s mnoha druhy a mnoha kanály“ (PDF). Journal of Physical Chemistry A. 104 (9): 1876–1889. Bibcode:2000JPCA..104,1876G. doi:10.1021 / jp993732q.
• Doob, Jacob L. (1942). „Témata v teorii Markoffových řetězců“. Transakce Americké matematické společnosti. 52 (1): 37–64. doi:10.1090 / S0002-9947-1942-0006633-7. JSTOR  1990152.
• Doob, Jacob L. (1945). "Markoffovy řetězy - nespočetné množství případů". Transakce Americké matematické společnosti. 58 (3): 455–473. doi:10.2307/1990339. JSTOR  1990339.
• Press, William H .; Teukolsky, Saul A .; Vetterling, William T .; Flannery, Brian P. (2007). „Oddíl 17.7. Stochastická simulace sítí chemických reakcí“. Numerické recepty: Umění vědecké práce na počítači (3. vyd.). New York, NY: Cambridge University Press. ISBN  978-0-521-88068-8.
• Kolmogorov, Andrey N. (1931). „Über die analytischen Methoden in der Wahrscheinlichkeitsrechnung“ [O analytických metodách v teorii pravděpodobnosti]. Mathematische Annalen. 104: 415–458. doi:10.1007 / BF01457949. S2CID  119439925.
• Feller, Willy (1940). „O integro-diferenciálních rovnicích čistě diskontinuálních Markoffových procesů“. Transakce Americké matematické společnosti. 48 (3): 4885–15. doi:10.2307/1990095. JSTOR  1970064.
• Kendall, David G. (1950). „Umělá realizace jednoduchého procesu„ narození a smrti “. Journal of the Royal Statistical Society, Series B. 12 (1): 116–119. JSTOR  2983837.
• Bartlett, Maurice S. (1953). "Stochastické procesy nebo statistika změn". Journal of the Royal Statistical Society, Series C. 2 (1): 44–64. JSTOR  2985327.
• Rathinam, Muruhan; Petzold, Linda R.; Cao, Yang; Gillespie, Daniel T. (2003). „Tuhost ve stochastických chemicky reagujících systémech: implicitní metoda tau-skákání“. Journal of Chemical Physics. 119 (24): 12784–12794. Bibcode:2003JChPh.11912784R. doi:10.1063/1.1627296.
• Sinitsyn, Nikolai A .; Hengartner, Nicolas; Nemenman, Ilya (2009). „Adiabatická hrubozrnnost a simulace stochastických biochemických sítí“ (PDF). Sborník Národní akademie věd Spojených států amerických. 106 (20): 10546–10551. Bibcode:2009PNAS..10610546S. doi:10.1073 / pnas.0809340106. PMC  2705573. PMID  19525397. Archivovány od originál (PDF) dne 14.07.2011.
• Salis, Howard; Kaznessis, Yiannis N. (2005). "Přesná hybridní stochastická simulace systému spojených chemických nebo biochemických reakcí". Journal of Chemical Physics. 122 (5): 054103. Bibcode:2005JChPh.122e4103S. doi:10.1063/1.1835951. PMID  15740306.
• (Slepoy Thompson Plimpton 2008): Slepoy, Alexander; Thompson, Aidan P .; Plimpton, Steven J. (2008). "Kinetický algoritmus Monte Carlo s konstantním časem pro simulaci velkých sítí biochemických reakcí". Journal of Chemical Physics. 128 (20): 205101. Bibcode:2008JChPh.128t5101S. doi:10.1063/1.2919546. PMID  18513044.
• (Bratsun et al. 2005): Bratsun, Dmitri; Volfson, Dmitri; Hasty, Jeff; Tsimring, Lev S. (2005). „Zpožděné stochastické oscilace v regulaci genů“. Sborník Národní akademie věd Spojených států amerických. 102 (41): 14593–8. Bibcode:2005PNAS..10214593B. doi:10.1073 / pnas.0503858102. PMC  1253555. PMID  16199522.
• (Barrio et al. 2006): Barrio, Manuel; Burrage, Kevin; Leier, André; Tian, ​​Tianhai (2006). "Oscilační regulace z hes1: Diskrétní stochastické modelování a simulace zpoždění ". PLOS výpočetní biologie. 2 (9): 1017. Bibcode:2006PLSCB ... 2..117B. doi:10.1371 / journal.pcbi.0020117. PMC  1560403. PMID  16965175.
• (Cai 2007): Cai, Xiaodong (2007). "Přesná stochastická simulace spojených chemických reakcí se zpožděním". Journal of Chemical Physics. 126 (12): 124108. Bibcode:2007JChPh.126l4108C. doi:10.1063/1.2710253. PMID  17411109.
• (Barnes Chu 2010): Barnes, David J .; Chu, Dominique (2010). Úvod do modelování pro biologické vědy. Springer Verlag.
• (Ramaswamy González-Segredo Sbalzarini 2009): Ramaswamy, Rajesh; González-Segredo, Nélido; Sbalzarini, Ivo F. (2009). "Nová třída vysoce efektivních přesných stochastických simulačních algoritmů pro sítě chemických reakcí". Journal of Chemical Physics. 130 (24): 244104. arXiv:0906.1992. Bibcode:2009JChPh.130x4104R. doi:10.1063/1.3154624. PMID  19566139. S2CID  4952205.
• (Ramaswamy Sbalzarini 2010): Ramaswamy, Rajesh; Sbalzarini, Ivo F. (2010). „Varianta částečného sklonu stochastického simulačního algoritmu odmítnutí složení pro sítě chemických reakcí“ (PDF). Journal of Chemical Physics. 132 (4): 044102. Bibcode:2010JChPh.132d4102R. doi:10.1063/1.3297948. PMID  20113014.
• (Indurkhya Beal 2010): Indurkhya, Sagar; Beal, Jacob S. (2005). Isalan, Mark (ed.). „Faktor reakce a bipartitní aktualizační grafy zrychlují Gillespieho algoritmus pro biochemické systémy velkého rozsahu“. PLOS ONE. 5 (1): e8125. Bibcode:2010PLoSO ... 5.8125I. doi:10.1371 / journal.pone.0008125. PMC  2798956. PMID  20066048.
• (Ramaswamy Sbalzarini 2011): Ramaswamy, Rajesh; Sbalzarini, Ivo F. (2011). „Formulace stochastického simulačního algoritmu pro sítě chemických reakcí s částečným sklonem se zpožděním“ (PDF). Journal of Chemical Physics. 134 (1): 014106. Bibcode:2011JChPh.134a4106R. doi:10.1063/1.3521496. PMID  21218996.
• (Yates Klingbeil 2013): Yates, Christian A .; Klingbeil, Guido (2013). "Recyklace náhodných čísel ve stochastickém simulačním algoritmu". Roční přehled fyzikální chemie. 58 (9): 094103. Bibcode:2013JChPh.138i4103Y. doi:10.1063/1.4792207. PMID  23485273.
• Gillespie, Daniel T. (2007). „Stochastická simulace chemické kinetiky“. The Journal of Physical Chemistry. 58: 35–55. doi:10.1146 / annurev.physchem.58.032806.104637.