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

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

1.1 Методология исследования количественных соотношений «cтруктура-свойство» / «cтруктура-активность»

Методы компьютерного прогнозирования свойств химических соединений на основе поиска количественных (QSAR, Quantitative Structure Activity Relationships, количественные соотношения «структура-активность»; QSPR, Quantitative Structure Property Relationships, количественные соотношения «структура-свойство») соотношений «структура — активность» / «структура-свойство» применяются при конструировании веществ с заранее заданными свойствами и прогнозировании свойств новых соединений. Применение этих методов позволяет сократить время и затраты экспериментальных процедур по синтезу и тестированию свойств новых соединений. В последние 20 лет в связи ростом числа новых материалов, разрабатываемых для нужд микроэлектроники, экологии, медицины, фармацевтики и др., востребованность методов QSAR/QSPR в самых разных областях науки и техники значительно возросла. Эти методы активно развиваются, и спектр объектов их применения расширяется, особенно в области прогнозирования свойств материалов и супрамолекулярных систем. Так, например, на сегодняшний день 29% публикаций в области QSPR материалов посвящено исследованию полимеров, 19% - катализаторам, 16% - наноматериалам, 12%, 10% и 2%, соответственно, суперкритическому СО₂, ионным жидкостям и керамике.

Методологии QSAR/QSPR и применению методов на её основе посвящено значительное число монографий, обзоров и работ [1−8]. Методология QSAR/QSPR основана на предположении, что изменения структуры молекулы влияет на изменение макроскопических свойств вещества. Связь между молекулярной структурой и макроскопическими свойствами (например, биологической активностью или физико-химическими свойствами) рассчитывают в виде функциональной зависимости с помощью статистического аппарата методов машинного обучения. Простейшая схема методологии QSAR/QSPR представлена на Рис. 1.

Рис. 1. Схема методологии QSAR/QSPR

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

Классификация методов QSAR/QSPR. В обзоре методы QSAR/QSPR классифицируют исходя из размерности описания структуры молекул или способа представления дескрипторов:

• 1D QSAR — поиск соотношений между активностью / свойством и общими молекулярными характеристиками: pK_a, logP и др. [10−14];

• 2D QSAR — построение корреляций на основе 2D (топологических) характеристик структуры, например, индексов связности, 2D фармакофоров и др. [15−18]. Эти методы имеют ограничения при конструировании новых молекул, т. к. не учитывают характеристик их пространственной структуры;

• 3D QSAR — построение корреляций на основе 3D геометрии молекул с использованием потенциалов молекулярных полей в качестве дескрипторов [1−7,19];

• 4D QSAR — в рамках 3D QSAR используют дополнительную информацию о конформациях ансамбля лигандов [20];

• 5D QSAR — в рамках 4D QSAR используют дополнительную информацию об изменении конформации лиганда при связывании с белком (induced-fit model) [21];

• 6D QSAR — в рамках 5D QSAR используют дополнительную информацию об эффекте растворителя.

Методы статистического анализа данных. Для построения QSAR/QSPR моделей используются различные методы машинного обучения. Их целью является нахождение зависимости между входными значениями (дескрипторами) и зависимыми параметрами (свойство, активность). Многие из этих методов рассмотрены в обзоре.

Авторы дифференцируют методы QSAR/QSPR по типам методов статистического анализа, используемых для построения моделей — линейные и нелинейные.

Классическими линейными методами статистического анализа, традиционно используемым для целей QSAR/QSPR, являются множественная линейная регрессия (Multiple Linear Regression, MLR), метод частичных наименьших квадратов (Partial Least Squares, PLS), регрессия на главных компонентах (Principal Components Regression, PCR), гребневая регрессия (Ridge Regression, RR).

К нелинейным методам относят искусственные нейронные сети (ANN) [24−25], метод ближайших соседей (kNN) и ряд других.

В последнее время в области QSAR/QSPR активно применяются ядерные методы машинного обучения [27], которые имеют большие перспективы в прогнозировании свойств химических соединений и материалов. Эти методы приспособлены для работы с математическим аппаратом ядер и позволяют представить произвольную нелинейную функцию в виде линейной комбинации нелинейных ядер. К ним относят машину опорных векторов (Support Vector Machine, SVM) [28], ядерную гребневую регрессию (Kernel Ridge Regression, KRR) [29], ядерные частичные наименьшие квадраты (Kernel Partial Least Squares, K-PLS) и ряд других методов.

Дескрипторы. На сегодняшний день разработано и описано более 5000 различных дескрипторов [8,31,32] для представления тех или иных свойств молекул и, соответственно программное обеспечение для их расчёта (DRAGON [33], CODESSA и др.). В обзоре приведена следующая классификация типов дескрипторов:

