HW9_Didkovskaya

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

In [16]:
import nglview
In [10]:
import IPython.display
import ipywidgets
from IPython.display import display,display_svg,SVG,Image, HTML
In [3]:
import __main__
from time import sleep
import pymol
In [2]:
import numpy as np
In [3]:
import matplotlib.pyplot as plt
%matplotlib inline
Скачиваем первичные файлы wget
In [43]:
%%bash
genconf -f dppc.gro -o b_64.gro -nbox 4 4 4
In [44]:
%%bash
editconf -f dppc.gro -o dppc.pdb
In [45]:
%%bash
editconf -f b_64.gro -o b_64.pdb
In [5]:
%%bash
ls

b_64.gro
b_64.pdb
b.top
dppc.gro
dppc.itp
dppc.pdb
em.mdp
HW9_Didkovskaya.ipynb
lipid.itp
md.mdp
pr.mdp
In [15]:
__main__.pymol_argv = ['pymol','-qc']
pymol.finish_launching()
pymol.cmd.load('b_64.pdb')
pymol.cmd.do('''
    png b_64.png
''')
In [16]:
Image(filename='b_64.png')
Out[16]:
In [18]:
__main__.pymol_argv = ['pymol','-qc']
pymol.finish_launching()
pymol.cmd.load('dppc.pdb')
pymol.cmd.do('''
    png dppc.png
''')
In [19]:
Image(filename='dppc.png')
Out[19]:
В текстовом редакторе в файле b.top установите правильное количество липидов в системе.

С помощью команды nano исправили количество липидов на 64

Сделаем небольшой отступ в ячейке от липидов, что бы добавить примерно 2500 молекул воды.
In [20]:
%%bash
editconf -f b_64.gro -o b_ec -d 0.5
Проведём оптимизацию геометрии системы, что бы удалить "плохие" контакты молекул.
In [21]:
%%bash
grompp -f em -c b_ec -p b -o b_em -maxwarn 2 
mdrun -deffnm b_em -v
Отметье в отчёте изменение максимальной силы в ходе оптимизации геометрии. Занесите начальное и конечное значение максимальной силы

Начальное значение максимальной силы - 4.37970e+05

Конечно значени максимальной силы - 6.16887e+02

Сила значительно уменьшилась

Добавим в ячейку молекулы воды типа spc.
In [22]:
%%bash
genbox -cp b_em -p b -cs spc216 -o b_s
Проведём "утряску" воды:
In [23]:
%%bash
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 изменеия в системах. Занесите наблюдения в отчёт.
In [24]:
%%bash
editconf -f b_pr.gro -o b_pr.pdb
editconf -f b_s.gro -o b_s.pdb

Визуализация

In [25]:
__main__.pymol_argv = ['pymol','-qc'] 
pymol.finish_launching()
pymol.cmd.load('b_pr.pdb')
pymol.cmd.do('''
    png b_pr.png
''')
In [27]:
Image(filename='b_pr.png')
Out[27]:
In [28]:
__main__.pymol_argv = ['pymol','-qc'] 
pymol.finish_launching()
pymol.cmd.load('b_s.pdb')
pymol.cmd.do('''
    png b_s.png
''')
In [30]:
Image(filename='b_s.png')
Out[30]:

Есть отличие в конформации и в положение воды. b_s более плотный

переходим на Ломоносов

scp -r ./hw9 lom:_scratch/chem/didkovskaya_hw9

ssh lom

cd _scratch/chem/didkovskaya_hw9/hw9

cp /home/users/golovin/progs/share/gromacs/top/residuetypes.dat .

grompp -f md -c b_pr -p b -o b_md -maxwarn 1

sbatch -n 4 -e error.log -o output.log -t 5 -p test impi /opt/ccoe/gromacs-5.0.4/build/bin/gmx_mpi mdrun -testverlet -deffnm b_md -v

Submitted batch job: 1644890

Запускаем основное моделирование на суперкомпьтере.

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

Submitted batch job: 1645769

Копируем получившиеся файлы обратно на сервер:

scp -r lom:_scratch/chem/didkovskaya_hw9/hw9/b_md.* ./

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

Любой анализ начинают с визуального анализа движений молекул. При вопросе о выводе групп выберите DPPC.
In [6]:
%%bash
trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_1.pdb -skip 20 -pbc mol
In [19]:
__main__.pymol_argv = ['pymol','-qc'] 
pymol.finish_launching()
pymol.cmd.load('b_pbc_1.pdb')
pymol.cmd.do('''
    png b_pbc_1.png
''')
In [20]:
Image(filename='b_pbc_1.png')
Out[20]:

