Молекулярная динамика биологических молекул в GROMACS

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

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

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

Типы файлов:

• *gro - файл с координатами системы.
• *top - файл с описанием ковалентных и нековалентных взаимодействий в молекулах.
• *mdp - файл с описанием параметров для работы молеклярно-механического движка.
• *tpr - файл для молеклярно-механического движка по сути есть объединение gro, top и mdp.
• *{trr,xtc} - файл с координатами после рассчёта.

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

Программы запускаются в консоли c префиксом gmx , флаги запуска программ начинаются с -, например -f. Как правило после флага следует либо имя файла либо значение параметра. Смотрите примеры ниже.

• editconf - манипуляция форматом координат и самими координатами. Пример:

%%bash 
gmx editconf -f my.gro -o my.pdb

• genbox - наполнение ячейки растворителем.Пример:

%%bash 
gmx genbox -cp my_protein.gro -cs mysolvent.gro -p my.top -o my_solvated_protein.gro

• grompp - объединение и проверка gro, top и mdp в tpr.

%%bash
grompp -f my.mdp -c my.gro -p my.top

• mdrun - молеклярно-механический движок. На входе принимает tpr файл

Объекты для практикума

На этом занятии Вам предлагается 6 различных систем для моделирования. Перейдите по ссылке для подробных инструкций по выполнен ию каждого задания.

• Моделирование самосборки липидного бислоя из случайной стартовой конформации.
• Моделирование поведения ДНК в формальдегиде.
• Моделирование перехода А-формы ДНК в В-форму в воде.
• Моделирование поведения короткого пептида в формальдегиде.
• Brand New !! QM/MM + метадинамика SN2
• Brand New !! Coarse Grain Martinni + метадинамика

---

Моделирование самосборки липидного бислоя из случайной стартовой конформации.

Создайте рабочую директорию mm-prac2

Вам даны файлы:

• дополнительной топологии для липида DPPC: /home/golovin/teach/prac2/dppc/dppc.itp
• параметры для липидов: /home/golovin/teach/prac2/dppc/lipid.itp
• координаты одного липида: /home/golovin/teach/prac2/dppc/dppc.gro
• Файл-заготовка тополгии системы: /home/golovin/teach/prac2/dppc/b.top
• файл праметров для минимизации энергии: /home/golovin/teach/prac2/dppc/em.mdp
• файл праметров для "утряски" воды pr.mdp: /home/golovin/teach/prac2/dppc/pr.mdp

• файл праметров для молекулярной динамики: /home/golovin/teach/prac2/dppc/md.mdp

  Скопируйте их в свою рабочую директорию

• На основе одного липида созадим ячейку с 64 липидами.

%%bash
gmx genconf -f dppc.gro -o b_64.gro -nbox 4 4 4

• Просмотрите результат в NGLview. Приведите результат в отчёте.

import nglview as nv
nv.show_structure_file('b_64.gro')

• В текстовом редакторе в файле b.top установите правильное количество липидов в системе. Это можно сделать прямо в Jupyter.
• Проведём оптимизацию геометрии системы, что бы удалить "плохие" контакты молекул.

%%bash
gmx editconf -f b_64 -o b_ec -d 1.5
gmx grompp -f em -c b_ec -p b -o b_em -maxwarn 2
gmx mdrun -deffnm b_em -v

• Отметье в отчёте изменение максимальной силы в ходе оптимизации геометрии. Занесите начальное и конечное значение максимальной силы.
• Добавим в ячейку молекулы воды типа spc.

%%bash
gmx solvate -cp b_em -p b -cs spc216 -o b_s

• Проведём утряску воды:

%%bash
gmx grompp -f pr -c b_s -p b -o b_pr -maxwarn 1
gmx mdrun -deffnm b_pr -v

• вероятно у вас будет взрыв системы тогда:

%%bash
gmx grompp -f em -c b_s -p b -o b_empr -maxwarn 1
gmx mdrun -deffnm b_empr -v

• и снова:

%%bash
gmx grompp -f pr -c b_empr -p b -o b_pr -maxwarn 1
gmx mdrun -deffnm b_pr -v

• Сравните визуально изменеия в системах. Занесите наблюдения в отчёт.
• Запускаем тестовое моделирование молекулярной динамики.

