Алгоритм Гиллеспи - Gillespie algorithm

В теории вероятностей , то Гиллеспи алгоритм (или изредка алгоритм Дуб-Gillespie ) генерирует статистически правильную траекторию (возможное решение) продукта в виде стохастической системы уравнений , для которых скорость реакции известна. Он был создан Джозефом Л. Дубом и другими (около 1945 г.), представлен Дэном Гиллеспи в 1976 г. и популяризирован в 1977 г. в статье, где он использует его для эффективного и точного моделирования химических или биохимических систем реакций с использованием ограниченной вычислительной мощности (см. стохастическое моделирование ). Поскольку компьютеры стали быстрее, этот алгоритм использовался для моделирования все более сложных систем. Алгоритм особенно полезен для моделирования реакций внутри клеток, где количество реагентов невелико, а отслеживание положения и поведения отдельных молекул возможно с вычислительной точки зрения. Математически это вариант динамического метода Монте-Карло, аналогичный кинетическим методам Монте-Карло . Он широко используется в вычислительной системной биологии .

История

Процесс, который привел к созданию алгоритма, включает несколько важных шагов. В 1931 году Андрей Колмогоров ввел дифференциальные уравнения, соответствующие эволюции во времени случайных скачкообразных процессов, которые сегодня известны как уравнения Колмогорова (процесс скачков Маркова) (упрощенная версия в естествознании известна как основное уравнение ). В 1940 году Уильям Феллер нашел условия, при которых уравнения Колмогорова допускали (собственные) вероятности в качестве решений. В своей теореме I (работа 1940 г.) он устанавливает, что время до следующего скачка было экспоненциально распределено, и вероятность следующего события пропорциональна скорости. Таким образом, он установил связь уравнений Колмогорова со случайными процессами . Позже Дуб (1942, 1945) расширил решения Феллера за пределы случая чисто скачкообразных процессов. Этот метод был реализован на компьютерах Дэвидом Джорджем Кендаллом (1950) с использованием компьютера Manchester Mark 1 и позже использован Морисом С. Бартлеттом (1953) в его исследованиях вспышек эпидемий. Гиллеспи (1977) получает алгоритм другим способом, используя физический аргумент.

Идея алгоритма

Традиционные непрерывные и детерминированные уравнения скорости биохимии не позволяют точно предсказать клеточные реакции, поскольку они полагаются на массовые реакции, которые требуют взаимодействия миллионов молекул. Обычно они моделируются как набор связанных обыкновенных дифференциальных уравнений. Напротив, алгоритм Гиллеспи позволяет дискретное и стохастическое моделирование системы с небольшим количеством реагентов, поскольку каждая реакция моделируется явно. Траектория, соответствующая единственному моделированию Гиллеспи, представляет собой точную выборку из функции массы вероятности, которая является решением основного уравнения .

Физическая основа алгоритма - столкновение молекул внутри реакционного сосуда. Предполагается, что столкновения часты, но столкновения с правильной ориентацией и энергией нечасты. Следовательно, все реакции в рамках концепции Гиллеспи должны включать не более двух молекул. Предполагается, что реакции с участием трех молекул крайне редки и моделируются как последовательность бинарных реакций. Также предполагается, что реакционная среда хорошо перемешана.

Алгоритм

В недавнем обзоре (Gillespie, 2007) представлены три различных, но эквивалентных формулировки; прямые методы, методы первой реакции и методы первого семейства, причем первые два являются частными случаями второго. Формулировка методов прямой и первой реакции сосредоточена на выполнении обычных шагов обращения Монте-Карло на так называемой «фундаментальной предпосылке стохастической химической кинетики», которая математически является функцией

${\ displaystyle p (\ tau, j | {\ boldsymbol {x}}, t) = a_ {j} ({\ boldsymbol {x}}) \ exp (\ tau \ sum _ {j} a_ {j} ( {\ boldsymbol {x}}))}$,

