Биоорганическая химия, 2022, T. 48, № 6, стр. 732-744

ВВЕДЕНИЕ

Нуклеозидфосфорилазы (в частности, тимидинспецифичная нуклеозидфосфорилаза КФ 2.4.2.4.) – белки-ферменты, играющие важную роль в синтезе нуклеозидов и азотистых оснований [1]. Нуклеозидфосфорилазы катализируют фосфоролитическое расщепление пуриновых и пиримидиновых нуклеозидов. Некоторые белки этого класса обнаружены во многих раковых клетках, что позволяет использовать их в качестве активаторов противоопухолевых препаратов [2, 3]. У некоторых классов прокариот присутствует широкоспецифичная пиримидиннуклеозидфосфорилаза (PyNP), способствующая расщеплению как тимидина, так и уридина примерно с равной каталитической активностью. При этом тимидинфосфорилаза и широкоспецифичная пиримидиннуклеозидфосфорилаза – единственные представители NP-II-семейства нуклеозидфосфорилаз. Показано, что инфицирование клеток некоторых тканей бактериями Mycoplasma hyorhinis препятствует нормальной фармокинетике фторпиримидинов. Это объясняется активностью фермента PyNP в этих бактериях в отношении 5-фтор-2'-дезоксиуридина, 5-трифтортимидина и 5-фтор-5'-дезоксиуридина [1, 2]. Препятствование ферментативой функции PyNP из M. hyorhinis [4] при фторпиримидиновой терапии обусловливает необходимость разработки конкурентных ингибиторов, специфичных к BsPyNP (пиримидиннуклеозидфосфорилаза из Bacillus subtilis) и неактивных в отношении hTP (тимидиннуклеозидфосфорилаза человека) [3].

Цель данного исследования – поиск предполагаемых ингибиторов тимидинфосфорилазы и широкоспецифичной пиримидиннуклеозидфосфорилазы и их сравнение по энергии связывания с белком-мишенью.

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

Молекулярный докинг и молекулярная динамика. В качестве стартовой трехмерной модели белка взята определенная нами ранее методом рентгеноструктурного анализа структура пиримидиннуклеозидфосфорилазы из Bacillus subtilis (BsPyNP) [5, 6], гомологичной пиримидиннуклеозидфосфорилазе из M. hyorhinis. Для поиска соединений, аффинных активному центру BsPyNP, среди существующих биологически активных соединений, удовлетворяющих правилам “пяти” Липински [7], были отобраны четыре соединения. Первые два взяты из работы [8]: 2',3'-дидегидро-3'-дезокситимидин (d4T) и 1-(2-дезокси-2-фтор-β-D-арабинофуранозил)-5-иодурацил (фиауридин, FIAU) (рис. 1). Поскольку в молекуле FIAU заместителем в пятом положении пиримидинового цикла является атом иода, который может создать дополнительные стерические ограничения при связывании вероятного ингибитора с активным центром фермента, было решено взять в качестве третьего лиганда то же соединение, но без заместителя – атома иода – 1-(2-дезокси-2-фтор-β-D-арабинофуранозил)-5-урацил (FAU). Моделирование связывания двух схожих молекул, но с разными по размеру заместителями по пиримидиновому кольцу, позволяет выяснить степень влияния атома иода на связывание лиганда с белком.

В работе [9] описано соединение 2-пиримидин-2-ил-1H-имидазол-4-карбоновая кислота (PIA) (рис. 1) – предполагаемый специфичный ингибитор BsPyNP, для которого в этой работе был проведен первый этап молекулярного докинга. Однако это соединение при разных значениях pH раствора присутствует в разных формах с разным общим электрическим зарядом: в растворах с нейтральной реакцией среды – в виде в целом не заряженного цвиттер-иона, а при кислой – в виде аниона. Различие общего заряда и состояния протонирования лиганда может кардинальным образом повлиять на его связывание с ферментом.

Результатом первого этапа стал набор конформаций для каждого лиганда. Разница в решениях докинга, полученных с помощью разных программ (Autodock и PLANTS), объясняется разницей в методике определения подходящих конформаций. Близкие друг к другу конформации объединены в кластеры по сходству расположения атомов и углов между связями атомов в конформере. Для d4T найден 1 кластер конформаций, для FIAU – 3 кластера (2 (Autodock) + 1 (PLANTS)), для FAU – два кластера (1 + 1) (рис. 2). Для каждой из конформаций был проведен анализ энергии взаимодействия и устойчивости комплекса (рис. 3) в зависимости от времени молекулярно-динамической симуляции частиц системы и получены, кроме стартовой, конформации с минимальной энергией взаимодействия по методу линейной интерполяции свободной энергии (LIE). Проведено сравнение конформаций лигандов с конформациями похожих лигандов в структурах экспериментально определенных комплексов тимидинфосфорилаз. По трехмерной структуре тимидин- и пиримидинфосфорилазы схожи, однако для последней не изучены структуры комплексов белок–лиганд.

