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

Моделирование плавления ДНК в формамиде

Цель данной работы: моделирование плавления ДНК в формамиде

Задачи и использованные данные на Kodomo: http://kodomo.fbb.msu.ru/wiki/Main/Modelling/BioNanoPrac/DnaMeltFBB

Все команды были выполнены через консоль.

In []: 
<a href="http://kodomo.fbb.msu.ru/~pankevich-ev/term8/pr8/"></a>

Для моделирования был использован готовый небольшой дуплекс с последовательностью GATCTA: dna1.pdb (в файле - структура дуплекса).

PDB-файл со структурой ДНК-дуплекса был приведен к формату, подходящему для GROMACS: dna.pdb. Для этого были удалены фосфаты на 5'-концах, название каждого нуклеотида было изменено (G на DG и т.п.), атом C5M тимина переименован в C7.

In [55]:
In []: 
%%bash
pdb2gmx -f dna.pdb -o dna -p dna -ff amber99sb -water tip3p #1
editconf -f dna.gro -o dna_ec -d 1.5                        #2
grompp -f em -c dna_ec -p dna -o dna_em -maxwarn 1          #3
mdrun -deffnm dna_em -v
  1. Были построены файл топологии системы в силовом поле amber99sb и файл с координатами в формате Gromacs.
  2. Был сделан небольшой отступ в ячейке от ДНК.
  3. Была проведена оптимизация геометрии системы, что бы удалить "плохие" контакты в молекуле.

Начальное значение максимальной силы:

4.11919e+03 (атом 68)

Конечное значение максимальной силы:

7.3891937e+02 (атом 176)

Значение максимальной силы уменьшилось. Ниже - график значений Fmax для каждого шага.

In [28]: 
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
In [64]: 
with open('dna_forces') as f:
    w = [float(x) for x in f.readline().split()]
    forces = [[float(x) for x in line.split()] for line in f]
steps = range(90)
plt.scatter(steps,forces)
plt.show()
In []: 
%%bash
genbox -cp dna_em -p dna -cs fam_em.gro -o dna_s            #4

4 Добавили в ячейку молекулы формамида (файл топологии для формамида: fam.itp). Всего добавленных молекул: 1288 (7728 атомов). ДНК среди среди формамида представлена на рисунках ниже - до и после утряски (см. ниже).

Изменили файл топологии dna.top.

После строчки: ''' ; Include forcefield parameters ''' добавили ''' #include "fam.itp" '''.

Добавили количество молекул формамида в запись [ molecules ]:

; Compound #mols

DNA 1

FAM 1288

In []: 
%%bash
grompp -f em -p dna -c dna_s -o dna_s    #5
genion -s dna_s -o dna_si -p dna -np 10  #6

Нейтрализовали заряд системы.

5: Узнали, что заряд системы -9.99999.

6: Добавили 10 положительных ионов.

In []: 
%%bash
grompp -f pr -c dna_si -p dna -o dna_pr -maxwarn 1 #7
mdrun -deffnm dna_pr -v

7: Произвели утряску воды.

In [2]: 
from xmlrpclib import ServerProxy
from IPython.display import Image
import os, sys, time

# pymol launching
import __main__
__main__.pymol_argv = [ 'pymol', '-cp' ]
# __main__.pymol_argv = [ 'pymol', '-x' ]  # for GUI

import pymol
from pymol import cmd
pymol.finish_launching()
In [3]: 
defaultImage = 'pymolimg.png'
def prepareImage(width=300, height=300, sleep=5, filename=defaultImage):
    ## To save the rendered image
    cmd.ray(width, height)
    cmd.png('pymolimg.png')
    time.sleep(sleep)
    
# Define some shortcuts
def focus(x):
    cmd.center(x)
    cmd.zoom(x)
In [4]: 
cmd.load("dna_si.pdb")
prepareImage()
Image(defaultImage)
Out[4]:

Можно вечно смотреть, как упорядочен формамид.

In [5]: 
cmd.load("dna_pr.pdb")
prepareImage()
Image(defaultImage)
Out[5]:

Формамид потрясен.

In []: 
%%bash
grompp -f md -c dna_pr -p dna -o dna_md -maxwarn 1 #8

8: Получили финальный фал для молдинамики

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

Перед финальным запуском:

In []: 
%%bash
module load cuda/6.5.14

Номер задачи: 1268785. Следить за счетом можно через просмотр последних файлов в директории на суперкомпьютере (ls -lt).

In []: