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

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

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

¡)

484(эиэо

Новиков Федор Николаевич

Алгоритм расчета конформационно-зависимых свойств белков для моделирования их координации с химическими соединениями

02.00.04 физическая химия 02.00.15 кинетика и катализ

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

1 2 МАЙ 2011

Москва 2011

4846056

Работа выполнена в отделе №58 Института органической химии им. Н.Д. Зелинского РАН

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

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

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

доктор химических наук, профессор Кустов Леон ид Модестович (Учреждение Российской академии наук Институт органической химии им. Н.Д. Зелинского РАН)

доктор химических наук Ефимов Александр Васильевич (Учреждение Российской академии наук Институт белка РАН)

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

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

Защита состоится 17 мая 2011 г. в 11:00 часов на заседании Диссертационного совета Д 02.222.02 при Институте органической химии им. Н.Д. Зелинского РАН по адресу: 119992 Москва, Ленинский пр-т, д. 47

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

Автореферат разослан 15 апреля 2011 г.

Ученый секретарь диссертационного совета, кандидат химических наук

Елисеев О.Л.

Общая характеристика работы

Актуальность темы.

Кислотно-основный катализ является одним из основных механизмов катализа, реализуемых ферментами. Кислотно-основные свойства каталитических остатков ферментов определяют рН-профиль активности соответствующих биокатализаторов. Известно, что кислотно-основные свойства аминокислотных остатков, принимающих участие в ферментативном катализе (таких как аспарагиновая кислота, глутаминовая кислота, гистидин, цистеин, лизин) могут сильно зависеть от пространственной структуры белка и соответствующего микроокружения остатка. Например, в случае промышленно значимого фермента пенициллинамидогидролазы, рКа каталитической Ы-концевой аминогруппы фермента смещено в кислую область на несколько единиц рН, что позволяет ферменту осуществлять ацилирование некоторых аминосоединений при термодинамическом оптимуме рН 3-4.

Основным экспериментальным методом для характеристики кислотно-основных свойств остатков белков является ядерно-магнтиный резонанс (ЯМР), применение которого на потоке сопряжено с рядом технических сложностей: наработкой достаточно большого количества высокоочищенного белкового препарата, использованием дорогостоящих ЯМР-спектрометров, а также трудоемкой расшифровкой экспериментальных данных.

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

Сложность решения задачи расчета ионизационных свойств белков обусловлена тем, что даже небольшой белок имеет десятки кислотных и основных остатков, а общее количество (конформационных и ионизационных) состояний макромолекулярной системы, которое надо учитывать для расчета соответствующих статистических сумм, выражается астрономическими числами (~1010°-101000). Существующие на сегодняшний день теоретические подходы ограничиваются неявным учетом множества возможных состояний фермента или учитывают лишь относительно небольшое количество различных состояний, получаемых с помощью методов молекулярной динамики или Монте Карло.

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

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

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

Таким образом, настоящее исследование посвящено разработке новых алгоритмов моделирования конформационно-зависимых свойств (ионизационных свойств и геометрий аминокислотных остатков) белков, характеристике точности разработанных вычислительных методов, а также практической проверке теоретических результатов при направленном поиске новых ингибиторов фермента поли (АДФ-рибозо)-полимеразы — перспективной терапевтической мишени для лечения онкологических и сердечно-сосудистых заболеваний. Цель работы.

В работе преследовались следующие основные цели:

1. Разработка нового метода моделирования кислотно-основных свойств белков.

2. Разработка нового метода моделирования трехмерной структуры белка -пространственного положения его аминокислотных остатков.

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

4. Апробация разработанных теоретических методов для направленного поиска новых ингибиторов фермента поли (АДФ-рибозо) полимеразы.

Научная новизна и практическая значимость работы.

Предложен принципиально новый граф-теоретический подход TSAR (Thermodynamic Sampling of Amino acid Residues) к моделированию свойств макромолекулярных систем, с учетом термодинамического усреднения по всевозможным состояниям системы. Алгоритм TSAR основан на представлении макромолекулярной системы в виде графа, узлы которого, соответствующие независимым молекулярным группам системы (например, аминокислотным остаткам, закрепленным на неподвижном остове полипептидной цепи), могут иметь множество (конформационных, ионизационных и т.д.) состояний, а ребра соединяют взаимодействующие между собой узлы. Основной идеей, реализованной в алгоритме

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

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

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

Апробация работы. Основные результаты работы были представлены на Московском международном конгрессе "БИОТЕХНОЛОГИЯ: состояние и перспективы развития" (Россия, Москва, ¡6-20 марта 2009 г.); конференции Американского химического общества ACS Fall 2010 (США, Бостон, 22 - 26 августа 2010 г.); IV молодежной конференции ИОХ РАН (Москва, 2010), конференции Американского химического общества ACS 2011 (США, Анахайм, 27 - 31 марта 2011 г.);

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

Структура и объем работы. Работа изложена на ÍSfyO страницах машинописного текста, содержит 26 рисунков, 47 таблиц и состоит из введения, обзора литературы, экспериментальной части, обсуждения полученных результатов, выводов и списка цитируемой литературы. Библиография насчитывает / Vl^ наименований.

Основное содержание работы Использованные методы и подходы.

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

Для конструирования полноатомных моделей белков-мишеней использовались экспериментальные данные о трёхмерных структурах из базы данных PDB [http://www.rcsb.org/]

При оценке точности алгоритмов предсказания ионизационных свойств остатков белка использован набор из 484 экспериментально определенных значений рКа (из 96 различных белков) с надежно измеренными ионизационными свойствами. Все экспериментальные данные были заимствованы из первичных источников и прошли экспертную аннотацию в соответствии с критериями наличия: прямого метода для определения состояния ионизации остатка; достаточно полной кривой титрования остатка, позволяющей сделать вывод о его значении рКа; трехмерной структуры белка.

Для исследования эффективности разработанных алгоритмов в задачах реконструкции боковых радикалов аминокислотных остатков белка был использован тестовый набор, состоящий из 65 белковых структур. Применение данного набора для оценки эффективности методов моделирования структуры белка описано в литературе [Protein Sei. - 2008. - 17. - 9. - P. 1576-85], что позволило провести прямое сравнение полученных данных с литературными.

При оценке влияния качества модели на точность молекулярного докинга использовался набор из 1300 комплексов белок-лиганд [Plewczynski, D., et al., JComputChem, 2010], содержащий 433 различных белка из 141 организма, с экспериментально определенной трехмерной структурой. В сформированный набор входили релевантные терапевтические мишени, такие как цАМФ -зависимая протеинкизана, фактор коагуляции X, глутаматный рецептор 2, фосфолипаза А2, CDK2, HIV-I протеаза и ряд других.

Для проведения виртуального скрининга использовали библиотеку STK, любезно предоставленную компанией Vitas-M Laboratory [http://vitasmIab.com/]. и библиотеку соединений компании Enamine [http://www.enamine.net/].

Методы молекулярного моделирования. Для моделирования координации химических соединений в активных центрах белков - терапевтических мишеней использовался метод молекулярного докинга, реализованный в программном комплексе Lead Finder [J. Chem. Inf. Model. - 2008. - 48. - 12. - P. 2371-85]. Суть метода состоит в том, что, зная трехмерную структуру белка и ее функционально значимый участок, (например, отвечающий за каталитическую функцию белка) можно рассчитать структуру и прочность комплекса потенциального лекарственного препарата с белком. В ходе подобных расчетов перебираются все возможные

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

Для расчета энергии взаимодействия в комплексах лиганд-терапевтическая мишень использовалось полуэмпирическое силовое поле из программы Lead Finder [J. Chem. Inf. Model. - 2008. - 48. - 12. - P. 2371-85]. Вклад каждого из молекулярно-механических слагаемых в общую энергию масштабировался с помощью трех наборов эмпирических коэффициентов, подобранных таким образом, чтобы а) максимально точно воспроизводить экспериментально измеренную энергию взаимодействия в системе химическое соединение — терапевтическая мишень; б) ранжировать различные положения связывания лиганда, таким образом, находя наиболее вероятное, и в) правильно ранжировать активные и неактивные молекулы в экспериментах по виртуальному скринингу. Использование по сути трех различных скоринговых функций в сочетании с модифицированным генетическим алгоритмом позволяет обеспечить высокую точность расчетов.

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

