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

Цель данного занятия - ознакомиться с возможностями моделирования молекулярной динамики. Команды выполнялись в командной строке, а основные вычисления - на суперкомпьютере Ломоносов-2. Здесь показаны результаты вычислений и выводы.

На основе одного липида создадим ячейку с 64 липидами (дипальмитоилфосфатидилхолин), c помощью editconf преобразуем dppc.gro и b_64.gro в pdb файлы, посмотрим на них в PyMol. Видим одну молекулу липида и 64 его копии.

In [ ]:
genconf -f dppc.gro -o b_64.gro -nbox 4 4 4
editconf -f dppc.gro -o dppc.pdb
editconf -f b_64.gro -o b_64.pdb
In [3]:
from IPython.display import Image
Image(filename='dppc.png')
Out[3]:
In [4]:
Image(filename='b_64.png')
Out[4]:

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

In [ ]:
editconf -f b_64.gro -o b_ec -d 0.5

Проведём оптимизацию геометрии системы, чтобы удалить "плохие" контакты молекул.

In [ ]:
grompp -f em -c b_ec -p b -o b_em -maxwarn 2
mdrun -deffnm b_em -v

Вот начальное и конечное значение максимальной силы:

Step= 0, Dmax= 2.0e-02 nm, Epot= 4.74007e+05 Fmax= 4.37970e+05, atom= 1842

Step= 51, Dmax= 1.9e-06 nm, Epot= 3.20911e+03 Fmax= 6.16887e+02, atom= 2765

Соответственно, изменение максимальной силы в ходе оптимизации геометрии равно 1,79.

Добавим в ячейку молекулы воды типа spc, проведём "утряску" воды:

In [ ]:
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 формат. И изучим визуально в PyMol изменения в файлах.

Появились молекулы воды, что и следовало ожидать. В b_s вода распределена более равномерно вокруг липидов, в b_pr ее меньше по бокам.

In [7]:
Image(filename='b_s.png')
Out[7]:
In [8]:
Image(filename='b_pr.png')
Out[8]:

Анализ молекулярной динамики

  • Силовое поле, используемое при построении топологии: lipid.itp
  • Заряд системы: 0 из-за амфифильности молекулы
  • Размер и форма ячейки: 6.26x4.44x5.78, параллелепипед
  • Минимизация энергии
    • Алгоритм минимизации энергии: integrator = steep
    • Алгоритм расчёта электростатики и Ван-дер-Ваальсовых взаимодействий: Cut-off
  • Модель, которой описывался растворитель: spc
  • Утряска растворителя
    • Число шагов: 1000
    • Длина шага: 0.002 ps
    • Алгоритм расчёта электростатики и Ван-дер-Ваальсовых взаимодействий: Cut-off
    • Алгоритмы термостата и баростата: Tcoupl = V-rescale, Pcoupl = no ; метод контроля температуры Berendsen (anisotropic)
  • Основной расчёт МД
    • Длина траектории: 50000
    • Число шагов: 10000000
    • Длина шага: 0.005
    • Алгоритм интегратора: md
    • Алгоритм расчёта электростатики и Ван-дер-Ваальсовых взаимодействий: pme и Cut-off
    • Алгоритмы термостата и баростата: Tcoupl = v-rescale, Pcoupl = Berendsen (semiisotropic)

Скопируем файлы на суперкомпьютер, запустим тестовое, а затем и основное моделирование. После того, как посчиталось, переделали .tpr файл, так как считалось оно по новой версии этого файла, а у нас - старый gromacs. При вопросе о выводе групп выбрали DPPC.

Скопировала b_pbc_1.pdb к себе с Ломоносова: Посмотрим, что получилось.

In [ ]:
scp b_pbc_1.pdb sunnymouse@kodomo.fbb.msu.ru:term8/Vasyutkina

Посмотрим, что получилось. Видим, что молекулы "вылезают" за границы ячейки, при этом они перерисовываются с другой стороны - это эффект pbc. Поэтому используем не первую, а вторую команду:

In [ ]:
trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_1.pdb -skip 20 
trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_1.pdb -skip 20 -pbc mol

Бислоя нет, мицелла. Так выглядит последнее состояние (101-е) системы:

In [10]:
Image(filename='b_pbc_1.png')
Out[10]:

