Работа 8

Материал из md
Перейти к: навигация, поиск
Моделирование белка в воде

Скачать файлы для выполнения работы: 08_lysozyme.tar.gz

Скачать описание работы в pdf-формате: 08_lysozyme.pdf

Цель

Познакомиться с процедурой моделирования сложной молекулы в водном растворе на примере молекулы белка лизоци́м (мурамидаза), взятого из банка данных PDB. Провести простейший анализ структуры молекулы белка в воде.

Задания

  1. Найти в базе PDB молекулу лизоци́ма (мурамидазы).
  2. Получить начальную конфигурацию системы (молекула белка в окружении воды).
  3. Провести релаксацию системы: (а) минимизация энергии, (б) МД-релаксация.
  4. Получить равновесную МД-модель (production run).
  5. Провести анализ молекулы белка в растворе.

Общие замечания

Данная работа показывает, что имея сравнительно небольшой опыт компьютерного моделирования, можно провести молекулярно-динамическое исследование весьма сложной молекулярной системы. В пакете Gromacs реализованы наработки многих исследователей, накопленные при моделировании биологических и органических молекул. Координаты атомов и связность для многих молекул можно найти в химических базах данных. Одной из них является банк данных белков PDB (protein data bank). Стандарт записи данных о структуре молекул в этом банке (PDB-формат) применяется сейчас для записи разных молекул, не только биологических, и читается практически всеми пакетами моделирования. Все это облегчает исследователю подготовительную работу по созданию исходной молекулярной конфигурации модели.

Fig8.1.png
Рис. 8.1. Лизоци́м (мурамидаза) —антибактериальный агент; фермент, разрушающий клеточные оболочки бактерий. В пищевой промышленности зарегистрирован в качестве разрешённой пищевой добавки E1105. Состоит из 129 аминокислотных остатков.

Имея готовую молекулу и стандартные поля сил, можно получить молекулярно-динамическую модель, например, белка в растворе, исследовать структуру и динамику в зависимости от времени, температуры, давления, состава растворителя. Если нас интересуют молекулы, которых нет в базах данных, или мы хотим модифицировать поле сил, то мы должны проводить подготовительную работу самостоятельно.

Данная работа основана на демонстрационной разработке Justin’a Lemkul’a, одного из соавторов пакета Gromacs:

http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/lysozyme/index.html

Указания к работе

Получение PDB-файла молекулы

Работа выполняется в директории 08_lysozyme/. Создадим в ней поддиректории preparing/, em/, equi1/, equi2/.

PDB файл данного белка можно скачать с сайта rcsb.org (в поисковой строке сайта введите код молекулы: 1AKI ). Нужный нам файл называется 1AKI.pdb, запишем его в директорию preparing. С помощью программы VMD можно посмотреть, как выглядит эта молекула. Используя меню программы VMD (Graphics: Representation: Drawing Method — NewCartoon, Coloring Method — Structure) можно выбрать различные представления белка, в том числе изобразить α-спирали и β-складки.

Ознакомьтесь с содержимым файла 1AKI.pdb. Вы найдёте там аминокислотные остатки (аргинин, лейцин, глицин, цистеин и т.д.), из которых состоит белок, 129 штук, а также кислороды воды, которые записаны под именем HOH.

Прежде всего, надо удалить «кристаллическую» воду из PDB-файла при помощи любого текстового редактора, т.е. удалить строчки, содержащие HOH. Напомним, что координаты атомов молекулы белка извлекают из экспериментов по рентгеновской дифракции, где используется закристаллизованный белок, в котором могут присутствовать «лишние» молекулы воды. Фрагмент файла 1AKI.pdb, содержащий воду:

......
ATOM   1000  CD2 LEU A 129      37.310  12.325   5.886  1.00 25.15           C
ATOM   1001  OXT LEU A 129      43.232  12.675   6.905  1.00 34.20           O
TER    1002      LEU A 129
HETATM 1003  O   HOH A 130      23.434  40.063  -6.661  1.00 19.48           O
HETATM 1004  O   HOH A 131      31.994  26.416  -6.047  0.90 22.43           O
HETATM 1005  O   HOH A 132      30.250  13.337   9.787  0.98 20.93           O
......
HETATM 1079  O   HOH A 206      37.255   9.634  10.002  0.46 18.39           O
HETATM 1080  O   HOH A 207      43.755  23.843   8.038  0.38 17.96           O
CONECT   48  981
......