%%bash
alias gmx=/home/golovin/progs/gmx/bin/gmx
gmx grompp -f md -c b_pr -p b -o b_md -maxwarn 1
sbatch -N1 --ntasks-per-node=1 -e error.log -o output.log -t 5 ompi /home/golovin/progs/gmx/bin/gmx mdrun  -deffnm b_md -v -nsteps 1000 -ntomp 4

Просмотреть ход счёта можно в файле output.log

• Запускаем основное моделирование

%%bash 
sbatch -N1 --ntasks-per-node=1 -e error.log -o output.log -t 600 ompi /home/golovin/progs/gmx/bin/gmx mdrun  -deffnm b_md -v -ntomp 4

---

Моделирование поведения ДНК в формальдегиде.

Создайте рабочую директорию. Вам даны файлы в директории :

/home/golovin/teach/prac2/dna-melt/

• Координаты дуплекса ДНК, dna.pdb
• Файл с ячейкой уравновешеных молекул формамида, fam_em.gro
• Файл дополнительной топологии для формамида, fam.itp
• файл праметров для минимизации энергии em.mdp
• файл праметров для "утряски" воды pr.mdp
• файл праметров для молекулярной динамики md.mdp

1. Теперь построим файл топологии системы в силовом поле amber99sb и файл с координатами в формате Gromacs. Предполагается, что структура дуплекса находится в файле dna.pdb.

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

• Если вы получили сообщение об ошибке, то удалите 5' фосфаты из структуры дуплекса. Их два на 5' конце кажой цепи. Вам также надо изменить имена нуклеотидов добвавив D к названию нуклеотида," Т"->"DT" (пример для vim: :%s/ ([GATC]) ([AB])/ D\1 \2/). И наконец замените имя атома "С5М" на " С7".

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

%%bash
gmx editconf -f dna.gro -o dna_ec -d 1.5

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

%%bash
grompp -f em -c dna_ec -p dna -o dna_em -maxwarn 1
mdrun -deffnm dna_em -v

• Отметье в отчёте изменение максимальной силы в ходе оптимизации геометрии. Знаесите начальное и конечное значение максимальной силы.

1. Добавим в ячейку молекулы формамида. Внимательно прочитайте вывод программы и обязательно запомните количество добавле нных молекул формамида.

%%bash
gmx genbox -cp dna_em -p dna -cs fam_em.gro -o dna_s

1. Теперь надо изменить в текстовом редакторе файл тополгии dna.top. После строчки:

; Include forcefield parameters

• добавим #include "fam.itp"'''

  ; Include forcefield parameters #include "fam.itp"

