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

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

История

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

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

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

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

Алгоритм

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

,

где каждый из члены являются функциями склонности элементарной реакции, аргумент которой , считается вектор видов. В параметр - время до следующей реакции (или время пребывания), а время конечно текущее. Перефразируя Гиллеспи, это выражение читается как «вероятность, заданная , что следующая реакция системы произойдет в бесконечно малом интервале времени , и будет иметь стехиометрию, соответствующую -я реакция ". Эта формулировка дает возможность использовать методы прямой и первой реакции, подразумевая - случайная величина с экспоненциальным распределением, а является "статистически независимой целочисленной случайной величиной с точечными вероятностями ".

Таким образом, метод генерации Монте-Карло состоит в том, чтобы просто нарисовать два псевдослучайных числа, и на , и вычислить

,

и

наименьшее целое число, удовлетворяющее .

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

1. Установите время  и состояние системы 2. С системой в состоянии  вовремя оцените все  и их сумма 3. Выполните следующую реакцию, заменив  и 4. Запись  по желанию. Вернитесь к шагу 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 молекула , а скорость реакции AB распад димера .

Если вовремя т есть одна молекула каждого типа, тогда скорость образования димера равна , а если есть молекулы типа А и молекулы типа B, скорость образования димера равна . Если есть димеров, то скорость диссоциации димеров равна .

Общая скорость реакции, , вовремя т тогда дается

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

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

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

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

вероятность того, что реакция

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

Теперь у нас обоих есть время, чтобы т + δт, и выполнил единственную реакцию. Алгоритм Гиллеспи просто повторяет эти два шага столько раз, сколько необходимо, чтобы моделировать систему, сколько мы хотим (то есть столько реакций). Результат моделирования Гиллеспи, который начинается с и в т= 0, а где и , показан справа. Для этих значений параметров в среднем 8 димеры и 2 из А и B но из-за небольшого числа молекул колебания вокруг этих значений велики. Алгоритм Гиллеспи часто используется для изучения систем, в которых эти колебания важны.

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

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

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

SIR graph.png

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

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

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

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


от 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], цвет="черный")пиплот.шоу()

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