где каждый из членов является функцией склонности элементарной реакции, аргументом которой является вектор количества видов. Параметром является время следующей реакции (или времени пребывания), и , конечно , в настоящее время. Перефразируя Гиллеспи, это выражение читается как «вероятность, при условии , что следующая реакция системы произойдет в бесконечно малом интервале времени и будет иметь стехиометрию, соответствующую th реакции». Эта формулировка предоставляет окно для методов прямой и первой реакции, подразумевая, что это экспоненциально распределенная случайная величина и «статистически независимая целочисленная случайная величина с точечными вероятностями ». ${\ displaystyle a}$${\ displaystyle {\ boldsymbol {x}}}$${\ Displaystyle \ тау}$${\ displaystyle t}$${\ displaystyle {\ boldsymbol {X}} (t) = {\ boldsymbol {x}}}$${\ Displaystyle [т + \ тау, т + \ тау + д \ тау]}$${\ displaystyle j}$${\ Displaystyle \ тау}$${\ displaystyle j}$${\ displaystyle a_ {j} ({\ boldsymbol {x}}) / \ sum _ {j} a_ {j} ({\ boldsymbol {x}})}$

Таким образом, метод генерации Монте-Карло состоит в том, чтобы просто нарисовать два псевдослучайных числа и далее и вычислить ${\ displaystyle r_ {1}}$${\ displaystyle r_ {2}}$${\ displaystyle [0,1]}$

${\ displaystyle \ tau = {\ frac {1} {\ sum _ {j} a_ {j} ({\ boldsymbol {x}})}} \ log \ left ({\ frac {1} {r_ {1}) }}\Правильно)}$,

а также

${\ displaystyle j =}$наименьшее целое удовлетворительное .${\ displaystyle \ sum _ {j '= 1} ^ {j} a_ {j'} ({\ boldsymbol {x}})> r_ {2} \ sum _ {j} a_ {j} ({\ boldsymbol { Икс}})}$

Используя этот метод генерации для времени пребывания и следующей реакции, алгоритм прямого метода сформулирован Гиллеспи как

1. Initialize the time ${\displaystyle t=t_{0}}$ and the system's state ${\displaystyle {\boldsymbol {x}}={\boldsymbol {x}}_{0}}$
2. With the system in state ${\displaystyle {\boldsymbol {x}}}$ at time ${\displaystyle t}$, evaluate all the ${\displaystyle a_{j}({\boldsymbol {x}})}$ and their sum ${\displaystyle \sum _{j}a_{j}({\boldsymbol {x}})}$
3. Effect the next reaction by replacing ${\displaystyle t\leftarrow t+\tau }$ and ${\displaystyle {\boldsymbol {x}}\leftarrow {\boldsymbol {x}}+\nu _{j}}$
4. Record ${\displaystyle ({\boldsymbol {x}},t)}$ as desired. Return to step 1, or else end the simulation.

Это семейство алгоритмов требует больших вычислительных ресурсов, поэтому существует множество модификаций и адаптаций, в том числе метод следующей реакции (Гибсон и Брук), тау-скачок , а также гибридные методы, в которых большое количество реагентов моделируется с детерминированным поведением. Адаптированные методы, как правило, ставят под угрозу точность теории, лежащей в основе алгоритма, поскольку он связан с основным уравнением, но предлагают разумные реализации для значительно улучшенных временных масштабов. Вычислительная стоимость точных версий алгоритма определяется классом связи реакционной сети. В слабосвязанных сетях количество реакций, на которые влияет любая другая реакция, ограничено небольшой константой. В сильно связанных сетях запуск одной реакции в принципе может повлиять на все другие реакции. Была разработана точная версия алгоритма с постоянным масштабированием для слабосвязанных сетей, позволяющая эффективно моделировать системы с очень большим количеством каналов реакции (Slepoy Thompson Plimpton 2008). Обобщенный алгоритм Гиллеспи, который учитывает немарковские свойства случайных биохимических событий с задержкой, был разработан Bratsun et al. 2005 и независимо Barrio et al. 2006 г., а также (Cai 2007). Подробности см. В цитируемых ниже статьях.

