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

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

История

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

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

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

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

Алгоритм

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

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

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

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

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

и

${displaystyle j=}$ наименьшее целое число, удовлетворяющее ${displaystyle sum _{j'=1}^{j}a_{j'}({oldsymbol {x}})>r_{2}sum _{j}a_{j}({oldsymbol {x}})}$.

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

1. Установите время ${displaystyle t=t_{0}}$ и состояние системы ${displaystyle {oldsymbol {x}}={oldsymbol {x}}_{0}}$2. С системой в состоянии ${displaystyle {oldsymbol {x}}}$ вовремя ${displaystyle t}$оцените все ${displaystyle a_{j}({oldsymbol {x}})}$ и их сумма ${displaystyle sum _{j}a_{j}({oldsymbol {x}})}$3. Выполните следующую реакцию, заменив ${displaystyle tleftarrow t+ au }$ и ${displaystyle {oldsymbol {x}}leftarrow {oldsymbol {x}}+u _{j}}$4. Запись ${displaystyle ({oldsymbol {x}},t)}$ по желанию. Вернитесь к шагу 1 или завершите симуляцию.

Это семейство алгоритмов требует больших вычислительных ресурсов, поэтому существует множество модификаций и адаптаций, включая следующий метод реакции (Гибсон и Брук), тау-прыжки, а также гибридные методы, в которых многочисленные реагенты моделируются с детерминированным поведением. Адаптированные методы, как правило, ставят под угрозу точность теории, лежащей в основе алгоритма, поскольку он связан с Основным уравнением, но предлагают разумные реализации для значительно улучшенных временных масштабов. Вычислительная стоимость точных версий алгоритма определяется классом связи реакционной сети. В слабосвязанных сетях количество реакций, на которые влияет любая другая реакция, ограничено небольшой константой. В сильно связанных сетях однократное срабатывание реакции в принципе может повлиять на все остальные реакции. Была разработана точная версия алгоритма с постоянным масштабированием для слабосвязанных сетей, позволяющая эффективно моделировать системы с очень большим количеством каналов реакции (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:

от математика импорт журналкласс SSA:    "" "Контейнер для SSA" ""    def __в этом__(я, модель, семя=1234):        "" "Инициализировать контейнер с моделью и генератором псевдослучайных чисел" ""        я.модель = модель        я.случайный = Мерсенн(семя=семя)    def непосредственный(я):        "" "Неопределенный генератор траекторий прямого метода" ""        в то время как Правда:            в то время как не я.модель.выход():                # оценить веса и разбиение                веса = [                    (rxn, сто, профи(я.модель))                    для (rxn, сто, профи) в я.модель.реакции                ]                раздел = сумма(ш[-1] для ш в веса)                # оценить время пребывания (шаг 1 MC)                пребывание = журнал(1.0 / я.случайный.плавающий()) / раздел                я.модель["время"].добавить(я.модель["время"][-1] + пребывание)                # оценить реакцию (шаг 2 MC)                раздел = раздел * я.случайный.плавающий()                в то время как раздел >= 0.0:                    rxn, сто, профи = веса.поп(0)                    раздел -= профи                для виды, дельта в сто.Предметы():                    я.модель[виды].добавить(я.модель[виды][-1] + дельта)                я.модель.священник()            Уступать я.модель            я.модель.сброс настроек()    def first_reaction(я):        "" "Неопределенный генератор траекторий 1-й реакции" ""        в то время как Правда:            в то время как не я.модель.выход():                # оценить время следующей реакции                раз = [                    (                        журнал(                            1.0 / я.случайный.плавающий()                        ) / профи(я.модель),                        сто                    )                    для (rxn, сто, профи) в я.модель.реакции                ]                раз.Сортировать()                # оценить время реакции                я.модель["время"].добавить(                    я.модель["время"][-1] + раз[0][0]                )                # оценить реакцию                для виды, дельта в раз[0][1].Предметы():                    я.модель[виды].добавить(                        я.модель[виды][-1] + дельта                    )                я.модель.священник()            Уступать я.модель            я.модель.сброс настроек()

Как видно, как и все методы Монте-Карло, SSA требуется семя для воспроизводимости, которое передается имени Мерсенн, генератор псевдослучайных чисел. Реализация Мерсенн можно найти в Мерсенн Твистер статья, но можно использовать и другую, например random.random. В непосредственный и first_reaction члены являются неопределенными генераторами, то есть они будут непрерывно создавать траектории, каждая из которых является полной симуляцией системы, в цикле, пока сигнал не разорвет этот цикл. Чтобы фактически реализовать такой цикл и дать некоторое количество траекторий для анализа, этот класс требует, чтобы модель была передана ему при создании экземпляра. Цель класса модели - избежать смешения кинетических свойств конкретного моделируемого процесса с логикой алгоритма Гиллеспи. Модель должна быть подклассом словаря с общедоступными членами. священник, выход, и сброс настроек что соответственно

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

Рассмотрим следующую реализацию модели, которая используется в приведенном ниже примере SIR:

класс SSAModel(диктовать):    "" "Контейнер для модели SSA" ""    def __в этом__(        я, первоначальные условия, склонности, стехиометрия    ):        """        Инициализировать модель со словарем начальных условий (каждый        """        супер().__в этом__(**первоначальные условия)        я.реакции = список()        я.excluded_reactions = список()        для реакция,склонность в склонности.Предметы():            если склонность(я) == 0.0:                я.excluded_reactions.добавить(                    (                        реакция,                        стехиометрия[реакция],                        склонность                    )                )            еще:                я.реакции.добавить(                    (                        реакция,                        стехиометрия[реакция],                        склонность                    )                )    def выход(я):        "" "Верните True, чтобы сойти с траектории" ""        # вернуть True, если реакции больше нет        если len(я.реакции) == 0: вернуть Правда        # вернуть False, если есть другие реакции        еще: вернуть Ложь    def священник(я):        "" "Подтвердить и аннулировать реакции модели" ""                # оценить возможные реакции        реакции = []        в то время как len(я.реакции) > 0:            реакция = я.реакции.поп()            если реакция[2](я) == 0:                я.excluded_reactions.добавить(реакция)            еще:                реакции.добавить(реакция)        реакции.Сортировать()        я.реакции = реакции        # оценить невозможные реакции        excluded_reactions = []        в то время как len(я.excluded_reactions) > 0:            реакция = я.excluded_reactions.поп()            если реакция[2](я) > 0:                я.реакции.добавить(реакция)            еще:                excluded_reactions.добавить(реакция)        excluded_reactions.Сортировать()        я.excluded_reactions = excluded_reactions    def сброс настроек(я):        "" "Очистить траекторию" ""        # сбросить виды к начальным условиям        для ключ в я: дель я[ключ][1:]        # сбросить реакцию на начальные условия        я.священник()

Примеры

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

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

Если вовремя т есть одна молекула каждого типа, тогда скорость образования димера равна ${displaystyle k_{mathrm {D} }}$, а если есть ${displaystyle n_{mathrm {A} }}$ молекулы типа А и ${displaystyle n_{mathrm {B} }}$ молекулы типа B, скорость образования димера равна ${displaystyle k_{mathrm {D} }n_{mathrm {A} }n_{mathrm {B} }}$. Если есть ${displaystyle n_{mathrm {AB} }}$ димеров, то скорость диссоциации димеров равна ${displaystyle k_{mathrm {B} }n_{mathrm {AB} }}$.

Общая скорость реакции, ${displaystyle R_{mathrm {TOT} }}$, вовремя т тогда дается

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

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

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

Сюжет номера А молекулы (черная кривая) и AB димеры как функция времени. Поскольку мы начали с 10 А и B молекулы во время т= 0, количество B молекул всегда равно количеству А молекул и так не показано.

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

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

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

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

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

Эпидемия 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}}-eta ight)n_{I}}$,

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

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

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

