Кристаллография, 2022, T. 67, № 3, стр. 411-418

ВВЕДЕНИЕ

В термодинамике свободной энергией называют количество внутренней энергии системы, которая может быть использована для выполнения работы. Изменение свободной энергии (ΔG) определяет направленность термодинамических процессов и управляет такими молекулярными процессами, как свертывание белков, молекулярные взаимодействия и химические превращения [1]. Таким образом, определение ΔG является одной из важных задач в биомолекулярных исследованиях белков и их комплексов с другими макромолекулами или низкомолекулярными лигандами (кофакторами, субстратами, ингибиторами). Поскольку экспериментальные измерения термодинамических свойств биомолекулярных систем часто являются дорогостоящими и трудоемкими, точные теоретические расчеты ΔG с помощью методов вычислительной биологии становятся все более востребованными при изучении макромолекулярных комплексов и рациональном дизайне лекарств. В настоящее время наиболее широко используемым инструментом для определения ΔG молекулярных систем являются подходы, основанные на моделировании методами молекулярной динамики (МД), которая позволяет не только моделировать физические процессы, происходящие в системе на атомарном уровне, но и выявлять термодинамические характеристики системы, определяющие ее устойчивость в различных состояниях. Наиболее популярными при вычислении ΔG являются метод, сочетающий использование молекулярной механики и решение уравнения Пуассона–Больцмана (ММ-PBSA), и метод, сочетающий использование молекулярной механики и решение обобщенного уравнения Борна (ММ-GBSA), поскольку они более точны, чем большинство скоринговых функций молекулярного докинга, но при этом не требуют значительных вычислительных мощностей [2].

Олигопептидазы В (OpB), представляющие собой сериновые пептидазы и обладающие трипсиноподобной субстратной специфичностью (КФ 3.4.21.83), найдены в бактериях, паразитических простейших и некоторых растениях [3]. OpB относятся к подсемейству S9А семейства пролил-олигопептидаз (РОР) [4]. Все представители семейства РОР состоят из С-концевого каталитического α/β-гидролазного домена и N-концевого β-пропеллерного регуляторного домена, ограничивающего доступ высокомолекулярных субстратов к каталитической триаде. Два домена соединены шарнирным участком, что позволяет ферментам переходить из каталитически неактивной (открытой) конформации, в которой домены и аминокислотные остатки каталитической триады разобщены, в каталитически активную (закрытую) конформацию, в которой сами домены и остатки каталитической триады сближены [5].

OpВ являются известными факторами патогенеза тяжелых паразитарных заболеваний и некоторых бактериальных инфекций [3, 6]. Данные пептидазы также участвуют в защите бактерий от антимикробных пептидов [7]. Эти факты, а также отсутствие ОрВ у млекопитающих делает эти ферменты перспективными мишенями для разработки противопаразитарных и антибактериальных препаратов. Моделирование низкомолекулярных ингибиторов требует знания пространственных структур ферментов-мишеней, имеющих каталитически активные конформации, и знания механизмов каталитической активации.

В пространственных структурах бактериальных ОpВ, полученных на сегодня, фермент находился в переходной конформации, в которой, несмотря на сближение доменов, каталитическая триада была разобщена [8, 9]. Однако структуры закрытых конформаций были получены для ферментов простейших, а также для бактериальных пролилэндопептидаз (РЕР), относящихся к тому же подсемейству S9А, что и ОрВ; оба типа ферментов были закристаллизованы в комплексах с ковалентно связанными ингибиторами – аналогами переходного состояния [5, 10–12]. Знание пространственных структур OpB из лейшмании Leishmania major (PDBID 2XE4, [10]) и трипаносомы Trypanosoma brucei (PDBID 4ВР9, [5]) в комплексах с антипаином, а также бактериальных РЕР из Myxococcus xanthus (PDBID 2BKL, [11]) и Aeromonas punctata (PDBID 3IVM, [12]) в комплексах с бензилоксикарбонил-Р/A-пролиналом позволяет моделировать каталитически активные конформации ОрВ бактерий. Однако механизм каталитической активации, предложенный в [5] для ОрВ простейших и бактериальных РЕР, не может быть автоматически перенесен на ОрВ бактерий, так как в последних необходимые для реализации данного механизма аминокислотные остатки границы между доменами не являются консервативными [13]. Ранее, комбинируя методы молекулярного моделирования и динамики с направленным мутагенезом, проведенным для подтверждения результатов вычислительных методов, были идентифицированы альтернативные аминокислотные детерминанты каталитической активации, находящиеся на границе между доменами ОрВ из бактерии Serratia proteamaculans (PSP), взятой в качестве модельного бактериального фермента [13, 14].