· Конституционные (характеризуют относительное число атомов различного типа);

· Топологические. В их основе лежит представление молекулы в виде молекулярного графа. Топологические дескрипторы подразделяют на топоструктурные, содержащие информацию о смежности и топологических расстояниях между атомами, и топохимические, которые, кроме этого, указывают на элементную принадлежность атомов и гибридизацию. Среди топологических дескрипторов особую роль играют фрагментные дескрипторы, которые показывают наличие или отсутствие тех или иных фрагментов в структуре молекулы [36]).

· Физико-химические (характеризующие растворимость, дипольные моменты, формальный заряд, липофильность и т. д.);

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

· Структурно-геометрические (характеризующие пространственную геометрию, форму и площадь поверхности молекулы, расстояние между функциональными группами);

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

Важным элементом методов QSAR/QSPR является представление химических структур в стандартных обменных форматах, которые обеспечивают возможность их хранения в базах данных и работы с ними с помощью широкого набора существующих компьютерных программ [2,4]. Для преобразования файловых форматов разработаны специальные программы (ChemAxon, OpenBabel, Avogadro и др.).

1.2 Методы 3D QSAR/QSPR

Несмотря на то, что по сравнению с подходами 3D QSАR/QSPR классические методы 2D QSАR/QSPR более просты и лучше приспособлены для анализа больших массивов данных, они имеют ограничения при конструировании новых соединений и, особенно, супрамолекулярных систем. В частности, они 1) не позволяют учитывать особенности пространственного строения молекул, и как следствие, различать стереоизомеры, 2) не позволяют детально описывать межмолекулярные взаимодействия лиганд-мишень, 3) в рамках этих методов невозможна наглядная интерпретация результатов путем рассмотрения пространственного строения комплексов лиганд-мишень.

Поскольку практически все свойства химических соединений, обусловленные образованием межмолекулярных комплексов, зависят от их пространственного строения, в настоящее время методы 3D QSAR являются ведущими при поиске новых биологически активных соединений, в частности, при создании лекарственных препаратов, а методы 3D QSPR представляют перспективный инструментарий для компьютерного прогнозирования свойств соединений и супрамолекулярных комплексов при конструировании новых материалов. Методам 3D QSAR и их применению для прогнозирования биологической активности соединений посвящено множество монографий и обзоров [1−7,38].

1.2.1 Методы 3D QSAR

К стандартными методам 3D QSAR можно отнести подходы, в основе которых лежит предположение о том, что биологическая активность лигандов обусловлена нековалентным взаимодействием с биологическими мишенями посредством молекулярных полей. В рамках этих методов для описания таких полей вычисляют энергию взаимодействия между атомами совмещенных в пространстве (выравненных) молекул и пробными атомами, помещенными в узлы воображаемой трёхмерной решётки. Такие энергии взаимодействия рассматривают как потенциалы молекулярных полей. На основе результатов расчёта формируют матрицу, каждая строка которой отвечает молекуле лиганда, а каждая колонка — энергии взаимодействия, рассчитанной на определенном узле решётки. Количественные соотношения между значениями энергий взаимодействия и значениями биологической активности получают с помощью статистического анализа на базе методов машинных обучения.

CoMFA. Исторически первым и до сих пор одним из наиболее распространённых методов 3D QSAR является CoMFA (Comparative Molecular Field Analysis, метод сравнительного анализа молекулярных полей), разработанный в 1988 г. Крамером. В рамках этого метода в качестве дескрипторов используются потенциалы электростатического и стерического полей, рассчитанные на узлах гипотетической трёхмерной решётки, по умолчанию имеющей шаг 2A и распространённая на 4A в каждом направлении от всех молекул (см. Рис. 2).

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

Электростатические и стерические поля обычно считаются достаточными для описания нековалентных взаимодействий между лигандом и биологической мишенью. Для расчёта потенциалов электростатического поля в узлы решётки помещают пробные атомы водорода с зарядом +1 (протон), а для расчёта потенциалов стерического поля — атомы углерода в sp³-гибридизации. Электростатические потенциалы рассчитывается по закону Кулона, а стерические — с использованием потенциала Леннард-Джонса 6−12. В качестве стандартного метода статистического анализа используется метод частичных наименьших квадратов PLS (Partial Least Squares). Метод CoMFA реализован в коммерческом программном продукте SYBYL [Sybyl]. Метод CoMFA описан в ряде монографий [5−7] и обзорных статей [40,41].