• Добавим количество молекул формамида в запись [ molecules ]'''

• было:

  [ molecules ] ; Compound #mols DNA 1

• стало:

  [ molecules ] ; Compound #mols DNA 1 FAM 2426

где 2426 надо заменить на количество из пункта 8

1. Нейтрализуем заряд системы. Это делаем в два шага: строим tpr и запускаем genion. В выводе grompp обратите внимание н а информацию о заряде системы.

%%bash 
gmx grompp -f em -p dna -c dna_s -o dna_s
gmx genion -s dna_s -o dna_si -p dna -np X

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

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

%%bash 
gmx grompp -f pr -c dna_si -p dna -o dna_pr -maxwarn 1
gmx mdrun -deffnm dna_pr -v

8. Cравните визуально изменеия в системах dna_pr.gro и dna_si.gro . Занесите наблюдения в отчёт. 

1. Запускаем тестовое моделирование молекулярной динамики.

%%bash
alias gmx=/home/golovin/progs/gmx/bin/gmx
gmx grompp -f md -c dna_pr -p dna -o dna_md -maxwarn 1
sbatch -N1 --ntasks-per-node=1 -e error.log -o output.log -t 5 ompi /home/golovin/progs/gmx/bin/gmx mdrun  -deffnm dna_md -v -nsteps 1000 -ntomp 4

Просмотреть ход счёта можно в файле output.log

• Запускаем основное моделирование

%%bash 
sbatch -N1 --ntasks-per-node=1 -e error.log -o output.log -t 600 ompi /home/golovin/progs/gmx/bin/gmx mdrun  -deffnm dna_md -v -ntomp 4

---

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

Создайте рабочую директорию. Вам даны файлы в директории :

/home/golovin/teach/prac2/dna-conf/

• Координаты дуплекса ДНК, dna.pdb
• файл праметров для минимизации энергии em.mdp
• файл праметров для "утряски" воды pr.mdp
• файл праметров для молекулярной динамики md.mdp

1. Теперь построим файл топологии системы в силовом поле amber99sb и файл с координатами в формате Gromacs. Предполагается, что структура дуплекса находится в файле dna.pdb.

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

• Если вы получили сообщение об ошибке, то удалите 5' фосфаты из структуры дуплекса. Их два на 5' конце кажой цепи. Вам также надо изменить имена нуклеотидов добвавив D к названию нуклеотида," Т"->"DT" (пример для vim: :%s/ ([GATC]) ([AB])/ D\1 \2/). И наконец замените имя атома "С5М" на " С7".

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

%%bash
gmx editconf -f dna.gro -o dna_ec -d 1.5

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

%%bash
gmx grompp -f em -c dna_ec -p dna -o dna_em -maxwarn 1
gmx mdrun -deffnm dna_em -v

• Отметье в отчёте изменение максимальной силы в ходе оптимизации геометрии. Знаесите начальное и конечное значение максимальной силы.

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

%%bash
gmx genbox -cp dna_em -p dna -cs  -o dna_s

1. Нейтрализуем заряд системы. Это делаем в два шага: строим tpr и запускаем genion. В выводе grompp обратите внимание на информацию о заряде системы.

%%bash 
gmx grompp -f em -p dna -c dna_s -o dna_s
gmx genion -s dna_s -o dna_si -p dna -np X

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

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

%%bash 
gmx grompp -f pr -c dna_si -p dna -o dna_pr -maxwarn 1
gmx mdrun -deffnm dna_pr -v

8. Cравните визуально изменеия в системах dna_pr.gro и dna_si.gro . Занесите наблюдения в отчёт. 

1. Запускаем тестовое моделирование молекулярной динамики.

%%bash
alias gmx=/home/golovin/progs/gmx/bin/gmx
gmx grompp -f md -c dna_pr -p dna -o dna_md -maxwarn 1
sbatch -N1 --ntasks-per-node=1 -e error.log -o output.log -t 5 ompi /home/golovin/progs/gmx/bin/gmx mdrun  -deffnm dna_md -v -nsteps 1000 -ntomp 4

Просмотреть ход счёта можно в файле output.log

• Запускаем основное моделирование

%%bash 
sbatch -N1 --ntasks-per-node=1 -e error.log -o output.log -t 5 ompi /home/golovin/progs/gmx/bin/gmx mdrun  -deffnm dna_md -v -ntomp 4

---

Моделирование поведения короткого пептида в формальдегиде.

Создайте рабочую директорию. Вам даны файлы в директории :

/home/golovin/teach/prac2/pep-melt/

• Координаты атомов пептида, 2lx1.pdb
• Файл с ячейкой уравновешеных молекул формамида, fam_em.gro
• Файл дополнительной топологии для формамида, fam.itp
• файл праметров для минимизации энергии em.mdp
• файл праметров для "утряски" воды pr.mdp
• файл праметров для молекулярной динамики md.mdp

1. Теперь построим файл топологии системы в силовом поле amber99sb и файл с координатами в формате Gromacs. Предполагается, что структура дуплекса находится в файле 2lx1.pdb.

%%bash
gmx pdb2gmx -f 2xl1.pdb -o pep -p pep -ff amber99sb -water tip3p

• Если вы получили сообщение об ошибке, то удалите 5' фосфаты из структуры дуплекса. Их два на 5' конце кажой цепи. Вам также надо изменить имена нуклеотидов добвавив D к названию нуклеотида," Т"->"DT" (пример для vim: :%s/ ([GATC]) ([AB])/ D\1 \2/). И наконец замените имя атома "С5М" на " С7".

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

%%bash
gmx editconf -f pep.gro -o pep_ec -d 1.5

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

%%bash
gmx grompp -f em -c pep_ec -p pep -o pep_em -maxwarn 1
gmx mdrun -deffnm pep_em -v

• Отметье в отчёте изменение максимальной силы в ходе оптимизации геометрии. Знаесите начальное и конечное значение максимальной силы.

1. Добавим в ячейку молекулы формамида. Внимательно прочитайте вывод программы и обязательно запомните количество добавле нных молекул формамида.

%%bash
gmx genbox -cp pep_em -p pep -cs pep_em.gro -o pep_s

1. Теперь надо изменить в текстовом редакторе файл тополгии dna.top. После строчки:

; Include forcefield parameters

• добавим #include "fam.itp"'''

  ; Include forcefield parameters #include "fam.itp"

