Моделирование самосборки липидного бислоя

In [1]:
import nglview
import IPython.display
import ipywidgets
from IPython.display import display,display_svg,SVG,Image, HTML
import __main__
from time import sleep
import pymol
import matplotlib.pyplot as plt
%matplotlib inline
In [1]:
%%bash

wget http://kodomo.cmm.msu.ru/~golovin/bilayer/dppc.itp
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/lipid.itp
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/dppc.gro
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/b.top  # отредактирован на 64 липида
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/em.mdp
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/pr.mdp
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/md.mdp
In [2]:
%%bash

# на основе липида DPPC созадим ячейку с 64 липидами
genconf -f dppc.gro -o b_64.gro -nbox 4 4 4
In [27]:
%%bash

# сохраним pdb-файлы
editconf -f dppc.gro -o dppc.pdb 
editconf -f b_64.gro -o b_64.pdb
In [1]:
# Launch PyMOL 
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]

import pymol
from pymol import cmd
pymol.finish_launching()

# Set image
import time
from IPython.display import Image
defaultImage = 'pymolimg.png'
def prepareImage(width=500, height=500, sleep=2, filename=defaultImage):
    cmd.ray(width, height)
    cmd.png(filename)
    time.sleep(sleep)

Смотрим на структуры одного DPPC и ячейки из 64 DPPC

In [34]:
cmd.delete('all')
cmd.load('dppc.pdb', 'dppc')
cmd.rotate('x', 45, 'dppc')
cmd.rotate('y', 70, 'dppc')
cmd.zoom('dppc',4)
prepareImage()
Image(defaultImage)
Out[34]:
In [36]:
cmd.delete('all')
cmd.load('b_64.pdb','b_64')
prepareImage()
Image(defaultImage)
Out[36]:
In [4]:
%%bash

# сделаем отступ в ячейке от липидов, что бы добавить примерно 2500 молекул воды
editconf -f b_64.gro -o b_ec -d 0.5
In [5]:
%%bash

# оптимизируем геометрию системы, что бы удалить "плохие" контакты молекул
grompp -f em -c b_ec -p b -o b_em -maxwarn 2
mdrun -deffnm b_em -v

Произошло уменьшение максимальной силы в ходе оптимизации:

Fmax (start) = 4.37970e+05

Fmax (end) = 6.16887e+02

In [6]:
%%bash

# добавим в ячейку молекулы воды типа spc
genbox -cp b_em -p b -cs spc216 -o b_s
In [7]:
%%bash

# проведём "утряску" воды
grompp -f pr -c b_s -p b -o b_pr -maxwarn 1
mdrun -deffnm b_pr -v
In [8]:
%%bash

# сохраним pdb
editconf -f b_pr.gro -o b_pr.pdb
editconf -f b_s.gro -o b_s.pdb
In [28]:
# Submitted batch job 1651987

Посмотрим на структуры 64-DPPC-ячейки после после добавления воды и после утряски

In [38]:
# после добавления воды
cmd.delete('all')
cmd.load('b_s.pdb', 'b_s')
prepareImage()
Image(defaultImage)
Out[38]:

Видно, что по краям ячейки ячейка заполнена молекулами воды, липид в центре.

In [39]:
# после утряски 
cmd.delete('all')
cmd.load('b_pr.pdb', 'b_pr')
prepareImage()
Image(defaultImage)
Out[39]:

После утряски молекулы воды расположены более упорядочено и плотно (хотя на глаз разница невелика).

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

Визуализируем движение молекул

In [40]:
%%bash

trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_1.pdb -skip 20 -pbc mol 
DPPC
In [4]:
# готовим .png из состояний липида для анимации
cmd.delete('all')
cmd.load('b_pbc_1.pdb', 'b_pbc_1')
cmd.mset('1-'+str(cmd.count_states('all')))
cmd.mpng('b_pbc_1/b_pbc_1')
In [ ]:
# %%bash
# ffmpeg -pattern_type glob -i "*.png" -c:v libx264 -pix_fmt yuv420p -vf "scale=trunc(iw/2)*2:trunc(ih/2)*2" ~/b_pbc_1.mp4
In [7]:
from IPython.display import HTML

def show_video(file):
    video = open(file, "rb").read().encode("base64")
    video_tag = '<video controls src="data:video/mp4;base64, %s">' % video
    display(HTML(data=video_tag))
In [8]:
show_video('b_pbc_1.mp4')

Упорядоченный бислой появляется около 65-го кадра (из pdb-файла: соответствует 30 сек) и далее сохраняет структуру:

In [10]:
Image('b_pbc_1/b_pbc_10065.png')
Out[10]:

Получить размеры ячейки (из траектории) и определим площадь, занимаемую одним липидом.

In [ ]:
%%bash
g_traj -f b_md.xtc -s b_md.tpr -ob box_1.xvg
System
In [4]:
import numpy as np
data = np.loadtxt('box_1.xvg', comments=['#', '@'])[:,0:4]
In [35]:
np.std(data[:,1:4], axis=0)
Out[35]:
array([ 0.88146623,  0.62561565,  0.8488636 ])

Стандартное отклонение наименьшее по оси Y, скорее всего, она является нормалью к поверхности бислоя.

Построим зависимость площади по осям X и Z от времени, нормируем на один липид в слое.

In [5]:
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
%matplotlib inline

time = data[:,0]
x_axis = data[:,1]
z_axis = data[:,3]

plt.plot(time, x_axis*z_axis/64, color='salmon', linewidth=2.)
plt.title('XoZ square')
plt.xlabel('Time')
plt.show()

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

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

In [ ]:
%%bash

g_sas -f b_md.xtc -s b_md.tpr -o sas_b.xvg -dt 100
In [11]:
data = np.loadtxt('sas_b.xvg', comments=['#', '@'])
time = data[:,0]
hydro_phob = data[:,1]
hydro_phil = data[:,2]

plt.plot(time, hydro_phob, color='orange')
plt.plot(time, hydro_phil, color='blue')
plt.title('Hydrophobic and hydrophilic surfaces')
plt.xlabel('Time')
phob_lab = mpatches.Patch(color='orange', label='Hydrophobic surface')
phil_lab = mpatches.Patch(color='blue', label='Hydrophilic surface')
plt.legend(handles=[phob_lab,phil_lab])
plt.show()

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

Посмотрим на изменение параметра порядка для оценки фазового состояния бислоя:

In [ ]:
%%bash

wget http://kodomo.cmm.msu.ru/~golovin/bilayer/sn1.ndx # индекс-файл для g_order
In [ ]:
%%bash

# начало траектории
g_order -s b_md -f b_md.xtc -o ord_start.xvg -n sn1.ndx -e 5000 -d Y
# конец траектории
g_order -s b_md -f b_md.xtc -o ord_end.xvg -n sn1.ndx -b 45000 -d Y
In [10]:
data_start = np.loadtxt('ord_start.xvg', comments=['#', '@'])
data_end = np.loadtxt('ord_end.xvg', comments=['#', '@'])
In [25]:
data_start = np.loadtxt('ord_start.xvg', comments=['#', '@'])
time_start = data_start[:,0]
order_start = data_start[:,3]
data_end = np.loadtxt('ord_end.xvg', comments=['#', '@'])
time_end = data_end[:,0]
order_end = data_end[:,3]

plt.plot(time, order_start, color='black')
plt.plot(time, order_end, color='red')
plt.title('Order')
plt.xlabel('Time')
start_lab = mpatches.Patch(color='black', label='Start of trajectory')
end_lab = mpatches.Patch(color='red', label='End of trajectory')
plt.legend(handles=[end_lab, start_lab])
plt.show()

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