Создание начальной конфигурации

Работаем в директории preparing/. Для молекул в PDB-формате Gromacs имеет специальную программу pdb2gmx (см. подсказку по команде pdb2gmx -h), которая помимо файла с координатами в формате .gro создает автоматически файл топологии:

pdb2gmx -f 1AKI.pdb -o 1AKI.gro -water spce

где ключ -water spce определяет в качестве растворителя воду SPC.

При запуске вам потребуется выбрать поле сил из списка:

 Select the Force Field:
 From '/usr/share/gromacs/top':
  1: AMBER03 protein, nucleic AMBER94 (Duan et al., J. Comp. Chem. 24, 1999-2012, 2003)
  2: AMBER94 force field (Cornell et al., JACS 117, 5179-5197, 1995)
  3: AMBER96 protein, nucleic AMBER94 (Kollman et al.,Acc.Chem.Res. 29, 461-469, 1996)
  4: AMBER99 protein, nucleic AMBER94 (Wang et al., J. Comp. Chem. 21, 1049-1074, 2000)
  5: AMBER99SB protein, nucleic AMBER94 (Hornak et al., Proteins 65, 712-725, 2006)
  6: AMBER99SB-ILDN protein, nucleic AMBER94 (Lindorff-Larsen et al., Pro-teins 78, 1950-58, 2010)
  7: AMBERGS force field (Garcia & Sanbonmatsu, PNAS 99, 2782-2787, 2002)
  8: CHARMM27 all-atom force field (with CMAP) – version 2.0
  9: GROMOS96 43a1 force field
 10: GROMOS96 43a2 force field (improved alkane dihedrals)
 11: GROMOS96 45a3 force field (Schuler JCC 2001 22 1205)
 12: GROMOS96 53a5 force field (JCC 2004 vol 25 pag 1656)
 13: GROMOS96 53a6 force field (JCC 2004 vol 25 pag 1656)
 14: OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)
 15: [DEPRECATED] Encad all-atom force field, using full solvent charges
 16: [DEPRECATED] Encad all-atom force field, using scaled-down vacuum charges
 17: [DEPRECATED] Gromacs force field (see manual)
 18: [DEPRECATED] Gromacs force field with hydrogens for NMR

Выберем поле сил GROMOS96 53a6 под номером 13. После названия поля указан литературный источник, в котором можно ознакомиться с деталями данного поля сил. Далее программа выдаёт подробный отчёт о выполняемых действиях:

; Include forcefield parameters
#include "gromos53a6.ff/forcefield.itp"

[ moleculetype ]
; Name            nrexcl
Protein_chain_A     3

[ atoms ]
;   nr       type  resnr residue  atom   cgnr     charge       mass  typeB    chargeB      massB
; residue   1 LYS rtp LYSH q +2.0
     1       NL    1    LYS      N      1     0.129    14.0067   ; qtot 0.129
     2        H    1    LYS     H1      1     0.248      1.008   ; qtot 0.377
     3        H    1    LYS     H2      1     0.248      1.008   ; qtot 0.625
     4        H    1    LYS     H3      1     0.248      1.008   ; qtot 0.873
     5      CH1    1    LYS     CA      2     0.127     13.019   ; qtot 1
......
 1314  1312  1316  1315     2    gi_1
 1316  1314  1321  1317     2    gi_2
 1317  1319  1320  1318     2    gi_2
 1321  1316  1323  1322     2    gi_1

; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif

; Include water topology
#include "gromos53a6.ff/spce.itp"

#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
;  i funct       fcx        fcy        fcz
   1    1       1000       1000       1000
#endif

; Include topology for ions
#include "gromos53a6.ff/ions.itp"

[ system ]
; Name
LYSOZYME in water

[ molecules ]
; Compound        #mols
Protein_chain_A     1

