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

Даны файлы:
дополнительной топологии для липида DPPC, dppc.itp;
параметры для липидов lipid.itp;
координаты одного липида dppc.gro;
файл-заготовка тополгии системы b.top;
файл параметров для минимизации энергии em.mdp;
файл параметров для "утряски" воды pr.mdp;
файл параметров для молекулярной динамики md.mdp.

На основе одного липида создадим ячейку с 64 липидами: 
genconf -f dppc.gro -o b_64.gro -nbox 4 4 4
С помощью editconf преобразуем dppc.gro и b_64.gro в pdb файлы: 
editconf -f dppc.gro -o dppc.pdb
editconf -f b_64.gro -o b_64.pdb
В одном из хвостов жирная кислота образует цикл - позднее этот "плохой" контакт будет исправлен оптимизацией геометрии. Установим в файле b.top правильное количество липидов в системе - 64.
Сделаем небольшой отступ в ячейке от липидов, чтобы добавить примерно 2500 молекул воды: 
editconf -f b_64.gro -o b_ec -d 0.5
Проведём оптимизацию геометрии системы, чтобы удалить "плохие" контакты молекул. 
grompp -f em -c b_ec -p b -o b_em -maxwarn 2
mdrun -deffnm b_em -v
Начальное значение максимальной силы = 4.37970e+05, конечное = 6.4541919e+02.
Добавим в ячейку молекулы воды типа spc: 
genbox -cp b_em -p b -cs spc216 -o b_s
Проведем "утряску" воды: 
grompp -f pr -c b_s -p b -o b_pr -maxwarn 1
mdrun -deffnm b_pr -v
Переформатируем b_pr.gro и b_s.gro в pdb-формат: b_pr.pdb и b_s.pdb. Изменение в системах - исчез лишний цикл, образованный жирной кислотой; "утрясённая" вода распределена по объему более хаотично.

Скопируем эти файлы на суперкомпьютер. 
ssh skif
mkdir Golikov
exit
cd md
scp -r md/* skif:Golikov/
Запустим тестовое моделирование на суперкомпьютере: 
ssh skif
cd Golikov
grompp -f md -c b_pr -p b -o b_md -maxwarn 1
mpirun  -np 16 -q test -maxtime 5  /home/golovin/progs/bin/mdrun_mpi -deffnm b_md -v
Файл не содержит ошибок.

Запускаем основное моделирование на суперкомпьютере: 
mpirun  -np 16 -maxtime 1200  /home/golovin/progs/bin/mdrun_mpi -deffnm b_md -v

Анализ результатов

  • Силовое поле используемое при построении топологии топологии - ffgmx
  • Размер и форму ячейки - параллелепипед (6.26000х4.44300х5.77800 нм)
  • Минимизация энергии:
    • Алогритм минимизации энергии integrator = l-bfgs
    • Алгоритм расчёта электростатики coulombtype = Cut-off
    • и Ван-дер-Ваальсовых взаимодействий vdw-type = Cut-off
  • Модель, которой описывался растворитель - flexspc
  • Утряска растворителя:
    • Число шагов - 10000
    • Длина шага - 0.001 пс
    • Алгоритм расчёта электростатики coulombtype = pme
    • и Ван-дер-Ваальсовых взаимодействий. vdw-type = Cut-off
    • Алгоритмы термостата Tcoupl = Berendsen
    • и баростата. Pcoupl = no
  • Основной расчёт МД:
    • Время моделирования 7 hours 15 minutes 03 seconds, количество процессоров - 16.
    • Длина траектории - 50нс
    • Число шагов - 10000000
    • Длина шага - 0.05 пс
    • Алгоритм интегратора md
    • Алгоритм расчёта электростатики coulombtype = pme
    • и Ван-дер-Ваальсовых взаимодействий. vdw-type = Cut-off
    • Алгоритмы термостата Tcoupl = v-rescale
    • и баростата Pcoupl = berendsen
1. Визуальный анализ движения молекул. 
trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_1.pdb -skip 20
b_pbc_1.pdb. Проблема - длинные связи, пересекающие ячейку от одного края до другого. Исправляем: 
trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_2.pdb -skip 20 -pbc mol
b_pbc_2.pdb. Теперь все нормально, на анимации видно, как из бислоя периодически вылезают два хвоста.



Бислой начинает образовываться через 10500 фс после начала моделирования.

2. Определим площадь, занимаемую одним липидом. Получим размеры ячейки из траектории. 
   g_traj -f b_md.xtc -s b_md.tpr -ob box_1.xvg
В файле box_1.xvg содержатся размеры ячейки. В первой колонке время, во 2-4 колонках - размер ячейки в виде длин по трём осям. Ось Х является нормалью к поверхности бислоя. Зависимость площади по осям Y и Z от времени, нормированнвя на один липид в слое (т.е. делим на 32), представлена на графике:


Площадь, занимаемая одним липидом до формирования бислоя - 0.8 нм, после формирования - чуть меньше 0.7 нм.

3. Определим изменение гидрофобной и гидрофильной поверхностей в ходе самосборки: 
g_sas -f b_md.xtc -s b_md.tpr -o sas_b.xvg
Зависимость изменения гидрофобной (синим цветом) гидрофильной (красной) поверхностей доступных растворителю от времени:


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

4. Традиционной мерой оценки фазового состояния бифильных молекул является мера порядка. Для анализа нам понадобится специальный индекс файл. Запустим анализ для конца траектории 
   g_order -s b_md -f b_md.xtc -o ord_end.xvg -n sn1.ndx -b 45000 -d P


и для начала траектории 
   g_order -s b_md -f b_md.xtc -o ord_start.xvg -n sn1.ndx -e 5000 -d P


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

Назад