CC BY
12
УДК 539.199
ОСОБЕННОСТИ РАСЧЕТОВ ПРОФИЛЕЙ СВОБОДНОЙ ЭНЕРГИИ ФЕРМЕНТАТИВНЫХ РЕАКЦИЙ: ГИДРОЛИЗ ГУАНОЗИНТРИФОСФАТА БЕЛКОВЫМ КОМПЛЕКСОМ Яаз-САР
В.А. Миронов*, Л.А. Лычко, М.Г. Хренова
(кафедра физической химии; *е-тай: [email protected])
Обсуждены способы статистической обработки результатов расчетов методом зонтичной выборки профилей свободной энергии для реакций ферментативного катализа на примере гидролиза гуанозинтрифосфата белковым комплексом Яаз-ОАР.
Ключевые слова: малые ГТФазы, свободная энергия реакции, молекулярная динамика.
Предсказание и направленное изменение действия ферментов - приоритетные задачи современной биохимии. При исследовании механизмов работы ферментов требуется информация о всех элементарных стадиях, в частности о барьерах активации. Традиционным подходом для теоретической оценки энергетических барьеров является расчет профиля минимальной энергии реакции. В основе данного подхода лежит поиск стационарных точек на поверхности потенциальной энергии (ППЭ) исследуемой системы. Реагентам и продуктам реакции отвечают минимумы на ППЭ, а переходному состоянию соответствует разделяющая их седловая точка.
Многие биологические системы обладают чрезвычайно большим числом степеней свободы; при этом реагентам, продуктам и переходным состояниям соответствует множество близких по энергии стационарных точек, что затрудняет проведение теоретического исследования. Одним из способов решения данной проблемы является расчет свободной энергии реакции, учитывающий плотность состояний исследуемой системы, например свободной энергии Гельмгольца. Для получения приемлемой точности этот подход требует трудоемких молекулярно-динамических (МД) расчетов с использованием квантово-ме-ханических потенциалов. Подобные вычисления стали доступны лишь относительно недавно с развитием суперкомпьютерных технологий и методов молекулярного моделирования. Для методов расчета свободной энергии с использованием МД точность расчета коррелирует с количеством набранных статистических данных (значения длины молекулярно-динамической траектории). В настоящее время для столь больших систем, как белковые макромолекулы, при использовании квантово-механических потенциалов время моде-
лирования редко превышает несколько десятков пикосекунд.
В данной работе проведен анализ влияния длины выбираемой траектории, а также метода статистической обработки данных на получаемые результаты для профилей свободной энергии. Для исследования выбрана реакция гидролиза гуанозинтрифосфата (ГТФ) белковым комплексом RasGAP. Выбор обусловлен не только биологической важностью данной реакции, но также и тем, что ее механизм до сих пор является предметом дискуссии. Белок Ras относится к классу малых ГТФаз, которые играют важную роль в передаче клеточных сигналов. Они являются молекулярными переключателями: активность передаваемого ими сигнала зависит от того, с какой молекулой (гуано-зинтрифосфатом или гуанозиндифосфатом (ГДФ)) они связаны. Гидролиз ГТФ до ГДФ в данных белках регулируется активирующими белками (GAP).
Гидролизу ГТФ белковым комплексом Ras-GAP посвящено много теоретических исследований [1-5]. Согласно одному из предполагаемых механизмов [2], реакция протекает посредством атаки атома кислорода реакционной молекулы воды на концевой атом фосфора молекулы ГТФ с одновременным разрывом связи между концевым атомом фосфора и мостиковым кислородом молекулы ГТФ (рис. 1). Затем происходит циклический перенос протона с молекулы воды с помощью остатка Gln61 с образованием продуктов реакции (ГДФ и неорганический фосфат, а также имидная форма остатка Gln61).
Близкий механизм предложен в работе [5]. Авторы использовали метод эмпирических валентных связей, откалиброванный по реакции гидролиза ГТФ в воде, чтобы получить профиль свободной энергии гидролиза ГТФ в белковом окружении. На втором шаге реакции были исследованы два пути,
Н
Mg о'" 'О
Л*0
О.....;—.Р— .Р:
н/ О I О О—GMP
о
Glnól
Н
н
Mgz
«s Я
>.....р ■
н/ АО" о
L Gln61
я
"•О
II,,о"
0х ^O-GMP
н Р'"
V // II ,*о
о—Р. р' ' i-"'о о' ^O-GMP ™ О
..Mg
о" -о
II II,* о р р
HO^l '""Ó Ó' ^O-GMP ОН
но
н /
«г \
01п61 н
Рис. 1. Схема механизма реакции гидролиза ГТФ в активном центре ферментативного комплекса КаБ-ОЛР
Gln61 Н
один из которых соответствовал описанному выше механизму, а другой включал участие двух молекул воды для передачи протона.
В данной работе использовались координаты тяжелых атомов модельной системы из кристаллической структуры комплекса Яая-САР, содержащего аналог ГТФ, конечная трифосфатная группа которого заменена на трифторид алюминия. Ко -ординаты тяжелых атомов брались из банка данных белковых структур РББ (РББ ГО ШР1) [6]. Группа АШ3 заменена на соответствующую группу Р03. Недостающие атомы водорода добавлены с помощью программного пакета УМБ [7]. Полученная полноатомная модель была погружена в куб воды, при этом расстояние от атомов белка до границы куба составляло не менее 10 А. Для нейтрализации заряда системы добавлены ионы На+ и С1 в физиологической концентрации (~0,15 моль/л). С помощью программы НАМБ [8] проведена процедура молекулярно-динамиче-ской релаксации полученной структуры в течение 20 нс. Для расчета использовали ланжевеновскую молекулярную динамику, поддерживали температуру 300 К и давление 1 бар. Шаг интегрирования МД составлял 1 фс. Параметры атомов белкового комплекса Яа8-ОАР, молекулы ГТФ и ионов взяты из силового поля СИАИММ [9]; для описания молекул воды использовали модель Т1Р3Р.
Молекулярно-динамическое моделирование в рамках приближения квантовой механики/молекулярной механики (КМ/ММ) проведено с помощью программы СР2К [10] при использовании гибридных базисов гауссовых функций и плоских волн. Квантово-механическая подсистема
2+
включала трифосфатную группу ГТФ, ион Mg , ближайшие к ГТФ молекулы воды, в том числе
каталитическую, и так называемый «аргинино-вый палец» белка GAP (GAP-Arg789). Кроме того, в подсистему входили остатки белка Ras (Gln61, Ser17 и Thr35), а также основные цепи остатков Gly12...Ala18 и Ala59...Gly60. ММ-подсистема описывалась в рамках силового поля CHARMM. Силы на атомах КМ-подсистемы рассчитывали с помощью теории функционала электронной плотности (функционал BLYP) с эмпирическим вкладом дисперсионной коррекции DFT-D3 [11]. Для расчетов использовали базис DZVP, оптимизированный под расчеты молекулярных систем [12] с псевдопотенциалом GTH. Расчеты МД проводили в ансамбле NVT (300 K) с использованием цепочки термостатов Нозе-Гувера.
Профили свободной энергии рассчитаны с помощью подхода зонтичной выборки (umbrella sampling, US). Данный метод позволяет получать статистические данные о высокоэнергетических областях конфигурационного пространства вдоль координаты реакции. Для этого к системе дополнительно прикладывается специально подобранный потенциал, компенсирующий изменение энергии в ходе реакции. Поскольку вклад данного потенциала в результирующий профиль аддитивный, то исходный профиль легко рассчитать из полученных при моделировании статистических данных. Для упрощения метода и ускорения сбора статистики реакционный путь обычно разделяют на несколько участков («окон»), а данные из всех окон объединяют с помощью специальных методов, о которых пойдет речь ниже. В каждом окне используется свой потенциал, обеспечивающий лучшую выборку на заданном участке; как правило, используется гармонический потенциал. В данной работе реакционный путь разбивали на 18 участков, на
Рис. 2. Сравнение результатов расчета профилей свободной энергии методами WHAM (1), UI (2) и TRAM (3)
каждом из которых проводили МД-расчет с внешним потенциалом. Длина траектории, полученной для каждого окна, составляла 15-20 пс. В качестве координаты реакции использовали расстояние от атома кислорода реакционной молекулы воды до атома фосфора концевого фосфата молекулы ГТФ.
Предложенный в настоящей работе механизм реакции согласуется с механизмом, описанным в работе [1]. Первая стадия реакции отвечает образованию интермедиата с присоединенной к ГТФ молекулой воды и является лимитирующей с барьером ~10 ккал/моль. На следующей стадии происходит циклический перенос протона с молекулы воды через остаток Gln61 на атом кислорода неорганического фосфата с образованием имидной формы глутамина.
Данные, полученные на разных участках реакционного пути, обрабатывали с помощью методов взвешенных гистограмм (WHAM [13]), интегрирования производной свободной энергии (umbrella integration, далее - UI [14]) и анализа матрицы переходов (TRAM, [15]). Как следует из названий этих методов, они существенно различаются по способу построения профиля свободной энергии. В методе WHAM комбинируются данные из разных окон и рассчитывается плотность вероятности нахождения системы вдоль реакционной координаты. Свободная энергия Гельмгольца пропорциональна логарифму плотности вероятности. В методе UI из предположения о нормальности распределения значения координаты реакции на каждом участке реакционного пути рассчитывается производная свободной энергии, а итоговый профиль свободной энергии получается численным интегрированием.
В методе TRAM используется совершенно другой подход: каждая МД-траектория рассматривается как марковская цепь, исследуются равновесия между некоторыми состояниями системы, на-
пример участками вдоль координаты реакции. Из полученной матрицы переходов рассчитывается плотность вероятности нахождения системы и ее свободная энергия в выбранных состояниях. Данный метод позволяет получить неплохие результаты даже тогда, когда выборка по конфигурационному пространству мала. Несмотря на различия в методологии, для данной системы все методы анализа результатов МД с добавлением внешнего потенциала привели к похожим результатам (рис. 2). Это свидетельствует о том, что при наличии достаточного количества статистических данных эти подходы эквивалентны.
Необходимо отметить, что для каждого окна в процессе МД должно устанавливаться локальное равновесие. Следовательно, нужно отбросить из начала траектории часть данных, соответствующих достижению этого равновесия. Для этого мы построили зависимости профилей свободной энергии от числа отброшенных кадров в каждом
Рис. 3. Сходимость профилей свободной энергии в зависимости от числа отброшенных данных в начале каждого окна метода US: а - UI, б - WHAM, в - TRAM
окне (рис. 3). Оказалось, что равновесие наступает через ~6 пс. С учетом того, что для каждого окна проводилось предварительное моделирование (2 пс), можно утверждать, что локальное равновесие для данной системы достигается за время порядка 8 пс. Следует ожидать, что сделанная оценка времени локальной релаксации будет верна и для многих других белков, в особенности для класса Работа поддержана грантом Российско
(проект №
ГТФаз и АТФаз, структура активного центра которых весьма консервативна (сходна для всего ряда белков).
Таким образом, можно заключить, что релаксация системы занимает существенную долю вычислительного времени (30-50%). При достаточной выборке все три рассмотренных подхода приводят к сходным результатам. го фонда фундаментальных исследований 14-03-31898).
СПИСОК ЛИТЕРАТУРЫ
1. Topol I.A., Cachau R.E., Nemukhin A. V., et al. Quantum chemical modeling of the GTP hydrolysis by the RASGAP protein complex. // Biochim. Biophys. Acta. 2004. Vol. 1700, № 1. P. 125.
2. Grigorenko B.L., Nemukhin A. V., Topol I.A., et al. QM/ MM modeling the Ras-GAP catalyzed hydrolysis of guanosine triphosphate. // Proteins. 2005. Vol. 60, № 3. P. 495.
3. Lu Q., Nassar N., Wang J. A mechanism of catalyzed GTP hydrolysis by Ras protein through magnesium ion // Chem. Phys. Lett. 2011. Vol. 516. N 4-6. P. 233.
4. Martín-García F., Mendieta-Moreno J.I., López-Viñas E., et al. The Role of Gln61 in HRas GTP hydrolysis: a quantum mechanics/molecular mechanics study. // Biophys. J. 2012. Vol. 102. N 1. P. 152.
5. Prasad B.R., Plotnikov N. V, Lameira J., et al. Quantitative exploration of the molecular origin of the activation of GTPase. // Proc. Natl. Acad. Sci. U. S. A. 2013. Vol. 110. N 51. P. 20509.
6. Berman H.M., Westbrook J., Feng Z., et al. The Protein Data Bank. // Nucleic Acids Res. 2000. Vol. 28. N 1. P. 235.
7. Humphrey W. VMD: Visual molecular dynamics // J. Mol. Graph. 1996. Vol. 14. N 1. P. 33.
8. Phillips J.C., Braun R., Wang W., et al. Scalable molecular dynamics with NAMD // J. Comput. Chem. 2005. Vol. 26. N 16. P. 1781.
9. MackerellA.D., FeigM., Brooks C.L. Extending the treatment of backbone energetics in protein force fields: limitations of gas-phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations. // J. Comput. Chem. 2004. Vol. 25. N 11. P. 1400.
10. CP2K developers group. CP2K: A general program to perform molecular dynamics simulations: 2.3. 2012.
11. Grimme S., Antony J., Ehrlich S., et al. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. // J. Chem. Phys. 2010. Vol. 132. N 15. P. 154104.
12. Vande Vondele J., Hutter J. Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases. // J. Chem. Phys. 2007. Vol. 127. N 11. P. 114105.
13. Roux B. The calculation of the potential of mean force using computer simulations // Comput. Phys. Commun. 1995. Vol. 91. N 1-3. P. 275.
14. Kästner J., Thiel W. Bridging the gap between thermody-namic integration and umbrella sampling provides a novel analysis method: "Umbrella integration". // J. Chem. Phys. 2005. Vol. 123. N 14. P. 144104.
15. Wu H., Mey A.S.J.S., Rosta E., et al. Statistically optimal analysis of state-discretized trajectory data from multiple thermodynamic states // J. Chem. Phys. 2014. Vol. 141. N 21. P. 214106.
Поступила в редакцию 12.02.16
METHODOLOGICAL ASPECTS OF THE CALCULATION OF THE FREE ENERGY PROFILE OF GUANOSINETRIPHOSPHATE HYDROLYSIS BY Ras-GAP PROTEIN COMPLEX
V.A. Mironov*, L.A. Lychko, M.G. Khrenova
(Division of Physical Chemistry; *e-mail: [email protected])
Applicability to enzyme catalytic reaction of methods of umbrella sampling data analysis was studied on the example of guanosinetriphosphate hydrolysis by the Ras-GAP protein complex.
Key words: small GTPases, free energy, molecular dynamics.
Сведения об авторах: Миронов Владимир Андреевич - ст. науч. сотр. лаборатории химической кибернетики, канд. физ.-матем. наук ([email protected]); Лычко Леонора Алексеевна - аспирант кафедры физической химии ([email protected]); Хренова Мария Григорьевна - вед. науч. сотр. лаборатории химической кибернетики, канд. физ.-матем. наук ([email protected]).