Исследование биологической активности. Определение ингибиторной способности предсказанных для поли-(АДФ-рибозо)-полимеразы (ПАРП) соединений было любезно проведено к.х.н. Александрой Леонидовной Захаренко (Институт химической биологии и фундаментальной медицины СО РАН). Для измерений использовалась человеческая ПАРП, экспрессированная в системе Escherichia coli. В экспериментах концентрация фермента составляла 500 нМ, концентрация активированной ДНК 2°ое/мл, степень активации 25 %.

Результаты и обсуждение

Алгоритм термодинамического семплирования.

Наблюдаемые макроскопические свойства белков — физико-химические или биохимические —определяются ансамблем различных состояний (конформационных, ионизационных и т.д.) макромолекулы, а не каким-то одним из ее конкретных состояний. Поэтому для полноценного моделирования свойств белковых молекул необходимо развитие новых, достаточно общих теоретических подходов к расчету функций распределения макромолекулярных систем, применимых для решения разнообразных задач моделирования, в которых необходимо статистическое усреднение по конфигурациям системы. Граф-теоретический алгоритм TSAR был разработан для расчета функций распределения состояний боковых радикалов аминокислотных остатков белка и применен при определении рН-зависимости состояний ионизации аминокислот в белках, а также для реконструкции боковых радикалов аминокислотных остатков.

В рамках предложенного алгоритма белок рассматривается как система независимых подвижных боковых радикалов, закрепленных на жестком полипептидном остове.

Белку сопоставляется граф, вершинами которого, являются структурно независимые группы белка (боковые радикалы, молекулы лигандов или воды). Каждой вершине присваивается набор состояний (конформационных или ионизационных), описывающих конфигурационное пространство группы белка. Каждому ребру графа сопоставляется таблица парных энергий взаимодействия состояний вершин, соединенных ребром.

Полученный граф белка представляет собой конструкцию, известную в математике как Байесовская сеть (belief network), в которой вероятность пребывания каждого узла в определенном состоянии зависит от состояний соседних вершин. Аналогия становится полной, если перевести энергии состояний в вероятности в соответствии с известным соотношением:

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

В терминах сконструированного графа задача определения ионизационного состояния белка формулируется как расчет сумм вероятностей нахождения каждой из вершин в состоянии определенного типа. Данная формулировка полностью созвучна задачам Байесовских сетей.

Отдельный научный интерес представляет также задача поиска состояния белка, отвечающего глобальному энергетическому минимуму. Решение данной задачи актуально для подготовки полноатомных моделей белка, а также для широкого круга других прикладных задач, таких как моделирование структур мутантных белков и т.д. Данный класс задач также можно свести к методам Байесовских сетей, если вместо условной вероятности, определяемой соотношением (1), использовать следующее соотношение:

Общая идея алгоритма.

В текущей реализации алгоритма TSAR нахождение функций распределения состояний индивидуальных вершин проводится следующим образом:

— Вначале выполняется триангуляция исходного графа, из которого далее

Р(А | В) = Р(В) ■

(1)

1 если Еа + ЕЛВ = min (Ел + Елв) О

(2)

выделяются клики (такие подмножества вершин неориентированного графа, что между каждой парой вершин этих подмножеств существует ребро), и строится дерево обхода клик (дерево сочленений).

■ Затем проводится распространение потенциала по дереву сочленений (передача сигнала) в обе стороны, и пересчет функций распределения клик, через которые передается сигнал.

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

распределения индивидуальных вершин находится простой сверткой.

©

А.........в.......Д........s.........А