Таким образом, в работе рассмотрены две конформации: start – начальный вариант после молекулярного докинга, min – варианты с наименьшей энергией связывания белок–лиганд, вычисленной по методу LIE. Для лигандов d4T, FAU и FIAU была подобрана одна конформация со стартовым и одна – с минимальным значением энергии связывания по методу LIE (табл. 1).

Для конформаций PIA, полученных в программах Autodock и PLANTS при pH < 6.5, были найдены два минимума энергии с близкими значениями энергии связывания (табл. 1), но сильно отличающимися положениями атомов, поэтому дальнейшие расчеты проводили для двух конформаций: min₁ и min₂. Аналогичная ситуация характерна для конформаций PIA, полученных с помощью программы PLANTS при pH 7.0. Для них также найдены два минимума со схожими значениями энергий (табл. 1), и для дальнейшей обработки отобраны три конформации (стартовая и две минимальные).

Молекулярная динамика и расчеты методом возмущения свободной энергии (FEP). Полученные величины ΔG_bind представлены в табл. 2–5. Для всех конформаций лигандов с ΔG_bind > 0 вероятность связывания низкая. Это связано с тем, что при большом значении энергии Гиббса реакция связывания лиганда с активным центром белка не спонтанная. Константа ингибирования в этом случает имеет слишком высокие значения, такие, что концентрация лиганда не является достижимой.

Молекулярное моделирование комплекса белок–лиганд FIAU. Присоединение лиганда осуществляется с помощью трех водородных связей между атомами О1 и О2 лиганда и атомом водорода остатка Gly85, между атомом О2 лиганда и атомом N остатка Phe207 (рис. 4). Величина ∆G_bind = = 16.0 кДж/моль (табл. 2) говорит об очень низкой аффинности связывания данного лиганда с активным центром белка.

Молекулярное моделирование для комплекса белок–лиганд d4T. Для d4T в конформации с минимальной энергией взаимодействия ∆G_bind = = –10.2 кДж/моль (K_i = 1.4 мкМ) при Т = 300 К (табл. 2). Отрицательное значение свободной энергии связывания говорит о возможности связывания лиганда d4T с активным центром BsPyNP, т.е. d4T – потенциальный ингибитор. В связи с этим проверили специфичность соединения, сравнив комплекс белок–лиганд для BsPyNP и для hTP, соответствующим катализатором ресинтеза азотистых оснований в теле человека. Для этой цели было проведено выравнивание аминокислотной последовательности BsPyNP и hTP (рис. 5). Выявлены аминокислотные остатки активного центра hTP, отличающиеся от соответствующих остатков BsPyNP. Далее был проведен анализ связывания d4T с активными центрами двух ферментов, d4T образует водородные связи с атомами аминокислотных остатков BsPyNP: атомы кислорода d4T связываются с атомами водорода Ser183 (2.2 Å), Val174 (2.1 Å), Arg168 (2.1, 1.9 и 1.7 Å) и Asp161 (1.8 Å) (рис. 6а). В то же время d4T образует связи в активном центре hTP с аминокислотными остатками Arg171 (3 Å), Ser186 (2.7 и 3.5 Å) и Lys109 (2.9 Å) (рис. 6б). Один из остатков – Arg168 (рис. 6а) и этот же остаток под номером Arg171 (рис. 6б) – встречается в обоих случаях.

Молекулярное моделирование комплекса белок–лиганд PIA. Молекула PIA была исследована в двух вариантах протонирования: в форме карбоксианиона (нейтральные значения pH) и цвиттер-иона (кислая среда).

При pН < 6.5 ΔG_bind = –6.2 кДж/моль для конформации min₂, что лучше с точки зрения энергии связывания с белком, чем стартовая конформация, полученная в результате молекулярного докинга, и конформация min₁ (табл. 3). Однако это значение намного превышает значение ΔG_bind для конформации анионной формы PIA. Для нее проведено моделирование не только с нативной BsPyNP, но и с BsPyNP с мутациями в активном центре, соответствующими активному центру hTP, для сравнения активности PIA в бактериальной клетке и клетках человека. Это сделано не только для определения различий в аффинности связывания лиганда, но и для выявления конкретных остатков активного центра BsPyNP, обусловливающих различие. При pН > 6.5 ΔG_bind = –65.5 кДж/моль для конформации min₁, что говорит о высокой аффинности связывания PIA с активным центром в этой конформации (табл. 4).