Составы с частичной предрасположенностью, независимо разработанные Рамасвами и соавт. (2009, 2010) и Indurkhya and Beal (2010) доступны для построения семейства точных версий алгоритма, вычислительная стоимость которого пропорциональна количеству химических соединений в сети, а не (большему) количеству реакций. Эти формулировки могут снизить вычислительные затраты до масштабирования с постоянным временем для слабосвязанных сетей и масштабирования не более чем линейно с числом видов для сильно связанных сетей. Также был предложен вариант частичной склонности обобщенного алгоритма Гиллеспи для реакций с задержками (Ramaswamy Sbalzarini 2011). Использование методов частичной склонности ограничивается элементарными химическими реакциями, т. Е. Реакциями максимум с двумя различными реагентами. Каждую неэлементарную химическую реакцию можно эквивалентно разложить на набор элементарных за счет линейного (в порядке реакции) увеличения размера сети.

Рассмотрим следующую объектно-ориентированную реализацию методов прямой и первой реакции, которые содержатся в классе Python 3:

from math import log
from random import random

class SSA:
    """Container for SSAs"""

    def __init__(self, model):
        """Initialize container with model and pseudorandom number generator"""
        self.model = model
        self.random = random

    def direct(self):
        """Indefinite generator of direct-method trajectories"""
        while True:
            while not self.model.exit():

                # evaluate weights and partition
                weights = [
                    (rxn, sto, pro(self.model))
                    for (rxn, sto, pro) in self.model.reactions
                ]
                partition = sum(w[-1] for w in weights)

                # evaluate sojourn time (MC step 1)
                sojourn = log(1.0 / self.random()) / partition
                self.model["time"].append(self.model["time"][-1] + sojourn)

                # evaluate the reaction (MC step 2)
                partition = partition * self.random()
                while partition >= 0.0:
                    rxn, sto, pro = weights.pop(0)
                    partition -= pro
                for species, delta in sto.items():
                    self.model[species].append(self.model[species][-1] + delta)

                self.model.curate()
            yield self.model
            self.model.reset()

    def first_reaction(self):
        """Indefinite generator of 1st-reaction trajectories"""
        while True:
            while not self.model.exit():

                # evaluate next reaction times
                times = [
                    (
                        log(
                            1.0 / self.random()
                        ) / pro(self.model),
                        sto
                    )
                    for (rxn, sto, pro) in self.model.reactions
                ]
                times.sort()

                # evaluate reaction time
                self.model["time"].append(
                    self.model["time"][-1] + times[0][0]
                )

                # evaluate reaction
                for species, delta in times[0][1].items():
                    self.model[species].append(
                        self.model[species][-1] + delta
                    )

                self.model.curate()
            yield self.model
            self.model.reset()

Члены directи first_reactionявляются неопределенными генераторами, то есть они будут непрерывно создавать траектории, каждая из которых является полной симуляцией системы, в цикле, пока сигнал не разорвет этот цикл. Важное предостережение для любой реализации этого алгоритма заключается в том, что он не предоставляет логики для различения возможных и невозможных элементарных событий. Таким образом, чтобы фактически использовать любую версию алгоритма в цикле и получить некоторое количество траекторий для анализа, этот класс требует, чтобы модель была передана ему при создании экземпляра. В дополнение к видам населения на жилье и функции склонности, модель имеет три метода , что Gillespie алгоритм делает вызовы: curate, exitи reset. Назначение этих методов указано ниже и по существу определяет, когда алгоритм завершается. Обоснование введения класса модели в основном состоит в том, чтобы избежать смешения кинетических свойств конкретного моделируемого процесса с внутренней логикой алгоритма Гиллеспи.

Ниже модель представляет собой словарь подкласса с общественными членами curate, exitи resetкоторые , соответственно ,

• определять, какие реакции допустимы, а какие нет в конце каждой итерации заданной траектории;
• вернуть, Trueесли больше нет возможных реакций;
• вернуть хеш-значения модели к их исходным значениям (то есть их начальным условиям) в конце заданной траектории.