/..A L~1 L_L—J ( . л

123 234 345

.......-о = = = =

Рисунок 1. Иллюстрация идеи алгоритма TSAR на модельном графе.

На рисунке 1 изображен граф модельной системы. Для начала, ребрами объединяются вершины, являющиеся соседними по отношению к одной из вершин. Это имеет простой физический смысл — если у вершин более одного соседа, то при исключении данной вершины возникают связи между ее соседями. Следует, однако, помнить, что на данном этапе добавляются не все возможные ребра, а лишь до тех пор, пока в получающемся графе не останется циклов длиной больше 3. В результате данной операции получается так называемый триангулированный граф. Триангуляция модельного графа приведена на Рисунке 1.

Далее, в триангулированном графе отыскиваются клики — кластеры вершин, каждая из которых соединена с каждой. По сути, клика представляет собой часть графа, исключение любой вершины которой имеет одинаковую сложность (а именно т", где п - число вершин в клике, am- число состояний вершины). А значит, клику можно рассматривать как некую элементарную составляющую графа, неделимую в смысле поиска оптимального состояния (или функций распределения вершин). Очевидно, что триангулированный граф можно разбить на клики разной степени таким образом, что весь граф будет представлять собой набор клик, на котором может быть задана структура дерева. Для построения такого дерева клики, как вершины гипер-графа, соединяются ребрами, которые в исходном графе отвечают общим вершинам соседних клик. При этом вершины, входящие в границу клик, соединенных ребром, называются сепаратором. Используя алгоритм последовательного объединения ребрами клик, можно однозначным образом получить дерево, которое в математике называется деревом сочленений (junction tree).

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

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

Анализ рассмотренного алгоритма приводит к следующим соображениям. Во-первых, вычислительная сложность задачи растет экспоненциально с основанием, равным количеству состояний вершины (аминокислотного остатка) и показателем, равным размеру максимальной клики. Таким образом, с какого-то момента точное решение системы становится практически неосуществимым, и необходимо применять различные упрощения системы. Во-вторых, триангуляция графа не является однозначной, а оптимальная триангуляция представляет собой NP-сложную задачу, т.е. есть дополнительные возможности оптимизации процедуры триангуляции для понижения сложности задачи. В процессе адаптации алгоритма TSAR под решение конкретных задач были разработаны различные методы упрощения систем и методы оптимизации процедуры триангуляции для понижения сложности задачи.

Моделирование ионизационных свойств остатков белка.

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

Применение алгоритма TSAR.

Возможность применения алгоритма TSAR определяется сложностью полученного графа, т.е. числом слагаемых в выражении для статистической суммы системы, соответствующей графу. При решении задач, связанных с расчетом ионизационных свойств остатков белка, сложность решения графа энергий (в зависимости от размера белка и количества ионизируемых остатков) находится в диапазоне 1О20-Ю250. Поэтому, нахождение функций распределения для подобных графов напрямую невозможно и необходимы упрощения — отбрасывание ребер

взаимодействий, которые не дают значительного вклада в статистическую сумму.

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

Из-за необходимости упрощения графа количество возможных состояний для каждого остатка должно быть сбалансировано. С одной стороны, оно не должно быть слишком большим, иначе задача будет нерешаемой. С другой стороны, оно не должно быть слишком маленьким, иначе расчет функций распределения перестанет иметь смысл. Было показано, что приемлемое количество конформаций — по 5 на каждую форму аминокислотного остатка, различающуюся состоянием протонирования. Результаты, приведенные в Таблице 1, показывают, что равномерное уменьшение количества форм приводит к ухудшению расчета рКа, в то же время равномерное увеличение количества форм, на практике возможное только за счет более интенсивного разрыва связей с соседями, также приводит к ухудшенным результатам.

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

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

превышает 0.05, то остаток считается нестабильным и количество его конформаций удваивается. Примеиение дополнительного семплирования для нестабильных групп заметно улучшает качество расчета рКа, что отражено в Таблице 1.

Скоринговая функция.

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

(3)

11--- ДЛ(л) -.У.(Л-)

N.

(4)

ijPp.

J*«» *

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

энергию и энергию перекрывания с окружением. Для вычислений по настоящему алгоритму может использоваться любое молекулярно-механическое силовое поле. В настоящей работе, для повышения скорости вычислений, было использовано упрощенное силовое поле докинговой программы Lead Finder.

Валидация алгоритма.

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

Рисунок 1. Корреляция расчетных и экспериментальных

-3.0 -1.5 0.0 1.5 3.0 г

значений рКа для набора из 484 аминокислотных остатков.

Существенное улучшение качества предсказаний в сравнении с нулевым приближением показано как по значениям RMSD и коэффициентов корреляций для каждого из типов остатков, так и по количествам остатков, для которых ошибка в расчете рКа превосходит 1 или 2 единицы рН.

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

. Asp

- G !u . His .Lys

- ТегС

• Г":]

• Туе

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

Таблица 1. Точность моделирования ионизационных свойств белков при помощи алгоритма TSAR.

| ,. Упрощенная скоринговая

3 viioin шапила гОиппнПЛп^ииА 4 1

Основная ¡ Уменьшенное семплирование

функция

Нулевое

модель 1ь ir Г IIя приближение

Asp ! 0.54 0.65 0.59 0.63 0.68 0.65

Giu | 0.70 0.70 0.68 0.71 0.93 0.70

His ; 0.64 0.70 0.69 0.72 0.90 0.70

RMSD1" Tyr í 0.47 0.46 0.47 0.49 1.06 0.46

Lys ; 0.46 0.45 0.44 0.62 0.47 0.45

TerC ! 0.19 0.23 0.19 0.47 1.01 0.23

TerN i 0.45 0.51 0.45 0.48 0.59 0.51

Asp i 0.69 0.55 0.64 0.57 0.42 -

Glu i 0.72 0.73 0.73 0.72 0.39 -

His i 0.79 0.76 0.77 0,71 0.59 -

R Tyr j 0.95 0.96 0.96 0.95 0.64

Lys ! 0.61 0.63 0.63 0.35 0.53 _

TerC í 0.93 0.89 0.92 0.87 0.86 -

TerN 0.91 0.88 0.91 0.90 0.93 -

Asp i 3 9 8 12 12 12

Glu i 15 15 16 15 15 17

d¡>l % His i 9 12 11 17 19 30

Tyr j 0 0 0 0 26 37

Lys i 0 2 3 10 2 8

TerC i 0 0 0 0 36 0

TerN 0 0 0 0 0 22

Asp 1 0 .................о....... 0............ 0 ........2 .......2

Glu ; 1 2 2 1 5 6

d">2 His i 1 0 1 1 1 4 5

% Tyr ¡ 0 0 0 16 16

Lys ; 0 0 0 0 0 0

TerC i 0 0 0 0 0 0

TerN ¡ 0 0 0 0 0 0

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

форму бокового радикала до 2 ' значения, полученные при отключении дополнительного еемплироваиия нестабильных остатков

г значения, полученные при отключении электростатических взаимодействий д значения, полученные при отключении Н-связей е значения для нулевой модели ж среднеквадратичное отклонение 3 процент структур с ошибкой больше 1 " процент структур с ошибкой больше 2

Реконструкция боковых радикалов остатков белка.

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

В рамках модифицированной версии алгоритма TSAR было проведено изменение метода генерации конформаций остатков за счет подключения библиотеки ротамеров. Использование библиотеки ротамеров позволило существенно увеличить радиус семплирования остатков, без значительного изменения количества конформаций, генерируемых для каждого из остатков. Для уменьшения количества конформаций состояний групп был использован метод отсечения тупиковых ветвей (dead-end elimination, DEE). Применение теорем DEE во многих случаях позволяет проводить нахождение оптимума без потери точности. Кроме того, были разработаны специфические методы упрощения с потерями точности, включающие: удаление «не влияющих» связей и «плохих» форм, а также упрощение по референсному решению.

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

Исследование эффективности разработанных алгоритмов в решении задач реконструкции боковых радикалов белка было проведено с использованием тестового набора, состоящего из 65 белковых структур. Точность реконструкции боковых радикалов оценивалась по количеству аминокислотных остатков, для которых первый двугранный угол предсказанного положения бокового радикала (%i) отличался от истинного не более чем на 40°. Также определялось количество аминокислотных остатков, для которых отклонение углов Xi и Хг не превышало 40°.

Было установлено, что точность реконструкции боковых радикалов всех аминокислотных остатков белка с помощью алгоритма TSAR составляет 88,3% по углу %i и 80,4% по углам Xi и %2, и сравнима с точностью лучших мировых аналогов

(88,9% и 79,0%) tProteinSci.-2008.- 17.-9.-P. 1576-85]. Среднее время, затраченное на моделирование одного мутанта, составило 2.3 минуты и не превышало 10 минут.

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

Тип остатка Кол-во ! остатков j Хь% Все Xi и Хг,% Поверхностные | Xi,% Х> и %2,% 1 Погруженные Хь% Xi и X2,%

Аргинин 616 81.5% 68.2% 74.6% 56.9% ! 86.2% 75.7%

Аспарагин 646 89.2% 80.2% 81.4% 71.0% ; 97.5% 88.0%

Аспарагиновая кислота 804 88.9% 77.9% 85.6% 73.4% i 94.3% 82.9%

Цистеин 248 | 95.6% 91.7% 95.6%

Глутамин 522 82.8% 68.8% 72.8% 49.4% 1 94.4% 87.2%

Глутаминовая кислота 651 : 73.3% 60.1% 65.2% 45.1% : 84.8% 80.0%

Гистидин 293 | 94.9% 86.0% 88.4% 84.1% j 98.6% 88.4%

Изолейцин 725 ; 97.0% 85.2% 97.3% 77.0% ! 97.3% 88.5%

Лейцин 1011 ! 93.5% 86.7% 89.4% 83.2% 96.0% 89.7%

Лизин 727 ! 81.8% 65.7% 73.6% 53.8% ; 91.4% 77.4%

Метионин 218 84.9% 72.0% 71.1% 50.0% 1 90.9% 80.3%

Фенилаланин 520 96.7% 93.8% 86.8% 84.2% j 97.6% 95.2%

Серин 874 70.3% 62.5% ! 78.9%

Треонин 884 | 91.5% 89.4% 93.8%

Триптофан 232 94.8% 94.8% 88.9% 88.9% ! 96.2% 96.2%

Тирозин 478 97.1% 94.1% 87.3% 82.5% j 99.0% 96.9%

Вылин 929 | 91.3% 81.5% | 93.5%

Таблица 3. Точность моделирования геометрий боковых радикалов для 65 структур из тестового набора в зависимости от погруженности остатков.

PDB Bee Погруженные PDB Bee Погруженные

ID" Xi,% Xl И X2,% Xi,% Xi и ID" Хь% Xl"X2,% Хь°/о Xi и X2,%

1531 91.9 83.8 95.7 96.1 lkoe 87.5 86.1 92.7 91.4

la7s 89.3 88.0 94.4 88.9 lmla 90.7 82.1 96.0 87.2

la8q 93.8 85.7 96.2 94.6 lmml 86.0 76.5 94.6 88.7

lagy 90.5 87.9 95.6 93.8 lnar 87.0 77.9 93.8 89.2

lako 83.8 75.9 89.8 86.8 lnls 89.7 81.6 97.6 84.2

lamm 88.0 79.5 97.1 92.5 lnoa 93.8 95.0 92.6 91.7

larb 94.1 89.0 96.7 93.7 lnpk 83.6 71.1 92.5 81.5

lb9o 91.1 74.7 93.3 81.8 1 pic 85.4 77.6 96.3 80.0

lbd8 81.7 76.7 91.3 93.5 lqj4 91.0 79.1 92.9 86.5

lbj7 88.1 70.5 93.8 76.9 lqlO 93.7 86.6 97.2 96.0

lbyi 88.1 80.3 96.0 95.9 lqlw 88.5 79.7 95.9 88.9

1с5е 91.1 88.6 98.1 97.1 lqnj 93.4 84.4 96.2 87.1

1с9о 81.1 74.4 91.1 80.0 lqq4 93.0 88.6 94.0 93.5

lcbn 97.3 95.2 100.0 100.0 lqtn 88.6 78.5 93.0 87.5

1сс7 87.9 76.1 100.0 87.5 lqtw 86.8 77.3 95.2 85.4

Icem 89.7 85.4 96.7 92.8 lqu9 91.1 84.8 93.5 89.8

lcex 90.4 86.8 97.0 95.7 lrcf 92.3 87.3 98.2 93.2

lchd 90.3 75.7 92.9 86.4 lthv 86.2 78.5 89.3 82.1

lcku 90.8 79.3 92.3 74.4 lthx 86.6 82.6 92.0 87.5

lctj 93.4 83.0 94.4 92.9 lvfy 84.1 79.1 71.4 83.3

lcz9 80.2 78.5 85.4 85.7 lvjs 84.1 76.6 89.4 84.7

lczb 86.5 81.3 85.0 88.9 lwhi 84.2 73.2 94.7 85.0

lczp 89.7 83.2 98.6 100.0 2baa 85.4 76.9 90.3 85.5

ld4t 90.9 83.1 88.2 90.9 2cpl 89.4 79.0 93.9 82.2

ldhn 87.6 79.5 96.2 93.8 2end 83.9 72.4 82.8 76.0

leca 88.0 77.8 100.0 92.3 2hvm 90.0 85.5 97.2 89.4

ledg 87.5 77.2 93.6 88.7 2pth 90.1 82.1 98.3 91.5

lgci 92.8 92.4 95.9 94.8 2rn2 85.8 72.3 93.2 89.7

lhcl 81.1 65.5 86.0 75.7 31zt 90.5 85.1 98.0 97.1

lic6 89.7 87.0 93.1 87.3 5p21 81.3 71.0 91.7 83.3

life 85.8 71.6 93.5 84.0 5pti 89.1 74.3 100.0 100.0

ligd 82.0 74.2 90.0 85.7 7rsa 84.4 77.6 90.5 83.3

lixh 84.9 78.4 87.0 80.4

а Идентификатор структуры в Protein Data Bank |"littp://www.rcsb.org/'|.

Анализируя данные по точности алгоритма TSAR в задачах реконструкции геометрий боковых радикалов, приведенные в таблице 2, можно отметить, что процент корректных предсказаний гидрофобных остатков существенно выше, чем у гидрофильных. По-видимому, это связано с тем, что гидрофобные остатки в большинстве случаев являются погруженными, и для них реализуется меньшее количество конформаций. Во многих случаях наблюдается корреляция между размером остатка и процентом корректных предсказаний: в общем случае, чем больше остаток, тем ниже процент правильных предсказаний (аспарагиновая кислота — 88.9%, глутаминовая кислота — 73.3%; аспарагин — 89.2%, глутамин — 82.8%).

Анализ причин, вызывающих ошибки реконструкции боковых радикалов, был проведен на основе доступных кристаллографических данных о стабильности конформаций боковых радикалов. Для каждого из белков, представленных в тестовом наборе, из базы данных PDB были отобраны структуры, содержащие тот же белок. Для оценки стабильности положения бокового радикала проводилось сравнение его конформации в структуре из тестового набора с конформациями этого остатка в известных кристаллографических структурах. Если во всех известных кристаллографических структурах первый двугранный угол отличался от референсного значения не более чем на 40° , то конформация бокового радикала

считалась стабильной по первому углу.

Анализ результатов показал, что в среднем структура белка содержит 88% аминокислотных остатков имеющих стабильную конформацию бокового радикала по углу xi и 79% остатков, имеющих стабильную конформацию по углам Xi и Хг■ При этом соответствующие значения для поверхностных остатков 79% и 66%, а для погруженных остатков 94% и 89%. Таким образом, точность моделирования конформаций боковых радикалов аминокислотных остатков сопоставима с точностью кристаллографических данных.

Структурная фильтрация результатов виртуального скрининга.

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

Связывание с белком активных лигандов характеризуется определенными специфическими взаимодействиями. В роли таких взаимодействий могут выступать водородные связи, координация металла, специфические стекинг взаимодействия и т.д. Можно ожидать, что такие же взаимодействия должны образовывать потенциально активные химические соединения; если же лиганды не образует специфических взаимодействий, они вероятнее всего не активны. Таким образом, выделяя определенные наборы взаимодействий в комплексе белок-лиганд и используя их для фильтрации результатов виртуального скрининга можно отсеять неверно локированные лиганды и уменьшить процент ложноположительных предсказаний. Структурная фильтрация была применена к набору из 9 белков: адренэргическому рецептору бета 2 (ADRB2), дигидрофолат редуктазе (DHFR), дегидрооротат дегидрогеназе (DHODH), 3-гидроксиметилглутарил-СоА редуктазе (HMG-CoAR), лейкотриен А4 гидролазе (LTA4H), метаботропному глутаматному рецепторуЗ (mGluR3), оротидин -монофосфат декарбоксилазе (OMPDC), рецептору активирующему полиферацию пероксисом (PPAR-g) и фосфотирозин фосфатазе 1В (РТР1В).

Для каждого белка список требуемых взаимодействий разрабатывался на основе анализа доступных PDB структур комплексов соответствующих белков с лигандами. Структурная фильтрация проводилась автоматически с помощью специально разработанного модуля "structure_fílter" программы Lead Finder. Входными данными для модуля структурной фильтрации служат структура белка и набор остатков (или атомов) белка, которые должны участвовать в образовании водородных связей с

17

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

Таблица 4. Результаты виртуального скрининга (I), и результаты виртуального скрининга и последующей структурной фильтрации (II). Результаты, соответствующие корректно локированным лигандам, приведены в скобках.

Белок Первый лиганда Последний лиганд6 ЕР40" 1 Отобрано I лигандовг

I II I И I II

АОЯВ2 117 4(4) 109475 177(177) | 25 2356 (2403) | 1476

одая 1(1) 20018 341 (291) | 564 10922(10922) | 3907

ОНШН 263 17(17) 28513 3422 (3422) I 139 3215 (3274) I 3622

НМС-СоАК 93 17(17) 19921 1372 (1200) 1 194 1876(1876) { 3520

ЬТА4Н 4 1(1) 29984 2341(1594) 191 1783 (1783) | 2574

шОиЯЗ 93 33(33) 2420 1488 (389) ! 524 1396(1396) 1488

ОМРОС 1 1(1) 55 33 (33) I 24009 24009 (24009) 1 5664

РРАИ^ 171 9(9) 9304 6718 (6718) | 39 259 (259) | 13572

РЬА2 34 6(6) 21013 7415(7415) ! 15 350 (350) | 10415

РТР1В 16 3(3) 1895 294 (294) 1 3001 9235 (9235) | 2043

а Позиция лучшего нативного лиганда в ранжированной библиотеке.

6 Позиция худшего нативного лиганда в ранжированной библиотеке. * Фактор обогащения на 40% проценте проскринированной библиотеки равен отношению процента найденных активных лигандов (в 40% библиотеки), к проценту проскринированной библиотеки. г Количество лигандов в проскринированной библиотеке, удовлетворяющих структурным критериям.

Из приведенных данных можно видеть, что применение структурной фильтрации увеличивает обогащение виртуального скрининга как минимум в несколько раз, а в отдельных случаях в несколько сотен раз. Значительное улучшение результатов виртуального скрининга можно наблюдать на примере адренэргического рецептора р2. Применение структурной фильтрации позволяет эффективно отделять действительно активные лиганды от больших, гибких неактивных лигандов. Учет специфических водородных связей также позволяет уменьшить процент ложноположительных предсказаний на таких белках, как НМО-СоАЯ, ЦГА4Н и РРАЯ-§, активный центр которых также представляет собой полость большого размера. Иллюстрация улучшения результатов виртуального скрининга приведена на рисунке 2.

1 | ! | | Докинг

I; 1 ЯЕ11111II lililí III III II11111 1 II j Докинг +- Структурная фильтрация

0,5 Доля проскринорованиой библиотеки, % Рисунок 2. Улучшение результатов виртуального скрининга ADRB2. По горизонтальной оси отложен процент проскршшроваиной библиотеки; позиции нативных лигандов представлены вертикальными линиями. Верхняя часть рисунка соответствует результатам, полученным в ходе виртуального скрининга, нижняя часть — соответствует результатам, полученным в ходе виртуального скрининга и последующей структурной фильтрации. Лиганды, предположительно корректно локированные (удовлетворяющие структурным фильтрам) изображены непрерывными линиями; остальные (предположительно, локированные не правильно) лиганды изображены пунктиром.

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

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

Влияние качества модели белка на точность молекулярного докинга.

Для оценки влияния качества модели белка на точность молекулярного докинга использовался набор из 1300 комплексов белок-лиганд, содержащий 433 различных белка из 141 организма, с экспериментально определенной трехмерной структурой. В сформированный набор входили релевантные терапевтические мишени, такие как цАМФ -зависимая протеинкизана, фактор коагуляции X, глутаматный рецептор 2, фосфолипаза А2, CDK2, HIV-I протеаза и многие другие. Для многих важных терапевтических мишеней в тестовом наборе присутствовало несколько структур комплексов белок-лиганд.

Для оценки влияния качества подготовки модели белка на точность докинга тестовый набор был подготовлен в двух вариантах: со стандартным протонированием остатков белка при рН=7.0 (базовая модель) и с корректным учетом протонирования всех остатков, выполненным с помощью алгоритма TSAR (модель TSAR).

Необходимо отметить, что вероятность успешного докинга нативных

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

Поскольку в тестовом наборе из 1300 структур присутствовало значительное количество комплексов, содержащих большие лиганды, которые не могут быть корректно локированы независимо от качества модели; а также белки (например, ядерные рецепторы), имеющие большие неполярные активные центры, в которых отсутствуют протонируемые остатки, помимо полного набора структур использовались дополнительные поднаборы. Для формирования поднаборов были использованы различные критерии правила пяти Липинского: например, для исключения белков, имеющих сильно гидрофобные центры, был сформирован набор комплексов, содержащих лиганды с logP<5, а для устранения больших лигандов, были отобраны комплексы с лигандами, имеющими менее 10 свободно вращаемых связей. Поскольку правило пяти Липинского позволяет отсеять химические соединения, для которых маловероятно наличие биологической активности, дополнительно исследовалось влияние качества модели на точность докинга соединений, полностью удовлетворяющих всем критериям Липинского.

Докинг проводился с помощью программы Lead Finder 1.1.14. Поскольку в программе Lead Finder реализован вероятностный алгоритм докинга, для устранения случайных ошибок для каждой пары белок-лиганд выполнялось 20 независимых запусков докинга. Решение, полученное в индивидуальном запуске докинга, считалось корректным, если RMSD между позой лиганда, найденной программой Lead Finder, и экспериментально известной позицией лиганда не превышало 2.0 Á. Для каждой пары белок-лиганд рассчитывалась вероятность корректного локирования по итогам 20 независимых запусков, при этом, если вероятность корректного докинга превышала 50%, воспроизведение геометрии комплекса белок-лиганд считалось успешным. Таким образом, для каждого из подготовленных тестовых наборов (базовая модель/модель TSAR, нативная/оптимизированная конформация лигандов) было проведено 26000 индивидуальных запусков докинга. Оценка точности докинга проводилась с использованием среднего процента корректных запусков докинга (average docking success rate, avDSR), вычисленного для каждого из тестовых наборов по всем 26000 запусков, и процента успешно локированных структур (structure docking success rate, stDSR).

Таблица 5. Влияние качества модели белка на точность молекулярного докинга.

, Нативная конформация Измененная конформация

Базовая модель

Кол-во структур

Модель TSAR

Базовая модель

Модель TSAR

: avDSR stDSR avDSR stDSR i avDSR stDSR avDSR stDSR

Набор из 1300 структур 1300 52.9 54.3 57.1 59.6 41.8 42.0 45.4 45.7

пНИП <=10а | 913 61.8 63.5 66.7 69.6 ! 51.1 51.1 54.3 55.7

пНВ акцепторы <= 10б 965 : 60.5 61.9 65.6 67.8 49.0 48.7 52.6 53.3

пНВ доноры <5° | 903 ! 58.3 59.3 63.3 67.1 46.1 46.1 49.3 50.3

1о§Р <5 1161 | 54.2 55.2 58.5 61.4 44.0 43.2 47.3 47.5

Мг < 500 971 61.3 62.5 66.7 69.2 51.9 51.9 55.9 57.5

Правило пяти Липинского 644 ; 60.4 61.5 66.3 70.0 ! 52.3 52.1 56.3 58.5

' Количество свободно вращаемых связей <=10 6 Количество акцепторов водородных связей <=10 " Количество доноров водородных связей <5

Результаты оценки влияния качества моделей на результаты молекулярного докинга приведены в таблице 5. Таким образом, можно утверждать, что корректная подготовка моделей белков повышает точность докинга независимо от выбора способа оценки. Использование корректных моделей белков позволяет повысить средний процент корректных запусков докинга на 4%, а количество успешно локированных структур более чем на 5%. Корректный учет протонирования особенно важен для докинга «лекарственно подобных» молекул, т.к. для химических соединений удовлетворяющих правилу пяти Липинского, количество успешно локированных структур увеличивается на 8.5%.

Поиск новых ингибиторов фермента поли-(АДФ-рибозо)-полимеразы.

Фермент поли-(АДФ-рибозо)-полимераза (ПАРП) участвует в большом количестве клеточных процессов, связанных в основном с репарацией ДНК и программируемой клеточной смертью. Важность поли-(АДФ-рибозо)-полимеразы в качестве терапевтической мишени обусловлена необходимостью блокировки репарации ДНК клеток при лечении раковых заболеваний методами химио и радиотерапии.

Подготовка полноатомной модели белка.

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

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

Так как для абсолютного большинства рассматриваемых ингибиторов отсутствовали кристаллографические данные о комплексах белок-лиганд, было необходимо определить критерии оценки корректности докиговых решений. Было замечено, что все известные ингибиторы ПАРП имеют общую структурную особенность — пару соседних атомов, один из которых является донором водородной связи, а другой акцептором. Эти атомы участвуют в формировании коррелированных водородных связей с карбонильной и амидной группой остатка Gly 863. Таким образом, формирование двух водородных связей было выбрано в качестве первого критерия оценки правильности докинга лигандов. Этот критерий позволял надежно определять корректность позы для небольших лигандов, однако его было недостаточно для определения правильности докинга больших, наномолярных ингибиторов, содержащих гидрофобный «хвост», присоединенный к фармакофорной группе. Поэтому в качестве второго критерия было выбрано положение «хвоста» лиганда в горлышке активного центра ПАРП (см. рисунок 3).

Оказалось, что модель на основе структуры lefy адекватно описывает связывания большинства ингибиторов, а ошибки докинга в основном были связаны с неправильным размещением «хвоста» лиганда. Невозможность правильного размещения «хвоста» лиганда в горлышке активного сайта указывает на наличие стерических затруднений. Дальнейшее изучение показало, что причиной стерических затруднений является положение остатка Arg 878. С использованием алгоритма TSAR была проведена реконструкция остатка Arg 878 и ближайших остатков активного центра. Повторная валидация показала адекватность модели для описания положения

Виртуальный скрининг и структурная фильтрация.

После подготовки модели белка, обеспечивающей корректный докинг большинства ингибиторов, и параметризации энергетической функции для точного предсказания энергии, был осуществлен поиск новых ингибиторов методом виртуального скрининга. Для проведения виртуального скрининга была использована библиотека из 300.000 химических соединений компании Vitas-M Laboratory.

связывания всех известных ингибиторов

ПАРП.

Рисунок 3. Структурные критерии, использованные для оценки корректности докинга ингибиторов ПАРП: коррелированные водородные связи с консервативным остатком Gly863 и положение «хвоста» лиганда в горлышке активного сайта. Красным цветом выделена позиция бокового радикала остатка Arg 878 из структуры 1 efy.

Коррелированные водородные связи лиганда с карбонильной и амидной группой остатка й1у 863, были использованы для проведения структурной фильтрации. В соответствии с идеей структурной фильтрации, лиганды, не образующие водородные связи с остатком в1у 863, считались не активными. Применение структурной

фильтрации показало, что только 1 526 из 300 ООО лигандов удовлетворяют указанному критерию. С практической точки зрения это означает, что фактор обогащения может быть увеличен в несколько раз, за счет уменьшения размера библиотеки.

Рисунок 4. Пример нового ингибитора ПАРП1, найденного в ходе виртуального скрининга.

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

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

Структура Л 0 нм^Ч N № о \\ .Лет н,с О и 0=( N—СН, И

О 6" Вг У Вг он

Расчет. 1С50 мкМ 10 97 84 94 31

Эксперимент Остаточная 3 70 64 69 43

активность, %

Мг 206 238 353 274 253

Поиск фрагментных ингибиторов.

Для проведения виртуального скрининга использовалась библиотека из 3)869 фрагментных соединений. После проведения виртуального скрининга и последующей структурной фильтрации из 31869 соединений было оставлено 487. Далее была проведена кластеризация оставшихся соединений по описанной ранее методике1.

На основе визуального анализа связывания из каждого кластера были отобраны соединения, наиболее хорошо связывающиеся с ферментом и дальнейшие эксперименты проводились именно с ними. Среди ингибиторов ПАРП для дальнейшей экспериментальной проверки отобрали 4 соединения. Три из них экспериментально показали значительное понижение каталитической активности фермента (рисунок 5). Для двух соединений определили значение 1С50. Стоит отметить, что скаффолды двух найденных активных ингибиторов не встречаются в известных ингибиторах ПАРП.

М\Л/166 М\ПП36

Ю50 8.7 мМ Ю50 2 мМ

Рисунок 5. Результаты экспериментального определения активности ингибиторов ПАРП. Приведены значения 1С50 и процент ингибирования фермента при заданной концентрации вещества.

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

1 Стройлов, B.C. Моделирование координации биологически активных соединений с терапевтическими мишенями : дис...канд. Хим. Наук : 02.004,020015 : защищена21.12.10 /Стройлов Виктор Сергеевич. — Москва, 2010. — 190 с.

Основные результаты и выводы

1. Разработан новый граф-теоретический алгоритм (TSAR), основанный на учете топологии графа макромолекулы (белка), для нахождения термодинамических функций (ионизационных, конформационных и т.д.) распределения ее состояний. Применение алгоритма к белку, имеющему mN состояний, где N -число аминокислотных остатков белка, а ш - число (ионизационных, конформационных и т.д.) состояний остатка, позволяет без потерь точности снизить сложность задачи до ш", где п - число вершин в максимальной клике графа, и п« N.

2. Расчет ионизационных свойств (рКа) аминокислотных остатков белков с помощью разработанного алгоритма TSAR характеризуется точностью в пределах 1 единицы рН для 94% остатков и точностью в пределах 2 единиц рН для 99,5% остатков по набору из 484 экспериментально определенных значений рКа. Средний коэффициент корреляции (R) между рассчитанными и экспериментально определенными значениями рКа равен 0.80 и находится в диапазоне от 0.95 (для остатков Туг) до 0.61 (для Lys).

3. Моделирование конформаций боковых радикалов белковых остатков (первого двугранного угла Xi и второго двугранного угла %2) с помощью разработанного алгоритма TSAR характеризуется точностью в пределах 40° для 88,3% остатков по углу Xi и 80,4% остатков по углам Xi и для всех остатков белка, и 93,5% по и 88,5%, соответственно, для погруженных остатков, что сопоставимо с точностью рентгеноструктурного анализа.

4. Применение алгоритма TSAR для построения полноатомных моделей структуры белка, исходя из данных рентгеноструктурного анализа, позволяет повысить точность методов докинга в моделировании структуры комплекса белок-лиганд в среднем на 8,5% на наборе 1300 комплексов, и в расчете свободной энергии образования комплекса белок-лиганд с 2.1 до 1.9 ккал/моль на наборе 343 экспериментально охарактеризованных комплексов.

5. Применение алгоритма TSAR для построения полноатомной модели фермента ПАРП1 позволило повысить точность моделирования структуры комплекса фермент-ингибитор с 82% до 99% на наборе 140 экспериментально охарактеризованных ингибиторов. Использование построенной модели фермента в виртуальном скрининге позволило найти 15 его новых ингибиторов с определенными экспериментально константами ингибирования, в микромолярном диапазоне концентраций.

Список публикаций

1. Novikov, F.N. Molecular docking: theoretical background, practical applications and perspectives / F.N. Novikov, G.G. Chilov // Mendeleev Commun. — 2009. —19. — P. 237-242

2. Stroylov, V.S. Novel fragment-like inhibitors of EphA2 obtained by experimental screening and modeling / V.S. Stroylov, T.V. Rakitina, F.N. Novikov, O.V. Stroganov, G.G. Chilov, A.V. Lipkin // Mendeleev Commun. — 2010. — 20. — P. 263-265.

3. Novikov, F.N. Improving performance of docking-based virtual screening by structural filtration / F.N. Novikov,V.S. Stroylov, O.V. Stroganov, G.G. Chilov // J. Mol. Model. — 2010. — 16. — 7. — P. 1223-30.

4. Stroganov, O.V. Lead finder: an approach to improve accuracy of protein-ligand docking, binding energy estimation, and virtual screening / O.V. Stroganov, F.N. Novikov, V.S. Stroylov, V. Kulkov, G.G. Chilov // J. Chem. Inf. Model. — 2008. — 48. — P. 2371-85.

5. Novikov, F.N. Developing novel approaches to improve binding energy estimation and virtual screening: a PARP case study / F.N. Novikov, V.S. Stroylov, O.V. Stroganov, V. Kulkov, G.G. Chilov //J. Mol. Model.

— 2009, — 15,—11, —P. 1337-47.

6. Novikov, F.N. Target Finder — база структур терапевтических мишеней для моделирования фармакологической активности химических соединений / F.N. Novikov, O.V. Stroganov, V.S. Stroylov, G.G. Chilov // Сборник тезисов докладов VII Московского международного конгресса "БИОТЕХНОЛОГИЯ: состояние и перспективы развития".— 2009.

— Москва, Россия. — Том 1. — С. 231.

7. Stroganov, O.V. Lead Finder — software for drug discovery applications / O.V. Stroganov, F.N. Novikov, V.S. Stroylov, G.G. Chilov // Сборник тезисов докладов VII Московского международного конгресса "БИОТЕХНОЛОГИЯ: состояние и перспективы развития".— 2009.

— Москва, Россия. — Том 1. — С. 236.

8. Стройлов, B.C. Lead Finder — программный продукт для поиска новых лекарственных препаратов / B.C. Стройлов, О.В. Строганов, Ф.Н. Новиков, Г.Г. Чилов // Сборник тезисов IV Российского симпозиума "Белки и пептиды". — 2009. — Казань, Россия. — С. 393.

9. Новиков, Ф.Н. Методы моделирования биологической активности органических молекул и их применение в направленном органическом синтезе. / Ф.Н. Новиков, О.В. Строганов // Тезисы докладов IV молодежной конференции ИОХ РАН. — 2010. — Москва, Россия.

10. Stroylov, V. Lead Finder in the CSAR scoring challenge / Victor Stroylov, Ghermes Chilov, Oleg Stroganov, Fedor Novikov, Val Kulkov // Тезисы докладов международной конференции «240-th ACS National Meeting». — 2010.

— Бостон, США. — COMP 145.

11. Chilov, G. Key directions of mastering docking and scoring approaches / Ghermes Chilov, Oleg Stroganov, Fedor Novikov,Victor Stroylov, Val Kulkov // Тезисы докладов международной конференции «240-th ACS National Meeting». — 2010. — Бостон, США. — COMP 24

12. Chilov, G. TSAR — a new graph-theoretical approach to computational modeling of ionization properties of proteins / Oleg Stroganov, Fedor Novikov,Victor Stroylov, Val Kulkov, Ghermes Chilov // Тезисы докладов международной конференции «241-st ACS National Meeting ». — 2011. — Анахайм, США. — COMP 36

13. Chilov, G. Lead Finder docking and scoring evaluation with a pharmaceutically relevant data set / Oleg Stroganov, Fedor Novikov,Victor Stroylov, Val Kulkov, Ghermes Chilov // Тезисы докладов международной конференции «241-st ACS National Meeting ». — 2011. — Анахайм, США. — COMP 94

Заказ № 214-1/04/2011 Подписано и печать 11.04.2011 Тираж 100 экз. Усл. п.л. 1,25

ООО "'Цифровичок", тел. (495) 649-83-30 www.cfr.nl; е-таП: info@cfr.ru

Список сокращении.

Определения.

Введение.

Глава 1. Литературный обзор.

§ 1. Методы моделирования ионизационных свойств остатков белка.

1.1. Ачгорипшы расчета ионизационных состояний остатков белка.

1.2. Силовые поля и скоринговые функции.

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

§ 2. Методы реконструкции боковых радикалов остатков белка.

2.1. Алгоритмы реконструкции боковых радикалов аминокислотных остатков.

2.2. Силовые поля и методы оценки энергии системы.

2.3. Программные продукты для моделирования структуры белка.

§ 3. Подходы к фильтрации результатов виртуального скрининга.

§ 4. Строение и функции фермента поли-(АДФ-рибоза)-полимеразы.

4.1. Общие сведения о ферменте.

4.2. Структура активного центра.

4.3. Ингибирование поли-(АДФ-рибозо)-пол1шеразы.

Глава 2. Экспериментальная часть.

§ 1. Пакеты программ, использованные в работе.

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

§ 3. Подготовка трехмерных структур лигандов.

§ 4. Докинг и виртуальный скрининг.

§ 5. Расчет свободной энергии образования комплексов.

§ 6. Фильтрация и расчет обогащения виртуального скрининга.

§ 7. Методика поиска фрагментных ингибиторов in silico.

§ 8. Экспериментальное определение параметров связывания.

Глава 3. Результаты и обсуждение.

§ 1. Предсказания ионизационных свойств остатков белка. /. Описание алгоритма TSAR.

1.2. Расчете ионизационных свойств рибонуклеазы Н.

1.3. Скоринговая функция.

1.4. Параметризация скоринговой функции.

1.5. Предсказание ионизационных свойств аминокислотных остатков.

1.6. Использование нескольких моделей белка для расчета рКа.

§ 2. Реконструкция боковых радикалов остатков белка.

2.1. Алгоритм TSAR для реконструкции боковых радикалов и моделирования структур мутантных форм ферментов.

2.2. Параметризация скоринговой функции.

2.3. Точность и вычислительная эффективность алгоритма TSAR в задачах реконструкции боковых радикачов.

§ 3. Структурная фильтрация и подготовка фокусных библиотек.

3.1. Виртуальный скрининг и структураная фильтарация.

3.2. Анализ результатов.

§ 4. Влияние качества модели на точности докинга и виртуального скрининга.

4.1. Влияние качества модели на точность молекулярного докинга.

4.2. Докинг нативных лигандов в H1V-I протеазу.

4.3. Влияние протонирования белка на точность виртуального скрининга.

4.4. Влияние точности докинга на качество виртуального скрининга.

4.5. Влияние качества модели на точность оценки энергии связывания.

§ 5. Поиск новых ингибиторов фермента поли (АДФ-рибозо) полимеразы.

5.1. Подготовка и оценка качества модели фермента.

5.2. Поиск новых ингибиторов фермента ПАРП1.

5.3. Поиск фрагментных ингибиторов фермента ПАРП1.

Выводы.

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

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

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

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

Цель работы. В работе преследовались следующие основные цели:

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

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

3. Направленный поиск новых соединений, способных к координации и ингибированию каталитической активности фермента поли (АДФ-рибозо) полимеразы

Научная новизна и практическая значимость работы.

Предложен алгоритм термодинамического семплирования аминокислотных остатков (Thermodynamic Sampling of Amino acid Residues, TSAR), позволяющий корректно описывать энергетические профили систем, содержащих большое количество степеней свободы. Граф-теоретический подход, реализованный в алгоритме TSAR, основан на анализе микроскопических состояний, характеризующих определенный аминокислотный остаток при фиксированных состояниях остальных остатков белка, и усреднении данного свойства по всевозможным состояниям системы с использованием теории Байесовских сетей. Подобный подход впервые применен для поиска оптимальных состояний и расчета термодинамических функций распределения индивидуальных групп белка (боковых радикалов, закрепленных на жестком полипептидном остове, молекул лигандов или воды). Алгоритм TSAR успешно применен для расчета ионизационных свойств аминокислотных остатков и решения задач реконструкции боковых радикалов белка. 7

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

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

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

1. Gilson, M. K. Energetics of charge-charge interactions in proteins / M. K. Gilson and B. H. Honig // Proteins. — 1988. — Vol. 3, № 1. — P. 32-52.

2. Bashford, D. pKa's of ionizable groups in proteins: atomic detail from a continuum electrostatic model / D. Bashford and M. Karplus // Biochemistry. — 1990. — Vol. 29, № 44, —P. 10219-25.

3. Antosiewicz, J. Prediction of pH-dependent properties of proteins / J. Antosiewicz, J. A. McCammon, and M. K. Gilson // J Mol Biol. — 1994. — Vol. 238, № 3. — P. 415-36.

4. Bashford, D. Generalized born models of macromolecular solvation effects / D. Bashford and D. A. Case //Annu Rev Phys Chem. — 2000. — Vol. 51, №. — P. 129-52.

5. Mehler, E. L. A self-consistent, microenvironment modulated screened coulomb potential approximation to calculate pH-dependent electrostatic effects in proteins / E. L. Mehler and F. Guarnieri // Biophys J. — 1999. — Vol. 77, № 1. — P. 3-22.

6. You, T. J. Conformation and hydrogen ion titration of proteins: a continuum electrostatic model with conformational flexibility / T. J. You and D. Bashford // Biophys J. — 1995.Vol. 69, № 5. —P. 1721-33.

7. Kieseritzky, G. Optimizing pKa computation in proteins with pH adapted conformations / G. Kieseritzky and E. W. Knapp // Proteins. — 2008. — Vol. 71, № 3. — P. 1335-48.

8. Koehl, P. Application of a self-consistent mean field theory to predict protein side-chains conformation and estimate their conformational entropy / P. Koehl and M. Delarue // J Mol Biol. — 1994. — Vol. 239, № 2. — P. 249-75.

9. Barth, P. Accurate, conformation-dependent predictions of solvent effects on protein ionization constants / P. Barth, T. Alber, and P. B. Harbury // Proc Natl Acad Sci USA.2007. — Vol. 104, № 12. —P. 4898-903.

10. Ji, C. Developing polarized protein-specific charges for protein dynamics: MD free energy calculation of pKa shifts for Asp26/Asp20 in thioredoxin / C. Ji, Y. Mei, and J. Z. Zhang // Biophys J. — 2008. — Vol. 95, № 3. — P. 1080-8.

11. Mongan, J. Constant pH molecular dynamics in generalized Born implicit solvent / J. Mongan, D. A. Case, and J. A. McCammon // J Comput Chem. — 2004. — Vol. 25, № 16. —P. 2038-48.

12. Thurlkill, R. L. pK values of the ionizable groups of proteins / R. L. Thurlkill, G. R. Grimsley, J. M. Scholtz, et al. // Protein Sci. — 2006. — Vol. 15, № 5. — P. 1214-8.

13. Grimsley, G. R. A summary of the measured pK. values of the ionizable groups in folded proteins / G. R. Grimsley, J. M. Scholtz, and C. N. Pace // Protein Sci. — 2009. — Vol. 18, № 1. —P. 247-51.

14. Georgescu, R. E. Combining Conformational Flexibility and Continuum Electrostatics for Calculating pKas in Proteins / R. E. Georgescu, E. G. Alexov, and M. R. Gunner // Biophysical journal. — 2002. — Vol. 83, № 4. —P. 1731-1748.

15. Song, Y. MCCE2: improving protein pKa calculations with extensive side chain rotamer sampling / Y. Song, J. Mao, and M. R. Gunner // J Comput Chem. — 2009. — Vol. 30, № 14, —P. 2231-47.

16. Milletti, F. Predicting protein pKa by environment similarity / F. Milletti, L. Storchi, and G. Cruciani // Proteins: Structure, Function, and Bioinformatics. — 2009. — Vol. 76, № 2. — P. 484-495.

17. Sandberg, L. A fast and simple method to calculate protonation states in proteins / L. Sandberg and O. Edholm // Proteins. — 1999. — Vol. 36, № 4. — P. 474-83.

18. Spassov, V. Z. A fast and accurate computational approach to protein ionization / V. Z. Spassov and L. Yan // Protein Sci. — 2008. — Vol. 17, № 11. — P. 1955-70.