На главную страницу
На страницу шестого семестра

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

Задание 1

Было выбрано задание моделирования самосборки липидного бислоя.
Даны файлы:
дополнительной топологии для липида 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 zhuravlka
exit
cd md
scp -r md/* skif:zhuravlka/
Запустим тестовое моделирование на суперкомпьютере: 
ssh skif
cd zhuravlka
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.60992х7.4551х3.410333 нм)
  • Минимизация энергии:
    • Алогритм минимизации энергии integrator = l-bfgs
    • Алгоритм расчёта электростатики coulombtype = Cut-off
    • и Ван-дер-Ваальсовых взаимодействий vdw-type = Cut-off
  • Модель, которой описывался растворитель - flexspc
  • Утряска растворителя:
    • Число шагов - 10000
    • Длина шага - 0.001 пс
    • Алгоритм расчёта электростатики coulombtype = pme
    • и Ван-дер-Ваальсовых взаимодействий. vdw-type = Cut-off
    • Алгоритмы термостата Tcoupl = Berendsen
    • и баростата. Pcoupl = no
  • Основной расчёт МД:
    • Время моделирования 6 hours 39 minutes 40 seconds, количество процессоров - 16.
    • Длина траектории - 35 200 пс
    • Число шагов - 7044900
    • Длина шага - 0.005 пс
    • Алгоритм интегратора md
    • Алгоритм расчёта электростатики coulombtype = pme
    • и Ван-дер-Ваальсовых взаимодействий. vdw-type = Cut-off
    • Алгоритмы термостата Tcoupl = v-rescale
    • и баростата Pcoupl = berendsen
Начинаем анализ с визуального движения молекул.
trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_1.pdb -skip 20
Cчитаем , что бислой начинает образовываться с 23 модели, т.е. с t= 11000.00000 пс.

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

В структуре преобладают гидрофильные поверхности. В ходе моделирования обе поверхности уменьшаются. В начале гидрофобная поверхность наибольшая, в это время молекулы липида равномерно распределены в объёме сольвента, дальнейшее уменьшение гидрофобной поверхности соответствует увеличению количества гидрофобных контактов между алифатическими цепями фосфолипида. Динамика гидрофобной поверхности в начале- резкое схлопывание. Фосфолипидный слой собирается.
Традиционной мерой оценки фазового состояния бифильных молекул является мера порядка. Создадим специальный индекс-файл для анализа.Для конца траектории:
g_order -s b_md -f b_md.xtc -o ord_end.xvg -n sn1.ndx -b 35200 -d x
И для начала траектории
g_order -s b_md -f b_md.xtc -o ord_start.xvg -n sn1.ndx -e 5000 -d x

Первый график- для конца, второй - для начала. Мы получили немного не характерные зависимости.Хвост должен двигаться свободней головки. Судя по графикам, у нас возможно ещё не закончилось моделирование сборки, и в конце процесс ещё продолжается. В начале головка даже подвижнее хвоста, да и в конце тоже подвижность, пусть и меньшая сохраняется.


© Zhuravleva Katya, 2009