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

Анализ результатов моделирование самосборки липидного бислоя.

Анализ началя с визуального анализа движений молекул. При вопросе о выводе групп выберали DPPC.

gmx trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_1.pdb -skip 20

Выполнив команду, которая указана выше, получили файл b_pbc_1.pdb и открыли его в PyMOL. Результат был не красивый поэтому выполнили команду:

gmx trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_1.pdb -skip 20 -pbc mol

In [2]:
import time
from IPython.display import Image
In [3]:
Image(filename="C:\\Users\\kseni\\Desktop\\головин 2019\\bpbc1.png")
Out[3]:

Бислой образовался на 21 модели. Для опредления соотвествия между номером модели и временем моделирования найшли в b_pbc_1.pdb выражение "MODEL 21" и двумя строчками выше время моделирования t = 10000.00000

Для определения площади занимаемой одним липидом, получили размеры ячейки из траектории

gmx traj -f b_md.xtc -s b_md.tpr -ob box_1.xvg

В полученном файле box_1.xvg содержатся размеры ячейки. Нормалью к поверхности является ось X.

Определим площадь, занимаемую одним липидом.

In [5]:
import pandas as pd
In [6]:
with open('box_1.xvg', 'r') as f:
    lines = [item.strip('\n') for item in f.readlines() 
             if (not item.startswith("@")) and (not item.startswith("#"))]
times = [item.split('\t')[1] for item in lines]
sq = [float(item.split('\t')[3]) * float(item.split('\t')[4]) / 32 for item in lines]
In [9]:
data = pd.DataFrame({'Squares':sq,'times':times})
data.head()
Out[9]:
Squares times
0 0.802239 0
1 0.805691 25
2 0.807099 50
3 0.810258 75
4 0.813023 100
In [12]:
import matplotlib.pyplot as plt
In [15]:
data.plot(x="times", y = "Squares", c='red')
plt.show()

Видно, что с течением времени площадь, занимаемая одним липидом уменьшается. Это из-за образования бислоя, так как липиды уплотняются и структурируются, занимаемая одним липидом площадь уменьшается.

Определили изменение гидрофобной и гидрофильной поверхности в ходе самосборки.

g_sas -f b_md.xtc -s b_md.tpr -o sas_b.xvg

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

In [16]:
with open("sas_b.xvg", 'r') as f:
    lines = [item.strip('\n') for item in f.readlines() 
             if (not item.startswith("@")) and (not item.startswith("#"))]
times = [int(item.split()[0]) for item in lines]
hphob = [float(item.split()[1]) for item in lines]
hphil = [float(item.split()[2]) for item in lines]
In [17]:
data = pd.DataFrame({'times':times, 'Hydrophobic':hphob, 'Hydrophilic':hphil})
In [18]:
data.head()
Out[18]:
times Hydrophobic Hydrophilic
0 0 206.446 135.480
1 25 202.647 134.835
2 50 200.633 136.411
3 75 197.463 144.996
4 100 194.358 149.460
In [54]:
data.plot(x="times", y = ["Hydrophobic", "Hydrophilic"])
Out[54]:
<matplotlib.axes._subplots.AxesSubplot at 0x25e0e02e710>
In [55]:
plt.show()

В начале липиды расположены отдельно, поэтому гидрофобная поверхность имеет максимальное значение. Но со временем липиды компактизуются сначала в подобную мицелле структуру, а затем в бислой,как в наиболее выгодную структуру. Где-то в районе 10000 пс кривая выходит на плато, это время образования бислоя. Гидрофильная поверхность на всем протяжении больше, чем гидрофобная. Однако она тоже немного снижается со временем, так как уменьшается площадь поверхности, доступная для растворителя, из-за того что липиды компактизуются.

Традиционной мерой оценки фазового состояния бифильных молекул является мера порядка. Для анализа скачали файл sn1.ndx. Запустили анализ для конца траектории:

gmx  order -s b_md -f b_md.xtc -o ord_end.xvg -n sn1.ndx -b 25000 -d X

и для начала траэктории:

gmx order -s b_md -f b_md.xtc -o ord_start.xvg -n sn1.ndx -e 5000 -d X
In [57]:
with open("ord_start.xvg", 'r') as f:
    lines = [item.strip('\n') for item in f.readlines() 
             if (not item.startswith("@")) and (not item.startswith("#"))]
xax = [int(item.split()[0]) for item in lines]
start = [float(item.split()[2]) for item in lines]

with open("ord_end.xvg", 'r') as f:
    lines = [item.strip('\n') for item in f.readlines() 
             if (not item.startswith("@")) and (not item.startswith("#"))]
end = [float(item.split()[2]) for item in lines]
In [67]:
data = data = pd.DataFrame({'End':end,'Start':start, 'xax':xax})
data.head()
Out[67]:
End Start xax
0 -0.080060 -0.070927 1
1 -0.095122 -0.073118 2
2 -0.097271 -0.089820 3
3 -0.059976 -0.086991 4
4 -0.122472 -0.080438 5
In [68]:
data.plot(x="xax", y = ["Start", "End"])
Out[68]:
<matplotlib.axes._subplots.AxesSubplot at 0x25e0e3a6ac8>
In [69]:
plt.show()

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

Информация о системе:

Силовое поле используемое при построении топологии топологии - ffgmx

Заряд системы - 0

Размер и форму ячейки 6.89978нм х 4.89708нм х 4.74459нм, прямоугольный паралиллепипед

Минимизация энергии:

Алгоритм минимизации энергии - l-bfgs
Алгоритм расчёта электростатики и Ван-дер-Ваальсовых взаимодействий -  Cut-off

Модель, которой описывался растворитель - flexspc

Утряска растворителя:

Для биополимеров, укажите параметр который обуславливает неподвижность биополимера.
Число шагов 1000
Длина шага 0.0002
Алгоритм расчёта электростатики и Ван-дер-Ваальсовых взаимодействий - pme, Cut-off
Алгоритмы термостата - Berendsen и баростата - no

Основной расчёт МД:

Время моделирования - 50000 ps 
Число шагов - 10000000
Длина шага - 0.005 ps
Алгоритм интегратора - md
Алгоритм расчёта электростатики и Ван-дер-Ваальсовых взаимодействий - pme, Cut-off
Алгоритмы термостата и баростата - v-rescale, Berendsen