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

In [29]:
import io
import base64
from IPython.display import HTML
from IPython.display import Image,display
import gromacs.formats
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

Начнем анализ с визуального анализа движений молекул при помощи команды trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_1.pdb -skip 20. При вопросе о выводк групп выберем DPPC.

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

In [11]:
video = io.open('b_pbc_1.mp4', 'r+b').read()
encoded = base64.b64encode(video)
display(HTML(data='''<video alt="test" controls>
                <source src="data:video/mp4;base64,{0}" type="video/mp4" />
             </video>'''.format(encoded.decode('ascii'))))

Выглядит хаотично, поэтому попробуем применить команду trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_1.pdb -skip 20 -pbc mol

In [12]:
video = io.open('b_pbc_1_new.mp4', 'r+b').read()
encoded = base64.b64encode(video)
display(HTML(data='''<video alt="test" controls>
                <source src="data:video/mp4;base64,{0}" type="video/mp4" />
             </video>'''.format(encoded.decode('ascii'))))

Стало намного лучше. Видно, что произошла самосборка билипидного бислоя.

Образование бислоя происходит в момент времени между 23 и 24 моделью (t= 11000.00000).

Силовое поле, используемое при построении топологии - ffgmx

Заряд системы - 0

Размер и форма ячейки - 6.26000 x 4.44300 x 5.77800 нм

Минимизация энергии (em.mdp):

  • алгоритм минимизации: integrator = l-bfgs
  • алгоритм расчёта электростатики: coulombtype = Cut-off
  • алгоритм расчёта Ван-дер-Ваальсовых взаимодействий: vdw-type = Cut-off
  • Модель, которой описывался растворитель - spc.itp

Утряска растворителя:

  • Число шагов: 10000
  • Длина шага: 0.001 пс
  • Алгоритм расчёта электростатики: coulombtype = pme
  • Алгоритм расчета Ван-дер-Ваальсовых взаимодействий: vdw-type = Cut-off
  • Алгоритм термостата: Tcoupl = Berendsen

Алгоритм баростата: Pcoupl = no

Основной расчёт МД (md.mdp):

  • Количество процессоров - 16
  • Длина траектории - 50000.0 ps
  • Число шагов - 10000000
  • Длина шага - 0.005
  • Алгоритм интегратора - integrator = md
  • Алгоритм расчёта электростатики - coulombtype = pme
  • Алгоритм расчёта Ван-дер-Ваальсовых взаимодействий - vdw-type = Cut-off
  • Алгоритм термостата - Tcoupl = v-rescale
  • Алгоритм баростата - Pcoupl = Berendsen

При помощи команды g_traj -f b_md.xtc -s b_md.tpr -ob box_1.xvg получим box_1.xvg файл, в котором содержатся 4 столбца. Первый стоблик - время, остальные - размер ячейки по оси. Ось Y является нормалью. Перемножим значения по осям, не являющимися нормалью (X,Z), а затем разделим на среднее количество молекул липидов в одном слое (32).

In [37]:
area, time = [], []

with open("../sem_11/out_lom/box_1.xvg") as f:
    for line in f:
        if '#' not in line and '@' not in line:
            cols = line.split()
            time.append(float(cols[0]))
            x = float(cols[1])
            z = float(cols[3])
            area.append(x * z / 32)

fig = plt.figure()
ax1 = fig.add_subplot(111)
# ax1.set_title("Solvent Accessible Surface")    
ax1.set_xlabel('Time (ps)')
ax1.set_ylabel('Area (nm\S2\N)')
ax1.plot(time, area, c='g')
leg = ax1.legend()
plt.show()

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

При помощи команды g_sas -f b_md.xtc -s b_md.tpr -o sas_b.xvg определим изменение гидрофобной и гидрофильной поверхности в ходе самосборки.

In [30]:
x, y, time = [], [], []

with open("../sem_11/out_lom/sas_b.xvg") as f:
    for line in f:
        if '#' not in line and '@' not in line:
            cols = line.split()
            time.append(float(cols[0]))
            x.append(float(cols[1]))
            y.append(float(cols[2]))


fig = plt.figure()
ax1 = fig.add_subplot(111)
ax1.set_title("Solvent Accessible Surface")    
ax1.set_xlabel('Time (ps)')
ax1.set_ylabel('Area (nm\S2\N)')
ax1.plot(time, x, c='g', label='Hydrophobic')
ax1.plot(time, y, c='b', label='Hydrophilic')
leg = ax1.legend()
plt.show()

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

Построим зависимость изменения гидрофобной гидрофильной поверхностей доступных растворителю от времени. Сделаем это при помощи команды 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 [49]:
end, start, time = [], [], []


with open("../sem_11/out_lom/ord_end.xvg") as f:
    for line in f:
        if '#' not in line and '@' not in line:
            cols = line.split()
            end.append(float(cols[2]))
            
with open("../sem_11/out_lom/ord_start.xvg") as f:
    for line in f:
        if '#' not in line and '@' not in line:
            cols = line.split()
            time.append(float(cols[0]))
            start.append(float(cols[2]))


fig = plt.figure()
ax1 = fig.add_subplot(111)
ax1.set_title("Order tensor diagonal elements")    
ax1.set_xlabel('Atom')
ax1.set_ylabel('S')
ax1.plot(time, end, c='g', label='end')
ax1.plot(time, start, c='b', label='start')
leg = ax1.legend()
plt.show()

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