Документ взят из кэша поисковой машины. Адрес оригинального документа : http://kodomo.fbb.msu.ru/~smivic/Projects/Term_6/molmod_.html
Дата изменения: Mon May 30 13:22:02 2011
Дата индексирования: Tue Oct 2 11:28:39 2012
Кодировка: Windows-1251
Молекулярная динамика и ее анализ - учебный сайт Смирновой Виктории

Учебный сайт Смирновой Виктории

Главная Проекты Семестры


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


Моделирование плавления алфа-спирального пептида в формамиде.


  1. Даны файлы:
    • Координаты пептида, 2xl1.pdb.
    • Файл с ячейкой уравновешеных молекул формамида, fam_em.gro
    • Файл дополнительной топологии для формамида, fam.itp.
    • Файл праметров для минимизации энергии em.mdp.
    • Файл праметров для "утряски" воды pr.mdp.
    • Файл праметров для молекулярной динамики md.mdp

  2. Построим файл топологии системы в силовом поле amber99sb и файл с координатами в формате Gromacs.
    pdb2gmx -f 2xl1.pdb -o pep -p pep -ff amber99sb -water tip3p
    Сделаем небольшой отступ в ячейке от ДНК.
    editconf -f pep.gro -o pep_ec -d 1.5 
    Проведем оптимизацию геометрии системы, что бы удалить "плохие" контакты в молекуле.
    grompp -f em -c pep_ec -p pep -o pep_em -maxwarn 1
    mdrun -deffnm pep_em -v
    Изменение максимальной силы в ходе оптимизации геометрии максимальная сила уменьшается (почти каждый раз на новом атоме), падая в итоге на 1 порядок: в начале F-max = 4.37039e+03 (atom 146), в конце Maximum force = 3.2095328e+02 (atom 30). Оптимизировать о значения Fmax <1, как видно, не получилось.
  3. Добавим в ячейку молекулы формамида. Количество добавленных молекул формамида - 902.
    genbox -cp pep_em -p pep -cs fam_em.gro -o pep_s

  4. Теперь надо изменить в текстовом редакторе файл тополгии pep.top. После строчки:
    ; Include forcefield parameters
    добавим
    #include "fam.itp"
    ; Include forcefield parameters
    #include "fam.itp"
    Добавим количество молекул формамида в запись [ molecules ]
    1 стало:
    [ molecules ]
    ; Compound        #mols
    Protein             1
    FAM         902

  5. Нейтрализуем заряд системы. Это делаем в два шага: строим tpr и запускаем genion. В выводе grompp сказано, что заряд системы -1, т. е. необходим 1 положительный ион для нейтрализации заряда системы (-np 1).
    grompp -f em -p pep -c pep_s -maxwarn 1 -o pep_s
    genion -s pep_s -o pep_si -p pep -np 1

  6. Проведем "утряску" формамида:
    grompp -f pr -c pep_si -p pep -o pep_pr -maxwarn 1
    mdrun -deffnm pep_pr -v
    Переформатируем pep_pr.gro и pep_si.gro в pdb формат.
    editconf -f pep_pr.gro -o  pep_pr.pdb
    editconf -f pep_si.gro -o pep_si.pdb
    
    Сравним визуально в PyMol изменеия в системах
    Изначально молекулы растворителя четко упорядочены. После утряски молекулы растворителя располагаются более хаотично.

  7. Скопируем файлы на суперкомпьтер. Запускаем тестовое моделирование на суперкомпьтере.
    grompp -f md -c pep_pr -p pep -o pep_md -maxwarn 1
    mpirun  -np 16 -maxtime 5 -q test  /home/golovin/progs/bin/mdrun_mpi -deffnm pep_md -v
    Можно записать номер задачи и проверить ход счета в файле mdrun_mpi.out-....
    При отсутствии ошибок переходим к основному моделированию. Запускаем основное моделирование на суперкомпьтере.
    mpirun  -np 16 -maxtime 1200 /home/golovin/progs/bin/mdrun_mpi -deffnm pep_md -v

    Анализ результатов.

    Параметры моделирования
    • Силовое поле используемое при построении топологии: amber99sb.
    • Заряд системы: -1.
    • Размер и форма ячейки: параллелепипед со сторонами (нм) 5.06500, 4.67000, 4.22100.
    • Минимизация энергии:
      • Алогритм минимизации энергии: l-bfgs
      • Алгоритм расчета электростатики и Ван-дер-Ваальсовых взаимодействий: Coulomb cut-off; LJ or Bcukingham cut-off
    • Модель, которой описывался растворитель: FAM (формальдегид)
    • Утряска растворителя:
      • Параметр который обуславливает неподвижность биополимера: define = -DPOSRES
      • Число шагов: 104.
      • Длина шага: 1 фс.
      • Алгоритм расчета электростатики - pme, Ван-дер-Ваальсовых взаимодействий - Cut-off.
      • Алгоритм термостата Berendsen, баростата - нет.
    • Основной расчет МД:
      • Время моделирования: 7 часов 12 минут 42 секунды количество процессоров: 16 эффективность масштабирования: 100%.
      • Длина траектории: 20 нс.
      • Число шагов: 10^7.
      • Длина шага: 2 фс.
      • Алгоритм интегратора: md (Leap-frog).
      • Алгоритм расчета электростатики - pme, Ван-дер-Ваальсовых взаимодействий - Cut-off.
      • Алгоритмы термостата v-rescale, баростата - Berendsen.

  8. Визуальный анализ движений молекул.
    trjconv -f pep_md.xtc -s pep_md.tpr -o pep_pbc_1.pdb -skip 20 -pbc mol
    Молекула все время перемещается по экрану, наблюдать за ней не удобно, поэтому применена следующая команда:
    trjconv -f pep_md.xtc -s pep_md.tpr -o pep_fit_1.pdb -skip 20 -fit rot+trans
    Наблюдения, которые можно сделать из полученного ролика:
    ? модели,
    время
    ОписаниеИзображение
    0, t=0Изначально пептид - альфа-спираль
    с 4-мя оборотами.
    5, t=1000Начинает плавиться N-конец.
    13, t=2600Начинает плавиться C-конец.
    33, t=660030-70 модели - наиболее расплавленные,
    количество витков 2-3
    58, t=1160030-70 модели - наиболее расплавленные,
    количество витков 2-3
    72, t=14400Пептид начинает обратно закручиваться
    в спираль.
    100, t=20000В конце моделирования пептид снова
    похож на спираль, хотя и не такую закрученную,
    как изначально (около 3-х витков).



  9. Определение средне-квадратичного отколнения в ходе моделирования.
    Так как у нас происходит конформационный переход, сначала рассчитаем отклонение в ходе всей симуляции относительно стартовой структуры.
    g_rms -f pep_md.xtc -s pep_md.tpr -o  rms_1
    - отклонение постепенно увеличивается примерно до 0.4, затем снова уменьшается до 0.2-0.3.

    И относительно каждой предидущей структуры на растоянии 400 кадров. Если ближе к концу закончился конформационный переход, то отколнение должно уменьшаться.
    g_rms -f pep_md.xtc -s pep_md.tpr -o rms_2 -prev 400
    - отклонение постепенно увеличивается примерно до 0.5, затем уменьшается и колеблется 0.2-0.3-0.4. (0.3 в конце)
  10. Определение изменение гидрофобной и гидрофильной поверхности в ходе конформационного перехода.
    g_sas -f pep_md.xtc -s pep_md.tpr -o sas_pep.xvg
    Зависимость изменения гидрофобной гидрофильной поверхностей доступных растворителю от времени (красный - гидрофобная, зеленый - гидрофильная):


    На ролике было видно, что наиболее расплавлен пептид в середине моделирования, когда, согласно графику, гидрофильная поверхность, доступная растворителю, максимальна. Но изменения в доступности гидрофильной поверхности не слишком велики, для гидрофобной их почти нет.
  11. Расчет колчества образуемых водородных связей.
    Если мы будем исследовать связи между пептидом и пептидом, то это будут водородные связи в пептиде. Для конца траектории:
    g_hbond  -f pep_md.xtc -s pep_md.tpr -num hbond_pep
    Не менее интересно будет изучить количество вдородных связей пептид-формамид.
     g_hbond  -f pep_md.xtc -s pep_md.tpr -num hbond_pep_sl 
    График зависимости (красный - в пептиде, зеленый - пептид-формамид):


    Количество водородных связей в пептиде почему-то почти не меняется, хотя на ролике было явно видно разрушение витков спирали при плавлении.
    Количество водородных связей с растворителем несколько варьирует в течение моделирования, но картинка, опять же, не коррелирует с увиденным в ролике.
  12. Зависимость вторичной структуры от времени:
     do_dssp  -f pep_md.xtc -s pep_md.tpr -o ss  
       # Для просмотра переведем xpm в eps
       xpm2ps -f ss.xpm -o ss.eps -by 10
    Изображение:


    Как и в ролике, на диаграмме видно, что в середине моделирования максимальные концевые участки пептида расплавились из альфа-спирали. Больше всего, как и в ролике, меняется С-конец.


© Smirnova Victoriya, 2011