Для начала рассмотрим общие положения.
Типы файлов:
- gro - файл с координатами системы.
- top - файл с описанием ковалентных и нековалентных взаимодействий в молекулах.
- mdp - файл с описанием параметров для работы молекулярно-механического движка.
- tpr - файл для молекулярно-механического движка, по сути это объединение gro, top и mdp.
- trr, xtc - файл с координатами после расчёта.
Вся работа проводилась на компьютере 172.16.0.140.
Попробуем смоделировать переход А-формы ДНК в В-форму в воде. Это и будет нашим заданием на семинар.
Задание 1
Нам даны файлы:- файл праметров для минимизации энергии em.mdp.
- файл праметров для "утряски" воды pr.mdp pr.mdp.
- файл праметров для молекулярной динамики
md.mdp.
Задание 2
С помощью программы fiber из пакета 3DNA построим небольшой дуплекс с последовательностью GATCTA. Чтобы fiber работал, надо задать путь и переменную среды:export X3DNA=/home/preps/golovin/progs/X3DNA
export PATH=/home/preps/golovin/progs/X3DNA/bin:${PATH}
Задание 3
Скопируем нашу рабочую директорию на 172.16.0.140. Зайдем с помощью Putty на 172.16.0.140.Теперь построим файл топологии системы в силовом поле amber99sb и файл с координатами в формате Gromacs.
pdb2gmx -f dna.pdb -o dna -p dna -ff amber99sb -water tip3p
Получив сообщение об ошибке, мы удалили 5'-фосфаты из структуры дуплекса. Их было два, на 5'- конце каждой из цепей. Также мы изменили имена нуклеотидов, добавив D к названию нуклеотида, т.е. " Т"->"DT". Еще надо не забыть заменить имя атома "С5М" на " С7".
Полученный после всего этого файл: dna.pdb
Задание 4
Сделаем небольшой отступ в ячейке от ДНК.editconf -f dna.gro -o dna_ec -d 1.5
Задание 5
Проведём оптимизацию геометрии системы, чтобы удалить "плохие" контакты в молекуле.grompp -f em -c dna_ec -p dna -o dna_em -maxwarn 1
mdrun -deffnm dna_em -v
Изменение максимальной силы в ходе оптимизации геометрии:
начальное значение максимальной силы = Step 0 Fmax = 4.789e+03;
конечное значение максимальной силы = Step 54 Fmax = 9.566e+02.
Задание 6
Добавим в ячейку молекулы воды.genbox -cp dna_em -p dna -cs -o dna_s
Задание 7
Нейтрализуем заряд системы. Это делаем в два шага: строим 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
X = 10 (это количество положительных ионов, необходимых для нейтрализации заряда системы).
Задание 8
Проведём "утряску" воды:grompp -f pr -c dna_si -p dna -o dna_pr -maxwarn 1
mdrun -deffnm dna_pr -v
Задание 9
Переформатируем dna_pr.gro и dna_si.gro в pdb формат:editconf -f my.gro -o my.pdb
Задание 10
Теперь надо скопировать наши файлы на суперкомпьтер. Сначала зайдем на суперкомпьтер и создадим папку с нашей (моей) фамилией и вернёмся на kodomo.ssh skif
mkdir Ivanov
exit
Задание 11
Копируем все файлы:cd ..
scp -r * skif:Ivanov/
Задание 12
Запускаем тестовое моделирование на суперкомпьтере:ssh skif
cd Ivanov
grompp -f md -c dna_pr -p dna -o dna_md -maxwarn 1
mpirun -np 16 -q test -maxtime 5 /home/golovin/progs/bin/mdrun_mpi -deffnm dna_md -v
Файл ошибок не содержит, значит, можем двигаться дальше.
Задание 13
Запускаем основное моделирование на суперкомпьтере:mpirun -np 16 -maxtime 1200 /home/golovin/progs/bin/mdrun_mpi -deffnm dna_md -v
Анализ результатов моделирования перехода А-формы ДНК в
В-форму в воде.
Параметры моделирования:
Силовое поле, используемое при построении топологии:amber99sb
Заряд системы. Причины этого значения:
0; поскольку мы нейтрализовали заряд -10 после добавления 10 противоионов, который был на отрицательнозаряженной ДНК. Возможно, это примерно объясняет количество пар в нашей спирали.
Размер и форма ячейки:
Кубическая ячейка с параметрами 2.11772 x 1.94562 x 2.49396
Минимизация энергии
Алгоритм минимизации энергии
integrator = l-bfgs
Алгоритм использования малого количества памяти Broyden-Fletcher-;Goldfarb-Shanno (квази-Ньютоновский), метод аппроксимирует Гауссиановскую матрицу из предыдущих конфигураций.
Алгоритм расчёта электростатики и Ван-дер-Ваальсовых взаимодействий:Алгоритм для расчёта электростатики - двойное обрезание (cut-off).
Ван-дер-Ваальсовых взаимодействий - двойное обрезание (cut-off).
Модель, которой описывался растворитель
implicit_solvent = No ;т.е. неявного растворителя нет. Явный растворитель - вода.
Утряска растворителя:
Для биополимеров укажите параметр, который обуславливает неподвижность биополимера.
В строке: define = -DPOSRES. Строка фиксирует все связи отрицательно заряженных молекул.
Число шагов
nsteps = 1000000
Длина шага
dt = 0.001 ; в pс
Алгоритм расчёта электростатики и Ван-дер-Ваальсовых взаимодействий
Алгоритм для расчёта электростатики - pme (суммирование по Эвальду).
Алгоритм для вычисления VdW - cut-off.
Алгоритмы термостата и баростата.
Для контроля температуры использовался метод Berendsen. Контроля давления нет.
Основной расчёт МД:
Время моделирования, количество процессоров, эффективность машстабирования.
4 hours 38 minutes 1 seconds; 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
Открыли полученный файл b_pbc_1.pdb в PyMol , включив анимацию. Результат визуализации нас явно не устраивает, поскольку в окошке PyMol полная неразбериха и месево.
Попробуем исправить ситуацию:
trjconv -f dna_md.xtc -s dna_md.tpr -o dna_fit_1.pdb -skip 20 -fit rot+trans
Теперь хотя бы можно смотреть только в центр окошка, там же и происходят все превращения. dna_pbc_1.pdb
Образование B-формы происходит где-то в районе 4 фрейма.
Для опредления соотвествия между номером модели и временем моделирования найдем в pdb файле выражение "MODEL 50" (цифра 50 это и есть номер модели, т.е. у нас это 4).
Двумя строчками выше можно найти упоминание о времени в моделировании. Это время составляет 800 пикосекунд, или 0,8 наносекунд. Т.е. на 0,8 наносекунде произошел переход А-формы ДНК в В-форму.
Пункт 2
Определим среднеквадратичное отклонение в ходе моделирования. Так как у нас происходит конформационный переход, сначала рассчитаем отклонение в ходе всей симуляции относительно стартовой структуры:g_rms -f dna_md.xtc -s dna_md.tpr -o rms_1
Как видно из графика, среднеквадратиное отклонение находится в районе 0,3 нанометров, или 3 ангстрем.
Рассмотрим среднеквадратичное отклонение относительно каждой предыдущей структуры на растоянии 400 кадров. Если ближе к концу закончился конформационный переход, то отколнение должно уменьшаться.
g_rms -f dna_md.xtc -s dna_md.tpr -o rms_2 -prev 400
В данном случае среднеквадратичное отклонение меньше, чем в первом случае, и составляет примерно 0,2 нанометра, или 2 ангстрема.
Пункт 3
Определим изменение гидрофобной и гидрофильной поверхностей в ходе конформационного перехода:g_sas -f dna_md.xtc -s dna_md.tpr -o sas_dna.xvg
Построим зависимость изменения гидрофобной гидрофильной поверхностей доступных растворителю от времени.
Гидрофобная поверхность:
Как видно из графика, площадь гидрофобной поверхности практически не меняется с течением времени.
Гидрофильная поверхность:
Разброс площади гидрофильной поверхности несколько больше, чем в случае гидрофобной, но, в целом, среднее значение также не меняется.
Т.е. конформационному переходу из одной формы ДНК в другую в этом плане ничего не мешает, ведь площади гидрофобной и гидрофильной поверхностей практически не меняются.
Пункт 4
Традиционным анализом для ДНК является расчёт количества образуемых водородных связей. Если мы будем исследовать связи между ДНК и ДНК, то это будут водородные связи между цепями ДНК. Для конца траектории:g_hbond -f dna_md.xtc -s dna_md.tpr -num hbond_dna
По каноническим представлениям в дуплексе между цепями должно быть 14 водородных связей. В общем-то, конечно, мы видим, что большинство точек графика расположено ровно на прямой со значением 14, но еще много точек лежат на прямых со значениями 13 и 15.Фактически, количество водородных связей не меняется с течением времени.
Не менее интересно будет изучить количество вдородных связей ДНК-Вода:
g_hbond -f dna_md.xtc -s dna_md.tpr -num hbond_dna_sol
Количество водородных связей ДНК и воды лежит в широких пределах от 100 до 120. С течением времени ситуация не сильно меняется. Так что, видимо, количество водородных связей также не может помешать ДНК перейти из одной формы в другую.