Моделирование конформационного равновесия и сольватации при связывании органических соединений с биологическими мишенями

Зейфман, Алексей Александрович

Моделирование конформационного равновесия и сольватации при связывании органических соединений с биологическими мишенями тема автореферата и диссертации по химии, 02.00.03 ВАК РФ

Зейфман, Алексей Александрович АВТОР

кандидата химических наук УЧЕНАЯ СТЕПЕНЬ

Москва МЕСТО ЗАЩИТЫ

2013 ГОД ЗАЩИТЫ

02.00.03 КОД ВАК РФ

Диссертация по химии на тему «Моделирование конформационного равновесия и сольватации при связывании органических соединений с биологическими мишенями»

Автореферат диссертации на тему "Моделирование конформационного равновесия и сольватации при связывании органических соединений с биологическими мишенями"

Федеральное государственное бюджетное учреждение науки Институт опганической химии им. Н.Д. Зелинского

;сийской академии наук

' на правах рукописи

Зейфман Алексей Александрович

Моделирование коиформационного равновесия и сольватации при связывании органических соединений с биологическими мишенями

02.00.03 — Органическая химия

АВТОРЕФЕРАТ

диссертации на соискание ученой степени кандидата химических наук

16 МАП Ш

Москва-2013

005059448

Работа выполнена в Учебно-научном отделе №58 Федерального государственного бюджетного учреждения науки Института органической химии им. Н.Д. Зелинского

Российской академии наук

Научный руководитель

Официальные оппоненты

Ведущая организация

Чилов Гермес Григорьевич

кандидат химических наук

Институт органической химии им. Н.Д. Зелинского Российской академии наук

Пивина Татьяна Степановна

доктор химических наук, профессор

Институт органической химии им. Н.Д. Зелинского

Российской академии наук

Осолодкин Дмитрий Иванович

кандидат химических наук

Московский государственный университет им. М. В. Ломоносова, химический факультет

Федеральное государственное бюджетное учреждение «Научно-исследовательский институт по изысканию новых антибиотиков им. Г.Ф.Гаузе» РАМН

Защита состоится 11 июня 2013 г. в 15.00 часов на заседании Диссертационного совета Д 02.222.01 при Институте органической химии им. Н.Д. Зелинского Российской академии наук по адресу: 119991, Москва, Ленинский пр-т, д. 47.

С диссертацией можно ознакомиться в библиотеке Института органической химии им. Н.Д. Зелинского Российской академии наук.

Автореферат разослан 8 мая 2013 г.

Ученый секретарь диссертационного совета, /?

доктор химических наук /7 Ш^Р^^рг Родиновская Л. А.

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Актуальность темы. Основной парадигмой современного дизайна лекарственных средств является направленное воздействие на определённую биологическую мишень в организме человека. Крайне актуальной и до конца не решённой проблемой остаётся моделирование связывания лигандов с белковыми мишенями известной структуры. На сегодняшний день наиболее точным методом расчёта прочности связывания лигандов с терапевтическими мишенями является метод возмущения свободной энергии (ВСЭ), основанный на молекулярно-механическом описании полноатомной системы белок-лиганд с явным учётом молекул растворителя. Данный подход является наиболее эффективным вычислительным методом оптимизации лекарственных кандидатов, поскольку позволяет учесть ряд факторов (например, подвижность белка, конформационное равновесие лиганда и сольватацию системы), которые игнорируются в методах (^БАЯ и молекулярного докинга, но зачастую определяют различие в связывании близких по структуре соединений.

Точность метода ВСЭ определяется корректностью молекулярно-механического силового поля и качеством сэмплирования системы. Предметом исследования настоящей работы являлись подходы к улучшению сэмплирования при использовании стандартных параметров силового поля.

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

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

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

1. Создание эффективного метода расчёта конформационного равновесия органических соединений с явным учётом растворителя;

2. Разработка метода выявления и учёта (в свободной энергии связывания) различных сольватационных состояний комплекса белок-лиганд, находящихся в динамическом равновесии;

3. Направленный дизайн новых ингибиторов 8ук-киназы и поли(АДФ-рибоза)-полимеразы (РАКР) с использованием разработанных методов;

4. Синтез новых ингибиторов Бук-киназы и сопоставление экспериментальных данных по активности с результатами моделирования.

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

предложенного метода показало, что заселённости связывающихся с белком кон-формеров в серии ингибиторов могут существенно (в 10-50 раз) отличаться друг от друга, и, таким образом, вносить значительный вклад в итоговое значение свободной энергии связывания. При этом относительные доли конформеров, рассчитанные в вакууме, существенно отличались от заселённостей в воде, что указывало на значительно влияние растворителя на конформационное равновесие.

Впервые предложен метод выявления альтернативных сольватационных состояний комплекса белок-лиганд, основанный на определении корреляции сольватационных карт системы. Практическое применение метода показало, что прочность комплексов PARP-ингибитор, рассчитанная по альтернативным сольватационным состояниям, может отличаться на 4-5 кДж/моль, и, таким образом, вносить существенный вклад в точность расчётов. Стоит отметить, что стандартные критерии (временная зависимость координат атомов белка или энергии системы) не позволили определить изменения структуры воды в ходе расчётов исследованных систем.

Практическая значимость. На основе результатов, полученных с применением разработанных методов, был проведен направленный синтез новых ингибиторов Syk-киназы 2,4-диариламиноазинового ряда. Синтезированные соединения проявили ингибирующую активность в наномолярном диапазоне концентраций и являются перспективными кандидатами для терапии некоторых заболеваний, таких как ревматоидный артрит и неходжкинская лимфома.

Апробация работы. Основные результаты работы были представлены на следующих конференциях: V Российский симпозиум «Белки и пептиды» (Петрозаводск, Россия, 8-12 августа 2011), XX Российский Национальный Конгресс «Человек и лекарство» (Москва, Россия, 15-19 апреля 2013)

Публикации. По теме диссертации опубликовано 4 статьи в международных научных журналах и 2 тезиса в сборниках докладов научных конференций.

Структура и объем работы. Диссертационная работа состоит из введения, обзора литературы, описания материалов и методов исследования, результатов и их обсуждения, заключения, выводов и списка цитируемой литературы. Работа изложена на страницах машинописного текста, содержит рисунка и таблицы. Список цитируемой литературы включает 106 наименований.

Сокращения, принятые в тексте. ВСЭ — возмущение свободной энергии, КФ - конформационная фокусировка, PARP - поли(АДФ-рибоза)-полимераза (poly-(ADP-riboso)-poIymerase), Syk - селезёночная тирозинкиназа (spleen tyrosine kinase).

СОДЕРЖАНИЕ РАБОТЫ.

Во введении дана общая характеристика диссертационной работы, изложена актуальность темы, сформулированы цели и задачи исследования.

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

Во второй главе представлены материалы и методы исследования.

Химический синтез ингибиторов Бук-киназы. Синтез новых ингибиторов осуществлялся по следующей схеме:

| СООМе

ин.

ЫаН. НМР

HNOi.Hj.SO4

2, 86%

02М

III. 10% Р4С

3, 80%

н,м

4, 84%

хлмлд

ЫМР

С1

м^м^Ч н

Х=Н (5а. 81%) Х=С1 (5Ь, 63%)

МеО

.ОМе ОзБ-

МеО.

рН=10

.ОМе

МеО.

ю% ын,. н2о

р-с6н4зо3-

.ОМе

8 (30%)

С1

N Н

5а (Х=Н) 5Ь (Х=С1)

ОМе

у 2 еч *

ОМе

(Х=С1>. КМР"

МеО

Х=Н.У=0Мс (МТ-030, 85%) Х=Н,У=Н (МТ-ОЗО-ОМе, 88%) Х=С1,У=ОМс (МТ-030-С1, 56%) Х-1 ГЛ'-О] I (МТ-ОЗО-ОН, 95%)

а ЫМР - Ы-метилпирролидон

*В случае Х=С1 использовали 1 экв. анилина и 1 экв. триэтиламина

Пакеты программ, использованные в работе. Молекулярно-динамические расчёты, анализ траекторий и расчёты изменения свободной энергии методом ВСЭ проводились с помощью программного пакета Gromacs 4.5.4 [J. Chem. Theory Comput., 2008, 4, 3, 435-447]. Для построения трёхмерных карт плотности молекул воды использовалась программа gridcount [Proc. Natl. Acad. Sei. U. S. A., 2003, 100, 12, 7063-8]. Для квантово-механической оптимизации трёхмерных структур лиган-дов и определения их топологических параметров применялись программы ACPYPE [ВМС Res Notes, 2012, 5, 367] и AmberTools 1.5 [J. Comput. Chem., 2004, 25, 9, 115774]. Молекулярный докинг лигандов проводился с помощью программы Lead Finder [J. Chem. Inf. Model., 2008, 48,12,2371-85].

Расчёт относительной свободной энергии связывания. Начальные координаты комплексов белок-лиганд рассчитывались на основе трёхмерных структур мишеней из базы данных PDB [J. Mol. Biol., 1977, 112, 3, 535-42] с помощью метода молекулярного докинга. Относительная свободная энергия связывания ингибиторов Syk-киназы и PARP определялась на основании результатов переходов соответствующих лигандов друг в друга в воде и в белковом окружении. Изменение кулоновских и Ван-дер-Ваальсовых параметров в методе ВСЭ проводилось через 10 промежуточных состояний, что являлось достаточным для обеспечения необходимого перекрывания фазовых пространств между соседними состояниями. Длина траектории при расчёте переходов в воде составляла 2 не, в белке — 5 не для каждого X. Ошибка оценивалась методом блок-среднего по 5 блокам.

Для оценки степени перекрывания фазовых пространств рассчитывалась функция П [The Journal of Chemical Physics, 2004, 121, 18, 8742-8747]:

— WL

-Ш-1У

-л/2^.

(1)

Пл->0 —

где WlQx) - W-функция Ламберта, являющаяся решением w уравнения х = wex, М - количество точек, используемых для расчёта свободной энергии, и sA и sB — относительные энтропии состояния А и В.

Расчёт конформационного равновесия лигандов Syk-киназы. Конформацион-ные переходы в расчётах ВСЭ проводился плавным изменением потенциала двугранного угла в соответствии с формулой:

Vrbivuki) = k{ 1 + cos[n<p(p - (1 - Х)ср? - Хер!]), (2)

где к - силовая константа, п,р — кратность двугранного угла, <р£ и <р$ — 180 и -180°, Я - шаг перехода. Для расчёта энергии конформационного перехода использовался

следующий набор значений h от 0 до 0,45 и от 0,6 до 1,0 - с шагом 0,0125, от 0,45 до

0.6 - с шагом 0,05. Усреднение проводилось по траектории длиной 2 не.

Анализ сольватации комплекса белок-лиганд. Анализ проводился для соль-ватационных карт комплекса белок-лиганд, полученных усреднением позиций атомов кислорода молекул воды по траектории длиной 500 пс. Корреляция двух карт i и j (далее «корреляция карт плотности») рассчитывалась по формуле [J. Phys. Chem. В, 2011, 115,2,319-28]:

«(»*)-, сМ .

Jc(Pi,Pi)C(Pj.Pj)

где - ковариация между картами плотности Pi и р7 :

C{pi,pj) = ^ (<р;> - Pi(x,y,z)) ■ ((pj) - pj(x,y,z)), (4)

x,y,z

где p(jx, у, z) - плотность молекул воды в ячейке x,y,z. Аналогичным образом рассчитывалась последовательная корреляция карт плотности растворителя, полученных по / и i+1 участку траектории, а также двумерная матрица корреляции (i,j). Двумерная матрица подвергалась последующей кластеризации с выделением участков с высокой совместной корреляцией в отдельные блоки.

Экспериментальное определение связывания ингибиторов с PARP и Syk-киназой. Экспериментальное определение ингибиторной способности соединений по отношению к Syk-киназе проводилось по методике с у-33Р меченым АТФ (суммарная концентраций АТФ 10 мкМ). Определение ингибиторной способности ингибиторов PARP было выполнено в Институте химической биологии и фундаментальной медицины СО РАН) в соответствии с методикой [Mendeleev Commun., 2012, 22,

1, 15-17].

