Предсказание аффинности прогестинов к ядерному рецептору прогестерона на основе данных с коррекцией RBA

А.В. Микурова^*, В.С. Скворцов

Научно-исследовательский институт биомедицинской химии имени В. Н. Ореховича, 119121, Москва, Погодинская ул. 10;^*e-mail: a.v.mikurova@ibmc.msk.ru

Ключевые слова: рецептор прогестерона; аффинность; молекулярный докинг; молекулярная динамика; MM-PBSA

DOI: 10.18097/BMCRM00080

ВВЕДЕНИЕ

В настоящее время опубликовано большое число работ, посвящённых оценке in silico афинности как стероидных, так и нестероидных прогестинов к рецептору прогестерона [1] Данные препараты широко применяются в контрацепции, заместительной и противоопухолевой терапии [2,3]. Однако большое число таких работ, демонстрирующих неплохие результаты, ограничиваются гомологичным рядом химических соединений [4-6] и чаще всего используют технологию CoMFA [7], a priory нацеленную на анализ высокогомологичных рядов. В отличие от CoMFA, технологии, использующие такие процедуры как молекулярный докинг, молекулярную динамику, MM-PBSA (MM-GBSA) [8] или другие способы анализа модельных комплексов, позволяют делать предсказания на моделях, обученных на одном гомологичном ряде химических соединений для рядов, отличных от использованных в обучении соединений.

Ранее мы уже успешно применяли процедуры молекулярных динамики и докинга для создания уравнений предсказания аффиности производных 16α,17α-циклоалканопрогестерона к рецептору прогестерона (RP) [6]. Однако тогда был использован ряд соединений, представляющий собой группу стероидов, имеющих несущественные структурные различия. В настоящей работе число стероидных соединений было расширено. Кроме того была добавлена небольшая группа соединений нестероидной природы.

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

Структуры химических соединений и активность в отношении ядерного рецептора прогестерона [6,9,10] представлены на рисунках 1-2 и в таблице 1. Следует отметить, что в различных источниках авторы используют неодинаковые параметры в качестве оценки афинности либо степени ингибирования: в двух выборках используется относительная связывающая активность (Relative Binding Affinity, RBA), а в третьей – IC₅₀. Кроме того, референсное соединение для разных выборок также различается. В нашем случае это либо прогестерон [6], либо мифепристон [10], либо соединение 74 [9]. RBA прогестерона по отношению к мифепристону ниже в 2 раза [6], а соединения 74 в 7 раз [9]. Эти соотношения были использованы для пересчета значений и приведения данных к единому референсному соединению. С учетом применения авторами [6,9,10] различных способов оценки биологической активности соединений, наиболее адекватным представляется использовать непосредственно величину RBA и расчётную величину RIA (Relative Inhibition Activity) относительно единого соединения (прогестерона). Для построения уравнений использовали десятичный логарифм RBA (RIA), что позволяет регуляризовать ошибку предсказания. Разумеется, напрямую объединять величины RBA и RIA в одну выборку нельзя, однако можно ожидать, что данные величины должны между собой хорошо коррелировать, так как в основе своей имеют одно и тоже явление – связывание с рецептором. Учитывая, что выборка с RIA мала, её использовали как тестовую.

Не все данные были получены по отношению к рецептору прогестерона человека. Данные, взятые из работы [6], получены на рецепторе кролика. Возможность использовать эти данные для анализа с использованием рецептора человека обсуждается в работе [6] и основана на высокой степени идентичности RP у этих видов и отдалённости имеющихся аминокислотных различий от сайта связывания лиганда.

Все структуры лигандов были построены и оптимизированы средствами программы Sybyl-X [11] с расстановкой частичных атомных зарядов по эмпирической схеме с использованием поля сил MMFF [12]. Эту же схему частичных атомных зарядов использовали для лигандов и на последующих шагах.

Закрыть окно ОТКРЫТЬ В ПОЛНОМ РАЗМЕРЕ

Таблица 1. Структура и экспериментальная оценка аффиности к рецептору прогестерона используемых соединений (см. рис. 1-2). Отбор варианта комплекса после докирования.

Моделирование комплексов осуществляли посредством докинга молекул потенциальных лигандов к участку связывания прогестерона в лиганд-связывающем домене (Ligand Binding Domain, LBD) RP. Структура LBD была взята из файла 1A28 [13] в Protein Data Bank.

