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

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

История

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

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

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

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

Алгоритм

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

$p(\tau ,j\mid {\boldsymbol {x}},t)=a_{j}({\boldsymbol {x}})\exp \left(-\tau \sum _{j}a_{j}({\boldsymbol {x}})\right),$

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

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

$\tau ={\frac {1}{\sum _{j}a_{j}({\boldsymbol {x}})}}\log \left({\frac {1}{r_{1}}}\right),$

и

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

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

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

где $\nu _{j}$ представляет собой добавление $j^{\text{th}}$ компонент данного вектора изменения состояния $\nu$. Это семейство алгоритмов требует больших вычислительных затрат, поэтому существует множество модификаций и адаптаций, включая метод следующей реакции (Гибсона и Брука), тау-прыжка , а также гибридные методы, в которых большое количество реагентов моделируются с детерминированным поведением. Адаптированные методы обычно ставят под угрозу точность теории, лежащей в основе алгоритма, поскольку она связана с основным уравнением, но предлагают разумные реализации для значительного сокращения временных рамок. Вычислительная стоимость точных версий алгоритма определяется классом связи реакционной сети. В слабосвязанных сетях количество реакций, на которые влияет любая другая реакция, ограничено небольшой константой. В сильно связанных сетях запуск одной реакции может в принципе повлиять на все остальные реакции. Была разработана точная версия алгоритма с масштабированием в постоянное время для слабосвязанных сетей, позволяющая эффективно моделировать системы с очень большим количеством каналов реакции (Слепой Томпсон Плимптон, 2008). Обобщенный алгоритм Гиллеспи, учитывающий немарковские свойства случайных биохимических событий с задержкой, был разработан Bratsun et al. 2005 г. и независимо Barrio et al. 2006 г., а также (Cai 2007). Подробности смотрите в статьях, цитируемых ниже.

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

Примеры

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

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

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

Суммарная скорость реакции, $R_{\mathrm {TOT} }$, в момент времени t тогда определяется выражением

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

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

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

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

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

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

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

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

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

Стохастическая самосборка

Модель Гарда описывает самосборку липидов в агрегаты. С помощью стохастического моделирования показано появление множества типов агрегатов и их эволюция.

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

