Моделирование

Моделирование производилось в поле Amber99sb с водой TIP3P. Ячейка триклиническая (по умолчанию), отсуп от системы 1.5 нм.

В ходе подготовки относительно текста задания пришлось сменить команды ...

Заряд системы составил -10: всего в системе шесть пар нуклеотидов, соответственно, по пять фосфатных групп с зарядом -1 на каждой цепи, сумме 10. Конечные фосфатные группы были удалены. Заряд был нейтрализован добавлением 10 катионов натрия.

Минимизация энергии производилась алгоритмом L-BFGS, для расчёта электростатических и ван-дер-ваальсовых взаимодействий использовалась простые отсечки (cut-off).

Утряска воды производилась обычным leap-frog интегратором (md). Расчёт ван-дер-ваальсовых взаимодействий производился с простыми отсечками, в то время как для электростатических использовался алгоритм Эвальда (PME). Продолжительность составила 1 000 шагов по 0.0002 пс, т.е. 0.2 пс. Термостат velocity rescale, фиксированный размер ячейки (баростат отсутсвует). ДНК была зафиксирована на месте при помощи параметра define = -DPOSRES, включающего в топологию файл posre.itp с ограничениями на позиции.

На изображениях ниже можно видеть, что до утряски молекулы растворителя были расположены периодически, а после они расположены хаотически.

Итоговое моделирование во многом аналогично утряске: leap-frog интегратор, электростатика по Эвальду, ван-дер-ваальсовы просто отсечками, термостат velocity rescale. Однако использовался баростат Берендсена, а позиционные ограничения на ДНК были сняты. Продолжительность симуляции 5 000 000 шагов по 0.002 пс, всего 10 нс

Реальное (wall) время симуляции в на двух ядрах составило 32.6 минуты, суммарное время процессора (core time) — 62.3 минуты. Разница практически ровно в два раза, наблюдается эффективное линейное масштабирование.

Анализ

Вместо устаревших команд g_rms, g_sas, g_hbond использовались обновлённые rms, sasa, hbond. Однако синтаксис выделения гидрофобной и гидрофильной поверхностей из справки к команде g_sas (ссылка) сработал и для sasa.

Визуализация

Для предотвращения "разрывания" двух цепей ДНК друг от друга траектория была преобразована с обработкой периодических граничных условий как nojump.

Переход А-В формы визуально не очень-то виден. Примерно на 2200 пс происходит выворачивание первного тимидина. Он остаётся в вывернутом состоянии где-то до 9800 пс.

RMSD

На графике среднеквадратичного отклонения от исходной структуры хорошо видно существование двух состояний: с RMSD 0–0.5 нм и c RMSD 2–2.5 нм, которые хотелось бы сопоставить с A и B формами ДНК, однако, вероятно, такое соответствие неверно (см. ниже), поэтому далее они будут называться состояниями 0 и 1. Можно видеть, что переходы между состояниями происходят почти мгновенно, однако сперва состояние-1 существует лишь краткое время и только после примерно 8 000 пс оно становится основной, хотя скачки в состояние-0 всё ещё наблюдаются весьма часто, т.е. можно сказать, что переход не успел полностью завершиться за 10 000 пс симуляции.

Для опеделения однородности состояния-1 (вдруг это смесь различных состояний с одинаковым RMSD относительно исходной структуры) был также посчитан RMSD относительно структуры определённое число кадров назад. Поскольку переход произошёл позже, чем за рекомендуемые в задании 400 кадров (4000 пс) до конца симуляции, при таком временном сдвиге сравнения со структурой того же состояния не будет, поэтому это бессмысленно. Для определения адекватного меньшего знаения я построила график автокорреляции RMSD относительно исходной структуры в зависимости от временного смещения.

На графике наблюдается два плато — первое на 500-600 пс, второе от 2000 пс — отражающие некие два уровня схожести близких по времени структур. В качестве временного сдвига я взяла конец первого плато, 600 пс (60 кадров). Можно видеть, что, за исключением периодических кратковременных выбросов, на 8000–8600 происходит переход, после которого наблюдется примерно константное RMSD около 0.5 нм, что пусть и немного выше, чем RMSD до перехода, вероятно всё же соответсвует флуктуациям одной структуры — т.е. состояние-1 действительно относительно однородно.

Состояния в симуляции vs. формы ДНК

Для определения, насколько ли два состояния соответствуют A и В формам ДНК я сравнила распределение параметров динуклеотидного шага, описывающих геометрию ДНК. Референсные значения для А и В форм взяты из работы (Lu&Olson, 2003)

Можно видеть, что какое-то нормальное соответствие есть только для Twist. Параметр Rise должен быть одинаковым для обеих форм, однако в состоянии-1 имеет гораздо более низкие значения. Roll и Slide практически одинаковы в конце и в начале. Единственный вывод, который можно сделать — значимый отход от исходной структуры и увеличение дисперсии параметров начинается не на 8000 пс, а гораздо раньше, на 2000 пс, что совпадает с временем выворачивания тимидина, при том, что отдельные события перехода в состояние-1 наблюдаются и до этого.

SASA

Можно видеть, что доступная поверхность колеблется не сильно и её изменения слабо коррелируют с состояниями по RMSD. Некоторое повышение гидрофильной поверхности с 2000 пс можно объяснить выворачиванием первого тимидина, однако непонятно, почему гидрофобная при этом падает.

Водородные связи

Колебания числа водородных связей не коррелирует с состоянием. Самое частое значение — 14 — соответсвует общему числу водородных связей при уотсон-криковском взаимодействии всех азотистых оснований (последовательность TAGATC, 2 пары G-C по 3 связи, 4 пары A-T по 2 связи), т.к. водородных связей с остовом ДНК в обычной конформации не происходит. Меньше связей может наблюдаться при нарушении спаривания, больше — в ситуациях сильной девормации структуры, когда конечные остатки изгибаются и начинают образовывать водородные связи с остовом и другой стороной азотистых оснований (также, возможно, в некоторых таких конфигурациях может идти завышение реального числа связей, если атом-донор с одним водородом, расположенный между двумя атомами-акцепторами, даёт как бы по водородной связи с каждым).

Водородные связи с водой колеблются вокруг приблизительно 110 сходным образом в течение всей симуляции. В целом потенциальных валентностей на водородные связи у данного фрагмента ДНК 136 (фосфаты — 5*2 штук — по 2 как акцептор от каждого из четырёх кислородов, конечные -OH — 4 шт — по 2 как акцептор и 1 как донор, -O- сахара — 6*2 шт — 2 как акцептор, плюс азотистые основания: три дополнительные у аденина (как акцептор две =N- и одна дполнительная как донор -NH2, вторая уотсон-криковская), четыре у гуанина (то же + 1 доп как акцептор от =O), две от тимона (от =O), две от цитозина (от =O и -NH2)). Однако, вероятно, расположение вод для одновременного образования всех этих связей стерически затруднено.

Итог

Пронаблюдать переход А формы ДНК в В форму не удалось, однако удалось пронаблюдать выворачивание первого тимидина и связанные с эти изменения.