Занятие 7&8. Молекулярная динамика биологических молекул в GROMACS

Цель: ознакомится с возможностями моделирования молекулярной динамики.

В процессе выполнения задания будем пользоваться пакетом молекулярной динамики Gromacs.

Общие положения:

Типы файлов:

ћ gro - файл с координатами системы.

ћ top - файл с описанием ковалентных и нековалентных взаимодействий в молекулах.

ћ mdp - файл с описанием параметров для работы молекулярно-механического движка.

ћ tpr - файл для молекулярно-механического движка по сути есть объединение gro, top и mdp.

ћ trr, xtc - файл с координатами после расчета.

Моделирование перехода ДНК из А в В форму:

1. Создаем рабочую директорию на диске Н: типа Tereshkova/md

2. Нам даны файлы, скачаем их в свою рабочую директорию:

ћ файл параметров для минимизации энергии em.mdp.

ћ файл параметров для "утряски" воды pr.mdp pr.mdp.

ћ файл параметров для молекулярной динамики md.mdp.

3. Зайдем на kodomo через Putty и перейдем в рабочую директорию:

 cd Tereshkova/md

4. С помощью программы fiber из пакета 3DNA построим небольшой дуплекс с последовательностью GATCTA. Для того, что бы fiber работал надо задать путь и переменную среды:

export X3DNA=/home/preps/golovin/progs/X3DNA

export PATH=/home/preps/golovin/progs/X3DNA/bin:${PATH}

5. Скопируем директорию Tereshkova на 172.16.0.140. Зайдем через Putty на 172.16.0.140. Теперь построим файл топологии системы в силовом поле amber99sb и файл с координатами в формате Gromacs. Предполагается, что структура дуплекса находится в файле dna.pdb.

pdb2gmx -f dna.pdb -o dna -p dna -ff amber99sb -water tip3p

Получили сообщение об ошибке. Для ее устранения следует удалить 5' фосфаты из структуры дуплекса. Их два на 5' конце каждой цепи.

Также надо изменить имена нуклеотидов добвавив D к названию нуклеотида," Т"->"DT" и заменить имя атома "С5М" на " С7".

6. Сделаем небольшой отступ в ячейке от ДНК:

editconf -f dna.gro -o dna_ec -d 1.5

7. Проведем оптимизацию геометрии системы, что бы удалить "плохие" контакты в молекуле:

grompp -f em -c dna_ec -p dna -o dna_em -maxwarn 1

mdrun -deffnm dna_em -v

Отметим изменение максимальной силы в ходе оптимизации геометрии:

начальное и конечное значение максимальной силы соответственно.

8. Добавим в ячейку молекулы воды:

genbox -cp dna_em -p dna -cs -o dna_s

9. Нейтрализуем заряд системы. Это делаем в два шага: строим tpr и запускаем genion.

В выводе grompp обратим внимание на информацию о заряде системы!

grompp -f em -p dna -c dna_s -o dna_s

genion -s dna_s -o dna_si -p dna -np X

, где Х это количество положительных ионов необходимых для нейтрализации заряда системы.

В данном случае х=10.

10. Проведем "утряску" воды:

grompp -f pr -c dna_si -p dna -o dna_pr -maxwarn 1

mdrun -deffnm dna_pr -v

11. Переформатируйте dna_pr.gro и dna_si.gro в pdb формат.

editconf -f my.gro -o my.pdb

И сравним визуально в PyMol изменения в системах:

Заметное различие между полученными изображениями. Слева - система до утряски воды (dna_pr.pdb), справа - после (dna_si.pdb).

Видно, что система справа упорядочена, в отличие от ее исходного состояния.

12. Следующий шаг: скопировать файлы на суперкомпьютер.

Для этого:

ћ зайдем на суперкомпьютер и создадим папку с фамилиейи вернемся на kodomo

ssh skif

mkdir Tereshkova

exit

ћ скопируем файлы:

scp -r * skif:Tereshkova/

13. Запускаем тестовое моделирование на суперкомпьютере:

ssh skif

cd Tereshkova

grompp -f md -c pep_pr -p pep -o pep_md -maxwarn 1

mpirun  -np 16 -maxtime 5 -q test  /home/golovin/progs/bin/mdrun_mpi -deffnm pep_md -v

Ошибок нет, поэтому можем перейти к основному моделированию:

mpirun  -np 16 -maxtime 1200 /home/golovin/progs/bin/mdrun_mpi -deffnm pep_md -v

Номер задачи: 240958.

Анализ результатов моделирования перехода А-формы ДНК в

В-форму в воде.

Параметры моделирования:

ћ Силовое поле, используемое при построении топологии:

amber99sb

ћ Заряд системы: 0.

Нейтрализовали заряд системы следющим образом:

genbox -cp dna_em -p dna -cs -o dna_s

grompp -f em -p dna -c dna_s -o dna_s

genion -s dna_s -o dna_si -p dna -np 10

