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

Цель занятия - ознакомиться с возможностями моделирования молекулярной динамики.
Из предложенных 4 задач практикума была выбрана следующая:

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

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

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

  3. На основе одного липида созадим ячейку с 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
    Геометрия системы неудачна: в одном из хвостов фосфолипида жирная кислота образует цикл.

  4. Установим в файле b.top правильное количество липидов в системе (64).

  5. Сделаем небольшой отступ в ячейке от липидов, чтобы добавить примерно 2500 молекул воды: 
    editconf -f b_64.gro -o b_ec -d 0.5
  6. Проведём оптимизацию геометрии системы, чтобы удалить "плохие" контакты молекул (циклы в остатках жирных кислот): 
    grompp -f em -c b_ec -p b -o b_em -maxwarn 2
mdrun -deffnm b_em -v
    Начальное значение максимальной силы равно 4.37970e+05, конечное (после оптимизации) - 6.4541919e+02.

  7. Добавим в ячейку молекулы воды типа spc: 
    genbox -cp b_em -p b -cs spc216 -o b_s
  8. Проведем "утряску" воды: 
    grompp -f pr -c b_s -p b -o b_pr -maxwarn 1
mdrun -deffnm b_pr -v
  9. Переформатируем b_pr.gro и b_s.gro в pdb-формат и сохраним результаты в файлах b_pr.pdb и b_s.pdb соответственно. Теперь геометрия системы оптимизирована: никаких лишних циклов не наблюдается. Кроме того, видна большая разница между системами без "утрясенной" воды (b_s.pdb) и с "утрясенной" водой (b_pr.pdb). Во втором случае вода распределена куда более хаотично и, как следствие, реалистично. Более того, она занимает большое пространство между молекулами липидов.

  10. Скопируем эти файлы на суперкомпьютер. Для этого сначала зайдем на него и создадим папку: 
    ssh skif
mkdir Tyshkovskiy
exit
  11. Теперь скопируем файлы: 
    cd Tyshkovskiy
scp -r md/* skif:Tyshkovskiy/
  12. Запустим тестовое моделирование на суперкомпьютере: 
    ssh skif
cd Tyshkovskiy
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
    Файл не содержит ошибок.

  13. Наконец, запускаем основное моделирование на суперкомпьютере 
    mpirun  -np 16 -maxtime 1200  /home/golovin/progs/bin/mdrun_mpi -deffnm b_md -v
    Номер моей задачи - 240913.

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

  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. В данном случае все выглядит вполне реалистично. Длинного бислоя не образуется, однако образуется достаточно вытянутая мицелла. Мицелла образуется примерно с 25 номера модели. Это соответствует времени 12500 фемтосекунд.

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

    Зависимость была построена в Excel.
    Таким образом, как видно, оптимальная площадь одного липида в мицелле равна около 0.8 квадратных нанометров. Как видно из графика, в середине моделирования молекулы липидов сгруппировались достаточно плотно (это соответствует примерно началу образования мицеллы). После этого мицелла стала вытягиваться и, соответственно, стала расти площадь, приходящаяся на каждую молекулу липида в отдельности.

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

    Таким образом, при образовании мицеллы (и бислоя) происходит уменьшение как гидрофобной, так и гидрофильной поверхностей. Такая компактизация приводит к уменьшения энергии системы в водном растворителе (что и является движущей силой образования мицеллы и бислоя).

  4. Традиционной мерой оценки фазового состояния бифильных молекул является мера порядка. Для анализа был скачен специальный индекс-файл sn1.ndx. После этого был запущен сам анализ. Для конца траектории: 
    g_order -s b_md -f b_md.xtc -o ord_end.xvg -n sn1.ndx -b 35000 -d X
    Для начала траектории: 
    g_order -s b_md -f b_md.xtc -o ord_start.xvg -n sn1.ndx -e 5000 -d X
    График с изображением меры порядка для разных атомов липида (с головки до хвоста) для начала траектории представлен ниже:

    Такой же график для конца траектории:

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

Назад