Водородные связи образуются между молекулой предполагаемого ингибитора – PIA – и остатками Lys108, His82, Lys188 и Ser110 белка BsPyNP (рис. 7). Различные найденные конформации PIA отличаются в основном положением карбоксильной группы, которая связывается водородными связями с Lys108 и His82. Эти аминокислотные остатки было решено заменить на соответствующие остатки hTP и провести расчет свободной энергии методом возмущения свободной энергии (FEP).

Созданы структурные модели двух мутантов: с заменой остатка Lys108 на Trp108 и с заменой остатка His82 на Gly82. В результате для получившихся моделей ∆G_bind варьируется от –25.1 до 4.4 кДж/моль (табл. 5), что существенно меньше по модулю энергии для лиганда в активном центре нативной BsPyNP при pH > 6.5.

После замены аминокислотных остатков в активном центре белка остатки, образующие водородные связи с лигандом, изменились (рис. 8). Водородные связи образуются между карбоксильной группой лиганда и атомами кислорода для Thr108 (2.9 Å), Lys188 (3.4 и 3.0 Å) и Ser110 (2.6 Å). Также образуется водородная связь между карбоксильной группой и атомом азота в случае Lys188 (3.1 Å). В целом наблюдаются большие расстояния между атомами при связывании PIA с мутантными формами PyNP, что подтверждается меньшей по модулю ∆G_bind.

ЭКСПЕРИМЕНТАЛЬНАЯ ЧАСТЬ

Молекулярный докинг. В качестве исходной структуры для молекулярного докинга была использована структура белка BsPyNP, полученная в ходе рентгеноструктурного анализа [5]. Для получения наиболее полного набора конфигураций использовали программы молекулярного докинга Autodock [10] и PLANTS [11], работающие с разными алгоритмами поиска глобального минимума: Аutodock использует алгоритм Монте–Карло [12], а PLANTS – алгоритм оптимизации подражания муравьиной колонии [13]. В обеих программах для параметризации атомов использовали полуэмпирическое силовое поле AMBER. В качестве активной области докинга был выбран куб с ребром 22 Å, в который были включены все атомы аминокислотных остатков активного центра белка BsPyNP. Шаг сетки докинга составлял 0.375 Å. Аминокислотные остатки считались неподвижными. В ходе докинга было отобрано по 10 конформаций для каждого лиганда с самыми высокими значениями оценочных функций, используемых по умолчанию в каждой из программ. Конформации объединены в кластеры по параметру среднеквадратичного отклонения координат атомов лиганда (рис. 2).

Молекулярное моделирование методом классической молекулярной динамики. Для исследования устойчивости связывания лигандов с PyNP проводили молекулярно-динамическое (МД) моделирование этого комплекса в пакете программ GROMACS (версия 2018.6) [13, 14]. Компоновку и параметризацию MД-систем комплексов белок–лиганд проводили в web-сервисе CHARMM-GUI [13–15] с использованием набора полноатомных силовых полей CHARMM [16, 17] версии C36m [15]. Модельные ячейки заполняли молекулами воды, которые описывались трехцентровой моделью TIP3P, с добавлением нейтрализующих ионов Na⁺ и Cl^– с концентрацией 0.1 M.

Для интегрирования с шагом 2 фс уравнений движения Ньютона использовали алгоритм “leapfrog”. Длина траектории молекулярной динамики составляла 20 нс. Расчет дальних электростатических взаимодействий проводили методом частица–сетка (PME) [18, 19]. Для описания ван-дер-ваальсовых взаимодействий использовали функцию сглаживания “force–switch” в интервале 10–12 Å. Давление в системе контролировали баростатом Паринелло–Рахмана [19] на уровне 1 бар. Температуру МД-системы поддерживали постоянной с помощью термостата V-rescale [20] на уровне 300 К. Процедуре молекулярной динамики предшествовал этап минимизации энергии системы с использованием выбранного силового поля в виртуальной ячейке методом градиентного спуска и уравновешивание системы путем моделирования с ограничениями, накладываемыми на положения неводородных атомов белка и лиганда в NVT- и NPT-ансамблях в течение 200 пс (использовали баростат Берендсена).