В третьей главе представлены результаты исследования и их обсуждение.

Учёт конформационного равновесия для ингибиторов Syk-киназы Теоретические основы учёта конформационного равновесия при образовании комплекса белок-лиганд. Для лиганда L, существующего в растворе в виде двух конформеров Li и L2, из которых с белком связывается только конформер (т. е. для связывания L необходима конформационная фокусировка (КФ) на конформере Lx), концентрация связывающегося с белком конформера L1 составляет:

[Ll] = [¿i][+k][L]tot = 77MlL]tot' (5)

1 + [¿i]

где [L]tot - общая концентрация лиганда в растворе, равная [LJ + [L2]- Если в ходе моделирования методом молекулярной механики конформеры и L2 не переходят друг в друга, то при расчёте методом ВСЭ (например, при расчёте абсолютной свободной энергии связывания) полученное значение AG соответствует константе равновесия связывания белка (Р) и конформера [¿г]:

[PLt]

ACfgJS = -ЯГ In;

(6)

'МП1 '"[р]^]'

в то время как экспериментально определяемому значению ДС^^ соответствует не равновесная концентрация [Ь^, а общая концентрация лиганда [¿]Со£:

KZ = -RTin-[m]t

(7)

Таким образом, поправка, вносимая фокусировкой на конформере в данном случае составляет:

AGb = AGebZ - AGbind = = RT In

[l2]

(»Eft

(8)

Величина -—- представляет собой относительную заселённость конформера Ь2 и определяется на основании расчётов методом ВСЭ:

(ДСЬ2-ДСЬ1)>

[L2] ( (AGL2-AGLI)\

й=еХР1--Rf-}

(9)

где ДС^2 — ДС£1 — свободная энергия перевода конформера Ь2 в рассчитываемая методом ВСЭ. Если лиганд Ь имеет более двух конформеров в растворе (обозначим их ¿1( 12, ...,Ьп), -то итоговое значение имеет вид:

Д G£f = Я7Чп 1 +

) . (10)

protein

Белок

Вода

Полученные поправки справедливы для расчётов абсолютной свободной энергии связывания соединения с белком. Тем не менее, при расчётах методом ВСЭ зачастую определяется не абсолютная, а относительная энергия связывания (т.е. разность свободных энергий связывания близких по структуре лигандов и V, см. рис. 1). Аналогичные выкладки в данном случае при-

Рис. 1. Термодинамический цикл для расчёта относительной свободной энергии связыва- водят к следующей формуле: ния с учётом КФ.

AGbind ~ Mbind = Krotei ~ + ^CF ~ ^tF . (11)

где AGprotein 11 1 ~ свободная энергий перевода конформера Lx в L\ в белке и

в воде, соответственно, а ДGcF и Д— свободные энергии конформационной фокусировки для L' и L.

Ситуация, в которой один из лигандов может связываться с белком более чем в одной конформации (обозначим такие конформации L\,L'2, — ,L'm где ш<п), анализируется аналогично, и приводит к выражению:

.ги l dti V ( AGbinding ~ AGbinding + AGcf \ ДСbinding - binding = ~RT IП ^ exp I--—- I. (12)

Данное выражение вырождается в формулу (11), если с белком связывается только один конформер L\.

Таким образом, для моделирования конформационного равновесия при связывании лиганда с белком необходимо определить заселённость всех конформаций L1,L2,—,Ln. Для решения данной задачи был разработан специальный протокол ВСЭ, описанный ниже.

Протокол расчёта свободной энергии конформационного перехода методом ВСЭ. При расчётах методом ВСЭ проводится модификация одного или нескольких параметров системы в соответствии с шагом X. Для конформационного перехода координатой реакции является значение двугранного угла (точнее, 4 смежных двугранных углов). В программе Gromacs изменение потенциала двугранного угла в зависимости от шага перехода X описывается формулой:

v((p) = kd{ 1 + cos[n,p(p - (1 - A)q>s ~ A<p?]), (13)

где kd — силовая константа на данном двугранном угле, п1р - кратность функции, <р$ и <Ps — значение равновесных двугранных углов для состояния А и В. Конформаци-онный переход достигается выбором <р$ = 180° и <pf = —180°. При данных параметрах начальный (1=0) и конечный (1=1) потенциалы на двугранном угле совпадают, однако при промежуточных значениях X значение равновесного угла изменяется, что обеспечивает переход от одного конформера к другому. Пример изменения потенциала двугранного угла, а также геометрии молекулы при переводе цис-изомера N-метилформамида в транс-изомер показаны на рис. 2. Переход осуществляется через невыгодное пирамидальное переходное состояние, однако, поскольку свободная энергия является функцией состояния, на корректность расчётов это не влияет.

Рис. 2. Схема перехода конформеров Ы-метилформамида (цис в транс). Показан схематический вид потенциальной функции двугранного угла при X = 0;1 (сплошная линия) и Х=0.5.

Валидация метода на примере кон-формационного равновесия N-метилформамида. Для валидации разработанной методики была выбрана модельная молекула N-метилформамида, для которой известны экспериментальные данные по соотношению цис:транс конформеров в воде. Результаты расчётов свидетельствуют об энергии перехода Д^транс-щис="6Л7±0,04 кДж/моль. Согласно экспериментальным данным, в воде отношение концентраций цис:транс конформеров составляет 8:92, что эквивалентно энергии перехода -6,05 кДж/моль при 25°С. Таким образом, полученные расчётные данные хорошо согласуются с данными эксперимента.

Поскольку моделирование конформационных переходов методом ВСЭ ранее описано не было, для определения корректности расчётов дополнительно изучались некоторые характеристики системы в процессе вычислений. AG перехода стабилизировалось на 500 пс траектории, что указывало на временную сходимость расчётов при использовании выбранного метода и на корректность итогового значения AG, рассчитанного по траектории длиной 2000 пс.

Зависимость изменения кумулятивной свободной энергии конформационного перехода N-метилформамида от координаты X представлена на рис. 3. Изменение свободной энергии на каждом шаге по X не превышало 1-2 кДж/моль, что косвенно указывало на достаточно большое перекрывание фазовых пространств между соседними состояниями. Отсутствие максимума при Х=0,5 (при «перпендикулярном» расположении заместителей) на данном графике объяснялось постоянным смещением минимума потенциала двугранного угла (см. рис. 2) от -180° к 180°, а при 1=0,5 равновесный угол составлял 90°. Выгодность пирамидальной конформации вследствие подобного изменения потенциала физически

Рис. 3. Изменение кумулятивной свободной энергии при переходе транс-Ы-метилформамида в цис-изомер в зависимости от шага X.

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

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

Конформационная фокусировка ингибиторов Бук-киназы. Исследуемые в данной работе ингибиторы Бук-киназы имеют 16 заторможенный конформаций (рис. 4). Для обозначения конформеров использовалась бинарная номенклатура, в соответствии с которой каждый из конформеров определялся 4 цифрами «0-1» в зависимости от значения конкретного двугранного угла. При этом «0» соответствовал конформации, реализующейся при связывании с белком, «1» - второй возможной конформации двугранного угла (конформер, связывающийся с белком, обозначался «0000»). Перпендикулярное расположение циклов, реализующееся в процессе перехода, обозначалось <А>. Пример подобного обозначения конформеров представлен на рис. 5 и 6.

Я406 МТ-030

МТ-050 МТ-037

Рис. 4. Структуры изученных ингибиторов 8ук-киназы и двугранных углов, вращение вокруг которых затруднено.

ф,1- {- Ф4

0 О

о/ К

н.мЛ А^н I/X Ал."

NN N /^М N<7 N N

(Ь 6 ' '

„/ р

оооо

0100

1010

Ф, = Ф2=Фз=Ф4=я Ф1 = Ф: = Ф3=Я, ф4=-л ф,= ф4= л, ф, = ф3=-л

Рис. 5. Обозначение конформеров ингибиторов Бук-киназы. «0000» соответствует связывающейся с белком форме.

Для расчёта итоговой свободной энергии конформационной фокусировки были определены значения свободной энергии перехода каждого из конформеров в кон-

0000

г? О7

1000

о о

31=0 31=0.5 /=1

Рис. 6. Конформационный переход «0000» в «1000» для ингибитора МТ-037. формер «0000». Изменение конфигурации одного из двух крайних двугранных углов (ф! или Ф4) молекулы (переходы ОхуО—>1ху0 или ОхуО—>0ху1) по отдельности оказалось невозможным ввиду стерических затруднений, возникавших из-за перпендикулярной ориентации колец при 1=0,5. Тем не менее, данные переходы могли быть проведены при одновременном изменении ф! или ф4, которое позволяет избежать невыгодной перпендикулярной ориентации колец. Таким образом, для перехода в кон-

формеры 1ху1 разумным являлся одновременный поворот обоих двугранных углов Ф1 и (р4, в то время как переход в конформеры 1ху0 и 0ху1 осуществлялся дополнительным вращением второго кольца сначала в одну (при Х<0,5), а затем в другую (при Х>0,5) сторону (рис. 6), что в итоге не влияло на его конечную конформацию.

Итоговые значения относительной заселённости конформеров исследованных соединений в воде приведены в таблице 1 (заселённость связывающегося конформе-ра 0000 принималась равной 1).

Таблица 1. Относительные заселённости конформеров исследованных ингибиторов Яук-кипазы в воде, а также энергии их конформационной фокусировки (кДж/моль). Жирным шрифтом и курсивом выделены группы конформеров, для которых наблюдалось наиболее существенное отличие в заселённости между различными соединениями.

Конформер Относительная заселённость

Я406 МТ-030 МТ-050 МТ-037а)

0000 1.00 1.00 1.00 1.00

0001 1.31 1.50 2.75 0.54

0010 0.73 0.61 5.97 12.4

ООН 0.95 0.63 11.67 3.64

0100 0.01 4.14 0.72 12.55

0101 0.00 3.26 0.41 2.49

0110 0.01 3.48 1.39 79.93

0111 0.00 4.75 1.11 34.12

1000 0.05 0.00 0.53 0.00

1001 0.05 0.00 0.46 0.00

1010 0.02 0.00 6.17 0.01

1011 0.02 0.00 7.46 0.00

1100 0.00 0.00 0.00 0.01

1101 0.00 0.00 0.00 0.00

1110 0.00 0.00 0.00 0.05

1111 0.00 0.00 0.00 0.01

ДвсГ, кДж/моль 1,8±0,5 5,6±0,5 7,4±1,1 12,4±0,6 (0) 14,0±0,7 (1)

' Соединение МТ-037 может связываться с белком в конформациях «0000» (0) и «0001»(1)

Среди полученных значений существенное отличие наблюдалось между заселённостями конформеров 01ху для соединения Я406 и остальных соединений (таблица 1, выделено жирным). Данная группа конформеров являлась крайне невыгодной (относительная концентрация <0.01) для 11406, в то время как для других ингибиторов её относительная заселён-

Рис. 7.

11406,0100

Стерическое отталкива-

ние между атомами Н и Р.

(?406,1000 Рис. 8. Внутримолекулярное электростатическое отталкивание атомов азота

ность значительно превышала 1. Предположительное объяснение данного отличия состояло в стерическом отталкивании атомов Н и F в конформациях 01хх соединения R406, которое отсутствовало в других соединениях из-за замены группы C-F на атом азота (рис. 7).

Другая группа конформеров Юху имела достаточно высокую заселённость для МТ-050, однако являлась практически не заселённой для соединений R406, МТ-030 и МТ-037 (таблица 1, выделено курсивом). Невыгодность данной группы конформеров для R406, МТ-030 и МТ-037 связана, по-видимому, с внутримолекулярным отталкиванием отрицательно заряженных атомов азота пиридинового и пиримидинового или триазинового циклов (рис. 8). В соединении МТ-050 центральный атом азота заменён на группу СН, что исключало данное отталкивание в конформере Юху и делало его достаточно выгодным.

Для оценки влияния растворителя на конформаци-онное равновесие проводились аналогичные расчёты конформационной фокусировки в вакууме. Наиболее значительное отличие от результатов в воде наблюдалось для соединения МТ-050: конформеры «0010», «0011», «0110», «0111», «1010», «1011» были достаточно выгодны в воде, однако оказались совершенно не заселены в вакууме. Данные конформеры отличаются от конформера «0000» поворотом вокруг угла фз, в результате которого исчезает выгодное внутримолекулярное взаимодействие между отрицательно заряженным атомом азота пиримидинового цикла и положительно заряженным полярным атомом водорода (рис. 9). В воде данное взаимодействие проявлялось в значительно меньшей степени из-за полярной сольватации как атома азота, так и NH-группы. Это приводило к меньшему

изменению энергии при пере-n'-=n - воде конформеров «0010»,

Н,.АА.. «0011», «0110», «0111»,

«1010», «1011» в конформер «0000» в воде и, как следствие, к их большей заселённости. Для остальных соединений подобного эффекта не наблюдалось, поскольку в них присут-

Рис. 9. Внутримолекулярное электростатическое притяжение атомов азота и водорода группы >Ш.

ствовал дополнительный атом азота в триазиновом цикле, который компенсировал потерю данного взаимодействия в конформере «0010».

Для валидации корректности проведённых расчётов исследовалась степень перекрывания фазовых пространств и зависимость величины свободной энергии переходов от шага X и от времени расчёта. Степень перекрывания фазовых пространств, рассчитанная на основе относительных энтропий вд и Бв, была достаточной, так как значение параметра П было существенно положительным. Зависимость свободной энергии от времени расчётов стабилизировалась примерно к 500 пс, что указывало на корректность выбранной длины траектории 2000 пс.

Исходя из структур исследованных ингибиторов (рис. 4), для всех соединений, кроме МТ-037, обе конформации правого цикла (угол <р4) были идентичны ввиду симметрии замещения бензольного кольца, т.е. для этих соединений конформеры хугО и хуг1 фактически совпадают. Согласно результатам проведённых расчётов, заселённости конформеров хугО и хуг1 для веществ 11406, МТ-030 и МТ-050 были достаточно близки, что соответствовало теоретическим ожиданиям.

Таким образом, в результате использования разработанной методологии была продемонстрирована корректность метода и возможность его применения для расчётов свободной энергии конформационной фокусировки ингибиторов 8ук-киназы.

Расчёт свободной энергии связывания конформера «0000» ингибиторов Бук-киназы. Для определения относительной свободной энергии связывания изученных соединений проводились их взаимные переходы в воде и в белке для ингибиторов Зук-киназы в конформации «0000». Полученные результаты представлены на рис. 10.

Как было отмечено ранее, для всех соединений, кроме МТ-037, конформеры «0000» и «0001» совпадали вследствие симметричности замещения бензольного кольца. В то же время, для соединения МТ-037 конформеры «0000» и «0001» отличались друг от друга и, поскольку эти отличия незначительны и достаточно далеки от критичных для связывания водородных связей, разумно предположить, что оба конформера могут связываться с ферментом. Поэтому для МТ-037 в белке рассчитывались два перехода: МТ-030 в конформер «0000» и МТ-030 в конформер «0001» соединения МТ-037.

Расчёт перекрывания фазовых пространств, а также анализ зависимости кумулятивной свободной энергии от времени и шага X свидетельствовали о релевантности полученных значений.

Расчёт относительной свободной энергии связывания ингибиторов Бук-киназы. Для расчёта относительной свободной энергии связывания соединений

11406, МТ-030 и МТ-050 с Яук-киназой использовалась формула (11), для расчёта энергии связывания МТ-037 - формула (12), а также результаты конформационной фокусировки (таблица 1) и энергий связывания конформера «0000». Итоговые значения относительной энергии связывания отражены на рис. 10.

Для перехода Ы406-+МТ-030 дополнительный учёт вклада КФ позволил существенно приблизить расчетное значение к экспериментальному (с 1,96 до 5,76 кДж/моль при экспериментальном значении 7,4 кДж/моль). Причина большого отклонения значения, рассчитанного без учёта КФ, состоит в меньшей заселённости конформера 0000 для МТ-030 по сравнению с Я406. По результатам расчётов без учёта КФ, данные соедине-

Рис. 10. Относительная свободная энергия связывания ния связываются с белком С ингибиторов вук-киназы, кДж/моль (Эксперименталь- близкой прочностью, однако ные данные, расчётные данные, расчётные данные с введение атома фтора снижаех учётом конформационнои фокусировки выделены курси-

вому заселённость конформеров 01ху,

увеличивая тем самым относительную концентрацию конформера «0000».

Связывание соединений МТ-050 и МТ-037 происходило слишком слабо, что не позволило определить экспериментальные значения ДДОьы- Тем не менее, возможна оценка нижней границы относительной энергии связывания этих веществ по сравнению с МТ-030 (она составляет 6,8 кДж/моль). Как для МТ-050, так и для МТ-037 оба рассчитанных значения относительной энергии связывания превышали данную оценку снизу, что не позволяло выявить влияние конформационной фокусировки для этих соединений.

Полученные значения свободной энергии указывают на то, что примерно половина молекул 11406 находится в связывающейся с белком конформации, в то время как для МТ-037 доля молекул в связывающейся конформации составляет лишь около 1%. Для соединения МТ-030 относительная концентрация связывающегося кон-

х хса.

О N N

ын ^ ын

Я406

О'

Д.. о.

7,4

1,96±0,93 5,76+1,4

О' N N N4 N N4 ' Н

мт-озо

11,0210,46. 13,2±1,5 / >6,8

15,37±0,33 22,1±1,2

Г1 '

О' N N N4 Н

МТ-050

О'

А о

(1 " N4 ' О

П г

О ' 'ГГ N N11 N NN

МТ-037

формера равняется 5%, и, таким образом, повышение его заселённости может являться предметом для дальнейшей оптимизации данного соединения.

Анализ временной стабильности сольватации в комплексах РАЫР с ингибиторами

Корреляция карт плотности воды для тестовых систем. Величина корреляции Я, рассчитываемая по формуле (3), может меняться от 0 (отсутствие корреляции) до 1 (при полной корреляции молекул воды). В качестве примера системы с Я~0 можно привести не центрированный куб воды, в котором небольшие локальные флуктуации плотности в ходе молекулярно-динамической траектории случайны и отличаются для разных участков траектории. В то же время, для центрированной по одной из молекул воды траектории карты плотности 4 её соседних молекул являются высоко скореллированными ввиду образования направленных водородных связей.

Таким образом, выбранная в качестве критерия временная корреляция карт

плотности растворителя являлась чувствительным параметром, отражавшим временную стабильность локального распределения молекул воды. Значения корреляции, близкие к 1, соответствовали наличию выраженных и стабильных во времени максимумов на картах плотности, в то время как низкие степени корреляции указывали на нестабильность положения максимумов, либо на их отсутствие. Корреляция является безразмерной величиной, находящейся в пределах от 0 до 1, и позволяет сравнивать степень схо-

Анализ карт плотности воды для комплексов РА ЯР с ингибиторами. Для

расчётов относительной свободной энергии связывания ингибиторов РАМ" методом ВСЭ использовалась схема переходов, отображённая на рис. 11. Поскольку в методе ВСЭ расчёты проводятся с предварительно уравновешенной структурой белка, для изучения динамики изменения сольватации по мере стабилизации структуры данный параметр определялся для подготовительной стадии (отжиг и подготовительная динамика), а также для продуктивной динамики расчётов (переход

Рис. 11. Относительные свободные энергии переходов исследованных ингибиторов РАЯР, кДж/моль. Жирным шрифтом выделены экспериментальные значения.

жести различных вариантов сольватации.

0,8

0,6 £

0,4 0,2

Отжиг Подготовительная РасчётВСЭ динамика

I I

PARPCl—»-PARPH, заряд). Зависимость корреляции от времени приведена на рис. 12.

В ходе предварительных расчётов корреляция колебалась в пределах 0,55-0,75, что указывало на наличие выраженных, однако нестабильных максимумов плотности. При переходе к расчётам методом ВСЭ (при которых, в частности, фиксируется остов полипептидной цепи) корреляция значительно возросла и превысила 0,9, что говорило о значительной стабилизации позиций ближайших к лиганду молекул растворителя в комплексе. Таким образом, выбранные параметры определения корреляции являлись корректными, и чувствительными к небольшим изменениям в окружающей лиганд сольватной оболочке. Стабилизация сольватации указывает на отсутствие низкочастотных колебаний в белке для данного перехода.

Для определения стабильности окружающих лиганд молекул растворителя в ходе расчётов методом ВСЭ временная корреляция определялась для всех переходов, отражённых на рис.11. Примеры корреляции для переходов PARP_C1—>PARP_H (Ван-дер-Ваальсовы параметры) и PARP_N3->PARP_C1 (куло-новские параметры) приведены на рис. 13.

! , PARP_CI-)PARP Н (Ван-дер-Ваальс) 1 - PARP N,->PARP CI (заряд)

О 2000 4000 6000 8000 10000 1, ПС

Рис. 12. Корреляция плотности ближайших 15 молекул воды (расстояние до ~4А) вокруг лиганда РАКР_Н в ходе подготовки (отжиг и препаративная динамика) и расчётов методом ВСЭ (первый шаг при модификации зарядов в переходе РАШМЧНг—^►РАЯР_С1). Усреднение производилось по 500 пс.

0 10 20 30 40 50 нс 0 10 20 30 40 50 i, нс

Рис. 13. Временная корреляция последовательных карт плотности, усреднённых по 500 пс. Изменение шага X показано точками. Для переходов РА11Р_С1->РА11Р_Н, РАИРМН2->РАКРС1 и РАИР_ОН—>РАЯР_С1 временная зависимость корреляции являлась достаточно рав-

номерной, а Я менялся в пределах 0,7-1, что говорит о стабильности сольватации в ходе расчётов и о её идентичности для разных шагов А.. КМ 8С) боковых радикалов изменялось незначительно и составляло порядка 1,2-1,6 А. Наибольшие значения К МЯО наблюдались между соседними шагами по А,, что явилось ожидаемым результатом.

В то же время, для перехода РАКР_Ыз—»РАКР_С1 корреляция была значительно менее равномерной, и в некоторых случаях составляла менее 0,5-0,6, что указывало на низкую временную стабильность позиций молекул растворителя. Тем не менее, зависимость КМЯГ) от времени оказалась достаточно гладкой и не имела значимых отличий от зависимостей для других переходов. Причина нестабильности сольватации для данных переходов исследовалась отдельно.

Исследование особенностей сольватации комплекса ингибитора РАЯР_N3 с белком. Сравнение участков с низкой корреляцией (11=0,39, 0-500 и 500-1000 пс траектории перехода РАММЧз—>РАКР_С1 с модификацией зарядовых параметров, шаг /.=0,2) показало, что вблизи цик-лопентанового и тиофенового колец карты плотности воды были практически идентичны, в то время как основное отличие наблюдалось у азидной группы (рис. 14). Анализ молекулярно-динамической траектории показал, что на протяжении всего участка 0-500 пс азидная группа находилась в одном положении (показано синим на рис. 14), в то время как на участке 500-1000 пс она переходила в другую конформа-цию (зелёный, рис. 14). Таким образом, в данном случае низкая степень корреляции карт плотности объяснялась медленным поворотом азидной группы, происходившим раз в несколько наносекунд и приводившим к существенному изменению позиций ближайших к лиганду молекул воды. Для остальных соединений подобного явления не наблюдалось, поскольку небольшой размер групп -СН2С1, -СНтМЬЬ и -СН1ОН позволял им свободно вращаться и не приводил к возникновению нескольких трудно переходящих друг в друга сольватационных состояний.

Из полученных данных последовательной корреляции карт плотности можно сделать вывод о наличии нескольких альтернативных сольватационных карт, однако

Рис. 14. Карты плотности молекул воды и положение азидной группы для участков 0-500 (синий) и 5001000 (зелёный) пс траектории перехода РАЮМЧз^РАЯР а по зарядам, шаг ^=0,2.

нельзя определить их количество. Для определения количества альтернативных вариантов сольватации строилась двухмерная матрица корреляции карт плотности. В результате последующей кластеризации была получена матрица, на которой достаточно чётко наблюдалось присутствие двух вариантов сольватации с существенным преобладанием (-10:1) одного из них (рис. 15, слева). Для сравнения можно привести аналогичную двухмерную матрицу корреляции для перехода

РА11Р_С1—>РАЯР_Н по зарядовым параметрам (рис. 15, справа). В данном случае корреляция между всеми кадрами всех траекторий достаточно высока, и очевидно, что сольватационная карта на каждом из участков в 500 пс одинакова.

Таким образом, для точного расчёта свободной энергии перехода РАЯР_^—>РАКР_С1 необходим учёт наличия альтернативных вариантов сольватации системы. Продолжительности полученных траекторий (5 не) являлись недостаточными, так как количество переходов сольватационных состояний было слишком низким и не соответствовало установлению равновесию. Для получения равновесного распределения молекул воды в комплексе продолжительность расчётов была увеличена до 40 не, что позволило добиться адекватного количества переходов между сольватационными состояниями. Результаты кластеризации полученных карт плотностей молекул воды приведены на рис. 16. В обоих случаях двухмерные матрицы сольватации разбивались на два кластера, как и в случае с более короткой траекторией. Разница в свободной энергии перехода между двумя различными состояниями сольватации комплекса составила более 2 кДж/моль для модификации зарядовых параметров и более 4 кДж/моль для Ван-дер-Ваальсовых параметров, что значительно превышало погрешность рассчитанных величин, и, таким образом, являлось статистически значимым отклонением.

Учёт альтернативных сольватационных состояний при расчётах методом ВСЭ. Проведённый анализ альтернативных сольватационных состояний комплекса

Рис. 15. Кластеризованная матрица корреляции карт плотности молекул воды при переходе Р АII Р_1М3 —► Р А К Р_С I (слева) и РАКР_С1->РАЯР_Н по заряду (справа).

Рис. 16. Кластеризованная матрица корреляции карт плотности молекул воды при переходе РАЯРЫз—»РАИРО по зарядовым (слева) и Ван-дер-Ваальсовым (справа) параметрам.

РАИР Ыз показал, что различное расположение молекул воды может вносить существенный вклад в итоговое значение свободной энергии перехода. Таким образом, для определения итогового значения ААО необходим учёт наличия обоих вариантов сольватации с помощью термодинамического цикла, приведённого на рис. 17 (ЕХ1 и ЕХ2 - альтернативные варианты сольватации комплекса, ЕХ - равновесное состояние в растворе, 01, а2, Р1 и р2 - заселённости соответствующих сольва-тационных состояний). Доли каждого из сольватационных состояний для лигандов А и В определяются на основании длительной молекулярной динамики, значения свободных энергий перехода для каждого из них рассчитываются

а, ДС"'-ГЕВ' р. ' ыпа *

ДСЕА2-ЕВ2 о

Ыпс!

Рис. 17. Термодинамический цикл для учёта наличия альтернативных вариантов сольватации

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

— ДСьАсг —

-ЯТ 1п ос! + ДС^ЕВ1 + ЯТ 1п рг

(14)

а=0,50 Р,=0,87

-7 ^п+п ая

Аналогичное выражение получается и на основании участков траектории со вторым типом сольватации. Замкнутость центрального термодинамического цикла на рис. 17, с одной стороны, является дополнительным критерием для оценки корректности расчётов, а с другой — позволяет рассчитать доли альтернативных сольватационных состояний в случае, если одно из них для одной из систем (ЕА или ЕВ) является крайне невыгодной и не присутствует на траектории. На основании проведённых расчётов для модификации зарядовых и Ван-дер-Ваальсовых параметров при переходе РАММчГ3—►РАЯРа были построены соответствующие термодинамические циклы (для Ван-дер-Ваальсовых параметров -на рис. 18). Анализ данных показал, что рассчитанные заселённости находились в хорошем соот-

-3,56±0,98 а2=0,50 Р2=0,13

Рис. 18. Дв (кДж/моль) и заселённости различных сольватационных состояний для модификации Ван-дер-Ваальсовых параметров при переходе РАЯР_Нз—>РАКР_С1. Стрелкой показана суммарная свободная энергия по циклу.

ветствии со значениями свободной энергии переходов, что выражалось в практически точном равенстве нулю свободной энергии по циклу (с учётом погрешности, энергия по циклу строго равна 0 кДж/моль). Итоговые свободные энергии переходов составили 2,53±0,17 кДж/моль для модификации зарядовых параметров и -6,12±0,98 кДж/моль для модификации Ван-дер-Ваальсовых параметров.

Дизайн новых иингибиторов Бук-киназы

В предыдущих разделах были предложены и валидированы способы дополнительного увеличения точности расчётов свободной энергии связывания методом ВСЭ за счёт явного учёта конформационной фокусировки лигандов в растворе и анализа сольватации комплекса белок-лиганд. Применение этих методов позволило рассчитать свободные энергии связывания новых, потенциально более эффективных ингибиторов Бук-киназы и выбрать направления дальнейшей модификации соединения МТ-030.

Выбор модификаций структуры МТ-030. На основании результатов молекулярного моделирования и анализа синтетической доступности был предложен ряд потенциально активных аналогов МТ-030: МТ-ОЗО-ОМе и МТ-ОЗО-ОН (замена 4-метоксигруппы на атом водорода или гидроксильную группу), МТ-030-Ме (введение метильной группы в 3-е положение пиридинового цикла), МТ-030-С1 (добавление атома хлора в 6-ое положение триазинового цикла). Опираясь на результаты молекулярного докинга, удаление метокси-группы (МТ-ОЗО-ОМе) не должно было существенно повлиять на прочность связывания, так как данная группа не образовывала прочных контактов с остатками белка. Усиление связывания для соединения МТ-030-С1 предполагалось за счёт образования дополнительной «галогенной» связи (не-ковалентное взаимодействие С=О...На1) между атомом хлора и атомом кислорода карбонильной группы С1и449. Наконец, замена 4-метоксигруппы на гидроксильную (МТ-ОЗО-ОН) была предложена для улучшения растворимости, и, по аналогии с МТ-ОЗО-ОМе, также не должна была существенно ухудшить связывание. МТ-030-Ме должен был обладать меньшей энергий конформационной фокусировки ввиду невыгодности конформера «1000».

Расчёт прочности связывания новых ингибиторов Бук-киназы методом ВСЭ. Проведённые приблизительные оценки прочности связывания соединения были проверены с помощью разработанного протокола ВСЭ. Предсказанные значения свободной энергии связывания (относительно МТ-030), а также свободные энергии конформационной фокусировки приведены на рис. 19. Наиболее перспективным на основании проведённых расчётов являлось соединение МТ-ОЗО-ОН. Для МТ-030-Ме рост заселённости конформации Юху ввиду отталкивания атомов водорода и ме-

МТ-030-С1 ЛС=2,б кДж/моль ДД0Ь1ак=12,6 кДж/моль ДЛСЬеч>9,2 кДж/моль

N ^

„АЛ^и

N ^

•ЧАЛи'»*

мт-озо

Дй =5,6 кДж/моль

МТ-030-Ме ДС;)=7,2 кДж/моль ДАй =12,6 кДж/моль

АЛ.,"

тильной группы привёл к повышению энергии конформационной фокусировки, и поэтому его синтез не проводился.

В случае соединения МТ-030-С1 и МТ-ОЗО-ОМе расчёт методом ВСЭ опроверг предложенные изначально структурные гипотезы. Так, атом хлора МТ-030-С1 находился слишком близко к атому кислорода белка, что приводило к существенному Ван-дер-Ваальсовому и электростатическому отталкиванию и ослаблению связывания. Удаление метокси-группы также привело к снижению прочности комплекса за счёт потери выгодных Ван-дер-Ваальсовых контактов с остатком Ьеи377. Тем не менее, замена метальной группы на атом водорода в соединении МТ-ОЗО-ОН, не изменяющая количества водородных связей, приводила к незначительному снижению прочности связывания, и, таким образом, могла рассматриваться как возможный способ увеличения растворимости соединения.

Экспериментальные значения активности синтезированных соединений согласовались с предсказанными (отличие энергии связывания составило 2 кДж/моль): МТ-ОЗО-ОН проявило несколько меньшую активность, чем МТ-030, в то время как МТ-ОЗО-ОМе и МТ-030-С1 оказались неактивны в изученных концентрациях, как и предсказывали расчёты. Результаты, полученные методом молекулярного докинга, не согласовывались с экспериментом. Соединение МТ-030, биологическая активность которого была исследована на клеточных линиях неходжкинской лимфомы КАМАЬ\\^А, РЗНЗ и РЗН11-1, эффективно ингибировало их пролиферацию в суб-микромолярной концентрации.

МТ-ОЗО-ОМе ДО^г.б кДж/моль

ДДСЬса|[=7,8 кДж/моль ^ >

Ь.схр

МТ-ОЗО-ОН ДСс[=1,4 кДж/моль ДДСЬык=5,6 кДж/моль ААй =3,7 кДж/моль

Рис. 19. Предсказанные и экспериментальные значения относительной свободной энергии связывания (относительно МТ-030) и свободной энергии конформационной фокусировки для новых ингибиторов 8ук-киназы.

ОСНОВНЫЕ РЕЗУЛЬТАТЫ И ВЫВОДЫ.

1. Впервые предложен способ расчёта энергии заторможенных конформацион-ных переходов в органических молекулах, основанный на использовании двугранного угла в качестве реакционной координаты в методе ВСЭ; метод позволяет с высокой точностью (~1 кДж/моль) рассчитывать конформационное равновесие биологически активных соединений с явным учётом растворителя;

2. Исследована заселённость конформеров впервые синтезированных биологически активных соединений, производных 2,4-диариламиноазинового ряда - ингибиторов Syk-киназы; показано, что в ряду изученных соединений заселённость конформеров отличается в 10-50 раз, что соответствует вкладу в свободную энергию связывания 4-6 кДж/моль;

3. Разработан метод учёта сольватации комплекса белок-лиганд, основанный на анализе временной корреляции карт плотности растворителя, полученных из мо-лекулярно-динамических траекторий, который позволяет выявлять альтернативные сольватационные состояния комплекса белок-лиганд, не идентифицируемые при помощи стандартных критериев (временная зависимость координат атомов белка или энергии системы);

4. Показано качественное изменение сольватации комплекса RARP-лиганд при переходе от R=-H, -СН2С1, -СН2ОН, -CH2NH2 к R=-CH2N3 в ряду синтезированных производных 3,5,6,7-тетрагидро-4#-циклопента[4,5]тиено[2,3-*/]пиримидин-4-онов с помощью разработанного метода анализа корреляции карт плотности растворителя; в случае производного R=-CH2N3 было обнаружено наличие альтернативных сольвата-ционных состояний комплекса, отличие в свободной энергии которых составляло 4 кДж/моль;

5. Проведён направленный синтез четырёх новых производных 2,4-диаминотриазина, ингибирующая активность которых в отношении Syk-киназы хорошо согласуется с разработанной моделью, включающей учёт конформационного равновесия лиганда и сольватацию комплекса белок-лиганд (отличие рассчитанной свободной энергии связывания от экспериментальной— 2 кДж/моль); показана биологическая активность соединения МТ-030 на клеточных линиях неходжкинской лим-фомы;

СПИСОК ПУБЛИКАЦИЙ

1. Zeifman, A.A., Stroylov, V.S., Novikov, F.N., Stroganov, O.V., Kulkov, V., Chilov, G.G. Alchemical FEP Calculations of Ligand Conformer Focusing in Explicit Solvent. J. Chem. Theory Comput., 2013,9 (2), pp 1093-1102

2. Rakitina, T.V., Zeifman, A.A., Titov, I.Y., Svitanko, I.V., Lipkin, A.V., Stroylov, V.S., Stroganov, O.V., Novikov, F.N., Chilov, G.G. Efficacy of novel Syk-kinase inhibitors MT-SYK-03 and MT-SYK-322 in cellular models of autoimmunity and cancer. Mendeleev Communications, 2012,22 (6), pp 287-287

3. Zeifman, A.A., Titov, I.Y., Svitanko, I.V., Rakitina, T.V., Lipkin, A.V., Stroylov, V.S., Stroganov, O.V., Novikov, F. N., Chilov, G.G. Rational design and synthesis of novel Syk-kinase inhibitors. Mendeleev Communications, 2012, 22 (2), pp 73-74

4. Romashov, L.V., Zeifman, A.A., Zakharenko, A.L., Novikov, F.N., Stroylov, V.S., Stroganov, O.V., Chilov, G.G., Khodyreva, S.N., Lavrik, O.I., Titov, I.Y., Svitanko, I.V. Rational design and synthesis of new PARP1 inhibitors. Mendeleev Communications, 2012, 22 (1), pp 15-17

5. Зейфман A.A., Стройлов B.C., Строганов O.B., Новиков Ф.Н., Чилов Г.Г. Метод явного учёта конформационной фокусировки лиганда при связывании с белком-мишенью. XX Российский Национальный Конгресс «Человек и лекарство» (Москва, Россия, 15-19 апреля 2013), с. 343.

6. Зейфман А.А., Новиков Ф.Н., Строганов О.В., Стройлов B.C., Чилов Г.Г. Методы явного учёта молекул воды для предсказания связывания низкомолекулярных лигандов с белками. V Российский симпозиум «Белки и пептиды» (Петрозаводск, Россия, 8-12 августа 2011), с. 321.

Заказ № 124-П/05/2013 Подписано в печать 29.04.2013 Тираж 120 экз. Усл. п.л. 1,2 "Цифровичок", тел. (495) 649-83-30 II www.cfr.ru ; e-mail: info@cfr.ru

Текст научной работы диссертации и автореферата по химии, кандидата химических наук, Зейфман, Алексей Александрович, Москва

Федеральное государственное бюджетное учреждение науки Институт органической химии им. Н.Д. Зелинского Российской академии наук

На правах рукописи

04201356919

Зейфман Алексей Александрович

Моделирование конформациоииого равновесия и сольватации при связывании органических соединений с биологическими мишенями

02.00.03 - Органическая химия

Диссертация на соискание ученой степени кандидата химических наук

Научный руководитель кандидат химических наук, Г.Г. Чилов

Москва-2013

ОГЛАВЛЕНИЕ

ВВЕДЕНИЕ........................................................................................................................................5

1. ЛИТЕРАТУРНЫЙ ОБЗОР.......................................................................................................8

1.1. Методы расчета свободной энергии связывания органических соединений с терапевтическими мишенями..................................................................................................................8

111 Теоретические основы 8

112 Практическая реализация метода 11

113 Расчет ошибок в методе ВСЭ 13

1.2. Методы анализа конформационного равновесия при связывании органических соединений с белками. ......................................................................................................................................17

12 1 Расчеты методом двухтопологического ВСЭ 18

12 2 Расчет потенциала средней силы 20

12 3 Методы «усиленного» сэмплирования 23

12 4 Анализ существующих методов моделирования конформационной равновесия 28

1.3. Методы учета растворителя в расчетах свободной энергии связывания органических соединений с белками............................................................................................................................29

13 1 Методы, основанные на неявном учете растворителя 30 132 Явный учет растворителя в расчетах свободной энергии 34 13 3 Анализ существующих методов учета сольватации при моделировании связывания белков с лигандами 41

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

2.1. Общие методы................................................................................................................................43

2 11 Пакеты программ, использованные в работе 43

2 12 Подготовка трехмерных структур и топологических файлов лигандов 43

2 13 Параметры молекулярно-динамических расчетов методом ВСЭ 44

2.2. Расчёт относительной свободной энергии связывания лигандов с Syk-киназой и PARP методом ВСЭ .........................................................................................................................................................44

2 2 1 Конструирование моделей трехмерной структуры белков 44

2 2 2 Расчет относительной свободной энергии различных лигандов Syk-киназы и PARP в растворе методом ВСЭ 45

2 2 3 Расчет относительной свободной энергии различных лигандов Syk-киназы и PARP в белке методом ВСЭ 46

2.3. Расчет относительной свободной энергии конформеров лигандов Syk-киназы методом ВО..........47

2 3 1 Подготовка топологических файлов для расчетов изменения свободной энергии конформеров 47

2 3 2 Расчет свободной энергии конформеров методом ВСЭ 49

2 3 3 Перебор конформаций лигандов Syk-киназы при расчете заселенности методом ВСЭ 51

2.4. Анализ корректности расчетов методом ВСЭ................................................................................52

2 4 1 Определение погрешностей расчетов методом ВСЭ 52

2 4 2 Определение RMSD атомов белка 53

2 4 3 Определение свободной энергии перехода при различной длине траектории 53

2 44 Апзлиз лоррелнции карт плотности воды 53

2.5. Экспериментальное определение эффективности связывания соединений с Syk-киназой и PARP.........................................................................................................................................................55

2.6. Синтез новых ингибиторов Syk-киназы.........................................................................................55

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

3.1. Конформационное равновесие лигандов Syk-киназы..................................................................61

3 11 Теоретические основы учета свободной энергии конформационной фокусировки в расчетах энергии связывания 62 3 12 Методология расчета свободной энергии конформационного перехода методом ВСЭ 64 3 13 Проверка методики расчета свободной энергии конформационной фокусировки на примере N-метилформамида 65 3 14 Особенности конформационной фокусировки ингибиторов Syk-киназы 67 3 15 Расчет конформационного равновесия ингибиторов Syk-киназы 71 3 16 Расчет свободной энергии связывания ингибиторов Syk-киназы методом ВСЭ 77 3 17 Расчет свободной энергии связывания ингибиторов Syk-киназы методом ВСЭ с учетом конформационного равновесия 81

3.2. Явный учёт сольватации комплекса белок-лиганд в расчётах методом ВСЭ..............................82

3.2.1. Расчёт свободной энергии связывания ингибиторов РАКР методом ВСЭ................................84

3.2.2. Расчёт корреляции карт плотности воды для определения степени сэмплирования молекул растворителя..............................................................................................................................................87

3.2.3. Временная корреляции карт плотности воды в белке в расчётах связывания ингибиторов РА1*Р .........................................................................................................................................................89

3.2.4. Изучение особенностей сольватации комплекса РАКР с лигандом РАКР_Ы3..........................93

3.2.5. Учёт альтернативных вариантов сольватации комплекса белок-лиганд в расчетах методом ВСЭ .........................................................................................................................................................95

3.3. Дизайн новых ингибиторов Бук-киназы с учётом конформационной фокусировки и сольватации комплекса белок-лиганд.................................................................................................100

3.3.1. Выбор модификаций структуры МТ-030...................................................................................101

3.3.2. Расчёты свободной энергии конформационной фокусировки предложенных ингибиторов Бук-киназы...............................................................................................................................................104

3.3.3. Расчёт относительной свободной энергии связывания новых ингибиторов Бук-киназы методом ВСЭ............................................................................................................................................108

3.3.4. Сопоставление результатов расчетов методом ВСЭ и предварительных структурных гипотез..................................................................................................................................113

3.3.5. Синтез и экспериментальная проверка активности новых ингибиторов Бук-киназы...........116

ЗАКЛЮЧЕНИЕ............................................................................................................................119

РЕЗУЛЬТАТЫ И ВЫВОДЫ.......................................................................................................122

СПИСОК ИСПОЛЬЗУЕМОЙ ЛИТЕРАТУРЫ............................................................................123

СПИСОК СОКРАЩЕНИЙ

ВСЭ - возмущение свободной энергии КФ - конформационная фокусировка

PARP - поли(АДФ-рибоза)-полимераза (poly-(ADP-riboso)-polymerase) Syk - селезёночная тирозинкиназа (spleen tyrosine kinase) МД - молекулярная динамика

ВВЕДЕНИЕ

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

Метод ВСЭ имеет строгое теоретическое обоснование, и при корректном проведении моделирования позволяет получить точное значение изменения свободной энергии при связывании лиганда с белком. Тем не менее, поскольку входные данные для расчётов ВСЭ формируются с помощью молекулярной динамики, итоговая точность определяется точностью молекулярно-динамических вычислений, которая в свою очередь зависит от корректности молекулярно-механического силового поля и качества сэмплирования системы. Предметом исследования настоящей работы являлись подходы к улучшению сэмплирования при использовании стандартных параметров силового поля.

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

Недостаточность сэмплирования может возникать и в ходе молекулярной динамики свободного лиганда в воде при наличии у него медленно переходящих друг в друга конформаций. Как и в случае с сэмплированием воды, пренебрежение конформацинонным равновесием лиганда может вносить в итоговую энергию связывания недопустимую ошибку порядка нескольких ккал/моль. Решению этой, а также схожей проблемы заторможенных конформаций белка посвящено большое количество теоретических работ, основанных на применении методов двухтопологического ВСЭ, потенциала средней силы и методов общего «ускорения» сэмплирования системы. Тем не менее, в большинстве случаев предложенные методы отличались низким быстродействием и сложностью практической реализации, хотя и позволяли существенно уменьшить ошибку расчётов свободной энергии связывания. Создание универсального и эффективного способа анализа конформационного равновесия для заторможенных степеней свободы лиганда и белка позволило бы не только повысить точность метода ВСЭ, но и выявить дополнительные направления модификации структуры лиганда для повышения заселённости связывающегося конформера.

Настоящее исследование направлено на решение данных задач направленного синтеза и посвящено разработке способов явного расчёта конформационного равновесия и сольватации комплекса белок-лиганд при моделировании связывания лигандов с терапевтическими мишенями методом ВСЭ. Эффективность и применимость предложенных подходов показана на примере направленного синтеза ингибиторов поли(АДФ-рибозо) полимеразы и 8ук-киназы.

1. Создание эффективного метода расчёта конформационного равновесия

органических соединений с явным учётом растворителя;

3. Направленный дизайн новых ингибиторов 8ук-киназы и поли(АДФ-рибоза)-полимеразы (РАЯР) с использованием разработанных методов;

1. ЛИТЕРАТУРНЫЙ ОБЗОР

1.1.Методы расчёта свободной энергии связывания органических соединений с терапевтическими мишенями.

Метод возмущения свободной энергии (ВСЭ) - метод статистической механики, позволяющий рассчитывать свободную энергию практически любых процессов на основании результатов молекулярно-динамических расчётов. На практике метод ВСЭ успешно применялся для дизайна ингибиторов ВИЧ-1 обратной транскриптазы [1; 2], секретируемой фосфолипазы А2 [3] и р38а MAP киназы [4], а также теоретических исследований ингибиторов катехол-О-метилтрансферазы inhibitors [5], нейраминидазы вируса гриппа [6], глюкозидазы [7] и т.д. Метод ВСЭ может быть реализован в классическом варианте, с использованием соотношения переходе Беннета либо в виде термодинамического интегрирования. Теоретические основы данных методов изложены ниже.

1.1.1. Теоретические основы

1.1.1.1. Метод возмущения свободной энергии

Метод возмущения свободной энергии (ВСЭ) был предложен Цванцигом в 1954 году для расчёта термодинамических функций газов [8]. Как известно, термодинамические свойства системы могут быть рассчитаны на основе её функции распределения Z. Свободная энергия Гельмгольца может быть рассчитана по следующей формуле:

А = U — TS = —kT\nZ = —RT In Z, (1.1)

где U - внутренняя энергия системы, S - энтропия, Т - температура в К, к - постоянная Больцмана 1,38*10"23 Дж/К. Функция распределения Z системы описывается следующим выражением:

( А\ [ [ ит Z = ехр — J = J - J е «< ал, (1.2)

где dX - фазовое пространство системы, U(X) - энергия системы в данной конфигурации X. Обозначим потенциальную энергию системы в невозмущенном состоянии U0(X), в возмущённом - иг{Х). Тогда отношение функций распределения Z-l/Zq будет иметь следующий вид:

-t/iW

Zx -(Л1-Л0) f... fe RT dX 0 S-fe-^dX

-uqW u0m

Умножим подынтегральное выражение в числителе на 1 = е rt х е rt :

-Ц^Х) -и0(Х) и0(Х) -[илп-иот] -ЫХ)

2Л \ ... \ е нт X е ят X е пт йх Г... Г е «г хе ИТ йх

_± _ I-1-_ I-1--л 4)

7 -и0(Х) -и0т ■ к }

Полученное соотношение есть не что иное, как усреднённая по состояниям невозмущённой системы разница потенциальной энергии возмущенной и невозмущенной систем:

-1 = е ИТ ={е кг )0, (1.5)

¿0

и разница в свободных энергия данных систем примет вид

-(Цх-Цр) п Лч

ДЛ = —НТ1п(е ИТ )0 (1-6)

Аналогичная формула может быть применена для расчёта свободной энергии Гиббса для М>Т-ансамбля [9]. Таким образом, итоговое значение свободной энергии перехода системы 0 в систему 1 получается термодинамическим усреднением разниц внутренних энергий системы 1 и 0 по состояниям системе 0. При экспериментальном расчёте свободной энергии методами молекулярного моделирования, как правило, проводится постепенный переход из состояния «0» в состояние «1» с несколькими шагами X, (см. ниже).

1.1.1.2. Метод термодинамического интегрирования

Метод термодинамического интегрирования был предложен в 1935 году Кирквудом [10]. Рассмотрим переход некой системы из состояния «0» в состояние «1», характеризующийся координатой перехода X. Разница свободных энергий Гиббса этих состояний по определению может быть записана в следующем виде (для ЫРТ-ансамбля):

дс = с1-с0 =/¿^«ы. (1.7)

Используя равенство (1.1), получаем:

/эс(Я)\ _ -кт /аг(Я)\ V ая )Р Т ~ г(Я) V ая )рт'

(1.8)

По определению функции распределения:

г (Я) = ¡...¡е-^йх. (1-9)

Тогда:

с-""

Подставляя данное выражение в формулу (1.8), получим:

гаясад ^ V ая )РТ г(Л)-> ■> ая

(1.11)

Данное выражение есть не что иное, как среднее значение дЕ(Х'Х) по

дЛ

конфигурационному пространству йХ. Таким образом, А С для данной системы можно записать в виде:

Г1 дЕ(Х,Х)

АС = I { дл )лйл- (1Л2)

На практике подобный расчёт проводят, постепенно модифицируя параметр X. в ходе одной молекулярно-динамической траектории и усредняя разницу в свободных энергия для каждого значения X.

1.1.1.3. Соотношение перехода Беннета

Данный подход к расчёту изменения свободной энергии системы был предложен Беннетом в 1976 году [11]. Рассмотрим две системы «О» и «1», взаимодействие в которых характеризуется различными потенциалами и0 и и1, действующими в одном и том же конфигурационном пространстве ... дп)}. Выведем разницу свободных энергий состояний «О» и «1» через средние значения функций Метрополиса М(х) = тт{1, ехр(—х)}. Функция Метрополиса обладает свойством М(х)М(—х) = ехр(—х), и поэтому используется в реализации алгоритма Монте-Карло для определения вероятности «пробных шагов», при которых изменение энергии АЦ происходит с вероятностью М(А11 Рассмотрим необычный вариант подобного «пробного шага», который полностью сохраняет конфигурацию системы (Дг—Ц-п), но переключает потенциал взаимодействия с и0 на иг или наоборот. Для каждой конфигурации будет справедл