Молекулярная динамика биологических молекул в GROMACS¶

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

*Силовое поле - gmx.ff, размер ячейки - 6.26000 4.44300 5.77800, растворитель - spc вода. При утряске растворителя: шаг 0.0002, макс число шагов 1000, термостат - velocity rescaling, баростат - нет. При МД шаг - 0.005, макс число шагов - 10000000, интегратор - md, термостат - velocity rescaling, баростат - по Берендсену.

Сначала на основе одного липида мы создали ячейку с 64 липидами, добавили 2500 молекул воды.

Потом провели оптимизацию геометрии системы, чтобы удалить "плохие" контакты молекул. При этом максимальная сила в начале была 4.37969e+05, а в конце 3.59368e+02.

Еще провели "утряску" воды. Можем сравнить визуально, что изменилось.

In [1]:
from IPython.display import Image

Состояние до:

In [2]:
Image('pr8_1.png', width = 800)
Out[2]:

Состояние после:

In [3]:
Image('pr8_2.png', width = 800)
Out[3]:

Видим, что положение воды действительно поменялось.

Далее провели уже собственно моделирование (результат можно найти на kodomo:/tmp/md/renata). Посмотрим на движение молекул.

In [8]:
display(Image('pr8_movie1.gif', format='png'))

Вышло что-то мягко говоря странное.. Причем интересно, что сначала все идет очень даже хорошо, головки повернулись к воду, хвостики друг к другу (уже где-то на 12 модели).

In [9]:
Image('pr8_4.png', width = 800)
Out[9]:
In [10]:
Image('pr8_5.png', width = 800)
Out[10]:

Но потом начало происходить что-то странное.

In [11]:
Image('pr8_6.png', width = 800)
Out[11]:

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

Определим площадь, занимаемую одним липидом. Зависимость нормированной площади (x*y/32) от времени:

In [14]:
Image('pr8_7.PNG', width = 800)
Out[14]:

Теперь посмотрим на зависимость поверхности, доступной растворителю, от времени (с помощью gmx sasa):

In [15]:
Image('pr8_8.PNG', width = 800)
Out[15]:

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

Традиционной мерой оценки фазового состояния бифильных молекул является мера порядка. Поэтому посмотрим и на нее.

В начале:

In [17]:
Image('pr8_9start.PNG', width = 800)
Out[17]:

В конце:

In [18]:
Image('pr8_9end.PNG', width = 800)
Out[18]:

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

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