от matplotlib импорт пиплотот тупой импорт внутреннее пространствоот scipy.integrate импорт odeint# начальный подсчет видов и время пребыванияinitital_conditions = {    "s": [480],    "я": [20],    "р": [0],    "время": [0.0],}# функции склонностисклонности = {    0: лямбда d: 2.0 * d["s"][-1] * d["я"][-1] / 500,    1: лямбда d: 1.0 * d["я"][-1],}# изменение видов для каждой склонностистехиометрия = {    0: {"s": -1, "я": 1, "р": 0},    1: {"s": 0, "я": -1, "р": 1},}# создать экземпляр контейнера модели SSA эпидемииэпидемия = SSAModel(    initital_conditions,    склонности,    стехиометрия)# создаем экземпляр контейнера SSA с модельюэпидемический_ генератор = SSA(эпидемия)# сделай красивую большую фигурупиплот.фигура(фиговый=(10,10), dpi=500)# создайте подзаговор для восприимчивых, инфицированных и выздоровевших людейaxes_s = пиплот.подсюжет(311)axes_s.set_ylabel("восприимчивые люди")axes_i = пиплот.подсюжет(312)axes_i.set_ylabel("инфицированные люди")axes_r = пиплот.подсюжет(313)axes_r.set_ylabel("выздоровевшие лица")axes_r.set_xlabel(«время (условные единицы)»)# моделируем и строим 30 траекторийтраектории = 0для траектория в эпидемический_ генератор.непосредственный():    axes_s.сюжет(траектория["время"], траектория["s"], цвет="оранжевый")    axes_i.сюжет(траектория["время"], траектория["я"], цвет="оранжевый")    axes_r.сюжет(траектория["время"], траектория["р"], цвет="оранжевый")    траектории += 1    если траектории == 30:        перерыв# численное решение с использованием решения обыкновенного дифференциального уравненият = внутреннее пространство(0, 14, число=200)y0 = (480, 20, 0)альфа = 2.0бета = 1.0def дифференциальный_SIR(n_SIR, т, альфа, бета):    dS_dt = -альфа * n_SIR[0] * n_SIR[1] / 500    dI_dt = ((альфа * n_SIR[0] / 500) - бета) * n_SIR[1]    dR_dt = бета * n_SIR[1]    вернуть dS_dt, dI_dt, dR_dtрешение = odeint(дифференциальный_SIR, y0, т, аргументы=(альфа, бета))решение = [[ряд[я] для ряд в решение] для я в ассортимент(3)]# построить численное решениеaxes_s.сюжет(т, решение[0], цвет="черный")axes_i.сюжет(т, решение[1], цвет="черный")axes_r.сюжет(т, решение[2], цвет="черный")пиплот.шоу()

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

