Таким образом, время наступления первой реакции рассчитывается методом обратной функции по формуле (15):
(15) (16) |
Номер реакции
определяется методом суперпозиции по формуле (16):(16) |
где
- независимые реализации базовой случайной величиныРисунок 2- Алгоритм прямого метода Гиллеспи
Схема алгоритма моделирования представлен на рисунке 2. Охарактеризуем каждый из блоков:
1. Задание начальных концентраций,скоростей и механизмов реакций, моделируемого объема, времени моделирования, и установление
.2. Вычисление
для всех и .3. Вычисление времени наступления реакции, используя
по формуле (15).4. Генерирование случайного
, определяющего номер произошедшей реакции. Определение номера произошедшей реакции по формуле (16).5. Замена количества молекул согласно произошедшей реакции,
устанавливаем
6. Переход к шагу 2.
Как уже писалось, этот алгоритм использует два случайных числа за итерацию, что занимает времени, пропорционально количеству реакций, для пересчета
и μ [12].Вместе с первым методом Гиллеспи был разработан алгоритм «первой реакции»,который генерирует предполагаемое время наступления каждой из реакции,- время, когда наступит данная реакция, если никакая другая не наступит раньше, потом полагается, что μ -реакция с наименьшим предполагаемым временем, τ переобозначается как
.Алгоритм :
1.Задание начальных концентраций,реакций и их скоростей, установление
;2.Вычисление
для всех i.3. Генерирование для каждого i с помощью метода обратных функций возможного времени
,4.Пусть μ-реакция с наименьшим временем
¸5.Полагаем τ=
¸.6.Меняем количество молекул согласно произошедшей реакции μ,устанавливаем
;7.Переход к шагу 2.
На первый взгляд, эти два метода кажутся очень разными, но они эквивалентны, так как распределения вероятности, используемые для вычисления τ и μ одинаковы. Как уже писалось, данный алгоритм также использует r случайных чисел за итерацию (r-количество реакций ), что затратит время, пропорционально r, чтоб пересчитать все
и столько же, чтоб вычислить наименьшее [15].Для построения имитационной модели полимеризации актина были сделаны следующие допущения:
· Процесс полимеризации моделируется в кубе,объемом V с периодическими граничными условиями.
· Количество молекул мономеров, например
,в начальный момент времени в объеме V высчитывается из начальных концентрации молекул.· Процесс диффузии происходит намного быстрее ,чем процессы полимеризации, поэтому при моделировании все вещества считаются равнораспределёнными по объему.
Моделировалась динамика роста отдельных филаментов и концентрация всех реагентов в системе. Предполагается, что в начальный момент времени в системе нет филаментов, поэтому начальное состояние задается концентрациями мономеров актина и протеинов в свободной форме.
На рисунке 3 приведена диаграмма, отражающая влияние на процесс полимеризации всех компонентов, реализованных в модели. На схеме представлены как реагенты, так и переходы между ними в результате реакций. Направление ребра в сторону или от эллипса, отображающего конкретную реакцию, соответствует связи с реагентом или продуктом реакции соответственно. Прерывистые линии используются для катализаторов, которые возвращаются в исходное состояние после окончания реакции. Пунктирные линии соответствуют реагентам, которые являются необязательными для протекания соответствующей реакции, однако их участие изменяет скорость реакции.
Рисунок 3-Основные реагенты и реакции в процессе актин-полимеризации
Для генерации событий были реализован вероятностный подход с использованием методики Монте-Карло, представленный в виде трех методов моделирования химических реакций: оригинальный метод Гиллеспи, метод первой реакции и оптимизированный прямой метод Гиллеспи. Метод τ-leap не использовался в данной работе, поскольку по результатам проведенных предварительных тестов доказал свою плохую сходимость при реализации в системах с малым числом молекул.
Реализация же метода Гибсон-Брук в данном случае не возможна, поскольку в соответствующем графе невозможно отразить распад коротких филаментов. Кроме того, из-за сложности системы построение графа становиться неоправданно ресурсоемким. Для филаментов и актинов использована объектная модель, остальные реагенты представлены только числом частиц. Следует отметить, что плюс и минус концы филаментов в различных формах (АДФ и АТФ), представлены как отдельные реагенты, также как и плюс концы филаментов с блокирующим протеином или формином.
4.1 Сравнение аналитической и имитационной моделей
Для проверки правильности построения имитационного алгоритма проведено сравнение результатов, полученных с использованием имитационной модели с аналитическим решением для упрощенной системы с меньшим числом реагентов и реакций, приведенным в статье [8]. Анализировалось две системы дифференциальных уравнений.
Аналитическая модель, которая сравнивалась с полученной имитационной дополнительно включает реакции блокировки ,разблокировки и нуклеации с участием блокирующего протеина[8].
Для тестирования использованы начальные концентрации и скорости реакций, представленные в таблице 1[8,16].Таблица 12. Начальные параметры, использованные для сравнения моделей (для неуказанных величин начальные значения равны нулю)
Обозначение | Описание | Значение |
скорость спонтанной нуклеации | ||
скорость элонгации плюс-конца | ||
скорость элонгации минус-конца | ||
скорость диссоциации плюс -конца | ||
скорость диссоциации минус -конца | ||
скорость блокировки плюс-конца | ||
скорость разблокировки плюс-конца | ||
скорость нуклеации с участием блокирующего протеина | ||
концентрация g-актина | ||
концентрация блокирующего протеина | ||
Результаты представлены на рисунке 4.