Добавив молекулы воды в систему, вычислили ее заряд, в соответствии с чем (заряд системы был равен "-10") было добавлено 10 катионов.

ћ Размер и форма ячейки:

Кубическая ячейка с параметрами

ћ Минимизация энергии :

- алгоритм минимизации энергии: integrator = l-bfgs

- алгоритм использования малого количества памяти Broyden-Fletcher-;Goldfarb-Shanno (квази-Ньютоновский),

метод аппроксимирует Гауссиановскую матрицу из предыдущих конфигураций.

- алгоритм расчета электростатики и Ван-дер-Ваальсовых взаимодействий: двойное обрезание (cut-off).

ћ Модель, которой описывался растворитель:

растворитель - вода. Модель воды: tip3p.

ћ Утряска растворителя:

Для биополимеров укажите параметр, который обуславливает неподвижность биополимера: define = -DPOSRES.

ћ Число шагов

nsteps = 1000000

ћ Длина шага

dt = 0.001 ; в pс

ћ Алгоритм расчета электростатики и Ван-дер-Ваальсовых взаимодействий:

Алгоритм для расчета электростатики - pme (суммирование по Эвальду).

Алгоритм для вычисления VdW - cut-off.

ћ Алгоритмы термостата и баростата.

Для контроля температуры использовался метод Berendsen. Контроля давления нет.

ћ Основной расчет МД:

- время моделирования: 4 часа 39 минут 7 секунд

- количество процессоров: 16

- эффективность машстабирования: 100%

- длину траектории : 10000 пс.

- число шагов: 5000000.

- длина шага: 0,002 пс.

- алгоритм интегратора: md

- алгоритм расчета электростатики и Ван-дер-Ваальсовых взаимодействий:

ћ алгоритм для расчета электростатики - pme (суммирование по Эвальду).

ћ Ван-дер-Ваальсовых взаимодействий - двойное обрезание (cut-off).

- алгоритмы термостата и баростата: для контроля температуры и давления использовался метод Berendsen.

АНАЛИЗ :

1. Любой анализ начинают с визуального анализа движений молекул. (при вопросе о выводе групп будем выбирать DNA).

 trjconv -f dna_md.xtc -s dna_md.tpr -o dna_pbc_1.pdb -skip 20 -pbc mol

Чтобы центрировать ячейку, и иметь возможно зрительно оценивать результат введем слеующую команду:

trjconv -f dna_md.xtc -s dna_md.tpr -o dna_fit_1.pdb -skip 20 -fit rot+trans

На мой взгляд, переход в В-форму происходит на 8 кадре:

dna_pbc_1.pdb , dna_fit_1.pdb .

Время соответствующего перехода записано в pdb-файле двумя строчками выше записи о 8 моделе.

Время перехода составляет 1600 пс.

2. Определим среднеквадратичное отклонение в ходе моделирования. Так как у нас происходит конформационный переход, сначала рассчитаем отклонение в ходе всей симуляции относительно стартовой структуры:

g_rms -f dna_md.xtc -s dna_md.tpr -o rms_1

Графическое представление: среднеквадратичное отклонение. rms_1

Из графика видно, что среднее квадратичное отклонение колеблется на уровне 0,25 нм.

Рассмотрим среднеквадратичное отклонение относительно каждой предыдущей структуры на растоянии 400 кадров. Если ближе к концу закончился конформационный переход, то отколнение должно уменьшаться:

g_rms -f dna_md.xtc -s dna_md.tpr -o rms_2 -prev 400

Графическое представление: rms_2

Действительно, к концу среднее квадратичное отклонение с 0,25 нм снизилось до 0,12-0,15 нм , что свидетельствует об окончании конформационного перехода.

3. Определим изменение гидрофобной и гидрофильной поверхности в ходе конформационного перехода:

g_sas -f dna_md.xtc -s dna_md.tpr -o sas_dna.xvg

Графическое представление: гидрофобная поверхность. sas_dna

Графическое представление: гидрофильная поверхность. sas_dna

Площадь гидрофобной и гидрофильной поверхности в процессе конформационного перехода остаются практически без изменений.

4. Традиционным анализом для ДНК является расчет колчества образуемых водородных связей.

Если мы будем исследовать связи между ДНК и ДНК, то это будут водородные связи между цепями ДНК.

Для конца траектории:

g_hbond  -f dna_md.xtc -s dna_md.tpr -num hbond_dna

Графическое представление: водородные связи ДНК-ДНК. hbond_dna

Количество водородных связей, образуемых с водой, колеблется в пределах от 100 до 120.

Не менее интересно будет изучить количество вдородных связей ДНК-Вода:

g_hbond  -f dna_md.xtc -s dna_md.tpr -num hbond_dna_sol

Графическое представление: водородные связи ДНК-вода. hbond_dna_sol

Известно, что водородных связей в дуплексе должно быть 14. Из графика видно, что, действительно, в процессе перехода количество водородных связей не сильно изменяется. Диапозон изменений:13-15 связей.

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