В настоящей работе с целью изучения прогностической значимости термодинамических расчетов методом MM-GBSA и их применимости для исследования механизма каталитической активации бактериальных ОрВ провели сравнительные исследования модельных фермент-субстратных комплексов, состоящих из PSP дикого типа или мутантных форм данного фермента с измененной активностью и вторичной субстратной специфичностью и тетрапептидов GRRG и GKRG. Рассчитанная методом MM-GBSA свободная энергия связывания пептидных субстратов была сопоставлена с эффективностью катализа, экспериментально измеренной для каждой фермент-субстратной пары. Выявлено, что увеличение ΔG связывания пептида в активном центре олигопептидаз, как правило, коррелирует с повышенной эффективностью катализа, что указывает на применимость метода MM-GBSA для прогнозирования каталитической эффективности бактериальных ОрВ в отношении различных субстратов, а также изучения влияния аминокислотных замен на каталитическую активность.

МАТЕРИАЛЫ И МЕТОДЫ

Построение начальных моделей комплекса. Модель PSP с закрытой конформацией построена как описано в [13], модели комплексов PSP с тетрапептидами GRRG и GKRG – как описано в [15]. Модели комплексов мутантных форм PSP получены с помощью программы PyMol, как описано в [12].

МД-моделирование. Полученные модели олигопептидаз в комплексах с тетрапептидами GRRG и GKRG использованы как стартовые для проведения МД-экспериментов. Для обеих систем моделирование проводили с помощью программы GROMACS 2020.3 [16]. В качестве силового выбрано поле ff14SB с улучшенной точностью моделирования боковых цепей аминокислотных остатков [17]. В качестве модели воды выбрана модель TIP3P, как наиболее подходящая для использования с силовым полем ff14SB. Молекула белка была помещена в ячейку в форме параллелепипеда таким образом, чтобы расстояние от молекулы белка до краев ячейки составляло 1.5 нм. Ячейка была заполнена молекулами воды. Кроме того, в каждую из систем было добавлено 0.15 мМ КCl, а также ионы K⁺ для уравновешивания заряда системы. На первом этапе моделирования для каждой из систем проведена минимизация энергии с шагом 1 фс до максимальной силы 1000.0 кДж/моль^–1 нм^–1. Далее системы были уравновешены при температуре 298 K и давлении 1 атм путем моделирования в NVT и NTP (100 пс в каждом) соответственно. Температуру и давление в системах контролировали с использованием модифицированного термостата Берендсена [18] и баростата Паринелло–Рахмана [19] с временными константами tau_t = = 0.1 и tau_p = 2 пс соответственно. Продуктивное 50 нс МД-моделирование для каждой системы проводили в изотермоизобарическом ансамбле с шагом в 2 фс. Алгоритм LINCS [20] использовали для ограничения водородных связей. Дальние электростатические взаимодействия рассчитывали с использованием схемы суммирования по Эвальду [21]. Ван-дер-Ваальсовы и кулоновские взаимодействия были усечены до 1.4 нм, что является оптимальным для используемого силового поля [17]. Поскольку при моделировании выполнялись периодические граничные условия, были выполнены перецентрирование белковых молекул и их возврат в ячейку моделирования с использованием команды gmxtrjconv программного пакета GROMACS с флагом-pbcnojump.