Построение моделей предсказания аффинности проводили на базе параметров комплексов, рассчитанных методами MM-PBSA и MM-GBSA [8]. Процедуры докирования и молекулярной динамики выполняли средствами программ Dock 6.7 [14] и Amber 16.0 [15] (поля сил FF14SB и GAFF2) по схеме, примененной ранее [16].

Для выполнения расчетов был использован гибридный высокопроизводительный вычислительный комплекс Федерального исследовательского центра «Информатика и управление» Российской Академии Наук (ФИЦ ИУ РАН) [17] IBM Minsky (2xPower8, 4хTesla P100). Основной объём вычислений проводили с использованием GPU ускорителей. Часть вычислений выполняли на CPU, так как использование GPU было либо неэффективно, либо не поддерживалось для данного типа вычислений. Среднее время расчётов на 1 комплекс составило 34 мин для вычислений на GPU и 5 ч на CPU.

Этапы вычислений:

• Формирование стартового набора ориентаций (поз) лиганда процедурой докирования. Всего отбирались 3 лучших по значению оценочной функции (Grid-Based Score) варианта после предварительной кластеризации с RMSD 2.0. Число вариантов ориентации 250000, шаги начальной оптимизации: 0.2 Å для смещения, 0.2° для вращения и 2° для торсионных углов. Программа Dock 6.7 не имеет штатной возможности компиляции для использования GPU ускорителей. Однако, в связи с тем, что докирование каждой пары лиганд/белок - независимая процедура, наилучшим решением является проведение данной процедуры сразу на всём массиве структур с использованием для каждой из пар отдельного потока общим числом по числу доступных CPU ядер. На IBM Minsky оптимальным было использование 64 из 128 ядер.
• Формирование набора файлов для запуска Amber модулем tleap, включающем сольватацию 3-х вариантов комплексов (либо меньше, если процедура докирования нашла менее 3-х решений) для каждой пары. Сольватация водой (модель TIP3P). Слой воды до границ прямоугольного бокса не менее 4 Å. Не требует использования GPU.
• Минимизация общей энергии системы (1000 шагов). Этот и оба предыдущих шага объединяли в один поток (шаг 1), запуская параллельно на 64 ядрах CPU. Каждую пару лиганд/белок обсчитывали независимо для всего набора лигандов (до 3-х вариантов комплексов).
• Разогрев системы от 0 K до 300 K. Использовали симуляцию молекулярной динамики на временном промежутке 20 пс с шагом в 2 фс в периодических граничных условиях (NTV ансамбль) с заторможенной структурой белка. Вычисления проводили на GPU ускорителе, запуская последовательно все варианты (шаг 2); время счёта 1-2 мин на каждый из вариантов.
• Выравнивание плотности системы при 300 K. Использовали симуляцию молекулярной динамики на временном промежутке 20 пс с шагом в 2 фс в периодических граничных условиях (NTP ансамбль) с заторможенной структурой белка. Данный этап характеризуется тем, что объём, форма и размеры бокса с системой изменяются, что не позволяет использовать GPU, и процедуру вынужденно проводили параллельно на доступных ядрах CPU (шаг 3). На используемой системе время счёта одного варианта составляло 30-35 мин, что сопоставимо со всем временем расчёта варианта на GPU.
• Уравновешивание системы при 300 K. Использовали симуляцию молекулярной динамики на временном промежутке 50 пс с шагом в 2 фс в периодических граничных условиях (NTV ансамбль).
• Симуляция молекулярной динамики (продуктивная динамика) при 300 K на временном промежутке 250 пс с шагом в 2 фс в периодических граничных условиях (NTV ансамбль), запись состояния системы каждые 10 пс. Оба последних этапа выполняли на GPU последовательно для всех вариантов (шаг 4).
• Запись состояний в процессе продуктивной динамики создавала набор из 25 точек для оценки энергии комплекса методом MM-PBSA. В существующем виде для данной процедуры использование GPU затруднено, так как требует модификации стандартных процедур, имеющихся в пакете Amber, однако процедура поддерживает экстенсивное распараллеливание, так как каждое состояние обрабатывается независимо, поэтому в нашем варианте использовали 25 ядер CPU одновременно, несмотря на то, что на используемой системе оптимальным было бы использование 16 ядер (после чего наблюдалось снижение производительности для конкретного варианта). Процедуру выполняли последовательно для всех вариантов комплексов (шаг 5).