• Гиллеспи, Дэниел Т. (1977). «Точное стохастическое моделирование связанных химических реакций». Журнал физической химии. 81 (25): 2340–2361. CiteSeerX  10.1.1.704.7634. Дои:10.1021 / j100540a008.
• Гиллеспи, Дэниел Т. (1976). «Общий метод численного моделирования стохастической временной эволюции связанных химических реакций». Журнал вычислительной физики. 22 (4): 403–434. Bibcode:1976JCoPh..22..403G. Дои:10.1016/0021-9991(76)90041-3.
• Гибсон, Майкл А .; Брук, Иошуа (2000). «Эффективное точное стохастическое моделирование химических систем с множеством видов и множеством каналов» (PDF). Журнал физической химии А. 104 (9): 1876–1889. Bibcode:2000JPCA..104.1876G. Дои:10.1021 / jp993732q.
• Дуб, Джейкоб Л. (1942). «Вопросы теории цепей Маркова». Труды Американского математического общества. 52 (1): 37–64. Дои:10.1090 / S0002-9947-1942-0006633-7. JSTOR  1990152.
• Дуб, Джейкоб Л. (1945). «Цепи Маркова - счетный случай». Труды Американского математического общества. 58 (3): 455–473. Дои: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. Дои:10.1007 / BF01457949. S2CID  119439925.
• Феллер, Вилли (1940). "Об интегродифференциальных уравнениях чисто разрывных марковских процессов". Труды Американского математического общества. 48 (3): 4885–15. Дои:10.2307/1990095. JSTOR  1970064.
• Кендалл, Дэвид Г. (1950). «Искусственная реализация простого процесса« рождения и смерти »». Журнал Королевского статистического общества, серия B. 12 (1): 116–119. JSTOR  2983837.
• Бартлетт, Морис С. (1953). «Случайные процессы или статистика изменений». Журнал Королевского статистического общества, серия C. 2 (1): 44–64. JSTOR  2985327.
• Ратинам, Мурухан; Петцольд, Линда Р.; Цао, Ян; Гиллеспи, Дэниел Т. (2003). «Жесткость в стохастических химически реагирующих системах: неявный метод тау-прыжка». Журнал химической физики. 119 (24): 12784–12794. Bibcode:2003ЖЧФ.11912784Р. Дои:10.1063/1.1627296.
• Синицын, Николай А .; Хенгартнер, Николас; Неменман, Илья (2009). «Адиабатическое крупнозернистое моделирование и моделирование стохастических биохимических сетей» (PDF). Труды Национальной академии наук Соединенных Штатов Америки. 106 (20): 10546–10551. Bibcode:2009PNAS..10610546S. Дои:10.1073 / pnas.0809340106. ЧВК  2705573. PMID  19525397. Архивировано из оригинал (PDF) на 2011-07-14.
• Салис, Ховард; Казнессис, Яннис Н. (2005). «Точное гибридное стохастическое моделирование системы связанных химических или биохимических реакций». Журнал химической физики. 122 (5): 054103. Bibcode:2005ЖЧФ.122э4103С. Дои:10.1063/1.1835951. PMID  15740306.
• (Слепой Томпсон Плимптон 2008): Слепой, Александр; Thompson, Aidan P .; Плимптон, Стивен Дж. (2008). «Кинетический алгоритм Монте-Карло с постоянным временем для моделирования больших сетей биохимических реакций». Журнал химической физики. 128 (20): 205101. Bibcode:2008ЖЧФ.128т5101С. Дои:10.1063/1.2919546. PMID  18513044.
• (Брацун и др., 2005): Брацун, Дмитрий; Вольфсон, Дмитрий; Поспешный, Джефф; Цимринг, Лев С. (2005). «Стохастические колебания, вызванные задержкой в ​​регуляции генов». Труды Национальной академии наук Соединенных Штатов Америки. 102 (41): 14593–8. Bibcode:2005ПНАС..10214593Б. Дои:10.1073 / pnas.0503858102. ЧВК  1253555. PMID  16199522.
• (Баррио и др., 2006): Баррио, Мануэль; Беррейдж, Кевин; Лейер, Андре; Тиан, Тяньхай (2006). "Колебательное регулирование hes1: Моделирование и симуляция дискретной стохастической задержки ». PLOS вычислительная биология. 2 (9): 1017. Bibcode:2006PLSCB ... 2..117B. Дои:10.1371 / journal.pcbi.0020117. ЧВК  1560403. PMID  16965175.
• (Цай 2007): Цай, Сяодун (2007). «Точное стохастическое моделирование сопряженных химических реакций с запаздыванием». Журнал химической физики. 126 (12): 124108. Bibcode:2007ЖЧФ.126л4108С. Дои:10.1063/1.2710253. PMID  17411109.
• (Барнс Чу 2010): Барнс, Дэвид Дж .; Чу, Доминик (2010). Введение в моделирование для биологических наук. Springer Verlag.
• (Рамасвами Гонсалес-Сегредо Сбальзарини, 2009 г.): Рамасвами, Раджеш; Гонсалес-Сегредо, Нелидо; Сбальзарини, Иво Ф. (2009). «Новый класс высокоэффективных алгоритмов точного стохастического моделирования для сетей химических реакций». Журнал химической физики. 130 (24): 244104. arXiv:0906.1992. Bibcode:2009ЖЧФ.130х4104Р. Дои:10.1063/1.3154624. PMID  19566139. S2CID  4952205.
• (Рамасвами Сбальзарини, 2010): Рамасвами, Раджеш; Сбальзарини, Иво Ф. (2010). "Вариант с частичной склонностью алгоритма стохастического моделирования отклонения состава для сетей химических реакций" (PDF). Журнал химической физики. 132 (4): 044102. Bibcode:2010ЖЧФ.132д4102Р. Дои:10.1063/1.3297948. PMID  20113014.
• (Индуркхья Бил 2010): Индуркхья, Сагар; Бил, Джейкоб С. (2005). Исалан, Марк (ред.). «Факторинг реакций и двухчастные графики обновления ускоряют алгоритм Гиллеспи для крупномасштабных биохимических систем». PLOS ONE. 5 (1): e8125. Bibcode:2010PLoSO ... 5.8125I. Дои:10.1371 / journal.pone.0008125. ЧВК  2798956. PMID  20066048.
• (Рамасвами Сбальзарини 2011): Рамасвами, Раджеш; Сбальзарини, Иво Ф. (2011). «Формулировка алгоритма стохастического моделирования для сетей химических реакций с задержками» (PDF). Журнал химической физики. 134 (1): 014106. Bibcode:2011ЖЧФ.134а4106Р. Дои:10.1063/1.3521496. PMID  21218996.
• (Йетс Клингбейл 2013): Yates, Christian A .; Клингбейл, Гвидо (2013). «Повторное использование случайных чисел в алгоритме стохастического моделирования». Ежегодный обзор физической химии. 58 (9): 094103. Bibcode:2013ЖЧФ.138и4103Г. Дои:10.1063/1.4792207. PMID  23485273.
• Гиллеспи, Дэниел Т. (2007). «Стохастическое моделирование химической кинетики». Журнал физической химии. 58: 35–55. Дои:10.1146 / annurev.physchem.58.032806.104637.