Практикум 11¶
import __main__
from IPython.display import Image
__main__.pymol_argv = [ 'pymol', '-x' ] # открыть pymol без использования внешнего GUI
import pymol
pymol.finish_launching()
from pymol import cmd,stored
import numpy as np
import time
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
Продемонстрирую последовательность команд, которую я выполняла чтобы получить необходимые для молекулярной динамики файлы.
Перейдем в директорию, в которой будем выполнять команды и загрузим пути к необходимым программам
%%bash
cd ../Perevoshchikova/
source /home/preps/golovin/progs/bin/GMXRC.bash
Создаем ячейку с 64 липидами
%%bash
gmx genconf -f dppc.gro -o b_64.gro -nbox 4 4 4
Создаем pdb файлы из dppc.gro и b_64.gro
%%bash
editconf -f ../Perevoshchikova/dppc.gro -o ../Perevoshchikova/dppc.pdb
editconf -f ../Perevoshchikova/b_64.gro -o ../Perevoshchikova/b_64.pdb
Теперь посмотрим на то, что у нас получилось
cmd.delete("all")
cmd.bg_color('white')
cmd.do('''load dppc.pdb
as cartoon, all
set stick_radius, 0.5
show stick, all
orient dppc''')
cmd.set("opaque_background", 'off')
cmd.set("ray_trace_mode", '4')
cmd.ray(600,400)
cmd.png('./task11_1.png')
time.sleep(2)
Image(filename='./task11_1.png')
Видим, что у меня истекла лицензия, а файл dppc содержит структуру фосфатидилхолина у которого одна из жирных кислот содержит циклопентановый элемент.
cmd.delete("all")
cmd.bg_color('white')
cmd.do('''load b_64.pdb
as cartoon, all
set stick_radius, 0.5
show stick, all
zoom b_64, complete = 1
turn x, 20
turn y, 20
turn z, 20''')
cmd.set("opaque_background", 'off')
cmd.set("ray_trace_mode", '0')
cmd.ray(600,400)
cmd.png('./task11_2.png')
time.sleep(2)
Image(filename='./task11_2.png', height = 7000, width = 7000)
Видим, что команда gmx genconf -f dppc.gro -o b_64.gro -nbox 4 4 4 создает ячейку 444, заполненную нашими молекулами липида
Необходимо отредактировать несколько файлов, прежде чем двигаться дальше.
(1) в файле b.top установить количество липидов всистеме 64
(2) в файле с параметрами для молекулярной динамики установим ша по времени 0.004 ps
Сделаем небольшой отступ в ячейке от липидов, чтобы осталось место воде
%%bash
gmx editconf -f b_64.gro -o b_ec -d 0.5
Проведем оптимизацию геометрии системы из 64 молекул липида, чтобы удалить плохие контакты воды
%%bash gmx grompp -f em -c b_ec -p b -o b_em -maxwarn 2
gmx mdrun -deffnm b_em -v
В ходе оптимизации согласно отчету GROMACS, записанному в файле b_em.log потенциальная энергия в ходе 68 шагов уменьшилась с 4.63573e+05 Кдж/моль до -3.11687e+04 Кдж/моль. Правда GROMACS отмечает, что программа "did not reach the requested Fmax < 1". Ну в прочем и ладно, эта наша система еще не окончательная.
Добавим в ячейку молекулы воды
%%bash
gmx solvate -cp b_em -p b -cs spc216 -o b_s
Произведем оптимизацию геометрии молекул воды. Часто такая оптимизация из за изначально некорректного положения молекум может привести к появлению больших, действующих на молекулы сил и водяному взрыву.
%%bash
gmx grompp -f pr -c b_s -p b -o b_pr -maxwarn 1
gmx mdrun -deffnm b_pr -v
#если водяной взрыв, то
gmx grompp -f em -c b_s -p b -o b_empr -maxwarn 1
gmx mdrun -deffnm b_empr -v
gmx grompp -f pr -c b_empr -p b -o b_pr -maxwarn 1
gmx mdrun -deffnm b_pr -v
Теперь переведем файл с добавленной водой b_s.gro и файл с оптимизированной водой b_pr.gro в pdb формат и посмотрим, как это выглядит
%%bash
editconf -f ../Perevoshchikova/b_s.gro -o ../Perevoshchikova/b_s.pdb
editconf -f ../Perevoshchikova/b_pr.gro -o ../Perevoshchikova/b_pr.pdb
Посмотрим как это выглядит в pymol
До оптимизации геометрии¶
cmd.do("reinitialize")
cmd.delete("all")
cmd.set("ray_trace_mode", '0')
cmd.bg_color('white')
cmd.do('''load b_s.pdb
as cartoon, all
set stick_radius, 0.5
show stick, all
show nonbonded
zoom b_s, complete = 1
turn x, 20
turn y, 20
turn z, 20''')
cmd.ray(1800,1800)
cmd.png('./task11_3.png')
time.sleep(2)
Image(filename='./task11_3.png',
width = 7000,
height = 7000)
cmd.do("reinitialize")
cmd.delete("all")
cmd.set("ray_trace_mode", '0')
cmd.bg_color('white')
cmd.do('''load b_s.pdb
as cartoon, all
set stick_radius, 0.5
show stick, all
select resn SOL
hide everything, sele
show nonbonded
zoom b_s, complete = 1
''')
cmd.ray(1800,1800)
cmd.png('./task11_3.png')
time.sleep(2)
Image(filename='./task11_3.png',
width = 7000,
height = 7000)
После оптимизации геометрии¶
cmd.do("reinitialize")
cmd.delete("all")
cmd.set("ray_trace_mode", '0')
cmd.bg_color('white')
cmd.do('''load b_pr.pdb
as cartoon, all
set stick_radius, 0.5
show stick, all
show nonbonded
zoom b_pr, complete = 1
turn x, 20
turn y, 20
turn z, 20''')
cmd.ray(1800,1800)
cmd.png('./task11_4.png')
time.sleep(2)
Image(filename='./task11_4.png',
width = 8000,
height = 7000)
cmd.do("reinitialize")
cmd.delete("all")
cmd.set("ray_trace_mode", '0')
cmd.bg_color('white')
cmd.do('''load b_pr.pdb
as cartoon, all
set stick_radius, 0.5
show stick, all
select resn SOL
hide everything, sele
show nonbonded
zoom b_pr, complete = 1
''')
cmd.ray(1800,1800)
cmd.png('./task11_5.png')
time.sleep(2)
Image(filename='./task11_5.png',
width = 7000,
height = 7000)
На самом деле видно, что относительно просто добавленной воды - файл b_s молекулы воды поменяли свое положение после оптимизации. На мой взгляд заметно, что в файле, где молекулы воды добавлены случайно (b_s) заметно ноличие странных незаполненных полостей - областей пространства не содержащих молекул. В то время как в структуре b_pr видимых полостей нет, молекулы воды расположены равномерно. К тому же на мой взгляд заметно, что в файле до оптимизации молекулы воды расположены случайно и найти водородную связь довольно сложно в то время как в файле после оптимизации молекулы воды сориентированы оптимально для образования водородной связи. Кроме тоо очевидно, что после оптимизации геометрии хвосты липидов стали лежать более свободно.
Далее заходим на суперкомпьютер и запускаем моделирование¶
%%bash
source /home/preps/golovin/progs/bin/GMXRC.bash
Сначала посмотрим как выглядит pdb файл с результатом без наложения граничных условий
%%bash
gmx trjconv -f ../Perevoshchikova/new/b_md.xtc -s ../Perevoshchikova/new/b_md.tpr -o b_pbc_1.pdb -skip 20
cmd.do("reinitialize")
cmd.delete("all")
cmd.set("ray_trace_mode", '0')
cmd.bg_color('white')
cmd.do('''load ./mmresults/b_pbc_1.pdb
as cartoon, all
set stick_radius, 0.5
show stick, all
''')
#покажем 50е состояние из 101
cmd.do('''zoom p_bc_1, complete = 1''')
cmd.ray(900,900)
cmd.png('./task11_7.png')
time.sleep(2)
Image(filename='./task11_7.png',
width = 7000,
height = 7000)
Видим странные связи, это все потому что мы не выставили pbc
%%bash
gmx trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_2.pdb -skip 20 -pbc mol
cmd.do("reinitialize")
cmd.delete("all")
cmd.set("ray_trace_mode", '0')
cmd.bg_color('white')
cmd.do('''load ./mmresults/b_pbc_2.pdb
as cartoon, all
set stick_radius, 0.5
show stick, all
''')
#покажем 50е состояние из 101
cmd.do('''zoom p_bc_2, complete = 1''')
cmd.ray(900,900)
cmd.png('./task11_8.png')
time.sleep(2)
Image(filename='./task11_8.png',
width = 7000,
height = 7000)
Добавим к хорошей структуре с граничными условиями молекулы воды
cmd.do('''load ./mmresults/b_pbc_2_sol.pdb
show stick, all
''')
Поищем момент образования бислоя. В принципе нечто довольно организованное, похожее на мицеллу образуется на шаге 15
#посмотрим на модель 14
cmd.ray(900,900)
cmd.png('./task11_9.png')
time.sleep(2)
Image(filename='./task11_9.png',
width = 7000,
height = 7000)
Хоть парочка молекул липидов торчит тут в другую сторону,я считаю что в целом формирование мицеллы произошло. Найдем время которое прошло до сборки мицеллы на шаге 14.
TITLE bilayer in water t= 0.00000 step= 0
Оценим площадь, приходящуюся на одну молекулу липида¶
#посмотрим на размер ячейки в состоянии 101
cmd.ray(900,900)
cmd.png('./task11_10.png')
time.sleep(2)
Image(filename='./task11_10.png',
width = 7000,
height = 7000)
Сгенерировали файл размера ячейки по трем измерениям.
gmx traj -f b_md.xtc -s b_md.tpr -ob box_1.xvg
Размер ячейки в конце 11,9 по Х 8,44 по У и 2,23 по Z, значит нормаль к поверхности бислоя - X, а Y и Z это те величины, которые нужно перемножить, чтобы оценить площадь бислоя
tim = []
square = []
with open("./mmresults/box_1.xvg", 'r') as infile:
lines = infile.readlines()
for line in lines:
if (not line.startswith("@")) and (not line.startswith("#")):
linelist = line.split("\t")
tim.append(linelist[1])
square.append(float(linelist[3])*float(linelist[4])/64)
square = [float(el) for el in square]
tim = [float(el) for el in tim]
#f, ax = plt.subplots(1,1,figsize=(7,7))
sns.lineplot(x = tim, y = square)
Мы видим, что площадь бислоя уменьшается со временем моделирования, что соответствует организации липидов в бислой и их уплотнению.
Определение гидрофобной и гидрофильной поверхности в ходе самосборки¶
К сожалению не смогла выполнить задание, потому что на ломоносове этой программы нет, а на кодомо она более старой версии. Но в целом понятно, что гидрофобная поверхность уменьшается в процессе сборки бислоя, движущая сила сборки бислоя - уменьшение площади контакта растворителя с липидом, потому что вода не может с липидом взаимодействовать.
Мера порядка¶
Оценим меру порядка на атом для углеродных хвостов в конце траектории командой :
gmx order -s b_md -f b_md.xtc -o ord_end.xvg -n sn1.ndx -b 35000 -d x
И в начале траектории:
gmx order -s b_md -f b_md.xtc -o ord_start.xvg -n sn1.ndx -b 5000 -d x
num = []
start1 = []
start2 = []
start3 = []
with open("./mmresults/ord_start.xvg", 'r') as f3:
lines = f3.readlines()
for line in lines:
if (not line.startswith("@")) and (not line.startswith("#")):
linelist = line.split()
num.append(linelist[0])
start1.append(linelist[1])
start2.append(linelist[2])
start3.append(linelist[3])
num_end = []
end1 = []
end2 = []
end3 = []
with open("./mmresults/ord_end.xvg", 'r') as f4:
lines = f4.readlines()
for line in lines:
if (not line.startswith("@")) and (not line.startswith("#")):
linelist = line.split()
num_end.append(linelist[0])
end1.append(linelist[1])
end2.append(linelist[2])
end3.append(linelist[3])
f, axs = plt.subplots(3,1, figsize=(10,30))
axs[0].plot(num, start1, label = "start order" )
axs[0].plot(num, end1, label = "end order" )
axs[0].legend()
axs[1].plot(num, start2, label = "start order" )
axs[1].plot(num, end2, label = "end order" )
axs[1].legend()
axs[2].plot(num, start3, label = "start order" )
axs[2].plot(num, end3, label = "end order" )
axs[2].legend()
gmx order оценивает три разных меры порядка, но вне зависимости от способа оценивания можно сделать 2 вывода.
(1) Мера порядка в конце моделирования для всех атомов больше, чем в начале моделирования, что логично, потому что возникает липидный бислой, предполагающий упорядочивание молекул.
(2)Мера порядка растет по мере увеличения индекса атома углерода в насыщенной цепи. То есть в середине бислоя упорядоченность выше, чем на ранице с растворителем.
Параметры запуска МД¶
- Силовое поле используемое при построении топологии топологии. fgmx
- Заряд системы. Причины этого значения. 0. Плюсы на холинах компенсируют свободные минусы фосфатов.
- Размер и форму ячейки. В конце моделирования - 11.89508 нм 8.44247 нм 2.23471 нм прямоугольный параллелепипед В начале моделирования 6.26000 нм 4.44300 нм 5.77800 нм прямоугольный параллелепипед
- Минимизация энергии:
em.mdp
- Алгоритм минимизации энергии.
steepest descent - Алгоритм расчёта электростатики и Ван-дер-Ваальсовых взаимодействий.
cutoff
- Алгоритм минимизации энергии.
- Модель, которой описывался растворитель
SPC - Утряска растворителя:
- Для биополимеров, укажите параметр который обуславливает неподвижность биополимера.
constraints = all-bonds. ограничение длин связей, но не углов - Число шагов.
1000 - Длина шага.
0.0002 ps - Алгоритм расчёта электростатики и Ван-дер-Ваальсовых взаимодействий.
PME, cutoff - Алгоритмы термостата и баростата.
v-rescale и no
- Для биополимеров, укажите параметр который обуславливает неподвижность биополимера.
- Основной расчёт МД:
- Время моделирования
16837.799 c - количество процессоров,
22 - Длину траектории
40 нс - Число шагов.
10000000 - Длина шага.
0.004 ps - Алгоритм интегратора.
md = leap-frog algorithm - Алгоритм расчёта электростатики и Ван-дер-Ваальсовых взаимодействий.
PME, cutoff - Алгоритмы термостата и баростата.
v-rescale и Berendsen
- Время моделирования