Кристаллография, 2021, T. 66, № 5, стр. 836-840

Pyfeprestr: программное дополнение к системе молекулярной графики PyMol для расчетов свободной энергии связывания лиганда с рецептором

А. А. Лашков 1*, И. В. Толмачев 2, П. А. Эйстрих-Геллер 1, С. В. Рубинский 1

1 Институт кристаллографии им. А.В. Шубникова ФНИЦ “Кристаллография и фотоника” РАН
Москва, Россия

2 Школа № 179
Москва, Россия

* E-mail: alashkov83@gmail.com

Поступила в редакцию 06.05.2020
После доработки 17.05.2020
Принята к публикации 19.05.2020

Полный текст (PDF)

Аннотация

Расчеты свободной энергии связывания лигандов с рецептором, чаще всего белковой природы, находят широкое применение в рациональном дизайне новых биологически активных соединений. Равновесные методы, основанные на применении молекулярной динамики и “алхимических” преобразований, такие как термодинамическое интегрирование и метод Беннета, являются “золотым стандартом” подобных вычислений. Технически эти методы требуют введения дополнительных ограничений на взаимное расположение атомов в системе лиганд–рецептор. Разработано дополнение к системе молекулярной графики PyMol, облегчающее введение требуемых ограничений и аналитический расчет поправки при вычислении свободной энергии связывания лиганда с рецептором (белком).

ВВЕДЕНИЕ

Создание новых лекарственных препаратов является одним из самых ресурсозатратных процессов в современном мире. Разработка биологически активных соединений – многостадийный процесс. На первом этапе нужно отобрать несколько веществ – кандидатов – для дальнейших доклинических и клинических испытаний. Если известна молекулярная мишень действия разрабатываемого биологически активного соединения, то проводят либо химический, либо виртуальный скрининг библиотеки соединений на степень аффинности к молекулярной мишени. Во втором случае крайне важными являются количественные оценки степени связывания молекулы биологически активного соединения (лиганда) с молекулярной мишенью – рецептором (чаще всего белковой молекулой).

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

Для количественной оценки аффинности связывания лиганда с рецептором необходимо применение более точных, но и более ресурсоемких методов, основанных на молекулярной динамике (МД): либо моделирование доставки лиганда к месту связывания, либо изменение характера связей лиганда с рецептором на основе так называемых “алхимических” преобразований [4]. На сегодняшний день среди методов, основанных на равновесных “алхимических” преобразованиях, широкое применение получили: термодинамическое интегрирование [4, 5], метод отношения вероятности принятия шага Беннета [6], метод множественного отношения вероятности принятия шага Беннета [7]. Технически эти методы требуют введения дополнительных ориентационных ограничений на взаимное расположение атомов в системе лиганд–рецептор для предотвращения при отключении взаимодействий лиганда с рецептором и изменении параметра связи Кирквуда [8] выталкивания лиганда из его сайта связывания. Лиганд оказывается в окружении, сильно отличающемся от окружения лиганда в сайте связывания. Это противоречит сути равновесных методов расчета свободной энергии, основанных на “алхимических” преобразованиях.

Для удобного введения ориентационных ограничений на положение лиганда и расчета поправки к свободной энергии связывания при введенных ограничениях авторами разработано программное дополнение к системе молекулярной графики PyMol [9]. Альтернативой является прямая правка файлов настроек МД-симуляции с использованием текстового редактора с поиском индексов атомов в файле структуры, что представляет собой довольно трудоемкий процесс, сопровождающийся повышением вероятности ошибок.

ТЕРМОДИНАМИЧЕСКИЙ ЦИКЛ

Для расчета свободной энергии связывания часто используется термодинамический цикл, описанный в [1]. Для расчета необходимо провести моделирование с последовательным изменением параметра связи [8] атомов лиганда для систем лиганд–рецептор и лиганд–растворитель. Параметр связи меняется последовательно сначала для кулоновских взаимодействий, потом для взаимодействий, описываемых потенциалом Леннард-Джонса. Ограничения на положение атомов лиганда в сайте связывания рецептора также необходимо вводить постепенно до начала изменения параметров связи. Это нужно для корректного расчета изменения энергии связывания, обусловленного введением ограничений в системе лиганд–рецептор в процессе моделирования.