Определим размер, который занимает один липид. Если открыть PDB файл в PyMOL, то изначально ось X слева направо, Y снизу вверх, Z на нас. Я считаю, что нормаль в моем случае - это ось Y. Тогда площадь - это X*Z/2. Возьмем X и Z из box_1.xvg. Размер box меняется. PyMOL показывает только ячейку стартового состояния. Надо построить зависимость площади от времени, нормированную на 1 липид. Это приблизительно, но сойдет для нас.

In [ ]:
g_traj -f b_md.xtc -s b_md.tpr -ob box_1.xvg -xvg none
In [13]:
%matplotlib inline

import matplotlib.pyplot as plt
import csv

t = []
s = []
with open('box_1.xvg', 'r') as box:
    box_reader = csv.reader((row for row in box if not row.startswith('#')), delimiter='\t')    
    for row in box_reader:
        lipid_time = row[1]
        t.append(int(lipid_time))
        x = float(row[2])
        z = float(row[4])
        s_one = x*z/32
        s.append(s_one)

plt.plot(t, s, "b-")
plt.title('Surface of one lipid during time') 
plt.xlabel('Time, ps')
plt.ylabel('SAS of one lipid, nm^2')
plt.show()

Действительно, при сворачивании липидного бислоя или мицеллы, площадь поверхности липида, доступная растворителю, стремится к минимуму. Уменьшается в итоге в 2 раза: от ~1,2 до 0,6 nm^2

Определим изменения в гидрофобной и гидрофильной поверхностях. Группы: calculation group, output group = 2 (DPPC)

In [ ]:
g_sas -f b_md.xtc -s b_md.tpr -o sas_b.xvg
In [23]:
%matplotlib inline

t = []
hf_list = []
hb_list = []

with open('sas_b.xvg', 'r') as box_reader:    
    for line in box_reader:
        row = line.split()
        lipid_time = row[0]        
        t.append(int(lipid_time))
        hb = float(row[1])
        hf = float(row[2]) 
        hb_list.append(hb)
        hf_list.append(hf)

plt.plot(t, hf_list, "r-",label="Hydrophylic")
plt.plot(t, hb_list, "b-",label="Hydrophobic")
plt.title('SAS during time') 
plt.xlabel('Time, ps')
plt.ylabel('Surface, nm^2')
plt.legend()
plt.show()

Мы видим, что при образовании бислоя происходит уменьшение как гидрофобной, так и гидрофильной поверхностей. Это приводит к уменьшению энергии системы в водном растворе. Площадь гидрофобной поверхности гораздо сильнее уменьшается. После 250 ps величина доступной растворителю гидрофобной поверхности становится меньше гидрофильной поверхности. За счет уменьшения гидрофобной поверхности происходит сборка липидов в мицеллу.

Традиционной мерой оценки фазового состояния бифильных молекул является мера порядка. Показывает, насколько упорядочены липидные хвосты. Некоторые части более упорядочены, чем другие. Как посчитать? Угол C-C относительно нормали бислоя. В среднем для постоянно изменяющегося угла <3/2 cos a - 1/2> или <cos a - 1> дают 0, и не 0, если упорядоченность. Для разных связей - разная упорядоченность.

Для анализа нам понадобится специальный индекс файл: http://kodomo.cmm.msu.ru/~golovin/bilayer/sn1.ndx

In [ ]:
g_order -s b_md -f b_md.xtc -o ord_end.xvg -n sn1.ndx -b 45000 -d Y   #для конца траектории
g_order -s b_md -f b_md.xtc -o ord_start.xvg -n sn1.ndx -e 5000 -d Y  #для начала траектории
In [37]:
%matplotlib inline

with open('ord_start.xvg','r') as order:
    mol_num = []
    s_start = []
    for line in order:
        if line.startswith(' '):
            line = line.strip(None).split()
            mol_num.append(int(line[0]))
            s_start.append(float(line[1]))
                   
                       
with open('ord_end.xvg','r') as order:                      
    s_end = []
    for line in order:
        if line.startswith(' '):
            line = line.strip(None).split()
            s_end.append(float(line[1]))
        
plt.plot(mol_num, s_start, "ro",label="start of modelling")
plt.plot(mol_num, s_end, "bo",label="end of modelling")
plt.title('Order of C-C') 
plt.xlabel('Atom number')
plt.ylabel('Order')
plt.legend(loc='center left')
plt.show()

В начале траектории значения меры порядка колеблются, что связано с хаотичностью движений молекул липидов.

В конце траектории видно, что головки липидов более ограничены в движении, чем хвосты. Когда мицелла приобрела четко сформированную структуру, подвижность молекул липидов падает, по сравнению с началом траектории.