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

Цель данного занятия ознакомится с возможностями моделирования молекулярной динамики. Используем пакет молекулярной динамики Gromacs. Была выбрана задача моделирования самосборки липидного бислоя. Сначала на основе одного липида созадим ячейку с 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

In [1]:
from IPython.display import Image
Image(filename='b_64.png')
Out[1]:

Рис.1. 64 липида
В текстовом редакторе в файле b.top установим правильное количество липидов в системе: DPPC 64. Сделаем небольшой отступ в ячейке от липидов, что бы добавить примерно 2600 молекул воды:

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

Добавим в ячейку молекулы воды типа 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 формат.

editconf -f b_pr.gro -o b_pr.pdb editconf -f b_s.gro -o b_s.pdb

In [2]:
from IPython.display import Image
Image(filename='b_pr-s.png')
Out[2]:

Рис.2. Один липид из 64. Все остальные липыда расположены идентично. Черный до утряски воды, желтый после.
После утряски воды произошли небольшие изменения в расположении липидов, также произошло изменение положения молекул воды (на рисунке вода не показана).
Последним шагом был запуск модулирования самосборки липидного бислоя на суперкомпьютере:

sbatch -N1 --ntasks-per-node=2 -e error-gpu.log -o output.log -t 350 -p gpu impi /opt/ccoe/gromacs-5.0.4/build/bin/gmx_mpi mdrun -testverlet -deffnm b_md -v

Теперь проанализируем результат моделирования.

Визуальный анализ изображения:

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

На 23 шаге уже сформировалась структура, похожая на бислой (11 нсек)

In [1]:
from IPython.display import Image
Image(filename='23step.png')
Out[1]:

Рис.3. Липидный бислой на 23 шаге моделирования
Определиь площадь занимаемую одним липидом. Для этого вам получим размеры ячейки из траектории.

g_traj -f b_md.xtc -s b_md.tpr -ob box_1.xvg

В полученном файле содержаться размеры ячейки. Поверхности бислоя перпендикурна ось Х, поэтому чтобы получить плоцадь поверхности бислоя необходимо перемножить Y и Z значения ячейки. Для 23 модели (рис.3) получили 23,5531 нм² площадь поверхности бислоя. Теперь построим зависимость площади от времени. И норимруем это значание на один липид в слое (с одной стороны бислоя у нас располагается 32).

In [5]:
with open("box_11.xvg","r") as box:
    time = []
    surface = []
    for line in box:
        line = line.strip().split()
        time = time+[int(line[0])]
        surface = surface+[float(line[2])*float(line[3])/32]
In [17]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

x=np.array(time)
y=np.array(surface)

plt.plot(x, y, "r-",alpha=0.5)
plt.ylabel('Surface area/32, nm^2')
plt.xlabel('Time, ps')
plt.show()
In [7]:
surface[-1]
Out[7]:
0.6203259329625

Согласно этому графику в процессе моделирования со временем бислой компактизуется. Примерно на 40нсек бислой площадь поверхности бислоя перестала изменяться. Конечная площадь поверхности бислоя составила 19,84 нм².
Теперь определим изменение гидрофобной и гидрофильной поверхности в ходе самосборки. И построим зависимость изменения гидрофобной и гидрофильной поверхностей доступных растворителю от времени.

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

In [10]:
with open("sas_b1.xvg","r") as box:
    time_2 = []
    hb_sas = []
    hf_sas = []
    for line in box:
        line = line.strip().split()
        time_2 = time_2+[int(line[0])]
        hb_sas = hb_sas+[float(line[1])]
        hf_sas = hf_sas+[float(line[2])]
In [13]:
x=np.array(time_2)
y_1=np.array(hb_sas)
y_2=np.array(hf_sas)

plt.plot(x, y_1, "r-", alpha = 0.5, label="Hydrophobic")
plt.plot(x, y_2,"y-", alpha=0.5, label="Hydrophylic")
plt.ylabel('Solvent Accessible Surface, nm^2')
plt.xlabel('Time, ps')
plt.legend()
plt.show()

На график видно, что гидрофобная и гидрофильная поверхности доступные растворителю с самого начала моделирования начинают стремительно уменьшаться. С определенного момента (320псек) гидрофобная поверхность становится меньше гидрофильной. Вероятно, именно уменьшение гидрофобной поверхности доступной растворителю (в данном случае воде) является причиной самосборки липидного бислоя.

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

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

In [15]:
with open("ord_start.xvg","r") as order:
    mol_num = []
    s_start = []
    for line in order:
        line = line.strip().split()
        mol_num = mol_num+[int(line[0])]
        s_start = s_start+[float(line[1])]

with open("ord_end.xvg","r") as order:
    s_end = []
    for line in order:
        line = line.strip().split()
        s_end = s_end+[float(line[1])]
In [21]:
x=np.array(mol_num)
y_1=np.array(s_start)
y_2=np.array(s_end)

plt.plot(x, y_1, "ro", alpha = 0.5, label="start")
plt.plot(x, y_2,"yo", alpha=0.5, label="end")
plt.ylabel('S')
plt.xlabel('DSSP atom number')
plt.legend()
plt.show()

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