Моделирование механизмов химических процессов в биомолекулярных системах методами молекулярной динамики
На правах рукописи
КАЛИМАН ИЛЬЯ АЛЕКСАНДРОВИЧ МОДЕЛИРОВАНИЕ МЕХАНИЗМОВ ХИМИЧЕСКИХ ПРОЦЕССОВ В БИОМОЛЕКУЛЯРНЫХ СИСТЕМАХ МЕТОДАМИ МОЛЕКУЛЯРНОЙ ДИНАМИКИ Специальность 02.00.17 – математическая и квантовая химия
АВТОРЕФЕРАТ
диссертации на соискание учёной степени кандидата физико-математических наук
Москва – 2011
Работа выполнена в лаборатории химической кибернетики на кафедре физической химии химического факультета Московского Государственного Университета имени М.В. Ломоносова
Научный консультант: доктор химических наук, профессор Немухин Александр Владимирович
Официальные оппоненты: доктор физико-математических наук, профессор Тихонов Александр Николаевич Физический факультет МГУ имени М.В. Ломоносова кандидат физико-математических наук Мазо Михаил Абрамович Институт химической физики имени Н.Н. Семенова РАН
Ведущая организация: Институт математических проблем биологии РАН
Защита состоится «8» декабря 2011 года в 15:00 в аудитории 446 Химического факультета МГУ на заседании диссертационного совета Д 501.001.50 при МГУ имени М.В. Ломоносова (119991, Москва, ГСП-1, Ленинские горы, д. 1, стр. 3, Химический факультет МГУ).
С диссертацией можно ознакомиться в научной библиотеке Химического факультета МГУ имени М.В. Ломоносова.
Автореферат диссертации размещен на сайте Химического Факультета МГУ имени М.В.
Ломоносова (www.chem.msu.ru) и на сайте Высшей Аттестационной Комиссии (vak.ed.gov.ru).
Автореферат разослан «2» ноября 2011 года.
Ученый секретарь диссертационного совета Д 501.001.50, кандидат химических наук Матушкина Н.Н.
Общая характеристика работы
Актуальность темы Метод молекулярной динамики является одним из важнейших теоретических методов изучения молекулярных систем. Применение молекулярно-динамического моделирования для исследования процессов, сопровождающихся разрывом и образованием химических связей, является особенно актуальным, и находит важное применение для расчёта профилей энергии Гиббса химических реакций. Знание активационных барьеров на поверхности энергии Гиббса позволяет непосредственно связывать данные теоретического моделирования с результатами кинетических измерений.
Диссертация посвящена развитию и применению новых подходов к описанию молекулярной динамики химических процессов в биомолекулярных системах. Реализованы процедуры расчёта градиентов энергии для метода потенциалов эффективных фрагментов, а также разработана оригинальная реализация метода молекулярной динамики с жесткими фрагментами. Эти методы позволяют проводить молекулярно-динамическое моделирование таких больших систем как белки, используя ресурсы современных суперкомпьютеров.
Молекулярные объекты, исследованные в данной работе, включают трансмембранный ионный канал грамицидина А, входящий в состав антибиотика грамицидина, и моторный белок миозин. Рассмотрена спонтанная изомеризация остатков аспарагина, которая играет важную роль в процессах деградации белков и является одной из начальных стадий в цепи деструктивных изменений в белках, приводящих к болезни Альцгеймера.
Использование неэмпирических и гибридных методов квантовой химии для вычисления сил в молекулярно-динамических расчётах из-за их больших требований к компьютерным ресурсам стало возможным только в последние несколько лет с появлением мощных суперкомпьютеров. Актуальность данной работы обусловлена возрастающим интересом к молекулярно-динамическим расчётам с использованием неэмпирических методов для вычисления сил. Эти методы позволяют рассчитывать профили энергии Гиббса для химических процессов в таких сложных системах как белки и растворы. Результаты компьютерных расчетов на основе современных методов молекулярного моделирования, а именно, молекулярной динамики, квантовой химии и комбинированных подходов квантовой и молекулярной механики (КМ/ММ), позволяют существенно дополнить экспериментальные исследования, добавить новую информацию об исследуемых объектах, а также обеспечить визуализацию процессов происходящих в молекулярных системах.
Цель работы заключалась в развитии и практическом применении методов молекулярной динамики с квантовыми силами для изучения механизмов важнейших биохимических реакций. В соответствии с этим были поставлены следующие задачи:
1. Реализовать процедуры расчёта градиентов энергии для метода потенциалов эффективных фрагментов и метод молекулярной динамики с жёсткими фрагментами.
2. Исследовать влияние выбора квантовой подсистемы на результаты теоретического моделирования методами КМ/ММ на примере ферментативной реакции гидролиза ацетилхолина.
3. Изучить реакцию образования сукцинимида из остатков аспарагина в водном окружении с использованием полностью квантового описания системы. Построить профили энергии Гиббса для всех стадий данного процесса.
4. Изучить транспорт протона в канале грамицидина А. Исследовать механизмы переноса протона и переноса отрицательного ионного дефекта.
5. Изучить реакцию гидролиза аденозинтрифосфата в миозине, а также последующий выход продуктов реакции из активного центра фермента.
Личный вклад диссертанта заключается в постановке и разработке путей решения поставленных задач, анализе литературных данных, разработке оригинального программного обеспечения для проведения расчётов методами квантовой химии и молекулярной динамики, проведении вычислений методами квантовой химии, комбинированными методами квантовой и молекулярной механики, методом классической молекулярной динамики, интерпретации результатов, подготовке публикаций и докладов по теме диссертационной работы.
Научная новизна результатов 1. Реализованы процедуры расчёта градиентов для метода потенциалов эффективных фрагментов в пакете молекулярного моделирования Q-Chem. Созданы оригинальные компьютерные программы, реализующие метод молекулярной динамики с жесткими фрагментами и метод анализа взвешенных гистограмм.
2. Комбинированным методом квантовой и молекулярной механики проведено моделирование реакции гидролиза ацетилхолина в активном центре ацетилхолинэстеразы. Проанализировано влияние выбора квантово-механической подсистемы на результаты теоретического моделирования.
3. Впервые рассчитаны профили энергии Гиббса с полным учётом водного окружения и полностью квантово-механическим описанием системы для реакции изомеризации аспарагина. Показано, что при моделировании данной реакции чрезвычайно важным является учет растворителя. Установлено, что скорость-лимитирующими стадиями являются стадии депротонирования и циклизации.
4. Исследованы два конкурирующих механизма протонного транспорта в канале грамицидина А с использованием метода молекулярной динамики с КМ/ММ потенциалами. Показано, что переориентация цепи из молекул воды является скорость определяющей стадией. Установлено, что стадии непосредственного переноса протона и переноса ионного дефекта являются безбарьерными.
5. Показано, что двух-водный механизм реакции гидролиза аденозинтрифосфата в миозине связан с процессом выхода продуктов гидролиза из активного центра фермента по механизму «чёрного хода».
Научная и практическая значимость работы заключается в том, что на примере различных классов биомолекул продемонстрированы возможности применения методов квантовой, КМ/ММ и классической молекулярной динамики для исследования таких систем.
Созданы новые программные продукты, реализующие ряд современных подходов к моделированию динамики биомолекулярных систем. В рамках разработки пакета молекулярного моделирования Q-Chem реализованы процедуры расчёта градиентов энергии для метода потенциалов эффективных фрагментов. Созданы оригинальные реализации метода молекулярной динамики с жёсткими фрагментами и метода анализа взвешенных гистограмм.
Результаты данного исследования помогают детализировать механизмы химических реакций и, следовательно, дают принципиальную возможность влиять на их ход. Например, понимание механизма изомеризации остатков аспарагина является важным при разработке лекарственных препаратов для лечения и профилактики болезни Альцгеймера. Результаты этой работы могут быть использованы в разработке биологически активных соединений и новых лекарственных препаратов.
Апробация работы и публикации. Материалы диссертации были представлены на международных конференциях «Ломоносов» (Москва 2007, 2008, 2009 и 2011), английской секции международной конференции «Ломоносов» (Москва 2009), IX Ежегодной международной молодежной конференции «Биохимическая физика ИБХФ РАН-ВУЗы» (Москва 2009), 6-ой Всероссийской Школе-Симпозиуме «Динамика и Структура в Химии и Биологии» (Москва 2008), XVI и XVIII международных конференциях «Математика. Компьютер.
Образование» (Пущино 2009 и 2011), международной конференции «Biocatalysis: Fundamentals & Applications» (Архангельск 2009), XXII симпозиуме «Современная Химическая Физика» (Туапсе 2010), Интернет-конференции «Информационно-вычислительные технологии в науке», ИВТН 2008.
Результаты опубликованы в 19 работах, в том числе, в 7 статьях в рецензируемых научных журналах, входящих в перечень журналов ВАК РФ.
Структура работы. Диссертация состоит из введения, 5 глав, заключения, выводов и списка цитируемой литературы из 144 наименований. Работа изложена на 129 страницах машинописного текста и включает 42 рисунка и 4 таблицы.
Содержание работы Глава 1 включает обзор литературы по основным методам расчёта энергетических характеристик и моделирования динамики молекулярных систем, а также рассматривает методы расчёта профилей энергии Гиббса для химических процессов. Подробно описан метод потенциалов эффективных фрагментов, для которого в рамках данной работы реализованы процедуры расчёта градиентов энергии в пакете молекулярного моделирования Q-Chem.
В первой части главы рассмотрены методы описания электронной структуры модельных систем, такие как метод Хартри-Фока и метод теории функционала электронной плотности.
Описаны недостатки и преимущества данных подходов при их использовании для расчёта сил в молекулярно-динамическом моделировании.
Другие традиционные методы квантовой химии, такие как метод теории возмущений, метод МКССП, а также методы связанных кластеров, ввиду их чрезвычайно высокой ресурсоемкости, не позволяют применять их в молекулярно-динамическом моделировании для расчёта сколько-нибудь значимых длин траекторий.
Важным классом методов описания поверхности потенциальной энергии молекулярных систем являются комбинированные методы квантовой и молекулярной механики (КМ/ММ). По принципам взаимодействия квантовой части с окружением методы комбинированной квантовой и молекулярной механики делятся на методы механического внедрения и методы электронного внедрения. Как отдельного представителя последнего класса следует выделить метод потенциалов эффективных фрагментов.
В методах механического внедрения взаимодействие квантово-механической и молекулярно-механической частей описывается на уровне молекулярной механики.
Молекулярно-механическая подсистема воздействует на квантовую часть на уровне сложения сил и не вносит поправок в электронный гамильтониан системы. Более корректное описание дают методы, основанные на электронном внедрении. Эти методы предполагают добавление вкладов от точечных зарядов атомов молекулярно-механической части в электронный гамильтониан системы. Метод потенциалов эффективных фрагментов обеспечивает наиболее точное описание взаимодействия между квантовой и молекулярно-механической подсистемами.
В рамках данного подхода влияние окружения на центральную подсистему описывается главными вкладами одночастичных операторов межмолекулярного взаимодействия в гамильтониан квантово-механической части.
Третий рассмотренный класс методов относится к варианту, использующему базисы плоских волн. Так как эти методы менее требовательны к компьютерным ресурсам по сравнению с методами, использующими гауссовы базисы, то они особенно популярны при молекулярно динамическом моделировании. Методы, использующие плоские волны, хорошо подходят для моделирования конденсированной фазы с периодическими граничными условиями.
Во второй части главы описаны методы расчета профилей энергии Гиббса: метод метадинамики и метод зонтичной выборки. Метод метадинамики1, или метод холмов, – это относительно новая разновидность метода молекулярной динамики. Этот подход направлен на улучшение сканирования отдельных областей конфигурационного пространства моделируемых систем. Он позволяет отобразить гиперповерхность энергии Гиббса как функцию небольшого числа обобщённых координат. Этот метод обладает достаточно высокой эффективностью, что позволяет применять его в комбинации с ab initio молекулярной динамикой, в частности, с Борн Оппенгеймеровской молекулярной динамикой и с методом Кара-Паринелло. Методы ab initio молекулярной динамики применяют для изучения процессов, в которых изменения в электронном строении играют ключевую роль, таких, например, как химические реакции.
Метод зонтичной выборки2 представляет собой метод неравновесной молекулярной динамики, позволяющий улучшить сканирование конфигурационного пространства в заданных областях. Основная идея данного подхода состоит в разбиении координаты реакции на отрезки, называемые окнами, в каждом из которых затем проводятся независимые молекулярно динамические расчёты. Для того чтобы система не выходила за пределы некоторого диапазона координаты реакции, на атомы, составляющие координату реакции, накладывают некоторый ограничивающий потенциал. Профиль энергии Гиббса можно получить, объединяя результаты моделирования с помощью метода анализа взвешенных гистограмм, подробно описанного в диссертации.
Laio A., Parrinello M. Escaping free-energy minima // Proc. Nat. Acad. Sci. 2002. V. 99. P. 12562-12566.
Roux B. The calculation of the potential of mean force using computer simulations // Comp. Phys. Comm. 1995. V. 91. P.
275-282.
Применение рассмотренных методов к конкретным системам для решения прикладных задач описано в главах 3, 4 и 5. Обзор литературных данных, касающихся этих систем, приведен в соответствующих главах.
В главе 2 описано исследование влияния выбора квантовой подсистемы на результаты моделирования с использованием комбинированных методов квантовой и молекулярной механики. Тщательный выбор квантовой подсистемы оптимального размера особенно важен при проведении молекулярно-динамических расчётов с квантовыми силами. Слишком большая квантовая часть может привести к тому, что расчёт станет чрезвычайно ресурсоёмким. С другой стороны, делая квантовую подсистему маленькой, исследователь рискует получить неадекватное описание моделируемой системы. В данной части работы на примере фермента ацетилхолинэстеразы изучена зависимость результатов моделирования реакции гидролиза ацетилхолина от выбора квантовой подсистемы.
Исследованная система включала часть фермента ацетилхолина, содержащую 3533 атома.
Для изучения влияния выбора квантово-механической подсистемы на энергетические характеристики реакции гидролиза был проведен ряд расчётов с последовательно увеличивающейся квантово-механической частью. В таблице 1 приведены аминокислотные остатки, входившие в квантовую подсистему восьми рассмотренных модельных систем.
Табл. 1. Сравниваемые квантово-механические подсистемы.
Число атомов в КМ Модель Состав КМ подсистемы подсистеме His447, Ser203, Субстрат Glu334, His447, Ser203, Субстрат Ser229, Glu334, His447, Ser203, Субстрат Ala204, Ser229, Glu334, His447, Ser203, Субстрат Gly120, Ala204, Ala207, Ser229, Glu334, His447, Ser203, Субстрат Gly121, Gly122, Gly120, Ala204, Ala207, Ser229, Glu334, His447, Ser203, Субстрат Gly206, Gly121, Gly122, Gly120, Ala204, Ala207, Ser229, Glu334, His447, Ser203, Субстрат Glu202, Gly206, Gly121, Gly122, Gly120, Ala204, Ala207, Ser229, Glu334, His447, Ser203, Субстрат Расчёты проводили с использованием комбинированного метода квантовой и молекулярной механики в варианте электронного внедрения. Для описания квантовой подсистемы применяли метод теории функционала электронной плотности с функционалами BB1K и B3LYP. Исследование показало, что обменно-корреляционные функционалы BB1K и B3LYP дают практически идентичные рассчитанные геометрические конфигурации, а отличия в энергиях составляют менее 1 ккал/моль. В дальнейшем приводятся результаты только для функционала BB1K.
Рассчитанные энергетические профили для восьми модельных систем представлены на Рис. 1. Светлой линией выделен профиль для системы с самой маленькой квантово-механической частью (модель 1), тёмной – с самой большой (модель 8). Анализ полученных данных показывает, что для корректного описания данной системы помимо субстрата, гистидина 447 и серина 203, в квантовую часть необходимо включать глютаминовую кислоту 334. Исключение аминокислотного остатка глютаминовой кислоты 334 из квантово-механической подсистемы (модель 1) ведёт к завышенному значению энергий тетраэдрического интермедиата и второго переходного состояния по сравнению с остальными модельными системами.
Рис. 1. Рассчитанные энергетические профили для реакции гидролиза ацетилхолина ацетилхолинэстеразой. ПС1 и ПС2 – первое и второе переходные состояния, ТИ – тетраэдрический интермедиат, АФ – ацилфермент, ФСК – фермент-субстратный комплекс.
Необычно высокая стабильность тетраэдрического интермедиата в реакциях гидролиза ацетилхолина холинэстеразами была подробно изучена экспериментально. Исследование вторичного изотопного эффекта показало, что в реакции гидролиза ацетилхолина ацетилхолинэстеразой идет накопление стабильного тетраэдрического интермедиата, и его гидролиз является лимитирующей стадией3. Полученные в данной работе результаты хорошо согласуются с этими экспериментальными данными.
Таким образом, на примере химических превращений в активном центре ацетилхолинэстеразы показано, что можно выделить минимально допустимый состав молекулярной подсистемы, относимой к КМ-части, который необходим для адекватных расчетов энергетических профилей элементарных стадий реакций ферментативного катализа. В соответствии с результатами анализа в данном случае КМ-подсистема должна включать все аминокислотные остатки каталитической триады – Ser203, His447 и Glu334, а также субстрат.
Невыполнение этого условия приводит к качественно неверному энергетическому профилю (модель 1 на Рис. 1). Во всех остальных случаях, рассмотренных в данной работе (Рис. 1, модели 2-8), качественно результаты близки. Сопоставление энергетических профилей для моделей 2- позволяет оценить возможный коридор погрешностей при расчетах методом КМ/ММ. Для величин энергий барьеров первого переходного состояния (относительно энергий реагентов) и второго переходного состояний (относительно энергий итермедиата) этот разброс составляет порядка 2 ккал/моль. Разброс для энергий интермедиата (относительно реагентов) составляет порядка 5 ккал/моль. Стоит отметить, что для таких методов, как комбинированные методы квантовой и молекулярной механики, общепринятыми являются величины ошибок порядка 1- ккал/моль, что позволяет утверждать о близости энергетических характеристик, полученных для моделей 2-8.
В главе 3 описано применение метода квантовой молекулярной динамики для моделирования реакции изомеризации аминокислотных остатков аспарагина в водном окружении. С помощью метода метадинамики проведено моделирование первой стадии изомеризации остатков аспарагина, а именно, образование циклического интермедиата сукцинимида. Предполагаемый механизм реакции, включающий стадии депротонирования, циклизации и деаминирования, представлен на Рис. 2.
Tormos J.R., Wiley K.L., Wang Y., Fournier D., Masson P., Nachon F., Quinn D.M. Accumulation of Tetrahedral Intermediates in Cholinesterase Catalysis: A Secondary Isotope Effect Study // J. Am. Chem. Soc. 2010. V. 132. P.
17751-17759.
Рис. 2. Механизм образования циклического имида из остатков аспарагина.
В первой части главы 3 представлен обзор литературных данных, касающихся данного процесса. Также описаны методики и детали молекулярно-динамического моделирования.
Описана процедура верификации применяемых методов, а также дан анализ возможных ошибок, вносимых используемым методом на основе плосковолнового базиса.
Для моделирования механизма образования циклического имида был применен метод метадинамики, подробно описанный в первой главе работы. Для описания электронной структуры использовали метод теории функционала электронной плотности с базисом плоских волн. Это позволило описать на квантовом уровне теории не только молекулу аспарагина, но и окружающие молекулы растворителя. Для того чтобы корректно описать поведение раствора использовали периодические граничные условия. Изображение элементарной ячейки модельной системы представлено на Рис. 3. Ещё раз следует подчеркнуть, что в данном случае вся система (элементарная ячейка содержит молекулу пептида и 104 молекулы воды) описана с использованием метода теории функционала электронной плотности.
Рис. 3. Периодическая ячейка модельной системы для исследования реакции изомеризации аминокислотных остатков аспарагина в водном окружении.
Во второй части главы изложены результаты теоретического моделирования.
Представлены рассчитанные профили энергии Гиббса для стадий депротонирования, циклизации и деаминирования.
Стадия депротонирования. На данной стадии происходит образование цвиттериона после ухода атома водорода пептидной группы в объёмную фазу. В качестве координаты реакции использовали расстояние от азота пептидной группы до уходящего протона. Из анализа молекулярно-динамической траектории установлено, что уходящий протон переходит на молекулу воды растворителя, а не непосредственно на аминогруппу. Протонирование аминогруппы является результатом переорганизации сети водородных связей, в результате которой происходит перескок другого протона с молекулы растворителя на данную аминогруппу. Важным наблюдением является то, что образующиеся на данной стадии цвиттерионные сущности стабильны только в окружении молекул воды. Попытки оптимизировать получившуюся структуру в отсутствие окружения приводили к тому, что протон с аминогруппы переходил обратно на азот пептидной связи. Из этого сделан вывод, что даже для качественного воспроизведения механизма данной реакции необходимо учитывать молекулы растворителя в явном виде. Рассчитанный профиль энергии Гиббса для стадии депротонирования представлен на левой половине Рис. 4. Из полученного графика видно, что процесс характеризуется энергетическим барьером порядка 20 ккал/моль.
Рис. 4. Потенциал средней силы для стадий депротонирования (слева) и циклизации (справа).
Стадия циклизации. На следующем этапе была рассмотрена стадия циклизации, которая соответствует образованию циклического тетраэдрического интермедиата (интермедиат 2 на Рис.
3). Формирование кольца происходит после атаки азотом пептидной группы углерода амидной группы. При рассмотрении данной стадии была использована координата реакции, соответствующая расстоянию между атомом азота пептидной группы и карбонильным углеродом. Рассчитанный профиль энергии Гиббса представлен на Рис. 4 справа. Как видно из представленного графика, энергетический барьер процесса также составляет порядка ккал/моль. Анализ молекулярно-динамической траектории показал, что изначально система находится в достаточно широкой потенциальной яме, что соответствует координате реакции между 3.2 и 4.8. Это объясняется свободным движением пептидной цепи. После того как система покидает эту потенциальную яму и проходит через переходное состояние, она попадает в локальный минимум, соответствующий образованию химической связи N-C (область в районе 1.5 на правом графике на Рис. 4), вследствие чего происходит образование циклического интермедиата. Завершение второй стадии процесса характеризуется переносом протона с аминогруппы на образующийся в ходе реакции отрицательно заряженный атом кислорода, что ведёт к формированию гидроксильной группы, образующей водородную связь с упомянутой аминогруппой.
Стадия деаминирования. Завершающей стадией исследуемой реакции является стадия деаминирования с образованием продукта – сукцинимида (см. Рис. 2). Обобщенной координатой реакции данного процесса является расстояние между атомом азота уходящей аминогруппы и смежным атомом углерода. Полученный профиль энергии Гиббса показан на Рис. 5.
Энергетический барьер для данной стадии составил порядка 15 ккал/моль (пик в районе 2 ).
Профиль, соответствующий данной стадии, характеризуется самым низким барьером, что говорит о том, что деаминирование проходит относительно быстро по сравнению с двумя предыдущими стадиями. Анализ молекулярно-динамической траектории показал, что уходящей группой является аммиак (NH3). Изначально существующая уходящая аминогруппа группа пептида протонируется атомом водорода из растворителя.
Рис. 5. Потенциал средней силы для стадии деаминирования.
Небольшой пик в районе значений координаты, соответствующих 4 по оси абсцисс, можно объяснить недостаточным сканированием тех областей, где образовавшийся аммиак отдалён от сформировавшегося сукцинимида. Однако данные области не представляют особого интереса, так как не влияют на величину энергетического барьера реакции.
Полученные в результате расчетов величины барьеров хорошо согласуются с экспериментальным значением общей эффективной энергии активации, составляющим порядка 22 ккал/моль4.
Сравнение полученных в настоящей работе данных с литературными источниками, в которых моделирование проводят без учёта водного окружения, показывает, что пренебрежение окружением даёт завышенные энергетические барьеры. Более того, в работе показано, что один из интермедиатов (интермедиат 1 на Рис. 2), образующийся в ходе данного процесса, устойчив только в присутствии водного окружения. Некорректный учёт окружения ведет к неправильному механизму процесса. Полученные данные показывают важность явного учета водного окружения для корректного моделирования данного процесса.
Глава 4 посвящена моделированию протонного транспорта в канале грамицидина А. В этой главе представлены результаты моделирования двух механизмов процесса переноса протона Patel K., Borchardt R.T. Chemical Pathways of Peptide Degradation. II. Kinetics of Deamidation of an Asparaginyl Residue in a Model Hexapeptide // Pharm. Res. 1990. V. 7. P. 703-711.
в канале грамицидина: механизма непосредственного переноса протона (механизм Гроттхуса) и механизма с участием переноса отрицательного ионного дефекта.
Грамицидин А представляет собой трансмембранный ионный канал, участвующий в процессах переноса однозарядных катионов щелочных металлов и протонов. Структура канала представлена на Рис. 6.
Рис. 6. Структура канала грамицидина А. Атомы водорода канала не показаны для наглядности представления.
Известно, что скорость переноса протонов по каналу превышает скорость переноса других однозарядных катионов более чем на порядок. Это объясняется существованием структурированной цепи из молекул воды в полости канала (см. Рис. 6). Эта цепь может служить переносчиком протонов по механизму Гроттхуса, одновременно затрудняя транспорт других катионов. Однако детали реализации этого механизма переноса протона в грамицидине до сих пор остаются неясными.
Новые экспериментальные данные позволили предложить другой вариант механизма протонного транспорта, заключающийся в переносе отрицательного ионного дефекта в сторону5.
противоположную переносу протона Подробное описание механизмов непосредственного переноса протона и механизма с участием переноса отрицательного ионного дефекта приведено в диссертации.
В данной части работы были рассчитаны потенциалы средней силы для стадий непосредственного переноса протона, переноса ионного дефекта, а также для стадии переориентации цепи из молекул воды, завершающей перенос в обоих случаях. На Рис. представлены профили энергии Гиббса для стадий непосредственного переноса протона (левый Rokitskaya T.I., Kotova E.A., Antonenko Y.N. Membrane dipole potential modulates proton conductance through gramicidin channel: movement of negative ionic defects inside the channel // Biophys. J. 2002. V. 82. P. 865-873.
график) и переноса ионного дефекта (правый график). Для их расчёта использовали оригинальную реализацию метода анализа взвешенных гистограмм.
Рис. 7. Потенциалы средней силы для стадий непосредственного переноса (слева), переноса ионного дефекта (справа), полученные с использованием модели PM6.
Для описания межатомных взаимодействий при моделировании был применен гибридный метод, сочетающий метод PM6 для описания водной цепи и метод молекулярной механики с силовым полем CHARMM для описания канала грамицидина. Как видно из графиков на Рис. 7, перенос, как ионного дефекта, так и непосредственный перенос протона, осуществляются безбарьерно. В то же время, стадия переориентации водной цепи характеризуется барьером в 5. ккал/моль (левый график на Рис. 8).
Поскольку проведенные расчеты показали, что стадия переориентации является скорость лимитирующей, в работе проведено дополнительное исследование этой стадии, используя гибридный метод теории функционала электронной плотности и молекулярной механики для расчёта молекулярно-динамических сил.
Рис. 8. Потенциал средней силы для стадии переориентации водной цепи по результатам моделирования с использованием модели PM6 (слева) и КМ/ММ моделирования (справа).
В данном случае для описания цепи из молекул воды, которая была выделена в квантовую часть, использовали метод теории функционала электронной плотности с гибридным градиент корректированным функционалом B3LYP и базисом 6-31G*. Для описания канала грамицидина использовали метод молекулярной механики и силовое поле AMBER. Рассчитанный профиль энергии Гиббса представлен на Рис. 8 справа. Барьер энергии Гиббса в данном случае составил 7.7 ккал/моль, что хорошо согласуется с экспериментальными данными, дающими значение эффективной энергии активации для всего процесса равное 6.5 ккал/моль 6. Также стоит отметить хорошее качественное согласие между профилями, полученными с использованием метода PM (левый график на Рис. 8) и метода КМ/ММ (правый график на Рис. 8) для расчёта молекулярно динамических сил.
Пятая глава работы посвящена изучению реакции гидролиза аденозинтрифосфата (АТФ) в миозине. Миозин – важнейший белок, являющийся движущей силой для работы мышечной ткани, получающий энергию в результате реакции гидролиза АТФ. В данной части работы проведено моделирование реакции гидролиза с использованием комбинированных методов квантовой и молекулярной механики, а также молекулярно-динамическое моделирование процесса открытия солевого мостика Arg238-Glu459. Полученные результаты подтверждают возможность выхода продуктов реакции из активного центра миозина через данный солевой мостик, то есть реализацию механизма «чёрного хода».
Рис. 9. Аденозинтрифосфат в активном центре миозина. Шарами показаны атомы, входящие в квантовую подсистему.
Chernyshev A., Cukierman S. Thermodynamic View of Activation Energies of Proton Transfer in Various Gramicidin A Channels // Biophys. J. 2002. V. 82. P. 182-192.
В первой части главы описано моделирование реакции гидролиза АТФ в миозине по так называемому двух-водному механизму. Суть данного механизма состоит в том, что в активном центре миозина, наряду с АТФ, присутствуют две молекулы воды, образующие водородную связь и участвующие в реакции гидролиза с образованием неорганического фосфата.
Схематичное изображение процессов, происходящих в активном центре миозина во время реализации данного механизма, представлено на Рис. 9.
Для исследования реакции гидролиза применяли гибридные методы теории функционала электронной плотности и молекулярной механики. При моделировании была использована квантовая часть, включающая все аминокислотные остатки активного центра фермента.
Используемая квантовая подсистема в окружении части атомов молекулярно-механической части приведена на Рис. 9.
Расчет энергетических профилей реакции гидролиза был проведен с использованием методики комбинированной квантовой и молекулярной механики в варианте электронного внедрения. Полученный активационный барьер для реакции гидролиза АТФ в миозине составил 10.4 ккал/моль, что отлично согласуется с экспериментальной оценкой барьера равной 10. ккал/моль7. Полученный энергетический профиль реакции представлен на Рис. 10.
Рис. 10. Энергетический профиль реакции гидролиза АТФ. ФС – фермент-субстратный комплекс, ПС – переходное состояние, ФП – комплекс фермент-продукт.
На втором этапе работы было исследовано то, как влияет протонирование остатка глютаминовой кислоты 459 в результате гидролиза молекулы АТФ на энергетику процесса открытия солевого мостика. В ходе исследования были рассчитаны профили энергии Гиббса для 7 Song Z., Parker K.J., Enoh I., Zhao H., Olubjao O. Myosin-catalyzed ATP hydrolysis elucidated by P NMR kinetic studies and H PFG-diffusion measurements // Anal. Bioanal. Chem. 2009. V. 395. P. 1453-1459.
процесса открытия солевого мостика Arg238-Glu459 с использованием структур фермент субстратного комплекса и комплекса фермент-продукт, полученных на предыдущем этапе с применением метода КМ/ММ (см. Рис. 11).
Рис. 11. Структура активного центра миозина для комплексов фермент-субстрат (сверху) и фермент-продукт (снизу).
Для расчета профилей энергии Гиббса был использован метод зонтичной выборки. Для объединения результатов применяли метод анализа взвешенных гистограмм. Полученные в результате молекулярно-динамических расчетов энергетические профили процесса открытия солевого мостика представлены на Рис. 12.
Рис. 12. Потенциалы средней силы для процесса открытия солевого мостика в комплексах фермент-субстрат (1) и фермент-продукт (2).
Анализ рассчитанных профилей показывает, что для открытия солевого мостика Arg238 Glu459 в случае комплекса фермент-продукт (кривая 2) необходимо порядка 16 ккал/моль, что почти вдвое меньше, чем для фермент-субстратного комплекса (кривая 1). Также обнаружено, что после протонирования остатка глютаминовой кислоты в ходе реакции гидролиза, расстояние между Arg238 и Glu459 увеличивается (кривая для комплекса фермент-продукт на Рис. смещается в сторону больших расстояний). Это можно объяснить уменьшением заряда на глютаминовой кислоте, а также создаваемыми лишним протоном стерическими затруднениями.
Опираясь на эти данные, можно сделать вывод, что в результате реакции гидролиза АТФ по рассмотренному механизму, включающему протонирование глютаминовой кислоты, облегчается открытие солевого мостика Arg238-Glu459. Это, в свою очередь, может способствовать выходу продуктов гидролиза АТФ из активного центра миозина по механизму «чёрного хода».
В заключительной части работы приведены рекомендации по использованию и дальнейшему развитию рассмотренных подходов моделирования динамики биомолекулярных систем, полученные на основании практического опыта при выполнении данной работы.
Выводы 1. Для интегрирования с методами молекулярной динамики реализованы процедуры расчёта градиентов энергии для метода потенциалов эффективных фрагментов в пакете молекулярного моделирования Q-Chem. Созданы оригинальные программные продукты, реализующие метод молекулярной динамики с жёсткими фрагментами и метод анализа взвешенных гистограмм.
2. Впервые с использованием полностью квантового описания проведено молекулярно динамическое моделирование реакции образования сукцинимида из остатков аспарагина в водном окружении. Показано, что при моделировании данной реакции чрезвычайно важным является учет растворителя, без которого не удается зафиксировать один из интермедиатов, а энергетический барьер реакции в газовой фазе вдвое превышает соответствующее значение в растворе. Установлено, что скорость-лимитирующими стадиями являются стадии депротонирования и циклизации.
3. С помощью теоретического моделирования процесса переноса протона в канале грамицидина А показано, что переориентация цепи из молекул воды является скорость определяющей стадией. Установлено, что стадии непосредственного переноса протона и переноса ионного дефекта являются безбарьерными.
4. Показано, что гидролиз аденозинтрифосфата в активном центре миозина по двух-водному механизму характеризуется барьером, не превосходящим 11 ккал/моль. Установлено, что данный механизм является инициирующей стадией для выхода продуктов реакции через так называемый «чёрный ход» в миозине.
5. На примере химических превращений в активном центре ацетилхолинэстеразы показано, что можно выделить минимально допустимый состав квантовой подсистемы, который необходим для адекватного описания данной системы методом КМ/ММ. В данном случае КМ-подсистема должна включать все аминокислотные остатки каталитической триады – серин, гистидин и глютаминовую кислоту, а также субстрат.
Основные публикации по теме диссертации 1. Московский А.А., Калиман И.А., Акимов А.В., Конюхов С.С., Григоренко Б.Л., Немухин А.В. Реализация метода молекулярной динамики с жесткими фрагментами для моделирования химических реакций в биомолекулярных системах // Вестник Моск. Ун та. Сер. Химия. 2007. Т. 48. С. 219-222.
2. Калиман И.А., Московский А.А., Конюхов С.С., Немухин А.В. Моделирование протонного транспорта в канале грамицидина А // Вестник Моск. Ун-та. Сер. Химия.
2008. Т. 49. С. 291-294.
3. Nemukhin A.V., Kaliman I.A., Moskovsky A.A. Modeling negative ion defect migration through the gramicidin A channel // J. Mol. Mod. 2009. V. 15. P. 1009-1012.
4. Kaliman I.A., Grigorenko B.L., Shadrina M.S., Nemukhin A.V. Opening the Arg-Glu salt bridge in myosin: computational study // Phys. Chem. Chem. Phys. 2009. V. 11. P. 4804-4807.
5. Конюхов С.С., Артемов Н.Н., Калиман И.А., Купченко И.В., Немухин А.В. Диффузия наноавтомобилей на основе фуллерена на поверхности кристалла золота // Вестник Моск. Ун-та. Сер. Химия. 2010. Т. 51. С. 267-269.
6. Kaliman I.A., Nemukhin A.V., Varfolomeev S.D. Free Energy Barriers for the N-Terminal Asparagine to Succinimide Conversion: Quantum Molecular Dynamics Simulations for the Fully Solvated Model // J. Chem. Theory Comput. 2010. V. 6. P. 184-189.
7. Grigorenko B.L., Kaliman I.A., Nemukhin A.V. Minimum energy reaction profiles for ATP hydrolysis in myosin // J. Mol. Graph. Model. 2011. V. 31. P. 1-4.
8. Калиман И.А. Моделирование процесса переноса протона в канале грамицидина А методом КМ/ММ молекулярной динамики // Международная конференция студентов, аспирантов и молодых ученых «Ломоносов». Россия, Москва 2007. Материалы конференции. С. 446.
9. Калиман И.А. Применение метода КМ/ММ молекулярной динамики для построения профилей свободной энергии // Международная конференция студентов, аспирантов и молодых ученых «Ломоносов». Россия, Москва 2008. Материалы конференции. С. 622.
10. Калиман И.А., Московский А.А. Реализация метода расчёта профилей свободной энергии химических реакций на основе гибридных методов КМ/ММ молекулярной динамики // Интернет-конференция «Информационно-вычислительные технологии в науке - ИВТН».
Россия, Москва 2008. Материалы конференции. С. 41.
11. Калиман И.А. Моделирование реакции образования сукцинимида из остатков аспарагина – одной из важнейших реакций деградации белков // Международная конференция студентов, аспирантов и молодых ученых «Ломоносов». Россия, Москва 2009. Электронный ресурс.
ISBN 978-5-317-02774-2.
12. Kaliman I.A. Modeling of Proton Transport in the Gramicidin A Channel // Международная конференция студентов, аспирантов и молодых ученых «Ломоносов». Секция английского языка. Россия, Москва 2009. Материалы конференции. C. 3.
13. Калиман И.А., Московский А.А., Немухин А.В. Применение методов квантовой механики – молекулярной динамики для моделирования биологических систем // Шестнадцатая Международная Конференция «Математика. Компьютер. Образование». Россия, Пущино 2009. Материалы конференции. С. 259.
14. Калиман И.А., Немухин А.В. Моделирование реакции гидролиза АТФ в миозине и открытия солевого мостика Arg238-Glu459 // IX Ежегоднпя международная молодежная конференция ИБХФ РАН-ВУЗы. Россия, Москва 2009. Материалы конференции. С. 106-108.
15. Kaliman I.A., Nemukhin A.V., Varfolomeev S.D. Computer simulation of the first stage of asparagine isomerization reaction in aqueous solution // International Conference “Biocatalysis:
Fundamentals & Applications”. Russia, Arkhangelsk 2009. Book of Abstracts. P. 28-29.
16. Калиман И.А., Григоренко Б.Л., Немухин А.В. Выход продуктов гидролиза АТФ через «чёрный ход» в миозине: новые свидетельства в поддержку механизма // XXII симпозиум «Современная Химическая Физика». Россия, Туапсе 2010. Материалы конференции. С. 165.
17. Калиман И.А., Немухин А.В. Применение методов квантовой молекулярной динамики для моделирования биомолекулярных систем в водном окружении // XVIII Международная Конференция «Математика. Компьютер. Образование». Россия, Пущино 2011. Материалы конференции. С. 68.
18. Калиман И.А. Методические вопросы применения метода комбинированной квантовой и молекулярной механики на примере реакции гидролиза ацетилхолина ацетилхолинэстеразой // Международная конференция студентов, аспирантов и молодых ученых «Ломоносов».
Россия, Москва 2011. Электронный ресурс. ISBN 978-5-317-03634-8.
19. Shao Y., Fusti-Molnar L., Jung Y., Kussmann J., Ochsenfeld C., Brown S.T., Gilbert A.T.B., Slipchenko L.V., Levchenko S.V., O'Neill D.P., DiStasio R.A. Jr., Lochan R.C., Wang T., Beran G.J.O., Besley N.A., Herbert J.M., Lin C.Y., Van Voorhis T., Chien S.H., Sodt A., Steele R.P., Rassolov V.A., Maslen P.E., Korambath P.P., Adamson R.D., Austin B., Baker J., Byrd E.F.C., Dachsel H., Doerksen R.J., Dreuw A., Dunietz B.D., Dutoi A.D., Furlani T.R., Gwaltney S.R., Heyden A., Hirata S., Hsu C.-P., Kedziora G., Khaliullin R.Z., Klunzinger P., Lee A.M., Lee M.S., Liang W., Lotan I., Nair N., Peters B., Proynov E.I., Pieniazek P.A., Rhee Y.M., Ritchie J., Rosta E., Sherrill C.D., Simmonett A.C., Subotnik J.E., Woodcock H.L. III, Zhang W., Bell A.T., Chakraborty A.K., Chipman D.M., Keil F.J., Warshel A., Hehre W.J., Schaefer H.F. III, Kong J., Krylov A.I., Gill P.M.W., Head-Gordon M., Gan Z., Zhao Y., Schultz N.E., Truhlar D., Epifanovsky E., Oana M., Baer R., Brooks B.R., Casanova D., Chai J.-D., Cheng C.-L., Cramer C., Crittenden D., Ghysels A., Hawkins G., Hohenstein E.G., Kelley C., Kurlancheek W., Liotard D., Livshits E., Manohar P., Marenich A., Neuhauser D., Olson R., Rohrdanz M.A., Thanthiriwatte K.S., Thom A.J.W., Vanovschi V., Williams C.F., Wu Q., You Z.-Q., Sundstrom E., Parkhill J., Lawler K., Lambrecht D., Goldey M., Olivares-Amaya R., Bernard Y., Vogt L., Watson M., Liu J., Yeganeh S., Kaduk B., Vydrov O., Xu X., Kaliman I.A., Zhang C., Russ N., Zhang I.Y., Goddard W.A. III, Besley N., Ghysels A., Landau A., Wormit M., Dreuw A., Diedenhofen M., Klamt A., Lange A.W., Ghosh D., Kosenkov D., Zuev D., Deng J., Mao S.P., Su Y.C., Small D. Q-Chem, Version 4.0.
http://www.q-chem.com/BetaQChem4DemoRequest.html. Q-Chem Inc. 2011. Pittsburgh, PA.