В случае цикла, в котором в отдельных ветвях реализуются переходы молекулы лигада из связанного состояния в несвязанное в системах лиганд–растворитель (“solv”) и лиганд–рецептор (“prot”), свободная энергия связывания рецептора с лигандом равна:

(1)
${{\Delta }}G_{{{\text{bind}}}}^{0} = {{\Delta }}G_{{{\text{coul}} + {\text{vdw}}}}^{{{\text{solv}}}} + {{\Delta }}G_{{{\text{restr}}\_{\text{on}}}}^{{{\text{solv}}}} - {{\Delta }}G_{{{\text{restr}} + {\text{coul}} + {\text{vdw}}}}^{{{\text{prot}}}},$
где ${{\Delta }}G_{{{\text{restr}} + {\text{coul}} + {\text{vdw}}}}^{{{\text{prot}}}}$ – изменение свободной энергии в системе лиганд–рецептор в процессе перехода лиганда из связанного состояния в несвязанное, но с введенными дополнительными ограничениями, ${{\Delta }}G_{{{\text{coul}} + {\text{vdw}}}}^{{{\text{solv}}}}$ – изменение свободной энергии в системе лиганд–растворитель в процессе перехода лиганда из взаимодействующего с растворителем состояния в невзаимодействующее, ${{\Delta }}G_{{{\text{restr\_on}}}}^{{{\text{solv}}}}$ – аналитически рассчитываемая поправка на то, что в системе лиганд–растворитель не вводятся соответствующие ограничения на положение атомов лиганда.

В случае заряженного лиганда при использовании методов, основанных на суммировании по Эвальду [10], для расчета кулоновских взаимодействий и периодических граничных условий необходимо в формулу (1) вводить дополнительные слагаемые, описанные в [11]. В случае применения схемы отсечки либо сдвигового потенциала для расчета дисперсионного члена потенциала Леннард-Джонса в формулу (1) желательно внести поправку, описанную в [12]. При использовании метода суммирования по Эвальду для расчета дисперсионного члена потенциала Леннард-Джонса [13] эта поправка не требуется.

ПРИМЕНЕНИЕ ОРИЕНТАЦИОННЫХ ОГРАНИЧЕНИЙ В ПАКЕТЕ GROMACS

Начиная с версии 2016 в пакете программ МД GROMACS [14] имеется возможность введения ограничений на взаимное положение атомов, даже относящихся к разным молекулам. Ограничения могут быть введены на расстояние между двумя атомами, угол между двумя лучами, выходящими из одного атома и проходящими через два других, двугранный угол между двумя плоскостями, которые проходят через четыре атома. Ограничения могут быть описаны разными функциями потенциальной энергии. Чаще всего используют гармоническое соотношение между отклонением параметра от равновесного значения (r0, θ0, φ0) и потенциальной энергией (U) при введенном ограничении:

(2)
$U(r) = \frac{{{{K}_{r}}}}{2}{{(r - {{r}_{0}})}^{2}},$
(3)
$U({{\theta }}) = \frac{{{{K}_{{{\theta }}}}}}{2}{{({{\theta }} - {{{{\theta }}}_{0}})}^{2}},$
(4)
$U(\phi ) = \frac{{{{K}_{\phi }}}}{2}{{(\phi - {{\phi }_{0}})}^{2}}.$

В [15] описана общая схема введения ориентационных ограничений на положение лиганда в системе лиганд–рецептор, включающая одно ограничение на расстояние между атомом рецептора (a) и лиганда (А), два ограничения на углы baA и aAB и три ограничения на двугранные углы cbaA, baAB, aABC (рис. 1). Там же даны рекомендации по выбору значений коэффициентов жесткости K. В большинстве случаев они должны лежать в диапазоне 5–50 ккал/(моль Å2) (ккал/(моль рад2) в случае ограничений по углам и двугранным углам). В [16] при вычислении свободной энергии связывания ингибитора и субстрата с уридинфосфорилазой из холерного вибриона использовали коэффициенты жесткости в диапазоне 10–20 ккал/(моль Å2) (ккал/(моль рад2)).