Конечное состояние:

In [21]:
Image(filename='end.png')#101 кадр
Out[21]:

В промежуточных кадрах происходит что-то страшное..

In [22]:
Image(filename='67.png')#67 кадр
Out[22]:

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

REMARK GENERATED BY TRJCONV

TITLE bilayer in water t= 33000.00000

CRYST1 78.069 55.409 37.098 90.00 90.00 90.00 P 1 1

MODEL 67

Определите площадь занимаемую одним липидом. Для этого вам придется получить размеры ячейки из траектории
In [32]:
%%bash
g_traj -f b_md.xtc -s b_md.tpr -ob box_1.xvg

С помощью SuperSym определяем, что нормалью будет ось Х

Постройте зависимость площади по соотвествующим осям от времени

Чтобы выудить координаты, отрежем первые 23 строки и запишем в новый файл "box_2.xvg"

In [35]:
box = 'box_2.xvg'
ax=np.loadtxt(box, skiprows=0)
In [50]:
ax[0][1:4]
Out[50]:
array([ 6.26 ,  4.443,  5.778])
In [48]:
time=[]
shape=[]
for i in range(len(ax)):
    time.append(ax[i][0]) # время 
    shape.append(ax[i][1:4])
shape=np.asarray(shape)  # координаты
In [49]:
S_xy=shape[:,0]*shape[:,1]
S_xz=shape[:,0]*shape[:,2]
S_yz=shape[:,1]*shape[:,2]
s_norm=S_yz/32
In [50]:
plt.plot(time, s_norm)
plt.xlabel('time')
plt.ylabel('Square of ZY')
Out[50]:
Text(0,0.5,u'Square of ZY')

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

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

Формируем файл и потом опять удаляем из него ненужные начальные строки. Важно: нулевой стобец-время, первый- гидрофобные, второй - гидрофильный

In [ ]:
%%bash
g_sas -f b_md.xtc -s b_md.tpr -o sas_b.xvg -dt 100
In [5]:
sas = 'sas_b.xvg'
sas1=np.loadtxt(sas, skiprows=0)
In [7]:
time=[]
fob=[]
fil=[]
for i in range(len(sas1)):
    time.append(sas1[i][0]) # время 
    fob.append(sas1[i][1]) # гидрофобный
    fil.append(sas1[i][2]) # гидрофильный
In [23]:
plt.plot(time, fil,'g',label='hydrofil')
plt.plot(time, fob,'r',label='hydrofob')
plt.xlabel('time')
plt.ylabel('gidro')
plt.legend()
Out[23]:
<matplotlib.legend.Legend at 0x7fa4a93b4d50>

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

In [ ]:
%%bash
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/sn1.ndx
Традиционной мерой оценки фазового состояния бифильных молекул является мера порядка. Для анализа нам понадобится специальный индекс файл. Теперь запустите сам анализ, где "P" это ось, которя является нормалью к поверхности бислоя (см. выше):

Для конца траектории:

In [51]:
%%bash
g_order -s b_md -f b_md.xtc -o ord_end.xvg -n sn1.ndx -b 45000 -d X

Для начала траектории:

In [52]:
%%bash
g_order -s b_md -f b_md.xtc -o ord_start.xvg -n sn1.ndx -e 5000 -d X

Снова с помощью функции sed отрезаем шапки файлов.

In [28]:
start = 'ord_start.xvg'
st=np.loadtxt(start, skiprows=0)
end = 'ord_end.xvg'
end=np.loadtxt(end, skiprows=0)
In [33]:
mol = st[:,0]
f_s = st[:,3]
f_e = end[:,3]
In [42]:
plt.plot(mol, f_s,'g',label='start',marker='o')
plt.plot(mol, f_e,'r',label='end',marker='o')
plt.xlabel('DSSP atom number')
plt.ylabel("Measure of order")
plt.legend()
Out[42]:
<matplotlib.legend.Legend at 0x7fa4a90e8b50>

Упорядоченность атомов одного из жирных хвостов липида в конце самосборки бислоя выше примерно в 4 раза, чем в начале самосборки бислоя, причем атомы, расположенные вблизи головки липида, упорядочены сильнее, чем атомы, расположенные на конце.