Расчет энергии взаимодействия тетрапептидных субстратов с ферментами. Для расчета ΔG взаимодействия тетрапептидов GRRG и GKRG с модельными ферментами использовали метод MM-GBSA [22]. Расчет осуществляли с использованием программы gmx_MMPBSA версия v1.4.0 (http://doi.org/10.5281/zenodo.4569307), которая с помощью программного пакета Amber Tools 2020 и программы MMPBSA.py [23] позволяет применять метод MM-GBSA для траекторий, полученных с использованием программного пакета GROMACS. Для расчета ΔG использовали 5000 структур для каждой из смоделированных систем. При вычислениях использовали модифицированную обобщенную модель Борна [24].

ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

Ранее были выявлены заряженные аминокислотные остатки PSP, участвующие в формировании границ между доменами. Замены этих остатков на аланин моделировали активность и вторичную субстратную специфичность фермента [13]. Данные о каталитической эффективности PSP дикого типа и ее мутантных вариантов в отношении п-нитроанилидных субстратов, содержащих остатки аргинина в P1-положениях и остатки аргинина или лизина в Р2-положениях, позволяли исследовать прогностический потенциал метода MM-GBSA, основанного на использовании МД-эксперимента, для определения изменения свободной энергии [2] и его применимость для исследования механизмов каталитической активации и поиска ингибиторов бактериальных ОрВ.

Метод MM-GBSA позволяет сравнить свободную энергию сольватированной молекулы в двух разных состояниях (связанном и свободном). Величина ΔG при связывании лиганда белком-акцептором может быть рассчитана по формуле

где ΔG_вак – изменение свободной энергии при переходе от несвязанного к связанному состоянию молекул в вакууме, ΔG_ком – изменение свободной энергии при переходе комплекса из вакуума в раствор, ΔG_лиг – изменение свободной энергии при переходе лиганда из вакуума в раствор, ΔG_рец – изменение свободной энергии при переходе акцептора из вакуума в раствор. Три последних параметра определяют энергию сольватации, изменение которой (ΔG_cолв) вычисляется путем решения обобщенного уравнения Борна [23, 24] для каждого из трех переходов: где ΔG_эл – электростатический, а ΔG_гидр – гидрофобный вклад в свободную энергию сольватации; где ΔG_мм рассчитывается с помощью молекулярной механики путем усреднения энергий взаимодействия лиганда с белком-акцептором на основе записанных снимков, а TΔS – энтропийный фактор, которым можно пренебречь, если требуется сравнить состояния с одинаковой энтропией, например связывание двух лигандов с одним и тем же белком-акцептором.

Впервые такой подход был применен в [22], где использовались расчеты на основе уравнения Пуассона–Больцмана (метод ММ-РBSA). Позднее для ускорения расчетов стали использовать обобщенное уравнение Борна [25, 26]. В настоящей работе применен усовершенствованный метод MM-GBSA, разработанный для макромолекул, в котором борновские радиусы масштабируются для учета межузельных пространств между атомными сферами [24, 27].

В качестве белков-акцепторов выбраны PSP дикого типа и мутантные формы фермента PSP-E125A и PSP-D649A, обладающие повышенной активностью и измененной вторичной субстратной специфичностью соответственно. Известно, что ОрВ, включая PSP, гидролизуют пептидную связь, образованную карбоксильной группой положительных аминокислотных остатков (R или К), находящихся в Р1-положении субстрата. Эффективность гидролиза увеличивается, если предшествующий остаток, занимающий P2-положение, тоже положительный. Все охарактеризованные ОрВ, включая PSP и PSP-E125A, предпочитают субстраты с остатком аргинина в P2-положении, тогда как PSP-D649A эффективнее гидролизует субстраты с остатком лизина в P2-положении. В качестве лигандов использованы тетрапептиды GRRG и GKRG, содержащие остаток аргинина в P1-положениях и остатки аргинина или лизина в Р2-положениях соответственно. Для всех шести систем проведено МД-моделирование продолжительностью по 50 нс.

Первоначальный анализ траекторий включал в себя расчет среднеквадратичных отклонений (СКО) Cα-атомов от начального положения (рис. 1), изменения радиусов инерции комплексов (рис. 2), а также анализ подвижности отдельных аминокислотных остатков в процессе МД-моделирования (рис. 3).

Из графика на рис. 1а видно, что в случае PSP дикого типа сначала (через 1.5 нс моделирования) величины СКО для обеих систем (комплексы с GRRG и GKRG) стабилизируются на уровне около 0.275 нм, затем (на 6-й нс) опять наблюдается увеличение СКО до 0.35 нм (комплекс с GKRG) и до 0.45 нм (комплекс с GRRG). В конце моделирования (начиная с 43 нс) значения СКО в обеих системах стабилизируются на уровне 0.33 нм. В случае PSP-E125A обе системы стабилизируются на 15 нс, при этом величины СКО составляют 0.29 и 0.42 нм для комплекса с GKRG и GRRG соответственно (рис. 1б). Комплекс PSP-D649A/ GKRG стабилизируется на 10 нс при значении СКО 0.31 нм, а комплекс PSP-D649A/GRRG – на 20 нс, СКО 0.36 нм (рис. 1в). Учитывая размер моделируемых комплексов (молекулярная масса PSP 76 кДа), можно сказать, что значения СКО свидетельствуют о высокой устойчивости систем, при этом комплексы с пептидом GKRG немного более стабильны при МД-моделировании.

Рисунки 2а и 2в показывают, что в случае PSP дикого типа и PSP-D649A радиусы инерции обеих систем практически не меняются, находясь в диапазоне от 2.6 до 2.65 нм, что указывает на то, что во время МД-моделирования компактность систем существенно не меняется. Это, в свою очередь, свидетельствует о том, что в процессе моделирования не происходит существенного движения доменов друг относительно друга. В то же время радиус инерции комплекса PSP-D649A/GRRG немного увеличивается, оказываясь в диапазоне 2.65–2.7 нм (рис. 2б), что коррелирует с увеличением СКО (рис. 1б) и указывает на незначительное расхождение доменов.

Из графиков на рис. 3 следует, что в случае PSP дикого типа и обоих мутантных белков значения среднеквадратичных флуктуаций (СКФ) Cα-атомов аминокислотных остатков каталитического домена PSP (С-концевая область белка, начиная с 410 аминокислоты) во всех системах не превышают 0.4 нм. Более высокие (>0.4 нм) значения RMSF Cα-атомов наблюдаются только в области петлевых участков пропеллерного домена.

Таким образом, все исследуемые комплексы являются стабильными и позволяют применить метод MM-GBSA для расчета энергии связывания пептидов. Рассчитанные с помощью MM-GBSA изменения свободной энергии при связывании тетрапептидов GRRG и GKRG молекулами PSP дикого типа и мутантными формами фермента приведены в табл. 1.

Из табл. 1 следует, что значение ΔG для комплекса PSP/GRRG примерно в 2 раза больше, чем для комплекса PSP/GKRG, что указывает на более сильное связывание фермента дикого типа с первым пептидом. Это коррелирует с ранее полученными данными о большем количестве межмолекулярных контактов, возникающих между PSP и тетрапептидом с остатком аргинина в Р2-положении [15]. Полученные в [13, 14] кинетические параметры гидролиза п-нитроанилидных субстратов: Z-RR-pNA и Z-KR-pNA (Z – бензилоксикарбонил), которые можно рассматривать как аналоги тетрапептидов GRRG и GKRG, подтверждают, что PSP гидролизует субстрат с аргинином в Р2-положении в 2.7 раза эффективнее, чем с лизином, что указывает на высокую точность расчетов методом MM-GBSA в случае фермента дикого типа.

Наиболее сильное связывание пептидов, исходя из расчетных данных, наблюдается в комплексах PSP-E125A, что коррелирует с более высокими значениями эффективности гидролиза обоих субстратов по сравнению с ферментом дикого типа (табл. 1). В случае субстрата, содержащего лизин в Р2-положении, 2.3-кратное усиление связывания в комплексе PSP-E125A/GKRG по сравнению с комплексом PSP/GKRG хорошо коррелирует с 1.8-кратным увеличением эффективности гидролиза в первой фермент-субстратной паре. В случае субстрата, содержащего аргинин в Р2-положении, увеличение связывания в комплексе PSP-E125A/GRRG по сравнению с комплексом PSP/GRRG составляет только 14%, тогда как каталитическая эффективность мутантного белка в отношении данного субстрата увеличивается в 2.2 раза.

Расчетные энергии связывания пептидов в комплексах PSP-D649A/GKRG и PSP-D649A/ GRRG увеличены в 2 раза и на 11% соответственно по сравнению с PSP-содержащими комплексами (табл. 1). В первом случае это хорошо коррелирует с кинетическими данными; каталитическая эффективность PSP-D649A в отношении Z-KR-pNA в 2.2 раза выше, чем у фермента дикого типа. Во втором случае расчетная величина ΔG оказалась выше, чем ожидалось согласно кинетическим данным. В этой связи нельзя не отметить, что для обоих мутантных белков точность предсказания вторичной субстратной специфичности была заметно более низкой, чем для фермента дикого типа. Можно предположить, что при введении мутаций стабилизация оптимальной конформации белка занимает определенное время, и расчет энергии по полной МД-траектории начинается еще до полной стабилизации системы, что и вносит ошибку. В то же время для большинства фермент-субстратных пар влияние мутации на величину свободной энергии связывания субстрата в активном центре фермента коррелировало с кинетическими данными.

Таким образом, результаты расчетов изменения свободной энергии при образовании фермент-субстратных комплексов в целом коррелируют с экспериментальными кинетическими данными. Это подтверждает целесообразность использования метода MM-GBSA для виртуальной оценки специфичности бактериальных ОрВ в отношении различных субстратов и изучения влияния аминокислотных замен на каталитическую эффективность ферментов, что в совокупности указывает на перспективность применения данного метода для решения таких задач, как изучение механизма каталитической активации бактериальных ОрВ и прогнозирование эффективности низкомолекулярных ингибиторов.

Работа выполнена при финансовой поддержке Российского научного фонда (проект № 21-74-20154) в части моделирования, МД-экспериментов и термодинамических расчетов. Сравнительный анализ методов ММ-РBSA и ММ-GBSA проведен при поддержке Министерства науки и высшего образования РФ в рамках выполнения работ по Государственному заданию ФНИЦ “Кристаллография и фотоника” РАН.