Рис. 1.

Схема, описанная в [15], для введения дополнительных ограничений на положение лиганда в сайте связывания рецептора (белка).

Для описания рассматриваемых в настоящей работе ограничений в программном пакете GROMACS в основном файле топологии, начиная с версии 2016, добавлена секция [intermolecular interaction]. Она должна быть последней секцией в этом файле. В этой секции необходимо описать три подсекции [bonds], [angles] и [dihedrals]. Каждое ограничение (из шести) описывается одной строкой в соответствующей подсекции. В строке задаются глобальные индексы атомов, на взаимное положение которых будет накладываться ограничение, тип функции потенциальной энергии (в соответствии с документацией GROMACS), равновесные значения параметра ограничения (в нанометрах для расстояния между атомами и в градусах для углов) и константы жесткости для двух состояний: A (λ = 0) и B (λ = 1). В задаче вычисления свободной энергии связывания лиганда с белком K (λ = 0) полагают равным нулю, а равновесные значения параметров ограничений для состояний A и B равными между собой. При моделировании реализуется переход из состояния A в состояние B через промежуточные состояния – число состояний и значения λ для них контролируют с помощью настроек моделирования.

Как показывает практика авторов настоящей и других работ [1], при использовании в методе МД алгоритма жестких ограничений LINCS [17] для описания валентных связей атомов водорода с другими атомами дополнительные ограничения на атомы водорода накладывать не следует.

ВВЕДЕНИЕ ОРИЕНТАЦИОННЫХ ОГРАНИЧЕНИЙ В ПАКЕТЕ NAMD

В пакете МД NAMD [18] также доступен расчет свободной энергии связывания лиганда с рецептором методами возмущения свободной энергии или термодинамического интегрирования. Возможность введения ориентационных ограничений присутствует за счет механизма коллективных переменных [19]. Однако в одном процессе в отличие от GROMACS невозможно совместить последовательное плавное включение ограничений и последующее отключение взаимодействий лиганда с рецептором. Тем не менее можно раздельно проводить эти преобразования, для каждого из которых нужен отдельный файл конфигурации, написанный на языке TCL (Tool Command Language) и содержащий описание коллективных переменных через индексы взаимодействующих атомов и параметры гармонического потенциала для них. Для составления шаблона файлов настроек ограничений в программном пакете NAMD также можно использовать описываемое в настоящей работе программное обеспечение.

ВЫЧИСЛЕНИЕ АНАЛИТИЧЕСКОЙ ПОПРАВКИ НА ВВЕДЕННЫЕ ОГРАНИЧЕНИЯ

Необходимую для вычисления ΔGbind поправку на энергию введенных ограничений для системы лиганд–растворитель с поправкой на стандартное состояние 1 M можно вычислить аналитически по формуле [15]:

(5)
$\begin{gathered} \Delta G_{{{\text{rest}}\_{\text{off}}}}^{{{\text{solv}}}} = - RT\ln \left[ {\frac{{8{{{{\pi }}}^{2}}{{V}_{0}}}}{{r_{{aA}}^{2}\sin {{{{\theta }}}_{a}}\sin {{{{\theta }}}_{A}}}}} \right. \times \\ \times \;\left. {\frac{{\surd {\kern 1pt} {{K}_{{{{r}_{{aA}}}}}}{{K}_{{{{{{\theta }}}_{a}}}}}{{K}_{{{{{{\theta }}}_{A}}}}}{{K}_{{{{{{\varphi }}}_{{ba}}}}}}{{K}_{{{{{{\varphi }}}_{{aA}}}}}}{{K}_{{{{{{\varphi }}}_{{AB}}}}}}}}{{{{{(2{{\pi }}RT)}}^{3}}}}} \right], \\ \end{gathered} $
где ${{r}_{{aA}}}$ – длина псевдосвязи [нм], ${{{{\theta }}}_{a}}$, ${{{{\theta }}}_{A}}$ – угловые ограничения, ${{K}_{{{{r}_{{aA}}}}}}$, ${{K}_{{{{{{\theta }}}_{a}}}}}$, ${{K}_{{{{{{\theta }}}_{A}}}}}$, ${{K}_{{{{{{\varphi }}}_{{ba}}}}}}$, ${{K}_{{{{{{\varphi }}}_{{aA}}}}}}$, ${{K}_{{{{{{\varphi }}}_{{AB}}}}}}$ – константы жесткости для ограничений на расстояние, углы и двугранные углы, ${{V}_{0}}$ – объем стандартного состояния (1.66 нм3), R – универсальная газовая постоянная, T – температура [К], $\Delta G_{{{\text{rest\_on}}}}^{{{\text{solv}}}} = - \Delta G_{{{\text{rest\_off}}}}^{{{\text{solv}}}}$.