Определение относительной аффинности связывания лигандов методом линейной интерполяции свободной энергии (LIE). Приблизительную оценку аффинности связывания лигандов методом LIE [21] осуществляли в программе gmx lie программного пакета GROMACS (версия 2018.6) [13, 14]. Метод использовали как для оценки энергии связывания лиганда в конформациях, полученных на этапе молекулярного докинга, так и для поиска дополнительных конформаций лиганда с помощью молекулярно-динамической симуляции. Для этого предварительно рассчитывались средние энергии силового поля ван-дер-ваальсовых и электростатических взаимодействий лигандов с водой в течение 20 нс.

Для расчета потенциальной энергии кулоновского взаимодействия лиганда с окружением траектории МД-симуляции белок–лиганд перерабатывали с использованием метода обобщенного поля реакции (Generalized RF) при отсечке 12 Å. Эмпирические константы α и β, определяющие вклад ван-дер-ваальсовых и электростатических взаимодействий в оценку энергии связывания лигандов, выбраны в соответствии со статьей [22]. Коэффициент сдвига γ был принят равным нулю в связи с невозможностью его подбора на данном этапе исследования. Коэффициент γ – дополнительный, он зависит от особенностей модели и гидрофобности активного центра и требуется для согласования данных LIE с экспериментальными данными.

Расчет свободной энергии связывания лиганда с ферментом BsPyNP методом возмущения свободной энергии (FEP). Расчет свободной энергии связывания лигандов с белком BsPyNP проводили, используя термодинамический цикл, описанный Aldeghi et al. [23]. Первый этап заключался в сборе отсчетов производной гамильтониана системы по параметру связи Кирквуда (λ-параметр, [24]) (для метода термодинамического интегрирования (TI) и разницы в потенциальной энергии системы между текущим состоянием и остальными (для методов на основе FEP)).

Для систем белок–лиганд и белок–вода проводили МД-симуляцию в NPT-ансамбле в течение 2.4 нс, используя алгоритм интегрирования “leapfrog” уравнений стохастической динамики для каждого набора λ-параметров. С целью удержания лиганда в области связывания при уменьшении силы взаимодействия с белком вводили дополнительные межмолекулярные ограничения – ориентационные рестрейны [23], используя программу PyFEPRestr [25]. Для системы белок–лиганд параметр связи λ менялись последовательно для λ_restr (0, 0.010, 0.025, 0.050, 0.075, 0.100, 0.200, 0.350, 0.500, 0.750, 1.000) и λ_coul (0, 0.25, 0.50, 0.75, 1.00); λ_VdW (0, 0.05, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, 1.00). В случае системы лиганд–вода параметры связи менялись лишь для кулоновских и ван-дер-ваальсовых взаимодействий лиганда с водой. Необходимую для вычисления свободной энергии связывания (ΔG_bind) поправку на энергию введенных ограничений для системы лиганд–вода (ΔG^solv_rest) считали аналитически по формуле, приведенной в работе [26]. Ван-дер-ваальсовые взаимодействия рассчитывали, используя потенциал “мягкого ядра” (soft core) Леннарда–Джонса. Для учета как дальних кулоновских, так и ван-дер-ваальсовых взаимодействий использовали метод частица–сетка (PME) [27].

До проведения МД-симуляций в дополнение к общему уравновешиванию системы при каждом значении параметра связи λ проводили процедуры оптимизации геометрии и уравновешивания системы в NVT- и NPT-ансамблях в течение 100 пс.

Обработку данных МД-симуляций при разных параметрах связи Кирквуда осуществляли в написанной в среде Jupyter Notebook программе, с использованием библиотек alchemlib и pymbar. Первый шаг – извлечение разницы редуцированных потенциалов (ΔU_ij) – выполняли с использованием библиотеки alchemlyb. Затем проводили отбор нескорелированных отсчетов ΔU_ijи отсечение неуравновешенной части траекторий с использованием функции detectEquilibration библиотеки pymbar [28]. В качестве входного параметра для расчета статистической неэффективности использовали значения ΔU_ijмежду МД-симуляциями с соседними λ. На третьем шаге статистически независимые отсчеты ΔU_ijобрабатывали с использованием метода множественного отношения вероятности принятия шага Беннетта (MBAR) [29]. Корректность расчетов контролировали по матрице перекрытия между λ-состояниями, а также сравнением полученных значений ΔG со значениями, полученными другими методами: EXP (метод экспоненциального усреднения) [28], BAR (метод отношения вероятности принятия шага Беннета) [29] и TI (метод термодинамического интегрирования).