• Гиллеспи, Дэниел Т. (1977). «Точное стохастическое моделирование связанных химических реакций». Журнал физической химии . 81 (25): 2340–2361. CiteSeerX   10.1.1.704.7634 . дои : 10.1021/j100540a008 . S2CID   2606191 .
• Гиллеспи, Дэниел Т. (1976). «Общий метод численного моделирования стохастической временной эволюции связанных химических реакций». Журнал вычислительной физики . 22 (4): 403–434. Бибкод : 1976JCoPh..22..403G . дои : 10.1016/0021-9991(76)90041-3 .
• Гибсон, Майкл А.; Брук, Иегошуа (2000). «Эффективное точное стохастическое моделирование химических систем со многими видами и множеством каналов» (PDF) . Журнал физической химии А. 104 (9): 1876–1889. Бибкод : 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 .
• Пресс, Уильям Х.; Теукольский, Саул А.; Веттерлинг, Уильям Т.; Фланнери, Брайан П. (2007). «Раздел 17.7. Стохастическое моделирование сетей химических реакций» . Численные рецепты: искусство научных вычислений (3-е изд.). Нью-Йорк, штат Нью-Йорк: Издательство Кембриджского университета. ISBN  978-0-521-88068-8 . Архивировано из оригинала 11 августа 2011 г. Проверено 17 августа 2011 г.
• Колмогоров, Андрей Николаевич (1931). «Об аналитических методах теории вероятностей». Математические летописи . 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. дои : 10.2307/2985327 . JSTOR   2985327 .
• Ратинам, Мурухан; Петцольд, Линда Р .; Цао, Ян; Гиллеспи, Дэниел Т. (2003). «Жесткость в стохастических химически реагирующих системах: неявный метод тау-прыжка». Журнал химической физики . 119 (24): 12784–12794. Бибкод : 2003JChPh.11912784R . дои : 10.1063/1.1627296 .
• Синицын, Николай А.; Хенгартнер, Николас; Неменман, Илья (2009). «Адиабатическое укрупнение и моделирование стохастических биохимических сетей» . Труды Национальной академии наук Соединенных Штатов Америки . 106 (20): 10546–10551. Бибкод : 2009PNAS..10610546S . дои : 10.1073/pnas.0809340106 . ПМК   2705573 . ПМИД   19525397 .
• Салис, Ховард; Казнессис, Яннис Н. (2005). «Точное гибридное стохастическое моделирование системы связанных химических или биохимических реакций». Журнал химической физики . 122 (5): 054103. Бибкод : 2005JChPh.122e4103S . дои : 10.1063/1.1835951 . ПМИД   15740306 .
• (Слепой Томпсон Плимптон 2008): Слепой, Александр; Томпсон, Эйдан П.; Плимптон, Стивен Дж. (2008). «Кинетический алгоритм Монте-Карло с постоянным временем для моделирования больших сетей биохимических реакций». Журнал химической физики . 128 (20): 205101. Бибкод : 2008JChPh.128t5101S . дои : 10.1063/1.2919546 . ПМИД   18513044 .
• (Брацун и др., 2005): Брацун, Дмитрий; Вольфсон, Дмитрий; Хэсти, Джефф; Цимринг, Лев С. (2005). «Вызванные задержкой стохастические колебания в регуляции генов» . Труды Национальной академии наук Соединенных Штатов Америки . 102 (41): 14593–8. Бибкод : 2005PNAS..10214593B . дои : 10.1073/pnas.0503858102 . ПМЦ   1253555 . ПМИД   16199522 .
• (Баррио и др., 2006 г.): Баррио, Мануэль; Беррейдж, Кевин; Лейер, Андре; Тянь, Тяньхай (2006). «Колебательное регулирование hes1 : моделирование и моделирование дискретной стохастической задержки» . PLOS Вычислительная биология . 2 (9): 1017. Бибкод : 2006PLSCB...2..117B . дои : 10.1371/journal.pcbi.0020117 . ПМК   1560403 . ПМИД   16965175 .
• (Кай 2007): Цай, Сяодун (2007). «Точное стохастическое моделирование связанных химических реакций с запаздыванием». Журнал химической физики . 126 (12): 124108. Бибкод : 2007JChPh.126l4108C . дои : 10.1063/1.2710253 . ПМИД   17411109 .
• (Барнс Чу 2010): Барнс, Дэвид Дж.; Чу, Доминик (2010). Введение в моделирование для биологических наук . Спрингер Верлаг. Бибкод : 2010itmf.book.....B .
• (Рамасвами Гонсалес-Сегредо Сбальцарини, 2009 г.): Рамасвами, Раджеш; Гонсалес-Сегредо, Нелидо; Сбалзарини, Иво Ф. (2009). «Новый класс высокоэффективных точных алгоритмов стохастического моделирования сетей химических реакций». Журнал химической физики . 130 (24): 244104. arXiv : 0906.1992 . Бибкод : 2009JChPh.130x4104R . дои : 10.1063/1.3154624 . ПМИД   19566139 . S2CID   4952205 .
• (Рамасвами Сбалзарини, 2010): Рамасвами, Раджеш; Сбальцарини, Иво Ф. (2010). «Вариант частичной склонности алгоритма стохастического моделирования с отклонением состава для сетей химических реакций» (PDF) . Журнал химической физики . 132 (4): 044102. Бибкод : 2010JChPh.132d4102R . дои : 10.1063/1.3297948 . ПМИД   20113014 .
• (Индурхья Бил 2010): Индурхья, Сагар; Бил, Джейкоб С. (2005). Исалан, Марк (ред.). «Факторинг реакций и двудольные графики обновления ускоряют алгоритм Гиллеспи для крупномасштабных биохимических систем» . ПЛОС ОДИН . 5 (1): е8125. Бибкод : 2010PLoSO...5.8125I . дои : 10.1371/journal.pone.0008125 . ПМЦ   2798956 . ПМИД   20066048 .
• (Рамасвами Сбалзарини, 2011): Рамасвами, Раджеш; Сбальцарини, Иво Ф. (2011). «Формулировка алгоритма стохастического моделирования с частичной склонностью для сетей химических реакций с задержками» (PDF) . Журнал химической физики . 134 (1): 014106. Бибкод : 2011JChPh.134a4106R . дои : 10.1063/1.3521496 . ПМИД   21218996 . S2CID   4949530 .
• (Йейтс Клингбейл 2013): Йейтс, Кристиан А.; Клингбейл, Гвидо (2013). «Переработка случайных чисел в алгоритме стохастического моделирования» . Ежегодный обзор физической химии . 58 (9): 094103. Бибкод : 2013JChPh.138i4103Y . дои : 10.1063/1.4792207 . ПМИД   23485273 .
• Гиллеспи, Дэниел Т. (2007). «Стохастическое моделирование химической кинетики». Ежегодный обзор физической химии . 58 : 35–55. Бибкод : 2007ARPC...58...35G . doi : 10.1146/annurev.physchem.58.032806.104637 . ПМИД   17037977 .