class SSAModel(dict):
    """Container for SSA model"""

    def __init__(
        self, initial_conditions, propensities, stoichiometry
    ):
        """Initialize model"""
        super().__init__(**initial_conditions)
        self.reactions = list()
        self.excluded_reactions = list()
        for reaction,propensity in propensities.items():
            if propensity(self) == 0.0:
                self.excluded_reactions.append(
                    (
                        reaction,
                        stoichiometry[reaction],
                        propensity
                    )
                )
            else:
                self.reactions.append(
                    (
                        reaction,
                        stoichiometry[reaction],
                        propensity
                    )
                )

    def exit(self):
        """Return True to break out of trajectory"""

        # return True if no more reactions
        if len(self.reactions) == 0: return True

        # return False if there are more reactions
        else: return False

    def curate(self):
        """Validate and invalidate model reactions"""
        
        # evaluate possible reactions
        reactions = []
        while len(self.reactions) > 0:
            reaction = self.reactions.pop()
            if reaction[2](self) == 0:
                self.excluded_reactions.append(reaction)
            else:
                reactions.append(reaction)
        self.reactions = reactions

        # evaluate impossible reactions
        excluded_reactions = []
        while len(self.excluded_reactions) > 0:
            reaction = self.excluded_reactions.pop()
            if reaction[2](self) > 0:
                self.reactions.append(reaction)
            else:
                excluded_reactions.append(reaction)
        self.excluded_reactions = excluded_reactions

    def reset(self):
        """Clear the trajectory"""

        # reset species to initial conditions
        for key in self: del self[key][1:]

        # reset reactions per initial conditions
        self.curate()

Примеры

Обратимое связывание A и B с образованием димеров AB

Простой пример может помочь объяснить, как работает алгоритм Гиллеспи. Рассмотрим систему молекул двух типов, A и B . В этой системе, и Б обратимо связываются вместе с образованием AB димеры , такие , что две реакции: либо А и В реагируют обратимо с образованием AB димер, или АВ димер диссоциирует на A и B . Константа скорости реакции для данной одиночной молекулы A, реагирующей с данной единственной молекулой B , равна , а скорость реакции разрушения димера AB равна . ${\ Displaystyle к _ {\ mathrm {D}}}$${\ Displaystyle к _ {\ mathrm {B}}}$

Если в момент времени t имеется одна молекула каждого типа, тогда скорость образования димера равна , тогда как если есть молекулы типа A и молекулы типа B , скорость образования димера равна . Если есть димеры, то скорость диссоциации димеров равна . ${\ Displaystyle к _ {\ mathrm {D}}}$${\ Displaystyle п _ {\ mathrm {A}}}$${\ Displaystyle п _ {\ mathrm {B}}}$${\ Displaystyle к _ {\ mathrm {D}} n _ {\ mathrm {A}} n _ {\ mathrm {B}}}$${\ Displaystyle п _ {\ mathrm {AB}}}$${\ Displaystyle к _ {\ mathrm {B}} п _ {\ mathrm {AB}}}$

Общая скорость реакции в момент времени t определяется выражением ${\ displaystyle R _ {\ mathrm {TOT}}}$

${\ displaystyle R _ {\ mathrm {TOT}} = k _ {\ mathrm {D}} n _ {\ mathrm {A}} n _ {\ mathrm {B}} + k _ {\ mathrm {B}} n _ {\ mathrm { AB}}}$

Итак, мы описали простую модель с двумя реакциями. Это определение не зависит от алгоритма Гиллеспи. Теперь мы опишем, как применить алгоритм Гиллеспи к этой системе.

В алгоритме мы продвигаемся вперед во времени в два этапа: вычисляем время до следующей реакции и определяем, какая из возможных реакций будет следующей. Предполагается, что реакции являются полностью случайными, поэтому, если скорость реакции в момент времени t равна , то время δ t , пока не произойдет следующая реакция, является случайным числом, полученным из экспоненциальной функции распределения со средним значением . Таким образом, мы продвигаем время от t к t + δ t . ${\ displaystyle R _ {\ mathrm {TOT}}}$${\ displaystyle 1 / R _ {\ mathrm {TOT}}}$