GRID. Программа GRID (Graphic Retrieval and Information Display) применяется в качестве альтернативы методу CoMFA. Программа GRID также рассчитывает взаимодействие между молекулой и пробными атомами, расположенными в узлах трёхмерной решётки, но имеет ряд преимуществ перед CoMFA: во-первых, вместо потенциалов Леннард-Джонса 6−12 используются более гладкие функции типа 6−4; во-вторых, для описания большого разнообразия типов межмолекулярного взаимодействия в методе GRID используют значительно большее число различных пробных атомов и даже групп атомов. В частности, в дополнение к электростатическим и стерическим потенциалам, программа вычисляет потенциалы водородной связи и гидрофобный потенциал. В работе силовое поле GRID использовали в сочетании с программой GOLPE (General Optimal Linear PLS Estimation) для исследования ингибиторов гликогенфосфорилазы b и получили хорошие статистические результаты.

CoMSIA. Другой метод, широко используемый в 3D QSAR, Метод сравнительного анализа индексов молекулярного подобия (Comparative Molecular Similarity Indices Analysis, CoMSIA) был разработан как развитие метода CoMFA. Подробное описание этого метода и его модификаций приведено в работах [46,47]. В рамках этого подхода рассчитываются индексы молекулярного подобия, которые используются в качестве дескрипторов. Расчёт проводят путём сравнения каждой молекулы базы с пробными атомами радиуса 1A с зарядом +1 и гидрофобностью +1, помещёнными в узлах решётки. Наиболее часто с помощью индексов молекулярного подобия описывают электростатические, стерические, гидрофобные поля, а также поля водородных связей. В отличие от CoMFA, для описания потенциалов в этом методе используются функции Гауссова типа, что позволяет избежать резких изменений при переходе из одной ячейки к другой и не требует введения ограничительных значений для потенциалов сверху. Кроме того, модели, полученные методом CoMSIA, легче интерпретировать визуально.

Описанные методы 3D QSAR предполагают выполнение определённого числа основных операций [3,19]:

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

Генерация 3D-геометрии. Для построения пространственной геометрии структур используют следующие подходы

· На основе экспериментальных данных. (Для многих молекул трёхмерная структура определена методом рентгеноструктурного анализа (РСА) и хранится в базах данных, доступных в электронном виде, например, Кембриджская кристаллографическая база структурных данных (Cambridge Crystallographic Structural Database) для малых молекул до 500 атомов или Банк данных белковых молекул (Protein Data Bank) для полипептидов и полисахаридов).

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

· Автоматическая конвертацией из 2D в 3D. Часто информация о строении молекулы хранится в одноили двумерном представлении, которое необходимо перевести в трёхмерную систему координат. Такие программы как CONCORD и CORINA[49] генерируют 3D-геометрию с учётом табулированных значений длин и углов связей.

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

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

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

1.2.2 Методы 3D QSPR

В различных областях науки и технологий (энергетики, микроэлектроники, фармацевтики и др.) важной научной и практической задачей является прогнозирования свойств, связанных с образованием супрамолекулярных комплексов. В последние 10 лет методы в литературе появился ряд публикаций, посвященных применению методов 3D QSAR, основанным на использовании гипотетической решетки, к объектам, отличным от традиционных для этих методов — межмолекулярных комплексов белок-лиганд. В этих работах исследованы возможности 3D QSPR прогнозирования свойств супрамолекулярных комплексов различной природы. Например, среди этих публикаций можно отметить работы по прогнозированию абсорбции красителей на целлюлозном волокне [51,52], и работы по прогнозированию каталитических свойств металлокомплексных катализаторов, которые обусловлены свойствами комплексов металлов с органическими лигандами [53,54].

Следует заметить, что в соответствии с терминологией, общепринятой в литературе [4], стандартным методами 3D QSAR принято считать методы, в рамках которых при построении 3D моделей дескрипторы рассчитывают на основе потенциалов молекулярных полей. Поэтому к методам 3D QSPR мы относим методы, в которых используют подход аналогичный стандартными методами 3D QSAR, в отличие от методов QSPR, в которых для расчёта 3D дескрипторов исследователи применяли геометрические или энергетические характеристики. Публикации, посвященные последним, подробно описаны в обзорах [9,32].

Полученные результаты демонстрируют перспективность применения методологии 3D QSAR/QSPR для целей прогнозирования свойств супрамолекулярных комплексов различной природы.

1.3 Концепция непрерывных молекулярных полей

Как уже отмечено выше, стандартные методы 3D QSAR основаны на аппроксимации молекулярных полей потенциалами, которые рассчитываются в узлах трёхмерной пространственной решётки и используются в последующем статистическом анализе в качестве дескрипторов [5−7,39]. Такой подход обладает рядом существенных недостатков. Во-первых, статистические модели оказываются чувствительным к пространственной ориентации и шагу решётки. Во-вторых, дискретные наборы дескрипторов не могут с достаточной точностью описать молекулярные поля, которые являются непрерывными физическими объектами. Уменьшения шага решётки возможно лишь до определённого предела, ниже которого статистические параметры модели вновь ухудшаются из-за появления очень большого числа дескрипторов, что ведёт к оверфиттингу (переподгонке). Увеличение шага решётки ведёт снижению числа дескрипторов, но при этом большой объём важной информации теряется.

