Занятие 8. Анализ результатов моделирование самосборки липидного бислоя.
Пакет программ Gromacs предоставляет много инструментов для анализа траекторий и свойств динамики. Суть любого анализа сводится к пониманию специфики динамики конкретной системы.
Результаты анализа выдаваемые GROMACS имеют расширение xvg ( программа GRACE), но формат самих файлов текстовой, так что построение графиков можно делать в gnuplot или excel.
Для gnuplot надо задать:
set datafile commentschars "#@&"
Часто в файлах есть более, чем два столбца построение первых двух столбцов в gnuplot тогда будет:
plot "file.xvg" using 1:2 with lines
Будем называть результаты анализа согласно общему шаблону:
Tool_system_param, где
Tool- это название программы которой проводили анализ
sytem- это либо b (бислой) либо dna (ДНК)
param- это некое дполнительное описание
Пример : g_rmsd_dna_1
Внимание ! Опции программы анализа вы можете узнать, набрав: имя_программы -h
Общая информация:
coulombtype = Cut-off.
Объединяется область отсечения со списком отсечения соседей (neighborlist cut-off) rlist и отсечения для кулоновских
взаимодействий (Coulomb cut-off) rcoulomb, где rcoulomb ? rlist.
Алгоритм расчета Ван-дер-Ваальсовых взаимодействий.
coulombtype = Cut-off
Объединяется область отсечения со списком отсечения соседей (neighborlist cut-off) rlist и отсечения
Ван-дер-Ваальсовых взаимодействий (VdW cut-off) rvdw, где rvdw ? rlist.
Модель, которой описывался растворитель
implicit_solvent = No
Нет растворителя.
Утряска растворителя:
pr.mdp
10000
Длина шага. 0.001
Алгоритм расчёта электростатики и Ван-дер-Ваальсовых взаимодействий.
coulombtype = pme
Быстрое вычисление электростатики по методу Эвальда Particle-Mesh.
vdw-type = Cut-off
Алгоритмы термостата и баростата.
Tcoupl = Berendsen
Контроль температуры термостатом Берендсена
Pcoupl = no
Без контроля давления. Это означает фиксированный размер ячейки.
Основной расчёт МД (md.mdp)
- Время моделирования, количество процессоров, эффективность маштабирования.
03:18:07 ;16; 21%
Если моделирование окончилось с ошибкой, указать ошибку.
The Z-size of the box (2.000160) times the triclinic skew factor (1.000000) is smaller than
the number of DD cells (2) times the smallest allowed cell size (1.000000)
Длину траектории
50000.0 ps
Число шагов
10000000 (прошло 8008000)
Длина шага
0.005
Алгоритм интегратора.
integrator = md
Алгоритм расчёта электростатики и Ван-дер-Ваальсовых взаимодействий.
coulombtype = pme
vdw-type = Cut-off
Алгоритмы термостата и баростата.
Tcoupl = v-rescale
Контроль температуры по методу "Velocity rescale"
Pcoupl = Berendsen
Контроль давления с помощью экспоненциальной релаксации. Ячейка пересчитывается каждый шаг
АНАЛИЗ ФАЙЛОВ:
- Любой анализ начинается с визуального анализа движений молекул (на вопрос о выводе групп выбераем DPPC).
trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_1.pdb -skip 20
Откроем b_pbc_1.pdb в PyMol , (не забываем включить анимацию). Не устроил результат визуализации - пробуем:
trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_2.pdb -skip 20 -pbc mol
Полученный файл. На мой взгляд образование бислоя наблюдается в моделе 80 (MODEL 80). Время моделирования t= 40000 фемтосекунд:

- Теперь определим площадь занимаемую одним липидом. Для этого получим размеры ячейки из траектории:
g_traj -f b_md.xtc -s b_md.tpr -ob box_1.xvg
В файле box_1.xvg содержатся размеры ячейки.
В первой колонке время, в следующих трёх то, что надо - 5,7639; 13,5574; 2,02057.
Определим какая ось является нормалью к поверхности бислоя.
Измерив параметры бислоя понимаем, что перпендикуляром к бислою является ось OX.
Построим зависимость длин осей OY и OZ ( не нормали к поверхности бислоя) от времени:
Получив изменение площади в зависимости от времени и нормировав это значание на один липид в слое, получаем следующие зависимости:
Плошадь вычислялась как произведение длин по осям, не являющимися нормалями к поверхности бислоя,
и деленным на 32 - молекулы в системе пополам.
Как видно средняя площадь для молекулы примерно 0.75 нм2. Это величина меняется со временем моделирования, скачок в конце связан с появлением группы липидов вне слоя.
- Определим изменение гидрофобной и гидрофильной поверхности в ходе самосборки:
g_sas -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 45000 -d X
Для начала траектории:
g_order -s b_md -f b_md.xtc -o ord_start.xvg -n sn1.ndx -e 5000 -d X
Графики имеют примерно одинаковый вид. Головы липидов подвижны, затем менее подвижная структура хвостов, образующих бислой, а потом - ближе к концам хвостов - атомы очень подвижны внутри самого бислоя. В конце траектории, когда бислой имеет практически сформированную структуру, подвижность атомов молекул липидов падает, по сравнению с началом траектории, когда структуры еще не сформированно.
©
Пискунова Юлия 2011