График зависимости числа молекул A (черная кривая) и димеров AB от времени. Поскольку мы начали с 10 молекулами A и B в момент времени t = 0, количество молекул B всегда равно количеству молекул A , поэтому оно не отображается.

Вероятность того, что эта реакция представляет собой связывание молекулы A с молекулой B , - это просто доля общей скорости, обусловленной этим типом реакции, т. Е.

вероятность того, что реакция ${\ displaystyle P ({\ ce {{A} + B -> AB}}) = k_ {D} n_ {A} n_ {B} / R _ {{\ ce {TOT}}}}$

Вероятность того, что следующей реакцией будет диссоциация димера AB, составляет всего 1 минус эта вероятность . Итак, с этими двумя вероятностями мы либо формируем димер, уменьшая и на единицу, и увеличивая на единицу, либо диссоциируем димер и увеличиваем и на единицу и уменьшаем на единицу. ${\ Displaystyle п _ {\ mathrm {A}}}$${\ Displaystyle п _ {\ mathrm {B}}}$${\ Displaystyle п _ {\ mathrm {AB}}}$${\ Displaystyle п _ {\ mathrm {A}}}$${\ Displaystyle п _ {\ mathrm {B}}}$${\ Displaystyle п _ {\ mathrm {AB}}}$

Теперь мы оба опередили время до t + δ t и выполнили одну реакцию. Алгоритм Гиллеспи просто повторяет эти два шага столько раз, сколько необходимо для моделирования системы, сколько мы хотим (т. Е. Столько реакций). Результат моделирования Гиллеспи, которое начинается с и в t = 0, где и , показан справа. Для этих значений параметров в среднем имеется 8 димеров и 2 димера A и B, но из-за небольшого числа молекул колебания вокруг этих значений велики. Алгоритм Гиллеспи часто используется для изучения систем, в которых эти колебания важны. ${\ Displaystyle п _ {\ mathrm {A}} = n _ {\ mathrm {B}} = 10}$${\ Displaystyle п _ {\ mathrm {AB}} = 0}$${\ Displaystyle к _ {\ mathrm {D}} = 2}$${\ Displaystyle к _ {\ mathrm {B}} = 1}$${\ Displaystyle п _ {\ mathrm {AB}}}$

Это был простой пример с двумя реакциями. Таким же образом обрабатываются более сложные системы с большим количеством реакций. Все скорости реакции должны быть рассчитаны на каждом временном шаге, и одна из них выбирается с вероятностью, равной ее дробному вкладу в скорость. Затем время продвигается вперед, как в этом примере.

Эпидемия SIR без жизненной динамики

Модель SIR - это классическое биологическое описание того, как определенные болезни проникают в популяцию фиксированного размера. В простейшей форме есть члены популяции, при этом каждый член может находиться в одном из трех состояний - восприимчивый, инфицированный или выздоровевший - в любой момент времени, и каждый такой член необратимо переходит через эти состояния в соответствии с ориентированным графом ниже. . Мы можем обозначить количество восприимчивых членов как , количество инфицированных членов как и количество вылеченных членов как . Поэтому на любой момент времени. ${\ displaystyle N}$${\ displaystyle n_ {S}}$${\ displaystyle n_ {I}}$${\ displaystyle n_ {R}}$${\ displaystyle N = n_ {S} + n_ {I} + n_ {R}}$

Кроме того, данный восприимчивый член переходит в инфицированное состояние, вступая в контакт с любым из инфицированных членов, поэтому заражение происходит быстро . Данный член зараженного состояния выздоравливает независимо от любого из трех состояний, что определяется показателем β . По этой базовой схеме можно построить следующую нелинейную систему. ${\ displaystyle n_ {I}}$${\ displaystyle \ alpha n_ {I}}$

Траектории, полученные прямым методом эпидемии SIR (оранжевые линии) с наложенным числовым решением (черные линии).

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