Поскольку суммарный заряд атомов PIA при pH > 6.5 равен –1e, необходимо было в этом случае вычислить поправки на суммирование по Эвальду и периодические граничные условия, их рассчитывали с использованием аналитической схемы [30]. Заряд белка был полностью нейтрализован противоионами, эффективный Q_P = 0. 1/4πε_o полагался равным 138.93545585 кДж нм e⁻² моль⁻¹, ε_s для TIP3P воды –82, ρ_300K = 997 кг/м³, γ_S = 0.0764 e нм². Остаточные интегральные потенциалы (RIP) молекулы белка и лиганда рассчитывали с использованием программы APBS [24].

Свободную энергию связывания белок–лиганд (ΔG_bind) рассчитывали по формуле:

где ΔG^coul – свободная энергия кулоновского взаимодействия; ΔG^vdw – свободная энергия ван-дер-ваальсового взаимодействия, определяемого как потенциал Ленарда–Джонса; ΔG^rest – свободная энергия рестрейнов (в случае системы лиганд–растворитель рассчитывается аналитически по формуле, описанной в работе [26].

В случае заряженного лиганда как к ΔG^prot, так и ΔG ^solv добавляются соответствующие поправки: ΔΔG_DSF – часть поправки, связанной с использованием явной (дискретной) модели растворителя, учитывающая конечные размеры ячейки моделирования; ΔΔG_NET – поправка на влияние периодических граничных условий на кулоновское взаимодействие бесконечного количества периодических копий точечного заряда, ΔΔG_USV – поправка на влияние эффекта неполной сольвации; ΔΔG_RIP – поправка на влияние остаточного интегрированного потенциала. Часть поправки, связанной с использованием явной (дискретной) модели растворителя, для случая системы бесконечного объема ячейки моделирования (DSI) не включена в выходные данные, т.к. ее значение (–74.11 кДж/моль) одинаково для всех рассматриваемых в статье систем белок–лиганд и лиганд–растворитель и обнуляется при расчете ΔG_bind [24, 30, 31].

ЗАКЛЮЧЕНИЕ

В ходе работы был проведен поиск специфичного ингибитора пиримидиннуклеозидфосфорилазы. В качестве возможных лигандов были изучены 2',3'-дидегидро-3'-дезокситимидин (d4T), 1-(2-дезокси-2-фтор-β-D-арабинофуранозил)-5-иодурацил (фиауридин, FIAU), 1-(2-дезокси-2-фтор-β-D-арабинофуранозил)-5-урацил (FAU) и 2-пиримидин-2-ил-1H-имидазол-4-карбоновая кислота (PIA) с помощью методов молекулярного моделирования комплексов белок–лиганд, таких как молекулярный докинг и молекулярная динамика, позволяющие найти энергию связывания белка с лигандом (ΔG_bind) с использованием методов LIE и FEP. По данным моделирования, а именно по величине свободной энергии связывания белок–лиганд, можно сделать вывод о том, что лиганды d4T и PIA – потенциальные ингибиторы бактериальных широкоспецифичных пиримидиннуклеозидфосфорилаз, т.к. их связывание с активным центром энергетически выгодно.

В настоящее время соединение d4T выпускается под названием Ставудин (ОАО Фармсинтез, Россия) и применяется для лечения ВИЧ [32, 33]. Возможный ингибитор PIA не проходил клинические исследования. Клинические исследования для лиганда FAU не проводились, а для лиганда FIAU в ходе эксперимента была обнаружена высокая токсичность [34].

Полученные в ходе расчетов результаты указывают на то, что оба соединения – PIA и d4T – связываются с активным центром BsPyNP с наибольшей аффинностью среди других предполагаемых субстратов. Причем при моделировании связывания PIA в средах с различным pH энергия принимает более высокие значения по модулю, чем при нейтральном pH. Это, вероятно, объясняется образованием шести водородных связей с 3 а.о. белка при pH > 6.5, в отличие от пяти водородных связей при pH < 6.5 в случае PIA. Дополнительная водородная связь образуется с Lys108. Особенную роль здесь играет значение pH, т.к. в бактериальных клетках значения pH приближены к 7, т.е. химические реакции проходят в нейтральной среде. Следствием специфического связывания белка BsPyNP с выбранными соединениями может стать снижение скорости процессов ресинтеза азотистых оснований, что, в свою очередь, может привести к гибели бактериальных клеток.

В ходе расчетов свободной энергии связывания для комплекса белок–лиганд с замененными аминокислотными остатками в активном центре, имитирующими активный центр тимидинфосфорилазы человека, было обнаружено: PIA хуже связывается с hTP чем с BsPyNP, это минимизирует возможные побочные эффекты применения этого соединения в терапевтических целях.