Альтернативный подход на основе непрерывных молекулярных полей был предложен в работах [56,57]. Суть его заключается в проведении статистического анализа молекулярных объектов, представленных не в виде набора дискетных дескрипторов, а в виде непрерывных гладких функций от пространственных координат (т.н. непрерывных молекулярных полей). Авторами было показано, что это может быть достигнуто путем конструирования специальных ядерных функций в Гильбертовом пространстве и их использования для построения регрессионных и классификационных моделей в рамках таких методов машинного обучения как регрессионный метод опорных векторов, ядерная гребневая регрессия и одноклассовая машина опорных векторов. Описание на основе непрерывных молекулярных полей точнее соответствует физической природе взаимодействий между мишенью и лигандом. Применение концепции непрерывных молекулярных полей на практике стало возможным благодаря развитию методов машинного обучения с использованием ядер (kernels) вместо векторов дескрипторов фиксированного размера. Такая статистическая модель содержит не дискретный набор конечного числа параметров, а непрерывное поле регрессионных коэффициентов.

Предложенный подход можно считать первым в мировой практике примером применения статистических методов анализа функциональных данных в хемоинформатике. В настоящее время в области хемоинформатики имеется очень ограниченное число публикаций эвристического характера, касающихся методов работы с непрерывными молекулярными полями. Так, в работах [59,60] описаны индексы молекулярного сходства Карбо, которые вычисляются как интеграл произведения функций электронной плотности для пары молекул. Их использование для целей QSAR можно рассматривать как частный случай применения непрерывных молекулярных полей без использования, однако, возможностей дуальных ядерных методов машинного обучения. В работе на основе концепции, внешне сходной со способом сравнения молекул при помощи непрерывных молекулярных полей, предложен метод проведения стерического и электронного совмещения структур молекул. Абсолютное же большинство работ, касающихся прогнозирования свойств соединений с использованием ядерных методов машинного обучения, например, мaшины опорных векторов основано на описании молекул с помощью векторов признаков ограниченного и фиксированного размера и не используют способность этих методов оперировать с функциональными данными, т. е. с фактически бесконечным числом переменных.

В группе хемоинформатики физического факультета МГУ под руководством И. И. Баскина предложенный подход реализован в рамках пакета программ CMF (Continuous Molecular Fields, метод непрерывных молекулярных полей). Математическое обоснование методологии непрерывных молекулярных полей описано в работе.

1.3.1 Методология непрерывных молекулярных полей

В основе метода непрерывных молекулярных полей лежит расчёт ядер (kernels) молекулярных полей. Ядро K (M_i, M_j) описывает сходство молекулярных полей молекул M_i и M_j и рассчитывается как линейная комбинация ядер, отвечающих каждому типу полей:

(1)

где h_f — коэффициент смешения для f-го типа поля, K_f(M_i, M_j) — ядро, описывающее сходство между полем f-го типа молекул i и j. Оно рассчитывается как сумма ядер всех пар атомов молекул i и j:

(2)

где ядро k_f(A_il, A_jm) описывает сходство между полями f-го типа l-го атома i-й молекулы и m-го атома j-й молекулы. Его значение вычисляют путём интегрирования произведения полей f-го типа для l-го атома i-й молекулы и m-го атома j-й молекулы по всему трёхмерному пространству:

(3)

где — значение поля f-го типа для l-го атома i-й молекулы в точке с радиус-вектором r; - значение для m-го атома j-й молекулы. Любое молекулярное поле может быть представлено с помощью одной функции Гаусса:

(4)

где w_fil — вес вклада l-го атома i-й молекулы в поле f-го типа; б_f — фактор аттенуации для поля f-го типа, показывает ширину кривой Гаусса; r_il — радиус-вектор l-го атома i-й молекулы. В методе CMF задаётся различная параметризация w_fil для разных типов полей, например, для электростатического поля w_fil представляет собой частичный заряд l-го атома на i-й молекуле, для стерического поля — величину:

(4.1)

где — ван-дер-ваальсова энергия; - ван-дер-ваальсов радиус. Эти параметры берутся из силового поля Tripos. Гидрофобное молекулярное поле задаётся коэффициентами w_fil, представляющими собой вклады атомов данного типа в величину гидрофобности молекулы. Поля кислотности и основности по отношению к образованию водородной связи задаются коэффициентами w_fil, представляющими собой вклады атомов данного типа в величины констант Абрахама A и B.

Аппроксимация молекулярного поля функцией Гаусса даёт возможность вычислить ядро k_f(A_il, A_jm) аналитически:

= (5)

В случае стерического поля выражение (5) несколько видоизменяется:

=

(5.1)

После вычисления ядра K (M_i, M_j) его можно использовать для построения модели при помощи одного из методов машинного обучения, приспособленных для работы с бесконечным числом переменных, таким как метод опорных векторов (Support Vector Machine, SVM), ядерная гребневая регрессия (Kernel Ridge Regression, KRR), ядерный вариант метода частичных наименьших квадратов (Kernel Partial Least Squares, KPLS) и др. В общем виде уравнение регрессии для свойства y_t может быть записано в виде:

(6)

Помимо параметров a_j и b в различных модификациях метода CMF могут использоваться дополнительные параметры. В ядерной гребневой регрессии (KRR) это коэффициент регуляризации г, а в методе опорных векторов? параметр v. Эти параметры оптимизируются для получения моделей с наилучшей прогнозирующей способностью. С этой же целью можно также проводить оптимизацию фактора аттенуации б_f и коэффициентов смешения полей h_f.

Схема графической визуализации моделей CMF показана на Рис. 3.

Рис. 3. Визуализация молекулярных полей и полей регрессионных коэффициентов Теоретические аспекты рассматриваемой методологии построения моделей SAR/QSAR/QSPR, в том числе методы построения ядер на основе молекулярных полей, визуализации молекулярных полей и соответствующих полей регрессионных коэффициентов, возможности применения методов статистического анализа функциональных данных для прогнозирования свойств химических соединений с использованием непрерывных полей при построении двухи многоклассовых классификационных моделей, а также перспективы применения этой методологии построения моделей для виртуального скрининга биологически активных соединений изложены также в работах [64−67].

1.3.2 Исследование биологической активности с использованием метода непрерывных молекулярных полей

Применению методологии непрерывных молекулярных полей для решения задач прогнозирования биологической активности соединений посвящен ряд работ. В работах [56,57] метод непрерывных молекулярных полей с использованием статистического аппарата регрессии на опорных векторах применяли для прогнозирования биологической активности ингибиторов ферментов. Показано преимущество этой методологии по сравнению со стандартными методами 3D QSAR в отношении прогнозирования ингибирующей активности производных 3-амидинофенилаланина по отношению к трём ферментам подгруппы сериновых протеаз: трипсину, тромбину и фактору Ха. Модели, построенные с помощью методологии CMF, показали более высокую прогнозирующую способность на скользящем контроле, чем модели, полученные стандартными методами CoMFA и CoMSIA.

В работах [68,69] исследована возможность применения одноклассовой классификации в сочетании с методом непрерывных молекулярных полей для проведения виртуального скрининга химических соединений. Апробация метода проводилась на базе DUD [70], для построения одноклассовых моделей по методу 1-SVM использовалась программа LIBSVM. Значения электростатических, гидрофобных и стерических ядер рассчитывали в соответствии с изложенной выше методологией. В качестве статистического ядра использовали суперпозицию электростатического, гидрофобного и стерического ядер. Оптимальные значения параметров одноклассового классификатора как для комбинации ядер, так и каждого ядра в отдельности находили путем максимизации площади под ROC-кривой. Неактивные аналоги применяли для оценки прогнозирующей способности моделей. Высокие значения площадей под ROC-кривыми построенных моделей свидетельствуют в пользу работоспособности предлагаемого метода. Сделан вывод, что использование непрерывных полей для сравнения структур позволяет применять построенные модели для поиска принципиально новых лекарственных препаратов с новыми мотивами расположения фармакофоров. Достоинством подхода являются его малая чувствительность к выбору контр-примеров, возможность объединения в рамках одной модели соединений, относящихся к разным структурным классам и поиска потенциально активных лигандов с различным расположением фармакофорных групп.