${\ displaystyle {\ frac {dn_ {I}} {dt}} = \ left ({\ frac {\ alpha n_ {S}} {N}} - \ beta \ right) n_ {I}}$,

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

Эта система не имеет аналитического решения. Один из подходов может заключаться в использовании алгоритма Гиллеспи, многократном моделировании системы и использовании метода регрессии, такого как метод наименьших квадратов, для подбора полинома по всем траекториям. По мере увеличения количества траекторий такая полиномиальная регрессия будет асимптотически вести себя как численное решение (черные линии). В дополнение к оценке решения трудноразрешимой проблемы, такой как эпидемия SIR, стохастический характер каждой траектории позволяет вычислять статистику, отличную от . При запуске приведенного ниже сценария иногда получаются реализации эпидемии SIR, которые резко отличаются от численного решения. Например, когда все люди излечиваются (случайно) очень рано или очень поздно. ${\ Displaystyle \ mathrm {E} [п | т]}$

Траектории, представленные на приведенном выше рисунке, были получены с помощью следующей реализации Python прямого метода вместе с классом модели, члены которого взаимодействуют с механизмом прямого метода для выполнения общего моделирования с теориями, лежащими в основе алгоритма Гиллеспи. Они представлены выше. Кроме того, решение ОДА из SciPy вызывается для получения численного решения системы дифференциальных уравнений, т.е. представления . ${\ Displaystyle \ mathrm {E} [п | т]}$

from matplotlib import pyplot
from numpy import linspace
from scipy.integrate import odeint

# initial species counts and sojourn times
initital_conditions = {
    "s": [480],
    "i": [20],
    "r": [0],
    "time": [0.0],
}

# propensity functions
propensities = {
    0: lambda d: 2.0 * d["s"][-1] * d["i"][-1] / 500,
    1: lambda d: 1.0 * d["i"][-1],
}

# change in species for each propensity
stoichiometry = {
    0: {"s": -1, "i": 1, "r": 0},
    1: {"s": 0, "i": -1, "r": 1},
}

# instantiate the epidemic SSA model container
epidemic = SSAModel(
    initital_conditions,
    propensities,
    stoichiometry
)

# instantiate the SSA container with model
epidemic_generator = SSA(epidemic)

# make a nice, big figure
pyplot.figure(figsize=(10,10), dpi=500)

# make a subplot for the susceptible, infected and recovered individuals
axes_s = pyplot.subplot(311)
axes_s.set_ylabel("susceptible individuals")

axes_i = pyplot.subplot(312)
axes_i.set_ylabel("infected individuals")

axes_r = pyplot.subplot(313)
axes_r.set_ylabel("recovered individuals")
axes_r.set_xlabel("time (arbitrary units)")

# simulate and plot 30 trajectories
trajectories = 0
for trajectory in epidemic_generator.direct():
    axes_s.plot(trajectory["time"], trajectory["s"], color="orange")
    axes_i.plot(trajectory["time"], trajectory["i"], color="orange")
    axes_r.plot(trajectory["time"], trajectory["r"], color="orange")
    trajectories += 1
    if trajectories == 30:
        break

# numerical solution using an ordinary differential equation solversir
t = linspace(0, 14, num=200)
y0 = (480, 20, 0)
alpha = 2.0
beta = 1.0

def differential_SIR(n_SIR, t, alpha, beta):
    dS_dt = -alpha * n_SIR[0] * n_SIR[1] / 500
    dI_dt = ((alpha * n_SIR[0] / 500) - beta) * n_SIR[1]
    dR_dt = beta * n_SIR[1]
    return dS_dt, dI_dt, dR_dt

solution = odeint(differential_SIR, y0, t, args=(alpha, beta))
solution = [[row[i] for row in solution] for i in range(3)]

# plot numerical solution
axes_s.plot(t, solution[0], color="black")
axes_i.plot(t, solution[1], color="black")
axes_r.plot(t, solution[2], color="black")

pyplot.show()

дальнейшее чтение