• Добавим количество молекул формамида в запись [ molecules ]

• было:

  [ molecules ]

  ; Compound #mols

  Protein 1

• стало:

  [ molecules ]

  ; Compound #mols

  Protein 1

  FAM 2426

где 2426 надо заменить на количество из пункта 8

1. Нейтрализуем заряд системы. Это делаем в два шага: строим tpr и запускаем genion. В выводе grompp обратите внимание н а информацию о заряде системы.

%%bash 
gmx grompp -f em -p pep -c pep_s -o pep_s
gmx genion -s pep_s -o pep_si -p pep -neutral

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

%%bash 
gmx grompp -f pr -c dna_si -p dna -o dna_pr -maxwarn 1
gmx mdrun -deffnm dna_pr -v

1. Cравните визуально изменеия в системах pep_pr.gro и pep_si.gro . Занесите наблюдения в отчёт.

1. Запускаем тестовое моделирование молекулярной динамики.

%%bash
alias gmx=/home/golovin/progs/gmx/bin/gmx
gmx grompp -f md -c pep_pr -p pep -o pep_md -maxwarn 1
sbatch -N1 --ntasks-per-node=1 -e error.log -o output.log -t 5 ompi /home/golovin/progs/gmx/bin/gmx mdrun  -deffnm pep_md -v -nsteps 1000 -ntomp 4

Просмотреть ход счёта можно в файле output.log

1. Запускаем основное моделирование

%%bash 
sbatch -N1 --ntasks-per-node=1 -e error.log -o output.log -t 600 ompi /home/golovin/progs/gmx/bin/gmx mdrun  -deffnm pep_md -v -ntomp 4

---

Brand New !! QM/MM + метадинамика SN2

Создайте рабочую директорию. Вам даны файлы в директории :

/home/golovin/teach/prac2/qmmm/

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

1. Создадим CH,,3,,Br

%%bash
echo "C[Cl] gg" > 1.smi
obgen 1.smi -ff UFF >1.mol
obabel -imol 1.mol -opdb 1.pdb

• Добавим ион F в mol файл. Достаточно скопировть строчку с CL в конец и изменить Х-координату на близкую к 0. Проверьте визуально

• Набросаем ручками топологию:

#include "./amber99sb-ildn.ff/forcefield.itp"
#include "./amber99sb-ildn.ff/tip3p.itp"
[ moleculetype ]
; molname       nrexcl
MOL             2
[ atoms ]
; id  at type     res nr  res name  at name  cg nr  charge    mass
  1   CT          1       MOL         C       1       0    12.00000
  2   Cl          1       MOL        Cl       2       0    35.90000
  3   H2          1       MOL        H1       3       0     1.00800
  4   H2          1       MOL        H1       4       0     1.00800
  5   H2          1       MOL        H1       5       0     1.00800
  6   F           1       MOL         F       6       0    19.00000
[ bonds ]
; i     j       funct   
1       2       5       
1       3       5       
1       4       5       
1       5       5       

[ system ]
; Name
Protein in water
[ molecules ]
; Compound     
MOL         1

• Сделаем gro:

%%bash
gmx editconf -f 1.pdb  -o 1.gro -d 1 -c

• Grompp и запуск молекулярной динамики:

%%bash
export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:/home/golovin/progs/lib/
export GMX_QMMM_VARIANT=1
export gmx /home/golovin/progs/gmx-dftbplus-012021/bin/gmx
gmx grompp -f mdqm.mdp  -c 1 -p 1 -o vac 
gmx mdrun  -deffnm vac -nt 1

