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

In [ ]:
# Любой анализ начинают с визуального анализа движений молекул. При вопросе о выводе групп выберите DPPC
trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_1.pdb -skip 20

Откроем полученный b_pbc_1.pdb в Pymol

In [3]:
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]

import pymol
from pymol import cmd,stored
import IPython
pymol.finish_launching()
In [4]:
cmd.bg_color("white")
cmd.load("b_pbc_1.pdb")

Посмотрим, что получилось

In [5]:
import io
import base64
from IPython.display import HTML
In [6]:
import pymol
from pymol import cmd,stored
import IPython
pymol.finish_launching()
In [7]:
cmd.set("ray_trace_mode", 0)
cmd.set("ray_trace_color", "black")
cmd.mpng('prac8')
In [8]:
video = io.open('movie8.mp4', 'r+b').read()
encoded = base64.b64encode(video)
HTML(data='''<video alt="test" controls>
                <source src="data:video/mp4;base64,{0}" type="video/mp4" />
             </video>'''.format(encoded.decode('ascii')))
Out[8]:

видео хаотично и непонятно. попробуем применить следующую команду:

In [ ]:
 trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_3.pdb -skip 20 -pbc mol
In [9]:
cmd.bg_color("white")
cmd.load("b_pbc_3.pdb")
In [10]:
cmd.set("ray_trace_mode", 0)
cmd.set("ray_trace_color", "black")
cmd.mpng('prac13')
In [11]:
import io
import base64
from IPython.display import HTML
In [12]:
video = io.open('movie123.mp4', 'r+b').read()
encoded = base64.b64encode(video)
HTML(data='''<video alt="test" controls>
                <source src="data:video/mp4;base64,{0}" type="video/mp4" />
             </video>'''.format(encoded.decode('ascii')))
Out[12]:

на видео видно, что что-то добавилось, но, в целом, видно плохо. покажем картинку.

In [13]:
import pymol
from pymol import cmd,stored
import IPython
pymol.finish_launching()
In [14]:
#cmd.reset()
#сделаем картинку (последний фрейм)
IPython.display.Image("prac130051.png", retina = True)
Out[14]:

Кажется, что-то со сборкой липидного бислоя получилось, но не до конца. Возможно, что-то пошло не так на этапе подсчета, видны артефакты. Выбрали DPPC. Выберем System - сделаем структуру с молекулами воды. И рассмотрим последнюю - 51ую модель. TITLE bilayer in water t= 25000.00000

In [15]:
cmd.bg_color("white")
cmd.load("b_pbc_4.pdb")
cmd.do('''
show spheres, solvent
''')
In [16]:
cmd.set("ray_trace_mode", 0)
cmd.set("ray_trace_color", "black")
cmd.mpng('prac14')
In [17]:
import pymol
from pymol import cmd,stored
import IPython
pymol.finish_launching()
In [18]:
IPython.display.Image("prac140051.png", retina = True)
Out[18]:
In [19]:
IPython.display.Image("prac140001.png", retina = True) #покажем самый первый фрейм
Out[19]:

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

In [20]:
import pymol
from pymol import cmd,stored
import IPython
pymol.finish_launching()
In [32]:
#видео,чтобы это понять
video = io.open('movie124.mp4', 'r+b').read()
encoded = base64.b64encode(video)
HTML(data='''<video alt="test" controls>
                <source src="data:video/mp4;base64,{0}" type="video/mp4" />
             </video>'''.format(encoded.decode('ascii')))
Out[32]:

видно, что вода на видео отступает

In [33]:
IPython.display.Image("prac140010.png", retina = True)
Out[33]:
In [34]:
IPython.display.Image("prac140011.png", retina = True)
Out[34]:

на 10 и 11 фреймах еще не совсем отступила. а вот между 22 и 23 фреймами - да

In [41]:
IPython.display.Image("prac140022.png", retina = True)
Out[41]:
In [42]:
IPython.display.Image("prac140023.png", retina = True)
Out[42]:

TITLE bilayer in water t= 11000.00000

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

In [ ]:
#размеры ячейки из траектории
g_traj -f b_md.xtc -s b_md.tpr -ob box_1.xvg
In [43]:
time1 = []
square = []

with open("box_1.xvg", 'r') as f1:
    lines = f1.readlines()
for line in lines:
    if (not line.startswith("@")) and (not line.startswith("#")):
        linelist = line.split("\t")
        time1.append(linelist[1])
        square.append(float(linelist[3])*float(linelist[4])/32)
In [44]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
In [45]:
plt.plot(time1, square, c='red',alpha=1)
plt.xlim()
plt.show()

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

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

In [ ]:
g_sas -f b_md.xtc -s b_md.tpr -o sas_b.xvg
In [54]:
import pandas as pd
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("#"))]
time = [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]
dt = pd.DataFrame({'time':time, 'Hydrophobic':hphob, 'Hydrophilic':hphil})
In [55]:
dt.plot(x="time", y = ["Hydrophobic", "Hydrophilic"])
plt.show()

В начале моделирования гидрофобная поверхность максимальна, весь липид доступен для воды. Но это состояние не выгодно в присутствии растворителя, и гидрофобная поверхность резко снижается - образуется мицеллоподобная структура. Выход на плато - образование бислоя (собственно, на 10 и 11 фрейме заметно довольно р гидрофобной поверхности).

Оценим меру порядка (традиционную меру оценки фазового состояния)

In [ ]:
#был скачан индекс-файл sn1.ndx
#был запущен анализ
#нормалью к поверхности бислоя является ось X
g_order -s b_md -f b_md.xtc -o ord_start.xvg -n sn1.ndx -e 5000 -d X  #для начала траектории
g_order -s b_md -f b_md.xtc -o ord_end.xvg -n sn1.ndx -b 45000 -d X #для конца траектории

Однако с параметром -b 45000 программа ломается. Судя по графикам, -b нужно поставить 25000 Получили файлы ord_start.xvg, ord_end.xvg

In [64]:
with open("ord_start.xvg", 'r') as f:
    lines = [i.strip('\n') for i in f.readlines() 
             if (not i.startswith("@")) and (not i.startswith("#"))]
x = [int(i.split()[0]) for i in lines]
start = [float(i.split()[2]) for i in lines]

with open("ord_end.xvg", 'r') as f:
    lines = [i.strip('\n') for i in f.readlines() 
             if (not i.startswith("@")) and (not i.startswith("#"))]
end = [float(i.split()[2]) for i in lines]
In [65]:
data = pd.DataFrame({'x':x, 'Start_order':start, 'End_order':end})
data.plot(x="x", y = ["Start_order", "End_order"])
plt.show()

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

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

  • алгоритм минимизации:steep
Утряска растворителя:
  • число шагов - 1000
  • длина шага - 0.0002 ps
  • алгоритм расчёта электростатики: pme
  • алгоритм термостата - Tcoupl = v-rescale
  • алгоритм баростата - Pcoupl = no;Berendsen
Основной расчет МД:
  • время моделирования - 6ч 14мин
  • количество процессоров - 16
  • эффективность масштабирования - 99,8%
  • алгоритм интегратора - md
  • число шагов - 10000000
  • длина шага - 0.005 ps
  • алгоритм термостата - Tcoupl = v-rescale
  • алгоритм баростата - Pcoupl = Berendsen
  • алгоритм расчета Ван-дер-Ваальсовых взаимодействий - Cut-off
  • алгоритм расчёта электростатики: pme