В работах [73,74] впервые в рамках одноклассового подхода с использованием метода опорных векторов получены серии регрессионных и одноклассовых SVM-моделей, демонстрирующие перспективность данного подхода для виртуального скринига потенциальных ингбиторов обратной транскриптазы ВИЧ-1. Для построения одноклассовой модели в этих работах использовали одноклассовый метод опорных векторов (1-SVM). База данных включала сведения о lg (EC₅₀) для 1045, 330, 330 и 128 ненуклеозидных ингибиторов обратной транскриптазы вируса ВИЧ-1 относительно 4-х штаммов — дикого и мутантных К103N, Y188L и IRLL98 — соответственно. Показано, что представление молекулярной структуры органических лигандов на основе непрерывных молекулярных полей позволяет получать классификационные модели более высокого качества по сравнению с подходами, базирующимися на использовании «молекулярных отпечатков», спектрофоров и фрагментных дескрипторов Кархарта. Наилучшие модели, построенные на основе непрерывных молекулярных полей, имеют статистические показатели, близкие к идеальному классификатору. Эти модели авторы рекомендовали для проведения широкомасштабного виртуального скрининга. Разработанное авторами сочетание методологии непрерывных молекулярных полей с методами одноклассовой классификации является совершенно уникальным инструментом, не имеющим близкого аналога в мировой практике. Примером дальних аналогов являются методы виртуального скрининга, осуществляемые при помощи рассмотрения близости молекулярных форм и фармакофорного описания. Недостатками подобных методов являются: использование только одного представителя активных соединений в виде шаблона, отсутствие стадии обучения и невозможность оптимизации метрики, определяющей близость молекул в химическом пространстве. Подход, сочетающий метод непрерывных молекулярных полей и метод одноклассовой классификации свободен от этих недостатков, что и объясняет значительно более высокие значения площади под ROC-кривой для этого сочетания. Ещё одним дальним аналогом является метод виртуального скрининга химических соединений на основе нейросетей Кохонена, решающих задачу одноклассовой классификации. Данная методология, однако, не позволяет использовать описание молекул при помощи молекулярных полей, не говоря уже о визуализации и интерпретации моделей.

1.4 Поиск количественных соотношений «структура-свойство» для комплексов органических лигандов с Am⁺³ и Eu⁺³

В настоящее время одной из практически значимых экологических проблем является проблема переработки ядерных отходов. На сегодняшний день в мире функционирует около 440 ядерных энергетических реакторов, в результате работы которых выделяется отработанное ядерное топливо с высокой остаточной радиоактивностью, обусловленной, в том числе, содержанием в нём минорных актинидов. Для их удаления наиболее перспективными являются экстракционные методы. Выделение минорных актинидов (например, Np, Am, Сm), однако, осложняется присутствием в топливе редкоземельных элементов, обладающих сходными химическими свойствами (Eu и др). Решение проблемы состоит в поиске высокоселективных лигандов, образующих устойчивые комплексы с разделяемыми элементами. Недостатками существующих лигандов являются их низкая селективность, неустойчивость к окислению и радиационному воздействию и т. д. Трудности при разработке высокоселективных лигандов связаны с наличием в растворе сложных внутрии межмолекулярных взаимодействий в комплексах металл — лиганд [81,82]. Поэтому особое значение для направленного синтеза лигандов имеет их конструирование с использованием теоретических методов.

В настоящее время теоретический анализ комплексов лиганд-металл осуществляется почти исключительно методами молекулярного моделирования, основанными на квантово-химических вычислениях и / или на применении метода молекулярной динамики с использованием эмпирических силовых полей (см., например, [84]). Оба подхода являются, однако, чрезвычайно трудоёмкими и поэтому не дают возможность осуществлять направленный дизайн таких супрамолекулярных систем с использованием методов виртуального скрининга либо de novo дизайна, требующих перебора очень большого числа возможностей. Этого недостатка лишены методы QSPR с применением современных методов машинного обучения [85,86]. В этом направлении в последние годы были достигнуты определённые успехи.

Следует заметить, моделирование экстрагирующей способности лигандов проводиться, как правило, раздельно для процессов комплексообразования и экстракции. Поскольку константа стабильности комплекса в воде (logK) и константа экстракции (logK_ex) могут быть связаны через коэффициенты распределения свободных лигандов и их комплексов в водной и органической фазах [87], на практике коэффициенты распределения часто бывают недоступны, поэтому связь между logK и logK_ex редко используется.

Первыми попытками прогнозирования способности лигандов к образованию комплексов с металлами было исследование эмпирических корреляционных соотношений. Были получены линейные корреляции способности к комплексообразованию со свободной энергией (Linear free energy relationships, LFER) и параметрические уравнения. LFER для серии металлов M_i при одинаковых лигандах L₁ и L₂ имело следующий вид:

logK (M_iL₁) = a * logK (M_iL₂) + b (7)

тогда как для серии лигандов L_i при одинаковых металлах M₁ и M₂:

logK (M₁L_i) = a * logK (M₂L_i) + b (8)

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

В работе предложено более сложное уравнение, представляющее logK (ML) как функцию от характеристик металла и лиганда:

log (K (M₁L) / K (M₂L)) = бE_n + вH (9)

где параметры б и E_n — мягкость, в и H — жёсткость для металла и лиганда соответственно.

Попытки прогнозирования устойчивости комплексов лигандов Ca⁺² с помощью 2D QSAR были предприняты в работе [91], в которой, наряду с топологическими и физико-химическими дескрипторами, применяли молекулярные фрагменты. Однако полученные модели не имели высокой прогнозирующей способности (R²_pred = 0.4).

