Содержание к диссертации
Введение
1. Моделирование строения, спектров и динамики молекул в инертных матрицах
1.1. Строение, спектры и динамика кластеров металлов в матрицах инертных газов 14
1.2. Малые молекулы и межмолекулярные комплексы в матрицах инертных газов 34
1.3. Моделирование спектров молекул в матрицах методами смешанной квантово-классической молекулярной динамики 45
2. Развитие метода двухатомных фрагментов в молекулах (ДФМ) для моделирования свойств межмолекулярных кластеров с водородными связями 64
2.1. Модель ДФМ для потециалыюй поверхности молекулы воды 67
2.2. Развитие модели ДФМ, и строение кластеров (Н20) 71
2.3. Развитие модели ДФМ, и строение кластеров (HF) 78
2.4. Комбинированный метод КМ/ДФМ: Диссоциация молекулы HF в кластерах (HF) 88
3. Развитие комбинированных методов квантовой и молекулярной механики (КМ/ММ) для моделирования химических реакций в сложном молекулярном окружении 106
3.1. Моделирование реакции ОН" + С02 —> НС03" в водных кластерах 110
3.2. Реакция глютатиона с гидроксиметильным радикалом в водных кластерах 113
3.3. Моделирование свойств и динамики НХеОН в водных кластерах 132
3.4. Метод КМ/ММ с конформационно-подвижными эффективными фрагментами. 135
3.5. Моделирование протонного транспорта по водным проводам 141
3.6. Моделирование механизмов ферментативных реакций 146
Основные результаты и выводы 165
- Малые молекулы и межмолекулярные комплексы в матрицах инертных газов
- Развитие модели ДФМ, и строение кластеров (Н20)
- Реакция глютатиона с гидроксиметильным радикалом в водных кластерах
- Моделирование протонного транспорта по водным проводам
Введение к работе
Роль молекулярного моделирования в химии достаточно велика, несмотря на явный приоритет экспериментальных исследований в этой области естествознания. Наиболее значимы такие теоретические результаты, которые невозможно, крайне трудно или слишком дорого получить экспериментальными средствами. Традиционно к задачам моделирования относят определение строения отдельных молекул, молекулярных ассоциатов или фрагментов твердых тел, а также описание механизмов химических реакций на молекулярном уровне. Очевидно, что методы квантовой механики позволяют давать достаточно полное представление о структуре и реакциях молекул. Область молекулярного моделирования, основанная на феноменологическом представлении энергии молекулярных систем, т.е. молекулярная механика, также испытывает серьезное влияние квантовой теории, поскольку очень часто параметры используемых потенциальных функций проще рассчитать методами квантовой химии, нежели подбирать их по большой совокупности экспериментальных данных.
Моделирование структуры и динамики химических процессов начинается с построения поверхностей потенциальной энергии (ППЭ) системы. Для расчетов молекулярных превращений в газовой фазе, как правило, достаточно знания ППЭ в основном и возбужденных электронных состояниях молекулярной системы. Основной особенностью при анализе превращений в конденсированных средах является необходимость учета вляния окружения (инертной или реакционной матрицы, растворителя, молекулярных групп белка и т.д.) на свойства, и соответственно, на
ППЭ, центральной подсистемы. Подразделение всей молекулярной системы, моделирующей фрагмент конденсированной фазы, на центральную и периферийную части зависит от постановки задачи. В задачах моделирования спектров частиц, изолированных в матрицах, центральная часть - это молекула внедрения, а периферийная часть - атомы матрицы. В задачах моделирования растворов на микроуровне центральная часть - это частицы растворенного вещества, а периферийная часть — молекулы растворителя. В задачах моделирования химических реакций в белках разделение менее очевидно, но, как правило, к центральной подсистеме можно отнести фрагменты в непосредственной близости от реакционного центра, а к периферийной - молекулярные группы в пределах десятков ангстрем от центра. Термин «периферийная подсистема» не означает, что атомные или молекулярные группы из этой подсистемы не оказывают влияние на свойства центральной части, это влияние может быть весьма существенным, и, в частности, учет подобных эффектов составляет предмет исследований, описанных в данной работе.
Известны два основных подхода к анализу молекулярных систем в конденсированных средах ~- на основе континуальных моделей среды и на основе кластерных моделей. Работы, составившие предмет данной диссертации, преимущественно относятся к направлению кластерных моделей, в рамках которого определенное число частиц среды явно включается в рассматриваемую молекулярную систему. В дальнейшем часто будет использоваться название «кластер» для фрагмента конденсированной фазы, выделенного для теоретического анализа.
Методы построения ППЭ для молекулярного моделирования
можно подразделить на полуэмпирические и неэмпирические. В
первом случае потенциалы взаимодействия между частицами
записываются в аналитическом виде, и параметры выбранных
функций подбираются по согласованию результатов моделирования
с выборкой экспериментальных данных. Наиболее популярны
потенциалы, используемые в методах молекулярной механики
(ММ), классической молекулярной динамики или потенциалы,
применяемые в физике твердого тела. В рамках подобного
подхода было выполнено огромное число работ по моде лир ованиго
свойств атомных кластеров, как гомогенных, так и гетерогенных.
Прежде всего, это расчеты равновесных геометрических
конфигураций кластеров, как совокупность минимумов на
многомерных ППЭ. Следует иметь в виду, что для подобных систем
характерно наличие большого числа таких минимумов, и
соответственно, наличие множества структурных модификаций. Как
правило, локальные минимумы разделены невысокими
потенциальными барьерами, и даже при небольших температурах
могут наблюдаться переходы между структурами. Следовательно,
требуется перечислить достаточно большое число структур, что
составляет достаточно сложную задачу. Далее, для найденных точек
рассчитываются энергетические характеристики и делаются
прогнозы о термодинамической стабильности системы. В частности,
используются молекулярно-динамические (МД) расчеты траекторий
частиц в кластерах при заданной температуре. Анализ траекторий,
построение различных функций распределения,
автокорреляционных функций позволяет характеризовать кластер, проанализировать зависимость свойств от размера кластера.
Наиболее существенным ограничением моделирования с эмпирическими потенциалами является то обстоятельство, что, оставаясь в рамках этого подхода, вообще говоря, невозможно рассматривать химические реакции, так как подобные потенциалы не описывают изменений электронной структуры.
Потенциалы взаимодействия, построенные неэмпирическими методами квантовой химии, естественно, более универсальны и принципиально позволяют решать все проблемы строения и химических реакций. Основные ограничения здесь связаны с размерами выделенной молекулярной системы. По-видимому, можно утверждать, что для частиц, включающих до двух десятков атомов, современные квантово-химические методы могут рутинно применяться и приводить к очень надежным результатам. Речь может идти о нахождении координат стационарных точек на поверхностях потенциальной энергии основных электронных состояний кластеров (т.е. точек локальных минимумов и барьеров на путях перегруппировок), расчетах относительных энергий в этих точках, а также энергетических профилей путей химических реакций в системе, прогнозах колебательных и электронных спектров, анализе деталей распределения электронной плотности. Стоимость подобных расчетов достаточно высока даже на уровне Хартри-Фоковского приближения. Кроме того, для молекулярных кластеров зачастую требуется применение методов с учетом эффектов электронной корреляции. Тем не менее, для систем с числом атомов от десяти до ста неэмпирические расчеты также технически осуществимы, но в этом случае придется значительно поступиться уровнем приближения. При этом равновесные геометрические конфигурации низколежащих минимумов на
потенциальных поверхностях можно определить относительно хорошо, равно как и оценить термодинамическую стабильность кластеров. Хотя прогноз спектров и реакционной способности кластерных частиц будет выполнен с меньшей точностью, но качественные тенденции будут воспроизведены правильно.
Существенный прогресс в развитии методов моделирования свойств больших молекулярных систем (для которых целесообразно подразделение на центральную и периферийную части) связан с т.н. комбинированными методами квантовой механики — молекулярной механики (КМ/ММ), интенсивно развивающимися в настоящее время. Основная идея таких подходов — использовать квантовое описание для активной части подсистемы, которая считается наиболее важной, и учесть строение периферийной части большой системы и ее влияние на центральную область посредством эмпирических или полуэмпирических потенциалов. Во многих случаях размер центральной части может быть выбран в пределах двух-пяти десятков атомов, и неэмпирические методы квантовой химии могут обеспечить хорошее количественное описание и строения, и реакций в выделенной подсистеме. В молекулярно-механическую подсистему могут быть включены сотни и тысячи атомов. Хотя не все принципиальные вопросы теории КМ/ММ решены, этот подход все более активно применяется при моделировании процессов в биомолекулярных системах и в материаловедении.
Все совеременные методы молекулярного моделирования ориентированы на использование компьютеров и специальных пакетов программ. Работы, составляющие предмет данной диссертации, преимущественно выполнены с применением
квантовохимического пакета программ GAMESS (The General Atomic and Molecular Electronic Structure System) и его версии PC GAMESS, ориентированной на компьютеры архитектуры Intel. Небольшая часть расчетов выполнена с использованием комплекса Gaussian, Задачи молекулярно-динамического моделирования решались с помощью компьютерных программ, разработанных в данной работе. Полностью оригинальными являются компьютерные программы реализации метода КМ/ММ, осуществляющие связку комплексов программ квантовой химии PC GAMESS и молекулярной механики TINKER.
Диссертация посвящена развитию новых методов расчета поверхностей потенциальной энергии молекулярных кластеров, моделирующих фрагменты конденсированных фаз, а также компьютерным реализациям новых подходов и исследованию конкретных задач матричной изоляции, задач теории микрорастворимости, и задач моделирования в биохимии.
Диссертация состоит из трех частей: Часть 1 «Моделирование строения, спектров и динамики молекул в инертных матрицах»; Часть 2 «Развитие метода двухатомных фрагментов в молекулах (ДФМ) для моделирования свойств межмолекулярных кластеров с водородными связями»; Часть 3 «Развитие комбинированных методов квантовой и молекулярной механики для моделирования химических реакций в сложном молекулярном окружении».
Основные результаты исследований представлены в 45 публикациях [1-45] в рецензируемых отечественных и международных научных журналах. Литературный обзор не выделен в отдельную главу, общую для всей работы, но в каждой части приводятся сведения о работах, представленных в литературе,
имеющих отношение к соответствующему разделу диссертации. Ссылки на общие вопросы квантовой механики, квантовой химии, молекулярной спектроскопии, изложенные, в частности, в учебниках и монографиях [46-57] и активно использованные в ходе исследований, как правило, не будут приводиться явно в тексте, чтобы не перегружать изложение.
Малые молекулы и межмолекулярные комплексы в матрицах инертных газов
Полностью неэмпирический анализ электронной задачи, следующий за молекулярно-динамическим моделированием, был использован при рассмотрении свойств молекулы NBr в аргоновых кластерах [4]. В данной работе ставилась цель интерпретировать данные спектроскопических экспериментов в условиях низкотемпературной матричной изоляции [96,97] и сравнить результаты различных теоретических подходов [4,65,98,99]. Для представления потенциальной поверхности использованы комбинации функций Леннард-Джонса с параметрами, приведенными в работе [98]. В рамках модели с числом атомов аргона до 170 оказалось возможно практически количественно воспроизвести индуцированный матрицей синий сдвиг (около +9 см"1) для колебательной полосы молекулы NBr. Колебания в системе анализировались методом молекулярной динамики, частоты колебаний определялись положениями пиков в спектре мощности по формулам (1.10). На Рис. 1.7 показан график спектра мощности для кластера NBr/Arno. На основании координат кластера NBr/Arno был выделен малый фрагмент NBr/Агц, показанный на Рис. 1.8, для которого были выполнены расчеты электронной структуры в основном Х32" и возбужденном b!S+ состояниях проводились методами квантовой химии. С этой целью рассчитывались электронные распределения в кластере NBr@Arn методами HF и МР2 с базисом 6-31G с использованием преобразования к натуральным связевым орбиталям. Сделан вывод, что следствием взаимодействия с матрицей могут быть достаточно заметные сдвиги в орбитальных энергиях связывающих орбиталей и орбиталей неподеленных пар, вплоть до 300 см Определенный прогресс в моделировании матричных эффектов был достигнут при использовании приближения двухатомных фрагментов в молекулах (ДФМ) [100, 101] для описания потенциалов взаимодействия молекул внедрения с сольватными оболочками атомов благородных газов. Метод ДФМ можно рассматривать и как полуэмпирический вариант одной из фундаментальных конструкций в теории электронного строения молекул, а именно, теории валентных схем, и как обобщение схемы парных потенциальных функций. Более полное обсуждение приближения ДФМ приведено в Части II; в данном разделе рассматривается только приложение к системам, содержащим атомы благородных газов.
В первых работах данного направления [6,14] рассчитывались сдвиги частоты колебаний Av молекулы HF при ее ассоциации с аргоновыми кластерами определенного размера Arn (n=l + 12, 62). Основной интерес был связан с моделированием зависимости свойства (в данном случае, частоты колебаний) молекулы внедрения от размера сольватной оболочки, что обычно рассматривается в теории микрорастворимости. Экспериментальные результаты микроволновой спектроскопии известны для начального участка ряда (n=l -j- 4) [102], а также известны спектральные данные для молекулы HF в аргоновой матрице [103]. Моделирование проводилось следующим образом. На построенной для гетерокластера каждого размера потенциальной поверхности находились точки глобального и локальных минимумов с использованием алгоритмов метода Монте Карло. Стартуя от конфигурации глобального минимума, далее рассчитывались классические молекулярно-динамические траектории при малых значениях температуры, определялись автокорреляционные функции скорости, фурье-преобразованием которых находились функции мощности (1.10). Пик спектра мощности в области высоких частот соотносился с частотой колебаний фрагмента HF в кластере HF»Arn. Разность с частотой колебаний молекулы HF в свободном состоянии давала оценку Av для данного значения п. Поверхности потенциальной энергии систем HF Arn строились как суперпозиции парных потенциалов взаимодействия атомов аргона V(ArjArj) и всевозможных троек V(AriHF) Для описания взаимодействия Аг-Аг был использован общепринятый многопараметрический потенциал Азиза-Чена [104], калиброванный по большой совокупности экспериментальных данных для аргона. Потенциал взаимодействия в тройной системе АгНР был построен по принципам приближения состояния 2 молекулы HF комбинируется с энергиями других пар, причем рассматриваются вклады как нейтральных ArF, АгН, так и ионных ArF", АгН4" электронных конфигураций. Смешивание нейтральных и ионных вкладов задается через параметр р\ Смешивание электронных конфигураций 2- и П-типа двухатомного фрагмента ArF зависит от угла у между направлениями F-H и F-Ar. На Рис. 1.9 приведена схема, иллюстрирующая детали построения поверхности потенциальной энергии трехатомной системы ArHF в основном синглетном состоянии А в рамках выбранной модели приближения ДФВМ. В качестве базиса выбраны три функции, комбинирующие атомные состояния атомов и ионов. Необходимые кривые потенциальной энергии двухатомных фрагментов были либо аппроксимированы по экспериментальным неэмпирическими методами квантовой химии (АгН, 2; АгН , 2). На Рис Л. 10 показаны потенциальные кривые нейтрального и ион-парного состояний молекулы HF, учет взаимного влияния которых, в существенной степени определяет успех построения поверхности потенциальной энергии ArHF.
Параметр смешивания нейтральных и ионных состояний р был выбран таким образом, чтобы для трехатомной системы ArHF, описываемой потенциальной поверхностью приближения ДФМ, правильно воспроизводилась энергия связи ( 200 см"1) по отношению к Ar + HF [105], а также сдвиг частоты колебаний (-10 см 1) молекулы HF при ее ассоциации с одним атомов аргона [102]. Коэффициент Р являлся единственным подгоночным параметром модели. На Рис. 1.11 приведены результаты расчетов Av, полученные в нашей модели (кривая 2), а также литературные данные квантовых расчетов с использованием известной эмпирической потенциальной поверхности Хатсона [106] для ArHF. Для некоторых кластеров показаны равновесные геометрические конфигурации основных изомеров. Результаты данной работы позволили сделать вывод, что потенциальные поверхности на основе метода ДФМ достаточно надежны при моделировании сольватных сдвигов в кластерах благородных газов. В наших расчетах получена такая же немонотонная зависимость Av от п, как и в более сложных подходах [107]. Важным результатом работы является то, что вычисленный сольватный сдвиг в кластере HF Are2 (-40 см"1) практически совпадает с экспериментальным сдвигом в аргоновой матрице (-41 см 1) [103]. Анализ структур кластеров HF Arn позволяет объяснить механизм роста сольватных оболочек: формирование HF Arn сходно с формированием кластеров Агп, и уже при п=12 молекула HF полностью входит в объем. Успехи моделирования структур и спектров кластеров HF Arn стимулировал работы по расчетам свойств димера фторида водорода (HF)2 в аргоновом окружении [5,10,11]. Строение и свойства простейшего комплекса с водородной связью, (HF)2, всегда были объектом пристального внимания и экспериментаторов, и теоретиков, При помещении этого комплекса в окружение атомов инертного газа взаимодействие солъватных оболочек с межмолекулярным комплексом хотя и меньше, но уже сопоставимо с взаимодействием мономерных единиц HF...HF в комплексе. Для (HF)2 известны результаты колебательной спектроскопии в газовой фазе [108], в объеме аргоновой матрицы [109], а также на поверхности аргоновых нанокластеров [ПО, 111]. Наиболее интересные и важные результаты моделирования относятся именно к различиям в спектрах комплекса в разном аргоновом окружении. Равновесная геометрическая конфигурация (HF)2 изображена на Рис.1.12. Указанные параметры отвечают точке минимума на эмпирической 6-мерной потенциальной поверхности V((HF)2), тщательно подобранной по большой совокупности экспериментальных данных для газофазной системы [112].
Развитие модели ДФМ, и строение кластеров (Н20)
Построение поверхностей потенциальной энергии ассоциатов полярных молекул HF, Н20 на основе приближения ДФМ потребовало существенного развития этой методологии [19,25,26,28]. Основное внимание было уделено расчетам энергии межмолекулярного взаимодействия для произвольных расположений в пространстве мономерных единиц. Для решения поставленной задачи применялась теория возмущений на основе метода ДФМ [70] в сочетании с модификацией методики ДФМ, учитывающей вклады от ионных электронных конфигураций мономеров. В приложении к кластерам воды (Н20)п предложенная схема сводится к следующему. Обычное представление матрицы гамильтониана переписывается таким образом чтобы выделить мономерные и димерные единицы и соответственно, внутримолекулярные и межмолекулярные вклады Если целью расчета является межмолекулярное взаимодействие, то задача сводится к выражению матричных элементов Н через потенциальные кривые пар атомных частиц, каждая из которых принадлежит разным мономерам выделенного димера. Матрица Н- е рассматривается как матрица оператора возмущения в базисе произведений функций мономеров Соответственно, Q H[ F) оценивает сумму энергий мономеров. После спецификации всех матриц преобразований R, Т, В (см. формулу (2.2)) оказывается возможным вывести явные выражения для каждого вклада в энергию взаимодействия для димера (ij), VyD1M). Эти вклады зависят от потенциальных кривых VaKb пары атомных частиц в определенном электронном состоянии к как функций расстояний между ядрами arty, а также от угловых переменных, определяемых геометрией всего кластера (Н20)п. Как принято в подходе ДФМ, относительные веса нейтральных и ионных электронных конфигураций в волновых функциях валентных схем мономерных молекул выражаются через параметр смешивания р (2.8) Значение параметра p можно выбрать по результатам квантовохимических расчетов молекулы Н20. Формула для Уу01М) разбивается на взвешенную сумму нейтральных-нейтральных (NN), нейтральных-ионных (N1) и ионных-ионных (II) составляющих (N относится к частицам О, Н; I относится к частицам О", ЬҐ) Весовые коэффициенты w t(p) определяются на основе выкладок с формулами (2.7, 2.8, 2.9).
В свою очередь, каждая составляющая Vy +) содержит вклады от перекрестных пар атомных частиц в димере Окончательное выражение для энергии взаимодействия п мономерньтх молекул воды получается при добавлении к / \ Vjj(DIM) электростатических слагаемых, зависящих от зарядов q (Н"1", О ), взимодействующих с нейтральными атомными частицами (Н, О) с поляризуемостями а, причем необходимо учесть масштабирующие факторы f((3), связанные с параметром смешивания; Как и в большинстве прежних приложений метода ДФМ, в данной работе часть двухатомных потенциалов хорошо известны (Н2, Н2+, некоторые состояния ОН), однако, возникала потребность рассчитать достаточно большое число потенциальных кривых двухатомных систем, недостаточно хорошо описанных в литературе. В частности, были необходимы энергетические кривые дублетных и квартетных состояний аниона 0{. Соответствующие кривые приведены в Приложении. После того, как все необходимые потенциальные кривые двухатомных фрагментов аппроксимированы подходящими аналитическими функциями и определены матричные элементы матриц вращения R и спиновой связи Т, методика позволяет рассчитывать потенциальные поверхности кластеров (Н20)п. Были рассмотрены молекулярные ассоциаты, состоящие из 2 -6 молекул воды. Основное внимание уделялось характеристикам систем в конфигурациях глобальных минимумов и низколежащих локальных минимумов (изомерных структур). Результаты сопоставлялись с литературными данными. Для малых кластеров воды известны результаты моделирования их структуры и свойств, выполненных в разных приближениях - от неэмпирической квантовой химии высокого уровня точности до молекулярно-механических моделей- Не всегда эти результаты согласуются между собой. Некоторые молекулярные постоянные кластеров (НгО) известны из экспериментальных спектральных исследований. На Рис.2.2 изображены некоторые структуры димера (а), тримера (b, с), тетрамера (d, е), пентамера (f) и гексамера (g, h, і, j). Наша модель ДФМ позволяет качественно правильно воспроизвести сетку водородных связей для всей совокупности кластеров и описать конфигурации основных изомеров. Определенное внимание уделялось анализу изомерных структур кластеров (Н20)п (п 2). Было показано, что использованная методика позволяет не только найти геометрические конфигурации низкоэнергетических изомеров, но правильно воспроизвести порядок следования изомеров на энергетической шкале и разности энергий между этими структурами. В таблицах 2.4 - 2.7 сопоставлены результаты нашего расчета и литературные данные для энергий нескольких изомеров для (Н20)п.
Результаты нашей модели ДФМ хорошо согласуются с результатами эталонных расчетов. Работы [19,25-27,29] были посвящены развитию метода ДФМ в приложении к моделированию свойств кластеров фторида водорода (HF)n. Формулы (2.3)-(2.7) общей теории ДФМ и ее модификации [28] для моделирования свойств кластеров воды при соответствующей коррекции относятся и к описанию комплексов (HF)n. В частности, волновая функция мономерной молекулы HF должна описываться уравнением Aj/ =(2P(F) 2S(H) cosP + IS(F) IS(H1 sinp (2.12) вместо (2.8). Хотя в случае фторида водорода аппарат теории несколько проще, чем при описании кластеров воды, но здесь также потребовалась большая аналитическая работа, прежде чем были получены рабочие формулы для вычислений потенциальных поверхностей (HF)n. Параметры линейных преобразований между спиновыми функциями определяются с помощью коэффициентов Клебша-Гордана [174], выписанных в формулах (2.22) — (2.26). В окончательном варианте, синглетные и триплетные состояния пространственной симметрии Л любой композиции FH, F2, Н2 входят в элементы двухатомных матриц в линейных комбинациях с весами 3/4 или 1Л, как показано в уравнении (2.27) Кластеры (HF)n часто рассматриваются как идеальные кандидаты для проверки методов конструирования поверхностей потенциальной энергии межмолекулярных комплексов с водородными связями. Для этих систем известны данные тщательных экспериментальных исследований спектральными методами высокого разрешения [175]. В частности, по этим данным известно, что основным изомерам кластеров (HF)n (по крайней мере до п 6) отвечают циклические структуры. Сравнение результатов нашего моделирования будет проводиться с экспериментальными структурами для циклических изомеров и с данными неэмпирических расчетов методами квантовой химии МР2/6-311+G(2d ,2р) для нециклических изомеров. Единственным подгоночным параметром модели ДФМ является параметр смешивания нейтральных и ионных состояний HF (3=38 в уравнении (2.12)), выбранный по значению энергии связи циклического кластера (HF)3.
Реакция глютатиона с гидроксиметильным радикалом в водных кластерах
Интерес к свойствам и реакциям глютатиона (у-глютамилцистеинилглицина) объясняется важностью этого биологического тиола как основного клеточного антиоксиданта [216]. Вследствие наличия высокореакционной сульфгидрильной группы SH от аминокислотного остатка цисте ина глютатион GSH легко вступает в реакции с углеродсодержащими свободными радикалами и превращается в глютатиильный радикал GS : Скорости радикальных реакций глютатиона с различными реагентами в воде могут быть достаточно высокими с небольшими величинами энергии активации [216]. При моделировании подобных реакций методами квантовой химии обычно рассматриваются существенно упрощенные молекулярные системы. К теме настоящей работы наиболее близки исследования Раука с сотр. [217-219], в которых рассчитывались реакции отщепления атомов водорода от тиолов углеродсодержащими радикалами и тиильными радикалами. Методами теории функционала электронной плотности рассчитывались энергии диссоциации связей, структуры переходных состояний и энергии активации для систем в вакууме, и затем эффекты растворителя учитывались в рамках континуальных моделей среды. В частности, для реакции гидроксиметильного радикала с метилтиолом, который часто рассматривается как прообраз глютатиона, по результатам расчетов методом B3LYP/6-311+G(d,p) была получена оценка (с учетом эмпирических поправок к вычисленным энергиям диссоциации реагентов и продуктов реакции) для энергии активации в газовой фазе Ea(g) = 9,4 кДж/моль. Применения различных вариантов континуальной модели приводили к величинам Ea(aq) в водном растворе, меняющимся от 10,0 до 19,3 кДж/моль [219]. Измерения констант скоростей реакции (3.7) в воде методом импульсного радиолиза позволили экспериментально оценить энергию активации Ea(aq) = 12,1 кДж/моль [219]. Мы рассматривали реакцию глютатиона с гидроксиметильным радикалом в водной среде в рамках супермолекулярной модели. Последнее означает, что все молекулярные группы трипептида, а также молекулы воды из гидратной оболочки явно включаются в систему. При физиологических значениях рН ионогенные группы должны быть ионизированы, поэтому рассматривается цвиттер-ионная форма пептида.
Основной целью моделирования является оценка энергии активации реакции (3.8) в воде, а также анализ влияния факторов среды, т.е. окружающих непосредственный реакционный центр пептидных групп, а также гидратной оболочки, на энергетический профиль реакции. Литературные результаты для реакции гидроксиметильного радикала с метилтиолом (3.7) рассматриваются как реперные данные. Для расчетов необходимых сечений потенциальной поверхности мы применили комбинированные методы квантовой и молекулярной механики (КМ/ММ) [189]. В этом подходе предполагается, что характеристики наиболее важной части молекулярной системы рассчитываются средствами квантовой химии, а периферийная часть системы описывается молекулярно-механическими силовыми полями. Используемая нами, оригинальная версия теории КМ/ММ основана на методе потенциалов эффективных фрагментов (ПЭФ) [190, 220, 221], ранее реализованном в квантово-химической программе GAMESS(US) [188]. В исходной формулировке метода ПЭФ, предназначенного для моделирования явлений сольватации, вся система подразделяется на квантовую часть, т.е. частицы растворенного вещества, и молекулы растворителя, рассматриваемые как эффективные фрагменты с фиксированными внутренними координатами. Влияние молекул растворителя на электронные состояния растворенного вещества и сответствующий отклик от квантовой подсистемы учитываются с помощью одноэлектронных потенциалов эффективных фрагментов, представленных в общем случае суперпозицией электростатических, поляризационных и обменных вкладов. Параметры ПЭФ могут быть определены в серии предварительных неэмпирических квантово-химических расчетов. В частности, для молекул воды эти параметры были тщательно подобраны ранее и включены в базу данных программы GAMESS(US). В предшествующих работах было показано, что использование метода ПЭФ позволяет успешно описать строение как малых кластеров воды, состоящих из 3-5 молекул [222], так и частиц (Н20)п среднего размера (п 20) [223]. Показано, что энергии связи в кластерах воды в приближении ПЭФ согласуются с результатами расчетов неэмпирическими методами квантовой химии в пределах 5-10 кДж/моль.
Особенно важна способность комбинированного метода КМ/ММ(ПЭФ) правильно описывать сетки водородных связей с участием молекул воды, относимых к ММ-подсистеме, и частиц растворенного вещества, составляющих КМ-подсистему. Это было показано на примере расчетов барьеров внутреннего вращения в молекуле формамида в водных кластерах [224], а также при сопоставлении структур нейтральной и цвиттер-ионной форм глютаминовой кислоты [225]. В последней работе [225] было найдено, что уже в кластере 10 молекул воды цвиттер-ионная форма становится более стабильной. Модель КМ/ММ(ПЭФ) была использована для построения потенциальных поверхностей возбужденных электронных состояний формамида в комплексе с тремя молекулами воды [224]. Структуры и колебательные спектры смешанных кластеров NaCl в окружении молекул воды были исследованы в работе [226]. Метод ПЭФ в исходной формулировке также применялся для анализа химических преобразований в пептидном окружении [227-234], однако, в этих приложениях жесткость эффективных фрагментов при оптимизации геометрических параметров системы представляла серьезное ограничение для моделирования. В работах [35,36,39] нами был предложен вариант метода КМ/ММ, снимающий это ограничение за счет того, что периферийная подсистема представляется совокупностью подвижных цепей малых эффективных фрагментов, взаимодействия между которыми описываются молекулярно-механическими силовыми полями. Концепция эффективных фрагментов также была использована для решения проблемы границы между КМ- и ММ-частями. Именно данная версия метода КМ/ММ применяется в данной работе для описания передачи атома водорода от глютатиона к гидроксиметильному радикалу при учете и перестроек в гидратнои оболочке, и изменений в конформациях трипетида в ходе химической реакции. При расчетах энергии для реакции (3.8) атомы гидроксиметильного радикала и сульфгидрильной группы GSH с примыкающим буферным фрагментом отнесены к КМ-подсистеме. На Рис.3.3 показано, каким образом можно осуществить разбиение молекулы глютатиона на совокупность эффективных фрагментов. Каждый такой фрагмент перемещается как единое целое в результате взаимодействий с квантовой подсистемой и с остальными эффективными фрагментами (включая и молекулы воды). Т.н. буферный фрагмент (здесь, -СН2-) на границе КМ- и ММ-областей относится к обеим подсистемам. Мы используем обычный для подходов КМ/ММ прием, дополняя при разрыве С-С связи квантовую область связующим ("link") атомом водорода.
Моделирование протонного транспорта по водным проводам
В другом приложении [38] нового метода КМ/ММ было выполнено моделирование реакций передачи протонов по ориентированным сетям водородных связей в белковом окружении. На Рис.3.13 показана одна из таких систем, описывающая участки канала мембранного белка М2 вируса гриппа. Передача протона осуществляется через последовательность молекул воды в полости канала (водные проводники), и очевидным препятствием на пути до концевой группы Asp является наличие аминокислотного остатка гистидина. Исходные координаты атомов модельной системы Рис.3.13 были получены на основании экспериментальных данных ЯМР и циркулярного дихроизма для белка М2, состоящего из четырех пептидных спиралей. Для наших расчетов была взята одна из таких спиралей, составленная из остатков глицина, общей длины 27А. Большинство частиц спирали были фиксированы в пространстве для симуляции стенки канала. Функциональные группы His и Asp, и ближайшие к ним 4 молекулы воды составляли квантовую подсистему. Буферные фрагменты и ближайшие к ним ЭФ из ММ-подсистемы были подвижны. Оптимизация геометрических параметров проводилась методом КМ/ММ в приближении RHF/6-31G/OPLSAA. Далее к системе был добавлен протон на крайнюю правую (Рис.3.13) молекулу воды, и проводилась минимизация энергии. Следуя координатам системы вдоль спуска по энергии, можно составить представление о движении протона. Весь маршрут подразделен на первую часть (до N5 остатка His) и вторую (от Ns His до Asp). Результаты моделирования методом КМ/ММ [14] показали, что вдоль ориентированных цепей молекул воды, прикрепленных к молекулярным стенкам, протон перемещается без потенциального барьера. Качественная картина перемещений протона вдоль молекул воды (с наблюдением промежуточных структур Н30+) похожа на картину Рис.3.8. С количественной точки зрения важно отметить, что рельеф поверхности потенциальной энергии оказывается достаточно пологим, и реориентация молекул воды позволяет избежать затрат энергии на перегруппировку химических связей. Единственный барьер на пути протона до концевой группы Asp возникает от необходимости разрывать связь Ne-H в His на второй части маршрута. Было показано, что в прицессе движения протона вдоль стенки спирали происходят достаточно заметные конформационные изменения, иллюстрированные на Рис.3.14. В частности, расстояние между атомами Су остатков His и Asp изменяется в широких пределах. В качестве одного из первых приложений новой методики КМ/ММ с конформационно-подвижными эффективными фрагментами был выбран каталитический цикл реакции гидролиза пептидных связей сериновыми протеазами [243].
Предшествующие теоретические подходы к моделированию отдельных участков данного цикла использовали неэмпирические и полуэмпирические методы квантовой химии [244-256], но полная картина превращений не была описана до нашей публикации [42]. Особенностью этих ферментов является наличие каталитической триады из аминокислотных остатков серина, гистидина и аспарагиновой кислоты. Следуя номенклатуре химотрипсина, в данной работе используются обозначения: Serl95, His57, Asp 102. Существенное значение имеют также группы NH от остатков Glyl93 и Serl95, способствующие образованию т.н. оксианионной дыры, позволяющей стабилизировать тетраэдрические интермедиаты и на стадии ацилирования, и на стадии деацилирования. На пути реакции гидролиза молекулярная система проходит через фермент-субстратный комплекс (ES), первый тетраэдрический интермедиат (INT1), фермент-ацильный комплекс (ЕА), второй тетраэдрический интермедиат (INT2) [243]. Каждой структуре отвечает область локального минимума на потенциальной поверхности. На Рис.3.16 пояснена схема разделения системы трипсин-субстрат на КМ- и ММ-части. В качестве начальных координат атомов был взят комплекс трипсин-ингибитор BPTI (код 2РТС из В квантовую область (Рис.3.16) были включены фрагменты каталитической триады и центральная часть модельного субстрата около разрываемой связи. ММ подсистема включала примерно 750 атомов в пределах 18 А от кислородного центра серина каталитической триады. Расчеты в квантовой области при локализации точек минимумов на потенциальной поверхности были проведены в приближении ССП с базисом 6-31G. Параметры молекулярно-механической части соответствовали библиотеке AMBER [258]. Энергии в найденных стационарных точках пересчитывались в квантово-химических приближениях более высокого уровня, т.е. по теории возмущений с расширенными базисами. Аминокислотный остаток His57 на двух стадиях действует как основание, которое отрывает протон от Serl95 на стадии ацилирования и от гидролитической воды на стадии деацилирования. Соответственно, протонированное имидазольное кольцо (His57) является донором протона, отдавая его субстрату на первой стадии и возвращая протон аминокислотному остатку Serl95 на второй. В качестве илюстрации на Рис.3.19(а,б) изображены структуры, относящиеся к лимитирующей стадии реакции, т.е. переходу от фермент-субстратного комплекса (ES) к первому тетраэдрическому интермедиату (INT1).
Показаны как химические (верхние части рисунков), так и молекулярные (нижние части) структуры. В таблице 3.3 приведены вычисленные в приближении МР2/6-31+G //.RHF/6-31G в квантовой части и при использовании силовых полей AMBER в молекулярно-механической части разности энергий между конфигурацией, отвечающей барьеру (В) при переходе ES — INTI, и исходной точкой (ES), а также между INT1 и В. Результаты для исходной системы (вторая строка) показывают, что энергетический эффект данной стадии оценивается как бликий к нулю, в то время как энергетический барьер составляет примерно 13 ккал/моль. В таблице 3,13 также указаны разности энергий, полученные для «мутированных» систем, в которых формальные заряды на остатках аргинина и лизина субстрата или на остатке аспарагиновой кислоты фермента полагались равными нулю. Эти результаты демонстрируют важность электростатических вкладов в энергетику ферментативных реакций. Таблица 3.3. Разности энергий (ккал/моль) между конфигурацией, отвечающей барьеру (В) при переходе ES - INT1, и исходной точкой (ES), а также между INT1 и В, вычисленные в приближении MP2/6-31+G //RHF/6-31G в квантовой части и при использовании силовых полей AMBER в молекулярно-механической части. Приведены данные для исходной системы и для «мутированных» систем, в которых формальные заряды на остатках аргинина и лизина субстрата или на остатке аспарагиновой кислоты фермента полагались равными нулю. Гидролиз гуанозинтрифосфата (GTP) белковым комплексом RAS-GAP Важнейшая ферментативная реакция гидролиза гуанозинтрифосфата (GTP), приводящая к гуанозиндифосфату (GDP) и неорганическому фосфату Pj является лимитирующей стадией всего цикла связывания/гидролиза GTP, ответственного за передачу сигналов на рост клеток. GTP/GDP-связывающие белки (G-белки) переходят из активной формы ("ON") в неактивную ("OFF") в результате реакции (3.11) [259,260]. Замена GDP на GTP восстанавливает активную форму. Эффективно гидролиз проходит в том случае, когда G-белок, содержащий GTP, образует комплекс с активирующим белком GAP. Сбой в работе этой молекулярной машины, в частности, замедление или прекращение реакции (3.11) при определенных мутациях одного из G-белков человека, p21ras (далее для краткости, RAS), имеет следствием развитие раковых опухолей. Экспериментальному и теоретическому изучению реакции гидролиза GTP посвящено огромное число публикаций [261-278], и тем не менее, четких и непротиворечивых ответов на вопросы о механизме реакции (3.11) в различных средах и о роли ключевых аминокислотных остатков не дано. Доминирующей точкой зрения является гипотеза ассоциативного механизма реакции гидролиза по схеме субстрат-ассистирующего катализа [265], однако, все попытки рассчитать энергетический профиль реакции по этому механизму в рамках различных молекулярных моделей приводят к нереально высоким активационньтм барьерам более 40 ккал/моль.
Похожие диссертации на Развитие методов моделирования молекулярных явлений в конденсированных средах