Аналитическая поправка, вычисленная по формуле (5), не учитывает возможное влияние ограничений на симметричные группы молекулы в том случае, когда параметры ограничивают переход из одной симметричной формы в другую [20] или имеется высокий потенциальный барьер между симметричными формами, не связанный с применением ограничений. Эта дополнительная поправка должна быть вычислена отдельно и добавлена в формулу (1). Так, если химическая группа преобразуется осью симметрии N-го порядка, то ${{\Delta }}{{G}_{{{\text{sym}}}}} = - RT\ln (N)$. Например, для фенильной группы (N = 2) ΔGsym = –0.413 ккал/моль.

ОПИСАНИЕ ПРОГРАММНОГО ДОПОЛНЕНИЯ К СИСТЕМЕ МОЛЕКУЛЯРНОЙ ГРАФИКИ PYMOL

Программное дополнение к системе молекулярной графики PyMol – PyFepRestr – предназначено для визуализации задания атомов, участвующих в ориентационных ограничениях, задания констант жесткости, автоматической конвертации используемых единиц измерения, аналитического расчета поправки на введенные ориентационные взаимодействия и записи в файл топологии программного пакета МД GROMACS настроек ограничений. Программное дополнение написано на языке программирования Python 2.7, но совместимо и с Python 3+. Для реализации графического интерфейса использован модуль Tkinter стандартной библиотеки подпрограмм. Для взаимодействия с PyMol в программе применяют интерфейс прикладного программирования, реализованный с помощью модулей pymol.cmd и pymol.wizard. Программа совместима с версиями PyMol, начиная с 1.5.0.6.

При запуске дополнения отображается первая форма графического интерфейса пользователя (рис. 1а), в которой необходимо указать ограничения констант жесткости и абсолютную температуру МД-моделирования. Значения констант можно задавать как в кДж/моль, так и в ккал/моль. После задания параметров управление передается главному окну PyMol (рис. 3), где в визуальном режиме на структуре комплекса лиганд–рецептор отмечают атомы, по которым строятся ограничения, удерживающие лиганд в сайте связывания в процессе МД-моделирования. После указания шести атомов отображается следующая форма графического интерфейса пользователя (рис. 2б) с рассчитанными значениями поправки $\Delta G_{{{\text{rest}}}}^{{{\text{solv}}}}$ на введенные ограничения. Также настройки ограничений для программного пакета GROMACS можно отобразить на экране и записать в основной файл топологии. Возможны отображение и запись шаблона файла настройки ограничений через механизм коллективных переменных для программного пакета NAMD. Перевод единиц измерения констант жесткости в принятые в программном пакете МД осуществляется автоматически.

Рис. 2.

Форма графического интерфейса пользователя программного дополнения, предназначенная для ввода входных параметров (а), для просмотра результирующего значения ${{\Delta }}G_{{{\text{restr}}}}^{{{\text{solv}}}}$ и записи в файл топологии ограничений, визуально описываемых на предыдущем этапе работы с программным дополнением (б).

Рис. 3.

Окно просмотра структуры системы молекулярной графики PyMol на этапе визуального описания вводимых ограничений. Шарами разного цвета обозначены атомы лиганда и рецептора, участвующие в задании ориентационных ограничений.

ЗАКЛЮЧЕНИЕ

