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

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

• 1 История
• 2 Идея алгоритма
• 3 Алгоритм
• 4 Примеры
  • 4.1 Обратимое связывание A и B с образованием Димеры AB
  • 4.2 Эпидемия SIR без жизненной динамики
• 5 Дополнительная литература

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

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

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

Гиллеспи разработал две разные, но эквивалентные формулировки; прямой метод и метод первой реакции. Ниже приводится краткое описание шагов для запуска алгоритма (математические данные опущены):

1. Инициализация : инициализация количества молекул в системе, констант реакций и генераторов случайных чисел.
2. Шаг Монте-Карло : Генерировать случайные числа для определения следующей реакции, а также временного интервала. Вероятность выбора данной реакции пропорциональна количеству молекул субстрата, временной интервал экспоненциально распределен со средним значением $1\; /\; RTOT\; \{\backslash \; displaystyle\; 1\; /\; R\; \_\; \{\backslash \; mathrm\; \{TOT\}\}\}$.
3. Обновить : увеличьте время на случайно сгенерированное время на шаге 2. Обновите количество молекул на основе произошедшей реакции.
4. Итерация : вернитесь к шагу 2, если количество реагентов не равно нулю или моделирование время превышено.

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

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

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

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

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

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

$RTOT\; =\; К\; D\; N\; A\; N\; B\; +\; К\; B\; N\; AB\; \{\backslash \; Displaystyle\; R\; \_\; \{\backslash \; mathrm\; \{TOT\}\}\; =\; k\; \_\; \{\backslash \; mathrm\; \{D\}\}\; n\; \_\; \{\backslash \; mathrm\; \{A\}\}\; n\; \_\; \{\backslash \; mathrm\; \{B\}\; \}\; +\; k\; \_\; \{\backslash \; mathrm\; \{B\}\}\; n\; \_\; \{\backslash \; mathrm\; \{AB\}\}\}$

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

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

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

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

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

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

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

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

Эпидемия SIR без динамики жизнедеятельности

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

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

Единственная реализация эпидемии SIR, полученная с помощью реализации алгоритма Гиллеспи и численного решения системы обыкновенных дифференциальных уравнений (пунктирная линия).

$dn\; S\; dt\; =\; -\; \alpha \; n\; SV\; n\; I\; \{\backslash \; displaystyle\; \{\; \backslash \; frac\; \{dn\_\; \{S\}\}\; \{dt\}\}\; =\; -\; \{\backslash \; frac\; \{\backslash \; alpha\; n\_\; \{S\}\}\; \{V\}\}\; n\_\; \{I\}\}$,

$dn\; I\; dt\; =\; (\alpha \; n\; SV\; -\; \beta )\; n\; I\; \{\backslash \; displaystyle\; \{\backslash \; frac\; \{dn\_\; \{I\}\}\; \{dt\}\}\; =\; \backslash \; left\; (\{\backslash \; frac\; \{\backslash \; alpha\; n\_\; \{S\}\}\; \{V\}\}\; -\; \backslash \; beta\; \backslash \; right)\; n\_\; \{I\}\}$,

$dn\; R\; dt\; =\; \beta \; n\; I\; \{\backslash \; displaystyle\; \{\backslash \; frac\; \{dn\_\; \{R\}\}\; \{dt\}\}\; =\; \backslash \; beta\; n\_\; \{I\}\}$.

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

Смоделированная траектория, представленная на рисунке выше, получена с помощью следующей реализации Python алгоритма Гиллеспи. Кроме того, для получения численного решения системы дифференциальных уравнений вызывается решатель ODE из SciPy.

.

import math import random # Входные параметры ################### # int; общая численность населения N = 350 # float; максимальное прошедшее время T = 100.0 # float; время начала t = 0.0 # float; пространственный параметр V = 100.0 # float; скорость заражения после контакта _alpha = 10.0 # float; скорость лечения _beta = 10.0 # int; исходная инфицированная популяция n_I = 1 ###################################### # Вычислительная чувствительность Население, установить восстановленное равным нулю n_S = N - n_I n_R = 0 # Инициализировать список результатов SIR_data = SIR_data.append ((t, n_S, n_I, n_R)) # Основной цикл while t < T: if n_I == 0: break w1 = _alpha * n_S * n_I / V w2 = _beta * n_I W = w1 + w2 # generate exponentially distributed random variable dt # using inverse transform sampling dt = -math.log(1 - random.uniform(0.0, 1.0)) / W t = t + dt if random.uniform(0.0, 1.0) < w1 / W: n_S = n_S - 1 n_I = n_I + 1 else: n_I = n_I - 1 n_R = n_R + 1 SIR_data.append((t, n_S, n_I, n_R)) with open("SIR_data.txt", "w+") as fp: fp.write("\n".join("%f %i %i %i" % x for x in SIR_data)) from scipy.integrate import odeint import numpy as np # Numerical solution using an ordinary differential equation solver ts = np.linspace(0, 2, num=200) initial_S_I_R = (N - 1, 1, 0) def differential_SIR(initial_S_I_R, t, _alpha, _beta, V): n_S, n_I, n_R = initial_S_I_R dS_dt = -_alpha * n_S / V * n_I dI_dt = (_alpha * n_S / V - _beta) * n_I dR_dt = _beta * n_I return dS_dt, dI_dt, dR_dt solution = odeint(differential_SIR, initial_S_I_R, ts, args=(_alpha, _beta, V))

• Gillespie, Daniel Т. (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.
• Gibson, Michael A.; Брук, Иошуа (2000). «Эффективное точное стохастическое моделирование химических систем с множеством видов и множеством каналов» (PDF). Журнал физической химии A. 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). «Искусственная реализация простого процесса« рождения и смерти »». Журнал Королевского статистического общества, серия 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). «Адиабатическое крупнозернистое моделирование и моделирование стохастических биохимических сетей» (PDF). Труды Национальной академии наук Соединенных Штатов Америки. 106 (20): 10546–10551. Bibcode : 2009PNAS..10610546S. doi : 10.1073 / pnas.0809340106. PMC 2705573. PMID 19525397. Архивировано из оригинального (PDF) 14.07.2011.
• Салис, Ховард; Казнессис, Яннис Н. (2005). «Точное гибридное стохастическое моделирование системы связанных химических или биохимических реакций». Журнал химической физики. 122 (5): 054103. Bibcode : 2005JChPh.122e4103S. doi : 10.1063 / 1.1835951. PMID 15740306.
• (Slepoy Thompson Plimpton 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.
• (Barrio et al. 2006): Barrio, Manuel; Беррейдж, Кевин; Лейер, Андре; Тиан, Тяньхай (2006). «Колебательное регулирование hes1: моделирование и симуляция дискретной стохастической задержки». PLOS Вычислительная биология. 2 (9): 1017. Bibcode : 2006PLSCB... 2..117B. doi : 10.1371 / journal.pcbi.0020117. PMC 1560403. PMID 16965175.
• (Cai 2007): Cai, Xiaodong (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. Бибкод : 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). «Повторное использование случайных чисел в алгоритме стохастического моделирования». Журнал химической физики. 138 (9): 094103. Bibcode : 2013JChPh.138i4103Y. doi : 10.1063 / 1.4792207. PMID 23485273.