• Опишите протекание реакции и оцените изменение энергии с помощью 'gmx energy'

• Добавим воды:

%%bash
gmx solvate -cp  1.gro -cs -p topol.top -o 1s

• Очень важно оптимизировать геометрию воды вокруг реакции. Для начала сделаем индекс файл с описанием группы атомов в реакции.

%%bash
gmx make_ndx -f 1s.gro

• И запретим им двигаться, так как нам нужно оптимизировать только воду. Укажите в отчете какие директивы в mdqm-s.mdp файле это делают.

%%bash
gmx grompp -f em -c 1s -p topol.top -o em
gmx mdrun -deffnm em -nt 1 -v

• Теперь можно запустить QM/MM молекулярную динамику

%%bash
alias gmx=/home/golovin/progs/gmx-dftbplus-012021/bin/gmx
gmx grompp -f mdqm-s.mdp -p topols -c em.gro -o mds  -n
gmx mdrun -deffnm mds -nt 1 -v

---

Brand New !! Coarse Grain Martinni + метадинамика

• Скачаем KALP пептид и утилиты

    /home/golovin/teach/prac2/cg

• Превратим его в CG описание ( топология и координаты )

%%bash 
./martinize.py -f kalp.pdb -o topo.top -x conf.pdb -n cg

• Скрипт генерирует старую ссылку на силовое поле, поможем ему:

%%bash
ln -s martini_v2.2.itp martini.itp

• Теперь создадим бокс вокруг пептида:

%%bash
gmx editconf -f conf.pdb -box 7 7 7 -c -o ec

• И добавим 128 липидов dppc:

%%bash
gmx insert-molecules -f ec -ci dppc.gro  -o s -nmol 128

• Добавить также придется и топологию. Догадайтесь как? это связано с директивой include
• Конечно, надо добавить воду

%%bash
gmx grompp -f minimization.mdp  -p -c solv.gro -o em
gmx mdrun -deffnm em -v

• Если оптимизация произошла успешно, то можно запустить молекулярную динамику

%%bash
gmx grompp -f dynamic.mdp  -c em.gro -p  -o md.tpr
gmx mdrun -deffnm md -nt 1 -v -nsteps 3000000

• Сравните ваш результат с http://cgmartini.nl/index.php/tutorials-general-introduction/proteins#membrane-protein , сделайте заключение об причинах, которые привели к разнице результатов.

Моделирование прохождения пептида сквозь бислой с помощью метадинамики

В результате предидущего задания вы получили бислой в воде, где пептид KALP лежит на интерфесе вода-липид

Давайте найдем стабильное положение пептида в этой системе с помощью метадинамики.

Для запуска достаточно использовать Gromacs с Plumed.

/home/golovin/progs/gmx-dftbplus-012021/bin/gmx

• Создадим файл с инструкциями для метадинамики. Представим, что две коллективные переменные определяют проникновение пептида в бислой. Это растояние между центрами масс бислоя и пептида и угол между центральным остатком, конечным остатком и нормалью к бислою.
• Итак сначала изучим как задать группу атомов https://plumed.github.io/doc-v2.7/user-doc/html/_group.html

• Запись будет иметь вид (ВНИМАНИЕ, что бы это работало читаем https://plumed.github.io/doc-v2.7/user-doc/html/_d_i_s_t_a_n_c_e.html):

  pep: COM ATOMS=1-x
  bl: COM ATOMS=y-100 
  dist: DISTANCE=pep,bl 
  ang:ANGLE ATOMS=bl,pep,1

• Теперь пропишем сами параметры метадиамики:

METAD ARG=dist SIGMA=0.2 HEIGHT=0.5 PACE=100 LABEL=metad

• И добавим запись переменных в файл текстовый COLVAR

PRINT ARG=dist,metad.bias   STRIDE=10 FILE=COLVAR

• Всё это запишем в один файл с расширением dat
• Создадим новый tpr

%%bash
gmx grompp -f dynamic.mdp  -c md.gro -p  -o meta.tpr

• Запустим, если молча упало, то смотрим в meta.log

%%bash
gmx mdrun -deffnm meta -nt 1 -v -plumed meta.dat -nsteps 3000000