В него включена «топология» воды (строка #include "gromos53a6.ff/spce.itp"), а также «топология» ионов (#include "gromos53a6.ff/ions.itp"), которые впоследствии нужно будет добавить в модель для нейтрализации заряда белка. В отчёте программы pdb2gmx выводится заряд белка (Total charge 8.000 e). Величину заряда нужно запомнить, чтобы позже ввести нужное число противоионов, чтобы полная система стала электронейтральной. Программа pdp2gmx создаёт также файл posre.itp, содержащий номера тяжёлых атомов белкового остова (без водородов), положение которых мы будем фиксировать в одном из последующих этапов релаксации.

Создание модельного бокса

Зададим кубический бокс для молекулы белка, используя команду editconf. Сделаем его такого размера, чтобы расстояние от атомов белка до краев бокса составляло минимум 1 нм:

editconf -f 1AKI.gro -c -d 1.0 -bt cubic -o 1AKI_box.gro

Здесь использованы ключи: -c (центрировать), -d 1.0 (минимальное расстояние от белка до края бокса), -bt cubic (кубический бокс). Результат записывается в файл 1AKI_box.gro.

Теперь заполним бокс молекулами воды командой genbox:

genbox -cp 1AKI_box.gro -cs spc216.gro -o 1AKI_solv.gro -p

Ключ -cp указывает, что координаты белка берутся из файла 1AKI_box.gro , который мы получили на предыдущем этапе подготовки, ключ -cs указывает, что координаты молекул растворителя лежат в файле spc216.gro , а ключ -o определяет имя выходного файла 1AKI_solv.gro. Ключ -p означает, что файл топологии, который используется программой, будет модифицирован. Если имя файла топологии не указано явно, то используется файл со стандартным именем topol.top.

Файл spc216.gro находится в стандартной директории /usr/share/gromacs/top/ и содержит 216 молекул воды, отрелаксированных при температуре 300К. Программа "покроет" ими всё пространство бокса, и удалит те молекулы воды, которые перекрываются с атомами белка.

Просмотрите отчёт программы о выполненных действиях. Обратим внимание, сколько в итоге было добавлено молекул воды. Убедитесь, что в topol.top добавилась строка SOL в раздел [ molecules ] в конце файла:

Protein_chain_A     1
SOL             10242

Добавление ионов

Для добавления ионов используется программа genion с файлом параметров ions.mdp, который находится в рабочей директории. Переместим ions.mdp в поддиректорию preparing/. В этом файле описывается, что при добавлении ионов будет использоваться минимизация энергии иона с учётом электростатического взаимодействия. Программа genion требует, чтобы ей подали параметры, прописанные в файле ions.mdp, в составе стартового файла topol.tpr. Таким образом, придётся сначала создать такой файл, воспользовавшись программой grompp:

grompp -f ions.mdp -c 1AKI_solv.gro

После того, как у нас появился файл topol.tpr, можно запускать genion:

genion -s -nname CL -nn 8 -p

Ключ -nname CL означает, что будет добавляться ион хлора, ключ -nn 8 указывает, что ионов будет восемь, а ключ -p означает, как и выше, что файл топологии topol.top будет модифицирован. Программа попросит выбрать группу растворителя. Надо выбрать SOL, т.е. воду, восемь молекул которой будут заменены ионами хлора:

Reading file topol.tpr, VERSION 4.5.5 (single precision)
Using a coulomb cut-off of 1 nm
Will try to add 0 NA ions and 8 CL ions.
Select a continuous group of solvent molecules
Group     0 (         System) has 11577 elements
Group     1 (        Protein) has  1323 elements
Group     2 (      Protein-H) has  1001 elements
Group     3 (        C-alpha) has   129 elements
Group     4 (       Backbone) has   387 elements
Group     5 (      MainChain) has   517 elements
Group     6 (   MainChain+Cb) has   634 elements
Group     7 (    MainChain+H) has   646 elements
Group     8 (      SideChain) has   677 elements
Group     9 (    SideChain-H) has   484 elements
Group    10 (    Prot-Masses) has  1323 elements
Group    11 (    non-Protein) has 10254 elements
Group    12 (          Water) has 10254 elements
Group    13 (            SOL) has 10254 elements
Group    14 (      non-Water) has  1323 elements
Select a group: 13
Selected 13: 'SOL'
Number of (3-atomic) solvent molecules: 3418

Processing topology
Replacing 8 solute molecules in topology file (topol.top)  by 0 NA and 8 CL ions.

Back Off! I just backed up topol.top to ./#topol.top.1#
Replacing solvent molecule 400 (atom 2523) with CL
Replacing solvent molecule 2589 (atom 9090) with CL
Replacing solvent molecule 1668 (atom 6327) with CL
Replacing solvent molecule 2891 (atom 9996) with CL
Replacing solvent molecule 1138 (atom 4737) with CL
Replacing solvent molecule 547 (atom 2964) with CL
Replacing solvent molecule 1319 (atom 5280) with CL
Replacing solvent molecule 3028 (atom 10407) with CL

В конце файла топологии появилась строка, где указано количество добавленных ионов:

CL               8

Скопируем этот файл топологии в основную директорию.

cp topol.top ../

Получившаяся конфигурация модели с ионами записана по умолчанию в файл out.gro.

Fig8.2.jpg
Рис. 8.2. Модель белка в воде с ионами хлора. Молекулы воды представлены точками (галочками), белок нарисован в виде вторичной структуры, ионы – в виде ван дер ваальсовских сфер.

С помощью VMD посмотрим на расположение вставленных ионов хлора. В основном меню VMD открыть окно Graphical Representations, создать дубликат системы, нажав Create Rep. В поле Drawing Method выбрать NewCartoon, в поле Coloring Method выбрать Secondary Structure. Вы увидите среди молекул воды, изображённых "галочками", молекулу белка в виде α-спиралей и β-складок с раскраской согласно вторичной структуре. Теперь отобразим ионы хлора. Создать дубликат системы, нажав Create Rep. В поле Selected atoms стереть надпись all, во вкладке Selections в поле Keyword выбрать имя остатка, кликнув дважды resname, тогда в правом окошке появятся названия, из которых дважды кликнуть на CL. В итоге в строке Selected atoms будет написана фраза resname CL, которую нужно «применить», нажав Apply. Затем перейти во вкладку Draw style и выбрать в поле Drawing method ван дер Ваальсовские сферы VDW.

Минимизация энергии

При приготовлении модели может оказаться, что некоторые атомы растворителя располагаются слишком близко к атомам белка и между ними действует большая сила, которая может привести к сбою программы. Поэтому прежде чем запускать релаксацию, следует сделать минимизацию энергии растворителя. Для этого предлагается дать возможность атомам смещаться до тех пор, пока максимальная сила между ними не станет меньше 1000 kJ mol-1 nm-1, см. файл параметров em.mdp. Обратим внимание, что в данном случае файл параметров моделирования имеет нестандартное имя (сокращение от «energy minimization», вместо стандартного grompp.mdp). Это сделано во избежание путаницы с файлами, однако при запуске программ придётся прописывать это имя в явном виде, например, см. ниже команду: grompp -f em.mdp.

Перейдем в поддиректорию em/ и скопируем туда файл em.mdp из основной директории. Сделаем ссылки на файл топологии и файл структуры приготовленной модели:

ln -s ../topol.top
ln -s ../preparing/out.gro conf.gro

Выполним препроцессинг с параметрами из файла em.mdp и запустим минимизацию энергии:

grompp -f em.mdp
mdrun -nt 1 -nice 11 -v

В конце минимизации на экран выводится максимальная сила, а также, потенциальная энергия системы.

Релаксация. Первый шаг

Теперь проведём обычную релаксацию растворителя. На этом шаге также предлагается зафиксировать молекулу белка, чтобы релаксация воды не вызывала нефизических изменений в белке, поскольку возникающие силы взаимодействия воды и белка все ещё могут быть большими. Для этого используется файл posre.itp, который был сгенерирован на этапе создания файла топологии программой pdb2gmx. Создаём ссылку на этот файл, а также на другие необходимые файлы (работаем в директории equi1/):

ln -s ../preparing/posre.itp
ln -s ../preparing/out.gro conf.gro
ln -s ../topol.top
ln -s ../equi1.mdp grompp.mdp

Обратите внимание на строчку в файле параметров моделирования, включающую режим удерживания атомов:

define = -DPOSRES ; position restrain the protein

Выполняем препроцессинг и моделирование:

grompp
mdrun -nt 8 -nice 11 -v

Убедимся, что релаксация прошла успешно, нарисовав поведение температуры со временем.

Fig8.3.png
Рис. 8.3. Зависимость температуры от времени в процессе релаксации.

Релаксация. Второй шаг

На первом шаге релаксации в NVT-ансамбле мы стабилизировали температуру системы. Теперь мы должны также стабилизировать давление и плотность. Работаем в NPT-ансамбле. Здесь также атомы белка удерживаются на своих местах.

Переходим в поддиректорию equi2/. Создаём все необходимые ссылки:

ln -s ../preparing/posre.itp
ln -s ../equi1/confout.gro conf.gro
ln -s ../equi2.mdp grompp.mdp
ln -s ../topol.top

проводим препроцессинг:

grompp

и запускаем релаксацию:

mdrun -nt 8 -nice 11 -v

Основное моделирование

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

ln -s equi2/confout.gro conf.gro
grompp -f prod_run.mdp
mdrun -nt 8 -nice 11 -v

Обратите внимание, что моделирование белка в течение заданного времени (1 нс) займёт приблизительно 1,5 часа машинного времени на 8 ядрах. Если время не позволяет провести расчёт до конца, то можно воспользоваться уже готовой траекторией, ссылку на которую можно получить у преподавателя.

Анализ белка

Среднеквадратичное отклонение

Первый вопрос, который возникает после помещения белка в раствор – насколько сильно молекула белка изменила свою форму? Для этого рассчитывается среднеквадратичное отклонение (RMSD) атомов белка в данной конфигурации (в текущий момент времени) относительно начальной конфигурации. Подразумевается, что текущая и начальная конфигурация белка совмещаются таким образом, чтобы сумма квадратов расстояний между соответствующими атомами была минимальна. (Заметим, что эта нетривиальная задача имеет красивое математическое решение, позволяющее находить искомое минимальное среднеквадратичное отклонение с помощью аналитических вычислений, т.е. без многократных попыток совмещения молекул в поисках наилучшего). Для этого используется команда g_rms (см. подсказку по команде g_rms -h):

g_rms -tu ns -w

ключ -tu ns означает, что мы хотим, чтобы в качестве единиц времени программа использовала наносекунды, поскольку сейчас мы занимаемся процессами, которые происходят в течение более долгого времени чем пикосекунды. Программа по умолчанию использует файл топологии topol.tpr. При запуске будет необходимо выбрать группу, для которой требуется рассчитать среднеквадратичное отклонение. Выбираем 1 (белок), дважды. Результат будет записан в файл rmsd.xvg и выведен на экран, см. Рис. 8.51.

Fig8.4.png
Рис. 8.4. Среднеквадратичное отклонение атомов белка от начальной конфигурации.

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

Радиус гирации

Радиус гирации (среднеквадратичное расстояние от атомов белка до его центра масс) характеризует размер и форму молекулы.

Fig8.5.png
Рис. 8.5. Радиус гирации белка и радиусы гирации, рассчитанные вокруг выделенных осей в зависимости от времени.

Gromacs рассчитывает не только радиус гирации R, но и радиусы гирации относительно выделенных осей, Rx, Ry, Rz. Для этого используется команда g_gyrate:

g_gyrate -w

Здесь также надо выбрать группу с атомами белка. На Рис. 8.5. показан результат для нашей молекулы. Видно, что за время моделирования форма белка существенно не изменилась.

Какую форму имеет молекула белка («вытянутую», «приплюснутую» или «сферическую»), судя радиусам гирации? Убедитесь в этом ещё раз, посмотрев молекулу белка в VMD.

Диаграмма Рамахандрана

Диаграмма Рамахандрана (Ramachandran Plot) используется биологами для изображения углов (φ, θ), отвечающих за скручивание аминокислотной цепи белка. Для выяснения, как интерпретировать различные значения углов φ, θ, воспользуйтесь, например, википедией (en.wikipedia.org/wiki/Ramachandran_plot ). Расчёт осуществляется командой g_rama (см. подсказку по команде g_rama -h)

g_rama

Результат записывается файл rama.xvg, Рис. 8.6.

Fig8.6.gif
Рис. 8.6. Диаграмма Рамачандрана белка лизоцим.


Наличие какого типа структур в белке демонстрируют нам данные углы?