Авторами разработано дополнение к системе молекулярной графики PyMol, облегчающее введение ориентационных ограничений на положение лиганда в сайте связывания и аналитический расчет поправки при вычислении свободной энергии связывания лиганда с рецептором – PyFepRestr. Описаны программа и используемые в ней алгоритмы, дано теоретическое обоснование. Программа опубликована на сервисе GitHub под лицензией MIT (https://github.com/tolmv/plugin_for_Pymol).

Работа выполнена с использованием оборудования ЦКП ФНИЦ “Кристаллография и фотоника” РАН при поддержке Министерства науки и высшего образования РФ (проект RFMEFI62119X0035) в рамках Государственного задания ФНИЦ “Кристаллография и фотоника” РАН.

Список литературы

  1. Aldeghi M., Heifetz A., Bodkin M.J. et al. // Chem. Sci. 2016. V. 7. P. 207. https://doi.org/10.1039/c5sc02678d

  2. Kitchen D.B., Decornez H., Furr J.R. et al. // Nat. Rev. Drug Discov. 2004. V. 3. P. 935. https://doi.org/10.1038/nrd1549

  3. Li J., Fu A., Zhang L. et al. // Interdiscipl. Sci. Comput. Life. Sci. 2019. V. 11. P. 320. https://doi.org/10.1007/s12539-019-00327-w

  4. Chipot C. // Wiley Interdiscipl. Rev.: Comput. Mol. Sci. 2014. V. 4 (1). P. 71. https://doi.org/10.1002/wcms.1157

  5. Klimovich P.V., Shirts M.R., Mobley D.L. // J. Comput. Aided Mol. Des. 2015. V. 29 (5). P. 397. https://doi.org/10.1007/s10822-015-9840-9

  6. Bennett C.H. // J. Comput. Phys. 1976. V. 22 (2). P. 245. https://doi.org/10.1016/0021-9991(76)90078-4

  7. Shirts M.R., Chodera J.D. // J. Chem. Phys. 2008. V. 129. P. 124105. https://doi.org/10.1063/1.2978177

  8. Kirkwood J. // J. Chem. Phys. 1935. V. 3. P. 300. https://doi.org/10.1063/1.1749657

  9. DeLano W.L., Lam J.W. // Abstr. Papers Am. Chem. Soc. 2005. V. 230. P. 1371.

  10. Hunenberger P.H., McCammon J.A. // Biophys. Chem. 1999. V. 78. P. 69. https://doi.org/10.1016/s0301-4622(99)00007-1

  11. Rocklin G.J., Mobley D.L., Dill K.A. et al. // J. Chem. Phys. 2013. V. 139. P. 184103. https://doi.org/10.1063/1.4826261

  12. Shirts M.R., Mobley D.L., Chodera J.D. et al. // J. Phys. Chem. 2007. V. 111 (45). P. 13052. https://doi.org/10.1021/jp0735987

  13. Wennberg C.L., Murtola T., Hess B. et al. // Chem. Theory Comput. J. 2013. V. 9. P. 3527. https://doi.org/10.1021/ct400140n

  14. Van Der Spoel D., Lindahl E., Hess B. et al. // J. Comput. Chem. 2005. V. 26. P. 1701. https://doi.org/10.1002/jcc.20291

  15. Boresch S., Tettinger F., Leitgeb M. et al. // J. Phys. Chem. B. 2003. V. 107. P. 9535. https://doi.org/10.1021/jp0217839

  16. Эйстрих-Геллер П.А., Рубинский С.В., Прокофьев И.И. и др. // Кристаллография. 2020. Т. 65. № 2. С. 271. https://doi.org/10.31857/S002347612002006x

  17. Hess B., Bekker H., Berendsen H.J.C. et al. // J. Comput. Chem. 1997. V. 18. P. 1463. https://doi.org/10.1002/(SICI)1096-987X(199709)18:123.0.CO;2-H

  18. Phillips J.C. // J. Comput. Chem. 2005. V. 26. P. 1781. https://doi.org/10.1002/jcc.20289

  19. Fiorin G., Klein M.L. // Mol. Phys. 2013. V. 111. P. 3345. https://doi.org/10.1080/00268976.2013.813594

  20. Mobley D.L., Chodera J.D., Dill K.A. // J. Chem. Phys. 2006. V. 125 (8). P. 084902. https://doi.org/10.1063/1.2221683

Дополнительные материалы отсутствуют.