• Гиллеспи, Дэниел Т. (1977). «Точное стохастическое моделирование связанных химических реакций». Журнал физической химии . 81 (25): 2340–2361. CiteSeerX  10.1.1.704.7634 . DOI : 10.1021 / j100540a008 .
• Гиллеспи, Дэниел Т. (1976). «Общий метод численного моделирования стохастической временной эволюции связанных химических реакций». Журнал вычислительной физики . 22 (4): 403–434. Bibcode : 1976JCoPh..22..403G . DOI : 10.1016 / 0021-9991 (76) 90041-3 .
• Гибсон, Майкл А .; Брук, Иошуа (2000). «Эффективное точное стохастическое моделирование химических систем со многими видами и многими каналами» (PDF) . Журнал физической химии . 104 (9): 1876–1889. Bibcode : 2000JPCA..104.1876G . DOI : 10.1021 / jp993732q .
• Дуб, Джейкоб Л. (1942). «Вопросы теории цепей Маркова» . Труды Американского математического общества . 52 (1): 37–64. DOI : 10.1090 / S0002-9947-1942-0006633-7 . JSTOR  1990152 .
• Дуб, Джейкоб Л. (1945). «Цепи Маркова - счетный случай». Труды Американского математического общества . 58 (3): 455–473. DOI : 10.2307 / 1990339 . JSTOR  1990339 .
• Press, William H .; Teukolsky, Saul A .; Веттерлинг, Уильям Т .; Фланнери, Брайан П. (2007). «Раздел 17.7. Стохастическое моделирование сетей химических реакций» . Числовые рецепты: искусство научных вычислений (3-е изд.). Нью-Йорк, Нью-Йорк: Издательство Кембриджского университета. ISBN 978-0-521-88068-8.
• Колмогоров, Андрей Н. (1931). "Über die analytischen Methoden in der Wahrscheinlichkeitsrechnung" [Об аналитических методах в теории вероятностей]. Mathematische Annalen . 104 : 415–458. DOI : 10.1007 / BF01457949 . S2CID  119439925 .
• Феллер, Вилли (1940). "Об интегродифференциальных уравнениях чисто разрывных марковских процессов" . Труды Американского математического общества . 48 (3): 4885–15. DOI : 10.2307 / 1990095 . JSTOR  1970064 .
• Кендалл, Дэвид Г. (1950). «Искусственная реализация простого процесса« рождения и смерти »». Журнал Королевского статистического общества, Series B . 12 (1): 116–119. JSTOR  2983837 .
• Бартлетт, Морис С. (1953). «Случайные процессы или статистика изменений». Журнал Королевского статистического общества, серия C . 2 (1): 44–64. JSTOR  2985327 .
• Ратинам, Мурухан; Петцольд, Линда Р .; Цао, Ян; Гиллеспи, Дэниел Т. (2003). «Жесткость в стохастических химически реагирующих системах: неявный метод тау-прыжка». Журнал химической физики . 119 (24): 12784–12794. Bibcode : 2003JChPh.11912784R . DOI : 10.1063 / 1.1627296 .
• Синицын, Николай А .; Хенгартнер, Николас; Неменман, Илья (2009). «Адиабатическое крупнозернистое моделирование и моделирование стохастических биохимических сетей» . Труды Национальной академии наук Соединенных Штатов Америки . 106 (20): 10546–10551. Bibcode : 2009PNAS..10610546S . DOI : 10.1073 / pnas.0809340106 . PMC  2705573 . PMID  19525397 .
• Салис, Ховард; Казнессис, Яннис Н. (2005). «Точное гибридное стохастическое моделирование системы связанных химических или биохимических реакций». Журнал химической физики . 122 (5): 054103. Bibcode : 2005JChPh.122e4103S . DOI : 10.1063 / 1.1835951 . PMID  15740306 .
• (Слепой Томпсон Плимптон 2008): Слепой, Александр; Thompson, Aidan P .; Плимптон, Стивен Дж. (2008). «Кинетический алгоритм Монте-Карло с постоянным временем для моделирования больших сетей биохимических реакций». Журнал химической физики . 128 (20): 205101. Bibcode : 2008JChPh.128t5101S . DOI : 10.1063 / 1.2919546 . PMID  18513044 .
• (Брацун и др., 2005): Брацун, Дмитрий; Вольфсон, Дмитрий; Поспешный, Джефф; Цимринг, Лев С. (2005). «Индуцированные задержкой стохастические колебания в регуляции генов» . Труды Национальной академии наук Соединенных Штатов Америки . 102 (41): 14593–8. Bibcode : 2005PNAS..10214593B . DOI : 10.1073 / pnas.0503858102 . PMC  1253555 . PMID  16199522 .
• (Баррио и др., 2006): Баррио, Мануэль; Беррейдж, Кевин; Лейер, Андре; Тиан, Тяньхай (2006). "Колебательное регулирование hes1 : дискретное стохастическое моделирование задержки и имитация" . PLOS вычислительная биология . 2 (9): 1017. Bibcode : 2006PLSCB ... 2..117B . DOI : 10.1371 / journal.pcbi.0020117 . PMC  1560403 . PMID  16965175 .
• (Цай 2007): Цай, Сяодун (2007). «Точное стохастическое моделирование сопряженных химических реакций с запаздыванием». Журнал химической физики . 126 (12): 124108. Bibcode : 2007JChPh.126l4108C . DOI : 10.1063 / 1.2710253 . PMID  17411109 .
• (Barnes Chu 2010): Barnes, David J .; Чу, Доминик (2010). Введение в моделирование для биологических наук . Springer Verlag.
• (Рамасвами Гонсалес-Сегредо Сбальзарини, 2009 г.): Рамасвами, Раджеш; Гонсалес-Сегредо, Нелидо; Сбальзарини, Иво Ф. (2009). «Новый класс высокоэффективных алгоритмов точного стохастического моделирования для сетей химических реакций». Журнал химической физики . 130 (24): 244104. arXiv : 0906.1992 . Bibcode : 2009JChPh.130x4104R . DOI : 10.1063 / 1.3154624 . PMID  19566139 . S2CID  4952205 .
• (Рамасвами Сбальзарини 2010): Рамасвами, Раджеш; Сбальзарини, Иво Ф. (2010). «Вариант с частичной склонностью алгоритма стохастического моделирования отклонения состава для сетей химических реакций» (PDF) . Журнал химической физики . 132 (4): 044102. Bibcode : 2010JChPh.132d4102R . DOI : 10.1063 / 1.3297948 . PMID  20113014 .
• (Индуркхья Бил 2010): Индуркхья, Сагар; Бил, Джейкоб С. (2005). Исалан, Марк (ред.). «Факторинг реакций и двухчастные графики обновления ускоряют алгоритм Гиллеспи для крупномасштабных биохимических систем» . PLOS ONE . 5 (1): e8125. Bibcode : 2010PLoSO ... 5.8125I . DOI : 10.1371 / journal.pone.0008125 . PMC  2798956 . PMID  20066048 .
• (Рамасвами Сбальзарини 2011): Рамасвами, Раджеш; Сбальзарини, Иво Ф. (2011). «Формулировка частичной склонности алгоритма стохастического моделирования для сетей химических реакций с задержками» (PDF) . Журнал химической физики . 134 (1): 014106. Bibcode : 2011JChPh.134a4106R . DOI : 10.1063 / 1.3521496 . PMID  21218996 .
• (Yates Klingbeil 2013): Yates, Christian A .; Клингбейл, Гвидо (2013). «Повторное использование случайных чисел в алгоритме стохастического моделирования» . Ежегодный обзор физической химии . 58 (9): 094103. Bibcode : 2013JChPh.138i4103Y . DOI : 10.1063 / 1.4792207 . PMID  23485273 .
• Гиллеспи, Дэниел Т. (2007). «Стохастическое моделирование химической кинетики». Журнал физической химии . 58 : 35–55. DOI : 10.1146 / annurev.physchem.58.032806.104637 . PMID  17037977 .

внешние ссылки

• Репозиторий, содержащий весь исходный код, использованный в этой статье.