Рассматривали два способа отбора вариантов. Первый – по усреднённому по 25 состояниям значению изменения свободной энергии комплекса, рассчитанному по методу MM-PBSA. Второй (для стероидов) – отбор решения, подобного расположению прогестерона в кристаллической структуре. Значения отдельных энергетических компонентов, рассчитанных методом MM-PBSA/MM-GBSA, использовали в качестве независимых переменных для уравнений предсказания аффинности. Были использованы: изменение величины электростатического взаимодействия (ELE); изменение величины ван-дер-ваальсовых взаимодействий (VDW); гидрофобный (PBSUR) и сольватационный (PBCAL) вклады в изменение свободной энергии, рассчитанной методом Пуассона-Больцмана (PB); аналогичные вклады, рассчитанные обобщенным методом Борна (GB) (GBSUR и GB); рассчитанные модулем NMODE [9] значения трансляционного (TSTRA), ротационного (TSROT) и колебательного (TSVIB) энтропийных вкладов; молекулярный вес лиганда. Уравнения линейной регрессии оценивали по величине Q² в процедуре скользящего контроля методом исключения по одному.

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

Результаты процедуры докирования отражены в таблице 1. Не для каждого случая программа Dock нашла решение. Кроме того, как видно из таблицы 1, не в каждом случае лучшим по изменению энергии рассчитанной методом MM-PBSA было положение, аналогичное положению прогестерона в кристалле. Для нестероидных соединений у нас не было априорного положения, поэтому в качестве основного была отобрана та ориентация лиганда, которая чаще всего была лучшей согласно энергетической оценке (рис. 2).

В таблице 2 представлены статистические характеристики ряда уравнений, полученных для двух отдельных выборок со стероидными соединениями и объединённой выборки, а также результаты тестирования на выборке нестероидных соединений. Наилучшие результаты в обучении показала выборка 1 (R²_{обучения} 0.73). Однако, учитывая, что в обучении использовали только 23 соединения, уравнение с 9 переменными может быть не совсем адекватным. Это видно на примере тестовой выборки (R²_{предсказания} 0.04) и выборки 2 (R²_{предсказания} 0.1). В данных случаях R²_{тестирования} значительно хуже, чем R²_{обучения}. Под R²_{тестирования} для тестовой выборки понимается величина корреляции (R²) между предсказанной величиной RBA и наблюдаемой величиной RBI. Выборка 2 показала весьма посредственный результат (R²_{обучения} 0.16). Соответственно, тестирование на выборке 1 и тестовой выборке также демонстрирует неадекватный результат. Объединение выборки позволяет, с одной стороны, увеличить статистическую значимость полученного уравнения (R²_{обучения} 0.44 – 0.47), а с другой – получить вполне удовлетворительные результаты (рис. 3) для предсказания величины RBI тестовой выборки (R²_{тестирования} 0.51). Модель позволяет уверено ранжировать соединения тестовой выборки по прочности связывания с RP.

Закрыть окно ОТКРЫТЬ В ПОЛНОМ РАЗМЕРЕ

Таблица 2. Параметры уравнений линейной регрессии, полученные при обучении, и результаты тестирования.

ЗАКЛЮЧЕНИЕ

Таким образом, можно заключить, что при использовании достаточно длительной процедуры молекулярной динамики параметры, полученные в результате применения процедуры MM-PBSA (MM-GBSA), можно использовать для создания уравнений предсказания аффинности, причём использование величины RBA позволяет объединить в одну обучающую выборку данные, полученные в различных условиях, а также оценить величины, отличные от RBA.

БЛАГОДАРНОСТИ

Работа выполнена в рамках Программы фундаментальных научных исследований государственных академий наук на 2013-2020 годы. Подготовка и отладка программного обеспечения выполнялась при поддержке гранта РФФИ № 18-29-03100.

ЛИТЕРАТУРА