Варнеком с соавторами был построен ряд 2D QSPR моделей, связывающих комплексообразующую способность комплексонов некоторых классов относительно металлов из групп лантанидов и актинидов со структурой органического лиганда. Все эти модели основаны на применении метода множественной линейной регрессии с отбором переменных в сочетании с использованием фрагментных дескрипторов. Построенные модели позволили осуществить направленный дизайн молекул комплексонов на основе процедуры виртуального скрининга сгенерированной виртуальной комбинаторной библиотеки молекул органического лиганда [94,95]. При этом, однако, были выявлены существенные недостатки этого подхода. Из-за чисто топологического характера используемых дескрипторов электронное влияние атомов и пространственные факторы оказываются учтенными лишь косвенно, что делает область применимость построенных моделей очень узкой и препятствует дизайну комплексонов с принципиально новой структурой.

В работах [96,97] методами 2D QSPR моделировали фактор разделения (Separation Factor, SF) америция и европия для набора, включающего 47 полиазагетероциклических лигандов. В работе получены модели на основе топологических дескрипторов с учётом некоторых характеристик атомов (электроотрицательность, поляризуемость).

В другой работе при построении 2D QSPR моделей использовали набор линейных и нелинейных методов в рамках трех программ — Cerius2, ISIDA (In SIlico Design and Data Analysis) и CODESSA-PRO (COmprehensive DEscriptors for Structural and Statistical Analysis). Для описания структур авторы использовали два вида дескрипторов: подструктурные молекулярные фрагменты (substructural molecular fragments, SMF), рассчитанные программой ISIDA, и молекулярные дескрипторы, рассчитанные программой CODESSA-PRO. В качестве фрагментных дескрипторов (ISIDA) использовали «последовательности», которые могли содержать атомы и связи, только атомы или только связи, и «расширенные атомы», представляющие отдельный атом с окружением. С помощью программы CODESSA-PRO были вычислены различные классы молекулярных дескрипторов: структурные, геометрические, топологические, электростатические, квантовые, химические и термодинамические. При построении статистически значимой модели проводился отбор дескрипторов. В программе ISIDA для этой цели применяли процедура t-test, а в CODESSA-PRO — «best multi-linear regression». Для нелинейного анализа моделей использовали алгоритмы Radial Basis Function Neural Networks (RBFNN) и Associated Neural Networks (ASNN).

Статистичеcкие характеристики 2D QSPR моделей для прогнозирования фактора разделения комплексов полиазагетероциклических лигандов с Eu⁺³ и Am⁺³ представлены в Табл. 1.

Рис. 4. Фактор разделения (logSF) америция и европия для t-Bu-hemi-BTP, рассчитанный: (1) по моделям (0.75 — 5.13); (2) по моделям (1.07 — 1.46). Экспериментальное значение logSF=1.0

Несмотря на хорошее качество некоторых моделей, они не показали удовлетворительную прогнозирующую способность в отношении синтезированного лиганда t-Bu-hemi-BTP (logSF = 1.0). В первом случае был получен разброс значений logSF в диапазоне 0.75 — 5.13, во втором более близкие значения 1.07 — 1.46.

Преодолеть недостатки подходов 2D QSPR можно только переходя к построению т.н. 3D QSPR (3D Quantitative Structure Property Relationships) моделей, основанных на явном рассмотрении геометрического строения молекул и полей, описывающих электронное влияние атомов. Эти методы, однако, до сих пор ни разу не применялись для прогнозирования свойств низкомолекулярных супрамолекулярных систем, таких как комплексы металлов с органическими лигандами. Более того, без существенных модификаций они и не могут быть использованы для этой цели.

Таким образом, в настоящее время при конструировании лекарственных препаратов и количественном прогнозировании биологической активности органических соединений большое значение имеют методы 3D QSAR. Эти методы основаны на анализе и сопоставлении пространственных структур молекул и поиске количественных соотношений между пространственными структурами молекул и проявляемой ими биологической активностью. Кроме того, аналогичные им методы 3D QSPR в последние годы начали использоваться для прогнозирования каталитических свойств в металлокомплексном катализе. Кроме того, имеются отдельные публикации, свидетельствующие о возможности применения методов 3D QSPR для прогнозирования адсорбционной способности красителей на целлюлозном волокне. Отметим, что во всех вышеупомянутых случаях речь идет о прогнозировании свойств, связанных с образованием и свойствами супрамолекулярных комплексов. В частности, биологическая активность молекул лекарств в большинстве случаев обусловлена образованием межмолекулярных комплексов белок-лиганд, каталитическая активность в металлокомплексном катализе связана со свойствами комплексом металлов с органическими лигандами, а абсорбция красителей на целлюлозном волокне также обусловлена образованием сложных межмолекулярных комплексов. Мы предполагаем, что это связано с тем, что при образовании супрамолекулярных комплексов большую роль играют не только электронные характеристики молекул, но и особенности их пространственного строения. Отсюда естественным образом вытекает идея о том, что методы 3D QSAR/QSPR с успехом могли бы быть применены при прогнозировании и других свойств, связанных с образованием комплексов. В частности, для нас представлял особый практический интерес распространение разрабатываемого метода непрерывных молекулярных полей, реализующего 3D QSAR/QSPR-анализ, на прогнозирование способности органических комплексонов проводить разделение ионов Am⁺³ и Eu⁺³, чему посвящена заключительная часть настоящей дипломной работы. Осуществление этого потребовало от нас разработки нового класса молекулярных полей как альтернативу либо возможное дополнение существующему набору физико-химических полей.

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

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

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

