Содержание к диссертации
Введение
Глава 1. Моделирование взаимодействия белка с лигандом 10
1.1 Теоретические основы для описания взаимодействия белка с лигандом 11
1.1.1 Основные понятия термодинамики и статистической физики 11
1.1.2 Выражение свободной энергии образования комплекса через конфигурационные интегралы 17
1.1.3 Потенциал усредненных сил 21
1.1.4 Неявный учет молекул растворителя 22
1.1.5 Пренебрежение гибкостью белка 24
1.1.6 Замена переменных в конфигурационных интегралах, «жесткие» степени свободы 1.2 Метод молекулярного докинга 30
1.3 Оценочные функции на основе усредненных межатомных взаимодействий
1.3.1 Статистические оценочные функции 37
1.3.2 Оценочные функции на основе усредненных потенциалов межатомного взаимодействия с дополнительными слагаемыми
1.4 Изменение внутренней энергии лиганда при образовании белок-лигандного комплекса 51
1.5 Тестовые наборы для оценки качества программ докинга 52
Глава 2. Методы и материалы 56
2.1 Используемая программа докинга и ее модификация 56
2.1.1 Базовая оценочная функция 56
2.2 Тестовый набор белок-лигандных комплексов 57
2.2.1 Подготовка тестового набора 58
2.2.2 Особые случаи среди комплексов тестового набора з
2.2.3 Используемые подмножества тестового набора 62
2.3 Критерии оценки качества докинга 65
Глава 3. Результаты и обсуждение 67
3.1 Анализ внутренних взаимодействий лигандов тестового набора 67
3.2 Модификация оценочной функции 76
3.3 Влияние учета внутрилигандных взаимодействий на точность докинга и оценки энергии связывания 81
3.4 Сравнение результатов докинга с другими программами 92
3.5 Лиганды белков, найденные с помощью расчетов, использующих предложенную модификацию оценочной функции 98
Заключение 101
Список литературы 104
Список рисунков
- Выражение свободной энергии образования комплекса через конфигурационные интегралы
- Статистические оценочные функции
- Тестовый набор белок-лигандных комплексов
- Влияние учета внутрилигандных взаимодействий на точность докинга и оценки энергии связывания
Выражение свободной энергии образования комплекса через конфигурационные интегралы
В обзоре «The Statisticalhermodynamic Basis for Computation of Binding Affinities: A Critical Review» [67] приводится детальный теоретический вывод свободной энергии образования комплекса белок-лиганд, на основе статистической термодинамики. Главная задача этого обзора — привести аккуратную теоретическую оценку изменения энтропии при образовании комплекса. Помимо основной задачи там решается также ряд дополнительных задач. Большинство результатов данного параграфа цитируются именно из этого обзора, а также из других работ Майкла Гильсона (Michael K. Gilson) в соавторстве [68; 69].
Рассмотрим реакцию образования комплекса белок-лиганд, описываемую химическим уравнением: Р + L PL (1.23) в природе эта реакция обычно протекает в растворителе. Условие равновесия реакции будет выражено уравнением: где /І - химический потенциал соответствующей частицы в растворителе, индекс «s» показывает наличие растворителя. Зависимость химического потенциала компонента раствора от концентрации компонента задается соотношением: Ms,« = Ms І + RTln fiCi/Co) fis j + RT\n(Ci/Co) (1.25) где jІ - коэффициент активности i-го компонента, а Со - значение концентрации компонента, выбранное в качестве стандартного. Последнее приближенное равенство верно при низких концентрациях компонента, а именно такой случай как правило имеет место при взаимодействии белков с лигандами и рассматривается в данной работе. Из равенства 1.25 с учетом определения 1.2 и условия равновесия смеси 1.24 можно получить связь между свободной энергией образования комплекса AGpj и константой диссоциации комплекса Kv\: AGp i = fis і — fis — fis і = —RTln(Kpi) (1.26) это равенство устанавливает связь между константой диссоциации комплекса Крі и свободной энергией образования комплекса AGP}1, определяемой как приращение свободной энергии при образовании белок-лигандного комплекса из свободно плавающих молекул белка и лиганда.
Чтобы равенство 1.26 имело смысл, логарифмируемое выражение должно быть безразмерной величиной, однако может создаться впечатление, что согласно определению 1.2, константа диссоциации имеет размерность обратной концентрации. Это кажущееся противоречие объясняется тем, что концентрации измеряются по отношению к стандартной концентрации Со, которую обычно принято брать равной 1М, то есть один моль на литр. В результате аккуратное выражение для связи между свободной энергией и константой диссоциации принимает вид: КР1 = e-AGp,i/RT = ( pl J (1.27) CpCi Из равенства 1.27 видно, что с ростом стандартной концентрации (при сохраняющихся равновесных концентрациях всех компонент смеси) свободная энергия образования комплекса убывает, а константа диссоциации возрастает - это наблюдение согласуется с законом действующих масс [68].
В свою очередь, стандартные химические потенциалы отдельных компонент смеси, фигурирующие в формуле 1.26, можно выразить через каноничеси-ке статсуммы для этих компонент в растворе для случая низкой концентрации растворяемых веществ [67] [64, с. 8,10]. Связь между стандартным химическим потенциалом молекулы А в растворителе /is и статистическими суммами для системы из N » 1 молекул растворителя с молекулой А и без молекулы А -QN,A(VN,A) и QN,O(VN,O) соответственно - задается формулой: О rtrrn I 4JN,A\VN.A) \ , nOт/ тг fJ,sA = Ш Iті гЧ + г ( VN А — у NO) (1.28) VN,AC QN,O{VN,O) тут VN,A и ЛГ,О объем системы с молекулой А и без нее в состоянии равновесия при давлении Р (в нашем случае это атмосферное давление). С0 - стандартная концентрация.
Последнее слагаемое в формуле 1.28, соответствующее работе против силы давления при внесении одной молекулы растворяемого вещества в N молекул растворителя, для рассматриваемого случая пренебрежимо мало при атмосферном давлении [61, с. 675], [67] [64, с. 11] и далее мы его писать не будем.
Логарифмируемое выражение в формуле 1.28 можно упростить, с учетом определения статсуммы как интеграла по фазовому пространству 1.10. Во первых, в классическом приближении вклад импульса в гамильтониан ограничивается кинетической энергией pf/2m,, поэтому интеграл по фазовому пространству распадается в произведение конфигурационного интеграла по координатам и итеграла по импульсам, причем интеграл по импульсам вычисляется и дает множитель вида (27гш ЛТ)3/2 на каждый атом системы. Q = С е ,x dpdx = С / е Ух + рі/ mi dpdx = = I [ (2irrriiRT) С е x dx (1.29) здесь U{x) - потенциальная энергия системы, ті - масса i-го атома, сумма и произведение берутся по всем атомам, образующим систему. Во вторых, если разделить координаты в конфигурационном пространстве системы на внешние координаты растворяемой молекулы А, ее внутренние координаты и координаты молекул растворителя, относительно молекулы А, то интегралы по внутренним координатам молекулы А и координатам растворителя не будут зависить от внешних координат молекулы А. Поэтому интеграл по внешним координатам молекулы А вычисляется независимо и дает множитель 8TT2VN,A перед конфигурационным интегралом [67; 70; 71] (для нелинейной молекулы). В результате в логарифмируемом выражении в формуле 1.28 сокращаются множители (2ттггцКГ)э /2, относящиеся к атомам растворителя, и множитель VN,A. В итоге формула для стандартного химического потенциала молекулы А в растворителе принимает вид:
Статистические оценочные функции
Можно отметить еще одну неточность в изложение Manfred J. Sippl: когда говорят об энергии всюду имеют ввиду сумму потенциальной энергии молекулы и ее энерги взаимодействия с растворителем. Кинетическая энергия нигде не учитывается, в частности — в статсумме и оценке вероятности микросостояния из больцмановского распределения. Однако, как отмечалось выше в разделе 1.1.2, итоговые оценки для свободной энергии системы действительно корректно записываются через конфигурационные интегралы, то есть без рассмотрения кинетической энергии системы.
Для оценки именно межатомных взаимодействия УПМВ был впервые введен тем же автором в работе 1993 года [87], как и ранее, УПМВ использовался для сравнения энергии различных конформаций пептида. В этой работе автор тоже использовал референсную систему. В качестве основной системы рассматривались пары атомов, входящих в пептид, состояние которых описывалось 6 переменными (a,b,c,d,k,r): a,b - типы аминокислотных остатков, которым принадлежат два атома, c,d - типы атомов, k - как и раньше, расстояние между аминокислотными остатками в штуках, r - геометрическое расстояние между атомами, дискретизированное на несколько интервалов. Референсная же система описывалась теми же параметрами, за исключением типов аминокислотных остатков: (c,d,k,r). В результате, аналогично формуле 1.53, эффективный УПМВ определяется по формуле: AEfbcdk(r) = -kTln[fabcdk(r)/fcdk(r)], где f - наблюдаемые частоты встречаемости систем в соответствующих состояниях.
Также в этой работе вводился еще один потенциал, оценивающий энергию взаимодействия отдельных атомов пептида с растворителем. Для введения этого поверхностного потенциала использовался такой же подход на основе обратного Больцмановского закона, как и для УПМВ. Для каждого атома пептида вычислялось количество других атомов пептида в сферической слое заданного радиуса, малое количество других атомов в сферической слое означает, что данный атом имеет большую площадь контакта с растворителем. Обозначая за s количество атомов белка в сферическом слое заданного радиуса, эффективный потенциал поверхностной энергии вычислялся через относительные частоты встречаемости системы в данном состоянии: AE s) = —kTln[f2C(s)/f2(s)], где j - частоты встречаемости основной и референсной системы в соответству 41 ющем состоянии. Итоговая оценка энергии конформации пептида вычисляется как сумма двух слагаемых, поверхностной энергия пептида (энергии взаимодействия с растворителем) и энергии межатомного взаимодействия: AE(S,C) = x s AE%c(si) + Шр У AEi c (rij) (1.54)
суммы берутся по атомам, входящим в пептид, весовые коэффициенты ш задают относительный вклад двух компонет в оценку энергии.
Для оценки взаимодействия лиганда с белком статистическая ОФ была впервые введена в работе Ingo Muegge и Yvonne C. Martin в 1999 году [7]. Новая ОФ называлась PMF, способ получения УПМВ базируется на описанном выше подходе Manfred J. Sippl, только анализируется распределение по расстоянию между атомами, один из которых принадлежит лиганду, а второй - белку. Используется 34 типа атомов лигандов и 16 типов атомов белка. Итоговая формула для УПМВ имеет вид: Aij(r) = —квТІп fJvol orr{r)j- (1.55) Pbulk здесь i,j - типы атомов в белке и в лиганде, соответственно; r - расстояние; ргд(г) - плотность распределения по расстояниям, вычисляемая как количество пар атомов с типами i,j на расстоянии от r до г + 5г, встречающихся в экспериментальных данных, поделенное на объем соответствующего сферического слоя. р%1а1к - отношение полного числа пар атомов с типами i,j на расстояние менее R к объему сферы радиуса R, где R = 12 A- радиус отсечки. fvol orr{r) - множитель корректировки объема, впервые введенный в этой работе, он учитывает для разных типов атомов лиганда средний объем, занимаемый другими атомами лиганда в сферическом слое радиуса r.
В работе отдельно отмечается, что полученный таким образом УПМВ является именно потенциалом свободной энергии Гельмгольца, что он неявным образом учитывает энергию взаимодействия с растворителем и часть энтропийной компоненты свободной энергии. Также в работе отмечается, что полученная оценочная функция не учитывает внутрилигандных взаимодействий. Для вычисления оценочной функции использовалась простейшая формула 1.51, то есть значение УПМВ просто суммировалось по всем парам атомов, один из которых принадлежит лиганду, а второй - белку.
В 2006 тот же автор опубликовал вторую версию ОФ - PMF04 [10]. На этот раз УПМВ выводился на существенно большем наборе белок-лигандных комплексов - 7125 штуки. Были введены потенциалы УПМВ для ионов металлов. Отмечается, что в новой версии ОФ УПМВ для галогенов приобрели статистическую значимость.
Статистическая ОФ ITScore, была разработана Sheng-You Huang с соавторами в 2006 году [12]. ITScore является классической статистической ОФ, с потенциалом УПМВ, вычисляемым по формуле 1.53, а итоговым значением ОФ, вычисляемым как сумма межатомных взаимодействий между лигандом и белком по простейшей формуле 1.51. Существенной особенностью ОФ ITScore является применение итеративной процедуры для вывода УПМВ, при которой частоты распределения межатомных расстояний на n-том шаге используются для вывода потенциала УПМВ на следующем шаге: +1() = () + (() - ()). Распределение расстояний между атомами белка и лиганда типов i и j на n-том шаге () получается с помощью докинга 786 белок-лигандных комплексов из тренировочного набора с потенциалом (). () - распределение расстояния между атомами в экспериментальных данных по тренировочному набору. Итеративный процесс завершался, когда для 99% комплексов из тренировочного набора среднеквадратичное отклонение позиции лиганда, полученной при докинге, и позиции лиганда, известной из эксперимента, становилось менее 2A . Процесс сходился за 20 итераций при значении параметра = 0.5. Эта итеративная процедура была предложена еще в 1996 году [95] для вывода потенциала для задачи сворачивания белков. Этот подход позволяет избежать проблем выбора референсной системы или множителя корректировки объема, упоминавшихся выше.
Тестовый набор белок-лигандных комплексов
Чтобы проверить гипотезу об участии внутрилигандных взаимодействий в формировании конфигурации белок-лигандного комплекса были проанализированы частоты контактов между атомами лиганда в 1300 белок-лигандных комплексах используемого тестового набора 2.2.
Под контактом подразумевается сближение двух атомов лиганда определенных типов на расстояние, менее заданного. Исследовалась частота встречаемости контактов двух типов, отражающих наиболее существенные нековалент-ные взаимодействия: водородные связи и гидрофобные взаимодействия. Нельзя сказать, что учитываемые контакты между атомами лиганда в точности совпадают с данными физическими взаимодействиями, тем не менее, количество найденных контактов качественно отражает частоту встречаемости соответствующих физических взаимодействий между атомами лиганда. Этот вопрос, а также точное определение рассматриваемых контактов и детальный анализ частот их встречаемости рассматриваются ниже в соответствующих подразделах, здесь же приводятся основные результаты.
Основные результаты анализа приведены в таблице 3.1, результаты приводятся для групп «all», «good», «gflex» и «gfflex», определенных в разделе 2.2.3. В таблице приводятся данные об общем количестве комплексов в группе, о доле комплексов, с внутрилигандными контактами каждого из рассматриваемых типов, а также о среднем количестве гидрофобных контактов среди лигандов рассматриваемой группы. Также в таблице приводится информация о доле комплексов, имеющих более чем 5, 10 или 15 внутрилигандных гидрофобных контактов Основной вывод из анализа таблицы 3.1 состоит в том, что внутрилиганд-ные водородные связи встречаются в белок-лигандных комплексах редко, в то время как внутрилигандные гидрофобные контакты встречаются часто, особенно в гибких лигандах из групп «gflex» и «gfflex». Так, в конфигурации лиганда при связывании с белком в 65% случаев присутствуют внутрилигандные гид 68 рофобные контакты, а при рассмотрении гибких лигандов из группы «gflex» (имеющих более 4 вращательных связей), доля лигандов с внутрилигандными контактами возрастает до 88%. Это наблюдение ясно показывает, что учет внут-рилигандных контактов в процедуре докинга может существенно влиять на результаты расчетов.
Также из таблицы видно, что все характеристики для групп «all» и «good» не сильно отличаются, это еще раз подтверждает, что удаление плохо подходящих для докинга комплексов из тестового набора не сильно меняет распределение различных свойств среди представленных в наборе лигандов.
Внутрилигандные водородные связи в комплексах из тестового набора
Водородные связи (ВС) представляют собой одно из ключевых нековалент-ных взаимодействий, поэтому отдельно исследовался вопрос о частоте встречаемости внутрилигандных ВС в белок-лигандных комплексах. Надо отметить, что нет строгого общепринятого определения понятия «водородная связь», это взаимодействие является комплексным и включает в себя электростатическую, ковалентную и Ван-дер-Ваальсову компоненты [118]. В частности, геометрические критерии взаимного расположения атомов донора и акцептора ВС, обеспечивающего возникновение ВС, тоже достаточно размыты, особенно для внутримолекулярных ВС [119; 120]. Поэтому нашей задачей было качественно оценить частоту встречаемости внутрилигандных ВС и оценить, насколько учет таких ВС может повлиять на результаты докинга, а не идентифицировать точно все эти ВС.
Для поиска внутрилигандных ВС сначала использовалась программа Chimera [121]. С используемыми по умолчанию значениями параметров программа обнаружила внутренние ВС в 87 лигандах, что составляет менее 7% от 1300 - общего количества рассматриваемых молекул (таблица 3.1). Однако, при внимательном рассмотрении оказалось, что даже для этих редких внутри-лигандных ВС, их наличие часто обусловлено не конформацией, которую лиганд принимает в комплексе с белком, а его внутренней химической структурой. То есть любая или почти любая конформация таких лигандов будет иметь внутренние ВС, поэтому учет таких ВС практически не будет влиять на результат оценки энергии конфигурации белок лигандного комплекса, а значит, и на результат докинга. Такие «вынужденные» внутрилигандные ВС возникают в двух случаях: когда структура ковалентных связей лиганда сближает атомы донора и акцептора ВС, не позволяя им двигаться друг относительно друга, или когда количество потенциальных доноров и акцепторов ВС в лиганде так велико, что в любой конформации лиганда будут образовываться ВС. Эти два случая «вынужденных» ВС проиллюстрированы на рисунке 3.1.
Влияние учета внутрилигандных взаимодействий на точность докинга и оценки энергии связывания
Помимо анализа эффекта от учета внутрилигандных взаимодействий на результаты докинга, было также проведено сравнение результатов докинга с использованием модифицированной ОФ Algo_score и результатов докинга, опубликованных в работе Плевзинского [27]. В работе Плевзинского проводилось тестирование программ докинга на том же самом тестовом наборе all, который используется в данной работе. Тестировались 7 программ докинга (Surflex, LigandFit, Glide, GOLD, FlexX, eHiTS и AutoDock) в сочетании с 5 способами для генерации конформаций лиганда. При этом, среди прочего, использовались критерии, аналогичные критериям dock_succ и correl, используемым в данной работе.
При сравнении полученных результатов докинга с результатами из публикации Плевзинского [27] возникает две сложности: во первых, в работе [27] не описано, что авторы делали с «плохими» комплексами, то есть с комплексами плохо пригодными для докинга, не попадающими в группу good (смотрите раздел 2.2.3). По этому поводу можно сказать, что результаты наших расчетов при включении «плохих» комплексов заметно ухудшаются (группа all по сравнению с группой good на графике 3.9), поэтому независимо от того, как именно обрабатывались эти комплексы в работе Плевзинского, маловероятно что наши расчеты получили какое-то несправедливое преимущество за счет включения «плохих» комплексов. Скорее возможно обратное, что при сравнении только по хорошим комплексам из группы good наша программа только еще больше бы выигрывала у результатов Плевзинского, но эти результаты Плевзинского нам неизвестны.
Вторая сложность связана с выбором параметров и , при которых следует сравнивать наши результаты с результатами Плевзинского. Самым простым способом было бы использовать параметры и, при которых получаются наилучшие результаты, но это приведет к несправедливому преимуществу наших расчетов по сравнению с расчетами Плевзинского. С точки зрения машинного обучения, это означает, что набор комплексов используется одновременно и как тренировочный для подбора параметров, и как тестовый набор. Однако, как отмечалось выше в разделе 3.3, принадлежность оптимальных значений параметров [0.5; 0.75], {4,5} уже не является случайным событием, если используется достаточно большой и разнообразный тестовый набор. Поэтому мы проводили сравнение результатов Плевзинского с двумя типами результатов, полученных нами с использованием модифицированной ОФ. Во первых, с наилучшими (завышенными) результатами, то есть с результатами при оптимальных для данного критерия значениях и , а во вторых с заниженными результатами, получаемыми как наихудший результат наших расчетов при условии [0.5; 0.75], {4,5}. Этот заниженный результат в таблице 3.5 обозначается как «результат с учетом внутрилигандных взаимодействий». Таблица 3.5 Сравнение результатов докинга всех 1300 комплексов с результатами из работы Плевзинского [27] Критерий Лучший Без учета внут- С учетом внут- Лучший с учетом резуль- рилигандных рилигандных вза- внутрилигандных тат в [27] взаимодействий имодействий взаимодействий dockjucc correl 62.5% 0.47 61.7% 0.5 68.5% 0.525 72.5% 0.539 Даже при = 0 в формуле 2.2 (внутрилигандные взаимодействия не учитываются, получается базовая ОФ Algo_score) результаты расчетов докин-га по всему тестовому набору оказываются довольно хорошими по сравнению с результатами других ОФ, приведёнными в работе [27]. Этот результат согласуется с публиковавшимися ранее результатами расчетов с использованием Algo_score [15], когда программа Алгокомб с ОФ Algo_score превзошла все остальные программы в тесте Рогнана [116] по поиску правильных конфигураций для 100 белок-лигандных комплексов.
По критерию dock_succ - доля комплексов, для которых докинг прошел успешно - при = 0 уже получается результат 61.7%, что превосходит все результаты из работы [27], за исключением комбинации программ GOLD + Omega_ten, для которой доля комплексов с правильно найденной конформацией лиганда составила 62.5% (рисунок 3.12). При = 0.5 доля успешных результатов докинга достигает 72.5% (график 3.9), что превосходит все результаты, приведенные в [27]. Рисунок 3.12 — копия рисунка 3 из работы Плевзинского [27]. Точки на диаграмме соответствуют результатам расчетов с помощью различных комбинаций программы докинга и способа генерации конформаций. Ордината точки показывает долю комплексов, для которых докинг прошел успешно, то есть в нашей терминологии результат докинга по критерию dock_succ. Функция Algo_score показала хорошие результаты и по корреляционному критерию correl. В работе [27] наилучшая корреляция Спирмана достигалась с помощью программы eHiTS и составила 0.47, тогда как немодифицированная ОФ Algo_score (при = 0) уже дает значение корреляции 0.5 по модулю. С учетом внутрилигандных взаимодействий значение корреляции увеличилось еще больше до 0.535 при = 0.5 и до 0.539 при = 0.75.
При анализе литературы было замечено, что набор из 195 комплексов, используемый для тестирования 16 оценочных функций в работе Ченга [36], является подмножеством используемого нами тестового набора. Это позволило нам вычислить результаты на наборе Ченга и сравнить их с опубликованными результатами. Используемый в работе Ченга критерий «docking power» соответ ствует критерию dock_succ в данной работе, также в работе Ченга вычислялась корреляция Спирмена между расчетной и экспериментальной оценками связывания, то есть критерий correl в терминах данной работы.
По набору Ченга доля комплексов, для которых докинг прошел успешно (критерий dock_succ), составила 62% без учета внутрилигандных взаимодействий, с учетом вунтрилигандных взаимодействий эта доля увеличилась до 68%, что соответствует среднему результату тестируемых в работе Ченга ОФ (рисунок 3.13). Однако, важно отметить, что в работе Ченга [36] отсутствовала процедура поиска в конформационном пространстве лиганда, вместо этого для каждого комплекса заранее генерировались 100 конформаций лиганда в активном сайте белка, для этих конформаций вычислялась оценка связывания и выбиралась конформация с наилучшей оценкой. Для гибких лигандов с очень большим конфигурационным пространством такая процедура может приводить к завышению результатов, по сравнению с нашими расчетами, так как она исключает возможность того, что в процессе перебора конфигураций программа докинга вообще не попала в правильную область конфигурационного пространства, так как среди 100 сгенерированных заранее конформаций лиганда заведомо присутствовала правильная конформация.