1. Lagarde, N., Delahaye, S., Jérémie, A., Ben Nasr, N., Guillemain, H., Empereurmot, C., ... & Zagury, J. F. (2017). Discriminating agonist from antagonist ligands of the nuclear receptors using different chemoinformatics approaches. Molecular Informatics, 36(10), 1700020. DOI
2. de Ziegler, D., & Fanchin, R. (2000). Progesterone and progestins: applications in gynecology. Steroids, 65(10-11), 671-679. DOI
3. Lundgren, S. (1992). Progestins in Breast Cancer Treatment: A Reveiw. Acta Oncologica, 31(7), 709-722.. DOI
4. Bursi, R., & Groen, M. B. (2000). Application of (quantitative) structure–activity relationships to progestagens: from serendipity to structure-based design. European journal of medicinal chemistry, 35(9), 787-796.Hillisch A., von Langen J., Menzenbach B., Droescher P., Kaufmann G., Schneider B., Elger W. (2003) Steroids, 68, 869-878. DOI
5. Hillisch, A., von Langen, J., Menzenbach, B., Droescher, P., Kaufmann, G., Schneider, B., & Elger, W. (2003). The significance of the 20-carbonyl group of progesterone in steroid receptor binding: a molecular dynamics and structure-based ligand design study. Steroids, 68(10-13), 869-878. DOI
6. Fedyushkina I.V., Skvortsov V.S., Romero Reyes I.V., Levina I.S. (2013). Molecular docking and 3D-QSAR on 16a,17a-cycloalkanoprogesterone analogues as progesterone receptor ligands, Biomeditsinskaya Khimiya, 2013, vol: 59(6), 622-635. DOI
7. Cramer, R. D., Patterson, D. E., & Bunce, J. D. (1988). Comparative molecular field analysis (CoMFA). 1. Effect of shape on binding of steroids to carrier proteins. Journal of the American Chemical Society, 110(18), 5959-5967. DOI
8. Massova, I., & Kollman, P. A. (2000). Combined molecular mechanical and continuum solvent approach (MM-PBSA/GBSA) to predict ligand binding. Perspectives in drug discovery and design, 18(1), 113-135.Van Helden, S. P., Hamersma, H., & Van Geerestein, V. J. (1996). In Genetic Algorithms in Molecular Modeling, 159-192. DOI
9. Van Helden, S. P., Hamersma, H., & Van Geerestein, V. J. (1996). Prediction of the progesterone receptor binding of steroids using a combination of genetic algorithms and neural networks. In Genetic Algorithms in Molecular Modeling (pp. 159-192). DOI
10. Söderholm, A. A., Lehtovuori, P. T., & Nyrönen, T. H. (2006). Docking and three-dimensional quantitative structure− activity relationship (3D QSAR) analyses of nonsteroidal progesterone receptor ligands. Journal of Medicinal Chemistry, 49(14), 4261-4268. DOI
11. SYBYL-X 2.1. Certara, Princeton, NJ, USA.
12. Halgren, T. A. (1996). Merck molecular force field. II. MMFF94 van der Waals and electrostatic parameters for intermolecular interactions. Journal of Computational Chemistry, 17(5–6), 520-552. DOI
13. Berman, H., Henrick, K., & Nakamura, H. (2003). Announcing the worldwide protein data bank. Nature Structural and Molecular Biology, 10(12), 980.
14. Kuntz, I. D., Blaney, J. M., Oatley, S. J., Langridge, R., & Ferrin, T. E. (1982). A geometric approach to macromolecule-ligand interactions. Journal of Molecular Biology, 161(2), 269-288. DOI
15. Case, D. A., Cheatham, T. E., Darden, T., Gohlke, H., Luo, R., Merz, K. M., ... & Woods, R. J. (2005). The Amber biomolecular simulation programs. Journal of Computational Chemistry, 26(16), 1668-1688. DOI
16. Mikurova, A. V., Skvortsov, V. S., & Raevsky, O. A. (2018). Computational Evaluation of Selectivity of Inhibition of Muscarinic Receptors M1-M4. Biomedical Chemistry: Research and Methods, 1(3), e00072. DOI
17. Federal Research Center Computer Science and Control of Russian Academy of Sciences [Electronic resource]: site. - Moscow: FRC CS RAS.- URL: http://hhpcc.frccsc.ru (application date: 09/12/2018)