В третьих, модели, построенные на основе физико-химических полей, не обеспечивают их структурно-химической интерпретации. Поэтому представляла интерес замена набора молекулярных полей физико-химической природы более универсальным набором функций («молекулярных полей»), позволяющих описывать произвольные свойства молекул, зависящие от различных типов взаимодействий, и позволяющих интерпретировать QSAR/QSPR модели на структурном уровне.

В задачу дипломной работы входила разработка такого набора функций и программного обеспечения на его основе в рамках метода непрерывных молекулярных полей (CMF). Представляла также интерес оценка применимости разработанного подхода для прогнозирования различных свойств, обусловленных образованием супрамолекулярных комплексов — биологической активности органических лигандов в отношении фармакологических мишеней, а также фактора разделения комплексов органических лигандов с катионами Am⁺³ и Eu⁺³.

2. Обсуждение результатов

2.1 Построение моделей 3D QSAR/QSPR на основе функций принадлежности точки атомным типам (непрерывных индикаторных полей)

В рамках выполнения задач дипломной работы в качестве универсального набора функций для описания свойств молекул нами были предложены и исследованы функции, определяющие меру принадлежности точки с заданными координатами атому определённого типа, или функции принадлежности атому. Эти функции представляют собой принципиально новый тип непрерывных молекулярных полей, для которого мы предложили название — непрерывные индикаторные поля (НИП, CIF — continuous indicator fields). Их дальним прототипом являются атомные индикаторные переменные, значения которых на узле решётки равно единице, если узел находится внутри какого-либо атома в молекуле, и нулю в противоположном случае. При этом рассматривается принадлежность расположенной на узле решетки точки любому атому независимо от его типа. Введённые таким способом индикаторные переменные были использованы в качестве дескрипторов при построении моделей 3D QSAR/QSPR. Замена такими индикаторными переменными 6−12 потенциалов стерического поля, рассчитанных с использованием потенциала Леннард-Джонса, позволила авторам повысить качество 3D QSAR моделей, полученных для 256 ингибиторов дигидрофолат — редуктазы, относительно стандартных моделей CoMFA.

В основе упомянутого подхода лежит предположение о том, что радиусы атомов имеют конечные размеры. В реальности, однако, атомы чётких границ не имеют, поэтому вместо бинарных (0 или 1) индикаторных переменных мы предлагаем использовать непрерывную функцию принадлежности точки с координатами r атому с номером i, которая может принимать значения в интервале от 0 до 1:

(10)

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

(11)

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

(12)

В рамках метода непрерывных молекулярных полей поле типа f для молекулы аппроксимируется следующим образом:

(13)

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

Значения коэффициентов w_fi могут быть табулированы для определённого молекулярно-механического типа атома:

(14)

где T_i — молекулярно-механический тип атома i. Подставляя (5) в (4), получаем:

(15)

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

(16)

Тогда:

. (17)

Таким образом, любое молекулярное поле X_f(r) может быть представлено как линейная комбинация функций принадлежности атомным типам, и введённые нами поля — функции принадлежности атомным типам — определяют универсальный набор полей, который может быть использован вместо стандартных физико-химических полей для построения моделей 3D QSAR/QSPR.

Иными словами, вместо бинарных индикаторных переменных мы предлагаем использовать непрерывную функцию, показывающую меру принадлежности (в интервале от 0 до 1) произвольной точки пространства атому определённого молекулярно-механического типа. Именно это и есть «функция принадлежности атомным типам» или НИП. Это функция от пространственных координат точки. Можно предположить, что преимущество использования НИП при построении моделей 3D QSPR будет проявляться тогда, когда стандартного набора полей недостаточно для адекватного описания межмолекулярных взаимодействий, ответственных за проявление моделируемого свойства. Кроме того, с помощью НИП можно описать моделируемую зависимость со структурной точки зрения. В частности, поля регрессионных коэффициентов таких моделей должны показывать, какие изменения на структурном уровне надо ввести для модификации проявляемого молекулой свойства.