Под молекулярно-механическим типом атома понимается тип, присвоенный атому в рамках определенного силового поля, например Tripos, Amber, MM3 и др [19]. В наших расчётах мы использовали силовое поле, разработанное фирмой Tripos (параметризация силовых полей Tripos приведена в Приложении, п. 1). Например, в силовом поле Tripos sp²-гибридизованному углероду соответствует тип C.2, а sp³-гибридизованному углероду - тип C.3. Для этих типов атомов в рамках силовых полей имеется набор табулированных значений: радиусы Ван-дер-Ваальса, длины связей с их участием, жёсткость связей с их участием, валентные углы, жёсткость валентных углов и др.

В рамках дипломной работы было создано программное обеспечение для построения 3D QSAR/QSPR моделей на основе НИП и ядерной гребневой регрессии. Описание программного обеспечения приведено в Экспериментальной части.

2.1.1 3D QSAR моделирование биологической активности органических лигандов в отношении фармакологических мишеней с использованием непрерывных индикаторных полей

С использованием созданного программного обеспечения были исследована активность органических лигандов супрамолекулярных комплексов по отношению к различным фармакологическим мишеням. Исследовали 8 выборок лигандов, традиционно используемых в исследованиях для сравнения возможностей различных методов 3D QSAR [100-102]. Наборы содержали данные о структурах органических лигандов, принадлежащих различным группам, а также сведения об их биологической активности. Данные о структурах и активностях соединений выборок были взяты из работы [102]. Наборы включали:

1) 114 ингибиторов ангиотензинпревращающего фермента (ACE), принадлежащих к 28 различным структурным типам [103];

2) 111 ингибиторов ацетилхолинэстеразы (AСhE) [104];

3) 163 лиганда бензодиазепинового рецептора (BZR) - производных 1,4 - бензодиазепин-2-она [105];

4) 322 ингибитора циклооксигеназы-2 (COX2) - производных 9 классов соединений: пирролы, имидазолы, циклопентены, бензолы, пиразолы, спирогептены, спирогептадиены, изоксазолы и тиофены [106];

5) 397 ингибиторов дигидрофолат-редуктазы (DHFR), содержащих 2,4 - диаминопиримидиновое кольцо[107];

6) 66 ингибиторов гликогенфосфорилазы b (GPB) [108];

7) 76 ингибиторов термолизина (THER) [109];

8) 88 ингибиторов тромбина (THR) [47].

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

Рис. 5. Структуры характерных органических лигандов для каждой из 8 баз: (А) эналаприл (ACE); (B) E2020 (AChE); (C) Ro14-5974 (BZR); (D) целекоксиб (COX2); (E) метотрексат (DHFR); (F) спирогидантоин глюкопираноза (GPB); (G) ZPLA (THER); (H) нафто-производное 4-TAPAP (THR)

При построении моделей на основе НИП для каждой из 8 баз формировали соответствующую обучающую выборку, которую использовали для построения моделей, и контрольную выборку, которая служила для независимой оценки их прогнозирующей способности. Разбиение на выборки было сделано согласно данными работы [102], в которой примерно 33% соединений были отобраны в качестве контрольных выборок, остальные соединения были включены в обучающие выборки. В Табл. 2 представлено разбиение лигандов на обучающие и контрольные выборки для восьми баз данных.

Данные о пространственной геометрии и выравненных пространственных структурах соединений были взяты из работы [102]. В этой работе структуры лигандов были построены с помощью программы SYBYL на основе фрагментов, минимизацию проводили с использованием силового поля TRIPOS.

В качестве моделируемых активностей для баз ACE, AChE, BZR, COX2 и DHFR исследовали отрицательный логарифм концентрации ингибитора, соответствующей 50% активности лиганда - log(IC₅₀), для баз GPB, THERM и THR - отрицательный логарифм константы связывания лиганда с биологической мишенью - log(K_i).

При построении моделей использовали четыре модификации метода непрерывных молекулярных полей: два варианта метода гребневой регрессии (KRR1, KRR2), описание которых представлено в Экспериментальной части и два режима расчета коэффициентов смешения полей h_f - с оптимизацией и без оптимизации. В последнем случае для всех типов молекулярных полей коэффициент смешения принимали равным единице (h_f=1). Хотя оптимизация коэффициентов смешения является более общим и гибким подходом, приводящим к меньшим ошибкам на обучающей выборке, однако, её использование приводит к увеличению числа настраиваемых параметров, и, следовательно, к большему оверфиттингу (overfitting, переподгонка) и большим ошибкам на контрольной выборке. Поэтому выбор между этими возможностями не всегда очевиден, и мы в дипломной работе рассматривали оба варианта. Кроме того, выбор между KRR1 и KRR2 также не является очевидным. Поэтому в рамках данной дипломной работы мы применяли оба этих метода. Из четырёх полученных моделей с различными параметрами выбирали модель с наилучшей прогнозирующей способностью.

Характеристики 3D QSAR моделей, полученных для восьми наборов лигандов с помощью методов KRR1 и KRR2 и непрерывных индикаторных полей. Формулы расчёта статистических параметров представлены в Экспериментальной части.

Как видно из Табл. 4, для пяти из восьми исследованных баз (AChE, BZR, COX2, GPB и THR) модели, полученные методом KRR1 (с оптимизацией коэффициентов смешения h_f), имеют сравнимые или более высокие значения параметра q², описывающем прогнозирующую способность в условиях скользящего контроля, чем модели, полученные методом KRR2 (с оптимизацией и без оптимизации коэффициентов смешения h_f) и KRR1 (без оптимизации h_f), соответственно - 0.67, 0.48, 0.56, 0.85, 0.77. Эти модели имеют более низкие или равные значения среднеквадратичной ошибки прогнозирования при скользящем контроле RMSE_cv, соответственно - 0.70, 0.48, 0.67, 0.42, 0.46. Для баз АCE, DHFR и THERM значения q² моделей, полученные методами KRR1 (с оптимизацией и без оптимизации h_f), немного ниже, чем для моделей, полученных методом KRR2 (с оптимизацией h_f), и, соответственно, равны 0.78 и 0.79, 0.71. и 0.72, а также 0.77 и 0.79. Соответственно, RMSE_cv для этих моделей - 1.10 и 1.26, 0.68 и 0.70, 0.91 и 1.14.

Наиболее высокой прогнозирующей способностью на независимой контрольной выборке характеризуются модели, полученные методом KRR2 (без оптимизации h_f) для 7 баз (ACE, BZR, COX2, GPB, DHFR, THERM и THR). Для этих баз значения R²_pred равны, соответственно - 0.56, 0.12, 0.13, 0.63, 0.58, 0.35, 0.35. Значения тех же параметров моделей, полученных методами KRR2 (с оптимизацией h_f) и KRR1 (как с оптимизацией h_f, так и без неё), значительно хуже. Только в случае базы AChE значения R²_pred модели (0.52), полученной методом KRR1 (с оптимизацией h_f и без нее), выше по сравнению с моделью, построенной методом KRR2 (с оптимизацией h_f). В целом же можно сделать вывод о том, что метод KRR2 без оптимизации коэффициентов смешения приводит к построению моделей с наилучшей прогнозирующей способностью на независимой выборке.

Полученные результаты можно интерпретировать следующим образом. При оптимизации коэффициентов смешения (h_f) индикаторных полей, общее число настраиваемых параметров модели значительно превышает число параметров модели, построенной без их оптимизации, когда значения коэффициентов смешения полей приняты постоянными. Эта разница равна числу типов атомов в параметризации Tripos, присутствующих в молекулах обучающей выборки. Поскольку эта разница в числе параметров весьма существенна, дополнительная ошибка прогнозирования, возникающая при оптимизации h_f вследствие оферфиттинга, в большинстве случаев может превышать ошибку прогнозирования, возникающую при отсутствии минимизации вследствие использования менее гибкой функции. Поэтому в большинстве случаев оптимизация параметров смешения, хотя и приводит к улучшению статистических характеристик, характеризующих прогнозирующую способность в условиях скользящего контроля, тем не менее, ведёт к ухудшению показателей, описывающих прогнозирующую способность на независимой выборке. С другой стороны, наличие лишь одного дополнительного настраиваемого параметра в KRR2 (свободного члена b) по сравнению с KRR1 ведёт лишь к незначительному увеличению возможности проявления оверфиттинга, что в большинстве случаев не может перевесить неблагоприятного эффекта, возникающего в KRR1 вследствие отсутствия настраиваемого свободного члена. Поэтому в среднем прогнозирующая способность моделей на контрольной выборке для моделей, построенных с помощью KRR2 выше, чем при использовании KRR1.

Решение проблемы оверфиттинга является первоочередной задачей будущего продолжения исследований.

Показатели моделей, полученных методами CoMFA, CoMSIA1 и CoMSIA2, были взяты из работы [102]. В методе CoMFA было использовано 2 типа молекулярных полей - электростатическое и стерическое. В CoMSIA1 использовали электростатическое и стерическое молекулярные поля, в CoMSIA2 к ним добавлены гидрофобное и поля водородных связей.

CMF лучшие модели, построенные на основе НИП для 3 баз данных (ACE AChE и DHFR) имеют значения статистических характеристик, описывающих прогнозирующую способность, как в условиях скользящего контроля, так и на независимой выборке, превышающие соответствующие показатели лучших моделей, полученных для этих же баз с помощью стандартных методов 3D-QSAR: CoMFA, CoMSIA1 и CoMSIA2 [102].

В случае наборов BZR и GPB CMF-модели имеют более высокие показатели на скользящем контроле по сравнению с вышеупомянутыми стандартными методами 3D QSAR, a прогнозирующая способность этих моделей выше, чем у соответствующих моделей CoMFA и CoMSIA1, и сравнима с аналогичными показателями для моделей CoMSIA2. Так, для BZR значения R²_pred моделей CMF и CoMSIA2 равны 0.12, а для GPB аналогичные показатели составляют, соответственно, 0.58 и 0.59.

Для базы COX2 при оценке предсказательной способности моделей метод CMF показал лучшие результаты (R²_pred = 0.13, RMSE_pred = 1.23) по сравнению с методом CoMSIA1 (R²_pred = 0.03, RMSE_pred = 1.44), но худшие по сравнению с методами CoMFA (R²_pred = 0.29, RMSE_pred = 1.24) и CoMSIA2 (R²_pred = 0.37, RMSE_pred = 1.17).

Для выборки THR модели, полученные с использованием метода CMF, имеют более низкие статистические показатели, описывающие прогнозирующую способность на независимой контрольной выборке, по сравнению с моделями CoMSIA2.

Как видно из Табл. 5 и 7, модели, построенные для исследованных баз с помощью метода CMF, практически во всех случаях имеют лучшую прогнозирующую способность на скользящем контроле (более высокие значения q²), чем соответствующие модели, полученные методами CoMFA, CoMSIA1, и CoMSIA2. Исключение составляют модели для COX2 и THR, полученные методом CoMSIA2. Особенно велика разность между значениями q² между моделями CMF и CoMFA для баз GPB (0.32) и DHFR (0.20). Существенная разница в значениях q² между моделями CMF и CoMSIA1 достигается для базы GPB - 0.31. Разность в значениях q² между моделями CMF и CoMSIA2 наиболее высока для баз GPB и THER (0.13) и DHFR (0.12).

Модели, построенные для исследованных баз с помощью метода CMF, имеют лучшую прогнозирующую способность на независимой контрольной выборке (более высокие значения R²_pred) для баз ACE, AChE, BZR и DHFR, чем соответствующие модели, полученные методами CoMFA, CoMSIA1, и CoMSIA2. Для базы GPB значения R²_pred CMF-моделей лучше, чем для моделей, полученных методами CoMFA и CoMSIA1, но хуже, чем для моделей, полученных методом CoMSIA2.

Таким образом, результаты анализа 3D QSAR моделей, полученных при помощи метода CMF с использованием НИП, свидетельствуют о том, что их прогнозирующая способность вполне сопоставима, а в ряде случаев и превышает прогнозирующую способность лучших моделей, построенных с помощью стандартных методов 3D QSAR и физико-химических полей. Это подтверждает выдвинутое выше утверждение о том, что НИП - это универсальный набор полей, который может быть использован вместо стандартных физико-химических полей для построения моделей 3D QSAR/QSPR. Это также свидетельствует о перспективности предлагаемого подхода и, как мы считаем, необходимости его дальнейшего развития.

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

Представлен пример визуализации перекрывания полей регрессионных коэффициентов CMF модели и непрерывных индикаторных полей, рассчитанных для двух типов атомов (ароматического углерода, C.ar и азота амидной группы, N.am) для наиболее и наименее активных соединений базы ACE, соответственно, соединения №54 (THIOL_9) и соединения №76 (SQ29852_2Q). Изо-поверхности полей регрессионных коэффициентов показаны как сплошные окрашенные поверхности. Красным цветом отмечены области предпочтительного для проявления активности расположения атомов данного типа, синим показаны области, в которых присутствие атомов данного типа нежелательно. Молекулярные поля для определённого типа атомов индивидуальных молекул представлены сеткой.

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

2.1.2 3D QSPR моделирование фактора разделения комплексов полиазагетероциклических лигандов с Am⁺³/Eu⁺³ с использованием непрерывных индикаторных полей

Для проверки возможности применения разработанного нами универсального набора непрерывных индикаторных полей для моделирования свойств комплексов органических лигандов с ионами металлов нами было предпринято построение 3D QSPR модели на основе методологии непрерывных молекулярных полей для прогнозирования фактора разделения Am⁺³/Eu⁺³ при образовании комплексов этих катионов с полиазагетероциклическими лигандам. Построение такой модели представляло для нас большой интерес как в теоретическом, так и в практическом плане. Теоретический интерес обусловлен тем, что из-за того, что строение внешней электронной оболочки атома Am (5f⁷7s²) аналогично строению внешней электронной оболочки атома Eu (4f⁷6s²), а значения ионных радиусов ионов Am⁺³ и Eu⁺³ различаются всего лишь на 0.012 нм, существование селективных лигандов к ним с трудом может быть объяснено «классическими» представлениями о геометрическом строении комплексов и «чистым» электростатическим взаимодействием иона металла с молекулами лиганда, лежащим в основе стандартного набора физико-химических полей. Поэтому можно было ожидать, что непрерывные индикаторные поля, благодаря своему универсальному характеру, окажутся способными в неявном виде аппроксимировать другие типы взаимодействия, недоступные для описания стандартным набором физико-химических полей, и, следовательно, привести к построению моделей 3D QSPR с более высокой прогнозирующей способностью. Также можно было надеяться, что в случае успеха дальнейший подробный анализ построенных моделей в сопряжении с высокоуровневыми квантово-химическими расчётами позволил бы определить ключевые факторы, необходимые для дизайна новых комплексонов с рекордно высокой разделяющей способностью по отношению к Am⁺³/Eu⁺³. Практический же интерес в построении такой модели состоит в том, что такая модель в дальнейшем могла бы быть использована для решения важнейшей задачи, связанной с утилизацией отходов атомной промышленности.

Для построения модели нами была использована база данных, на основе которой авторами работы [97] ранее были получены количественные модели «структура-свойство» с помощью фрагментных 2D дескрипторов. Исходная база включала сведения о 48 производных полиазагетероциклических соединений, принадлежащих к 6-ти структурным группам: L1: 5,6 - диалкил-3 - (4,6 - диалкилпиридин-2-ил) - 1,2,4 - триазины; L2: 5,6 - диалкил-3 - (пиразин-2-ил) - 1,2,4 - триазины; L3: 5,6 - диалкил-3 - (5,6 - диалкил - 1,2,4 - триазин-3-ил) - 1,2,4 - триазины; L4: 3 - (4,6 - диалкилпиридин-2-ил) - 5,6 - ди (пиридин-2-ил) - 1,2,4 - триазины; L5: 2,6 - бис (5,6 - диалкил - 1,2,4 - триазин-3-ил) пиридины; L6: 5,6 - диалкил-3 - (4-алкил-6 - (4-алкилпиридин-2-ил) пиридин-2-ил) - 1,2,4 - триазины.

В качестве моделируемого свойства нами был взят фактор разделения (SF, separation factor) ионов Am⁺³/Eu⁺³, представляющий собой отношение коэффициентов распределения между органической и водной фазами для этих металлов (D_am/D_eu). Значения SF были получены путём экстракции катионов металлов (Am⁺³, Eu⁺³) из водного раствора азотной кислоты в 1,1,2,2 - тетрахлороэтановую фазу, содержащую альфа-бромкаприновую кислоту в качестве соэстрагента [110]. Концентрации европия (¹⁵²Eu) и америция (²⁴¹Am) в жидких фазах были измерены методом гамма-спектрометрии. Исходная концентрация азотной кислоты для разных соединений варьировалась в пределах ([HNO3] = 0.012 - 0.045 M).

База 2D-структур лигандов была сформирована нами с помощью программы ChemAxon Instant JChem версии 5.12.0 и экспортирована в файл в формате Tripos Mol2. Так как структура соединения №47 ((5,6 - ди (пиридин-2-ил) - 1,2,4 - триазин) - (5,6 - ди (пиридин-2-ил) - 1,2,4 - триазин)) в значительной степени отличалась от скелета шаблона, выбранного нами для процедуры выравнивания при построении трехмерной геометрии лигандов, мы исключили его из выборки. Таким образом, исследованная нами база производных полиазагетероциклических соединений включала 47 соединений. База данных приведена в Приложении, п. 4.

Для понимания того, какие атомы лиганда участвуют в образовании координационных связей с ионом металла, мы использовали взятые из работы [111] данные рентгеноструктурного анализа комплекса 5,6 - диэтил-3-пиридин - 1,2,4 - триазина с Eu⁺³, имеющего код FIYFID в Кембриджской базе данных. Структура этих лигандов наиболее близка к структурам соединений исследуемой нами выборки. На основе анализа структуры комплекса мы уточнили, что связывание лиганда с катионом Eu⁺³ происходит при участии 2-го атома азота 1,2,4 - триазинового кольца.

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

Для построения 3D геометрии лигандов в комплексе с ионами Am⁺³/Eu⁺³ мы воспользовались входящей в состав комплекса ChemAxon программой MarvinSketch 5.11.3. С её помощью мы в интерактивном режиме строили и редактировали 3D cтруктуры комплексов органических лигандов с Eu⁺³ при соотношении 2:1. Конформацию лигандов в связанном с металлом состоянии определяли на основе реализованного в комплексе ChemAxon подхода, который включает в себя методы многомерной метрической геометрии и использование библиотеки фрагментов [112]. Далее, чтобы упростить процедуру пространственного выравнивания структур, один из лигандов удаляли и результаты сохраняли в файл.

Далее базу лигандов, связанных с Eu⁺³, выравнивали с использованием структурного шаблона, представленного на Рис. 11. Он является наибольшим общим фрагментом, содержащимся во всех соединениях базы, и поэтому он был использован для ее выравнивания. Поиск наибольшего общего фрагмента проводили с помощью специального приложения Library MCS, входящего в состав программного комплекса ChemAxon 5.11.3.

Выравнивание проводили с использованием алгоритма совмещения структур с использованием шаблона по методу наименьших квадратов [113]. Алгоритм был реализован нами в виде скрипта на языке R (mcmf-do-alignment.R) в рамках комплекса программ CMF https://sites.google.com/site/conmolfields. Для создания шаблона использовали скрипт mcmf-do-make-template.R.

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

С использованием этой базы выравненных пространственных структур лигандов в рамках дипломной работы с помощью метода CMF с использованием предложенных нами непрерывных индикаторных полей нами были построены 3D QSPR модели для фактора разделения для комплексов набора полиазагетероциклических лигандов с Am⁺³ и Eu⁺³. Число индикаторных полей определяли по числу присутствующих в базе данных молекулярно-механических типов атомов (в рамках параметризации силового поля Tripos).

Статистические характеристики 3D QSPR моделей, полученных методом CMF на основе непрерывных индикаторных полей.

Статистические показатели построенных 3D QSPR моделей практически не зависят от выбора модификации метода KRR. Лучший результат по прогнозированию на независимой контрольной выборке (R²_pred=0.76, RMSE_pred=0.20) был получен для трёх модификаций метода: с гребневой регрессией KRR2 без оптимизации коэффициентов смешения h_f и в случае KRR1 вне зависимости от того, оптимизируются ли коэффициенты смешения.

Следует отметить, что 3D QSPR модели, полученные на основе НИП и метода CMF, дают лучшие результаты прогноза, чем 2D QSPR модели, полученные, согласно данным из литературных источников, с использованием программы Cerius2, а также молекулярных дескрипторов, рассчитанных программой CODESSA-PRO, при построении статистической модели методом множественной линейной регрессии [97]. Параметры 2D QSPR моделей, полученных вышеупомянутыми методами: R²_pred = 0.44, RMSE_pred = 0.34 и R²_pred = 0.73, RMSE_pred = 0.24 [97]. Кроме того, прогнозирующая способность построенных нами 3D QSPR моделей сопоставима с таковой для 2D QSPR модели, полученной с помощью программы ISIDA с использованием подструктурных молекулярных фрагментов (substructural molecular fragments, SMF), и статистического метода множественной линейной регрессии MLR - R²_pred = 0.76, RMSE_pred = 0.22 [97].

Как уже было отмечено выше, использование непрерывных индикаторных полей в рамках метода CMF и для построения 3D QSPR моделей позволяют оценивать вклады атомов определенного типа в активность (свойства) молекул. Представлена визуализация перекрывания полей регрессионных коэффициентов 3D QSPR модели (KRR2, без оптимизации h_f) и молекулярных индикаторных полей для атомов, принадлежащих типам C.2 (углерод гибридизации sp²) и C.3 (углерод гибридизации sp³), для наиболее и наименее активных соединений. Красным цветом отмечены области предпочтительного для проявления активности расположения атомов данного типа, синим показаны области, в которых присутствие атомов данного типа нежелательно. Молекулярные поля для соответствующих типов атомов представлены сеткой.

В случае наиболее активного соединения области непрерывных индикаторных полей атомов обоих типов перекрываются с областями положительных значений (красного цвета) регрессионных коэффициентов модели. В частности, из рассмотрения перекрывания полей регрессионных коэффициентов, соответствующих типу C.2 (sp²-гибридизованному углероду) с соответствующими индикаторными полями, следует важное значение, которое имеет для высокого значения коэффициента разделения Am⁺³/Eu⁺³ наличие дополнительного пиридинового кольца, атом азота которого способен координироваться с атомом металла. Это кольцо присутствует в наиболее активном соединение (слева), но отсутствует в наименее активном (справа). С другой стороны, в случае наименее активного соединения (справа) перекрывание полей положительного значения (красного цвета) регрессионных коэффициентов модели с соответствующими непрерывными индикаторными полями атомов обоих типов либо вообще отсутствует (как в случае типа C.2), либо, наоборот, присутствует перекрывание непрерывных индикаторных полей с областями отрицательного значения (синего цвета) соответствующих полей регрессионных коэффициентов (как в случае типа C.3). В частности, в наименее активном соединении присутствует отрицательное перекрывание непрерывного индикаторного поля sp³-гибридизованного атома углерода, входящего в состав присоединённой к пиридиновому кольцу метильной группы, с областью отрицательных регрессионных коэффициентов (синего цвета), свидетельствующее о нежелательности присутствия в ней атомов такого типа. Подобное отрицательное перекрывание присутствует в наименее активной структуре, но отсутствует в наиболее активной. Следует также отметить, что рассмотренная выше область, неблагоприятная для sp³-гибридизованного атома углерода, совпадает в пространстве с областью, благоприятной для наличия sp²-гибридизованного атома углерода, которая, как следует из приведенного до этого анализа, свидетельствует о благоприятном для повышения коэффициента разделения присутствии дополнительного пиридинового кольца, образующего координационную связь с атомом металла. Отсюда следует, что для повышения активности метильная группа в соединении с низкой активностью, находящаяся внутри соответствующей зоны отрицательных значений поля регрессионных коэффициентов для типа C.3 должна быть заменена на присоединённое через орто-положение пиридиновое кольцо. Также из рассмотрения приведённых на Рис. 10 полей регрессионных коэффициентов для типа C.3 следует, что одним из возможных способов повышения активности самого активного соединения могла бы стать замена метоксифенильной группы на метильную.

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

2.2 3D QSAR/QSPR моделирование на основе непрерывных молекулярных полей физико-химической природы

Для того чтобы оценить сравнительные преимущества и недостатки использования в рамках метода CMF непрерывных индикаторных и физико-химических полей, мы осуществили построение 3D QSAR/QSPR моделей для всех рассмотренных выше баз данных (т.е. для связывания с 8 фармакологическими мишенями и для разделения Am⁺³/Eu⁺³) тем же методом, но c заменой индикаторных полей на физико-химические.

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

Все данные, использованные для моделирования биологической активности в отношении 8 фармакологических мишеней, включая наборы структур и активностей, разделение на обучающие и контрольные выборки, пространственную геометрию и выравненные структуры описаны в разделе III.1.1. Значения частичных зарядов на атомах были взяты из приложения к статье [102]. При построении моделей CMF мы использовали 5 типов физико-химических полей - электростатическое, стерическое, гидрофобное, а также поля кислотности и основности по отношению к образованию водородных связей. Параметризация молекулярных полей описана в Экспериментальной части.

Статистические характеристики 3D QSAR моделей, полученных с использованием обучающих и контрольных выборок для восьми наборов лигандов с помощью метода CMF и физико-химических непрерывных молекулярных полей, приведены в Табл. 9. В таблице представлены четыре статистических показателя: q² и RMSE_cv, характеризующие прогнозирующую способность моделей, определённую на основе 10-кратного скользящего контроля, и R²_pred и RMSE_pred, характеризующие прогнозирующую способность моделей на независимой контрольной выборке.

Как видно из таблицы, для баз ACE и THERM модели с наилучшей прогнозирующей способностью на независимой выборке получены методом KRR1 (без оптимизации h_f), для баз AChE, GPB и DHFR - методом KRR2 (без оптимизации h_f), для базы BZR - KRR1 (с оптимизацией h_f) и для баз COX2 и THR - KRR2 (с оптимизацией h_f).

Разности в значениях q² (характеризующих прогнозирующую способность моделей, рассчитанную на основе скользящего контроля) и R²_pred (характеризующих прогнозирующую способность моделей, рассчитанную на основе независимой контрольной выборки) для CMF и моделей, описанных в литературе на основе стандартных методов 3D QSAR, представлены в Табл. 12.

Модели, построенные для исследованных баз с помощью метода CMF, практически во всех случаях имеют лучшую прогнозирующую способность на скользящем контроле (более высокие значения q²), чем соответствующие модели, полученные методами CoMFA, CoMSIA1 и CoMSIA2. Исключение составляют модели для BZR, полученные методами CoMSIA1 и CoMSIA2 [102].

Особенно велика разность между значениями q2 между моделями CMF и CoMFA для баз GPB (0.27). Существенная разница в значениях q² между моделями CMF и CoMSIA1 достигается для баз GPB (0.26) и DHFR (0.24). Разность между моделями CMF и CoMSIA2 лучшая для базы THER (0.12).

Модели, построенные для исследованных баз с помощью метода CMF, имеют лучшую прогнозирующую способность на независимой контрольной выборке (более высокие значения R²_pred) для баз ACE, AChE, BZR, DHFR и THR, чем соответствующие модели, полученные методами CoMFA, CoMSIA1, и CoMSIA2. Для базы GPB значения параметров R²_pred CMF-моделей лучше, чем для моделей, полученных методами CoMFA и CoMSIA1, но хуже, чем для моделей, полученных методом CoMSIA2.

Следует отметить небольшое преимущество в прогнозирующей способности (параметры R²_pred и RMSE_pred, рассчитанные на основе контрольных выборок) для CMF моделей относительно CoMFA и CoMSIA1 моделей для всех выборок, также CoMSIA2 моделей для 6 выборок (ACE, AСhE, BZR, DHFR, GPB и THR).

Представляет интерес сопоставить результаты, полученные методом СMF на основе НИП и непрерывных физико-химических молекулярных полей. Результаты сравнения значений q² и R²_pred лучших моделей, полученных для 8 выборок лигандов с использованием НИП и непрерывных полей физико-химического типа.

Параметр q², описывающий прогнозирующую способность моделей в условиях внутреннего скользящего контроля, для моделей, полученных на основе НИП для 4-х выборок (BZR, DHFR, GPB и THERM), выше по сравнению с q² соответствующих моделей, полученных на основе непрерывных полей физико-химической природы. В то же время, для 3-х выборок (AСhE, COX2 и THR) значения q² моделей на основе физико-химических полей несколько выше. Прогнозирующая способность моделей на независимой контрольной выборке, полученных с использованием физико-химических полей, несколько выше, чем для моделей на основе непрерывных индикаторных полей для 7 выборок (АСЕ, AСhE, BZR, COX2 DHFR, THERM и THR), за исключением GPB.

Таким образом, прогнозирующая способность 3D QSAR моделей, полученных с помощью метода CMF с использованием непрерывных индикаторных полей, вполне сопоставима с прогнозирующей способностью моделей, построенных с помощью традиционных физико-химических полей. Это ещё раз подтверждает выдвинутое ранее утверждение о том, что НИП - это универсальный набор полей, который может быть использован вместо стандартных физико-химических полей для построения моделей 3D QSAR/QSPR. Следовательно, в рамках метода CMF НИП могут рассматриваться как разумная альтернатива и полезное дополнение физико-химическим полям.

2.2.2 3D QSPR моделирование фактора разделения комплексов полиазагетероциклических лигандов с Am⁺³/Eu⁺³ с использованием непрерывных молекулярных полей

Представляло интерес провести сравнение прогнозирующей способности 3D QSPR моделей, построенных с помощью метода CMF на основе непрерывных индикаторных полей и полей физико-химической природы для базы полиазагетероциклических комплексонов для разделения Am⁺³ и Eu⁺³. Для построения 3D QSPR моделей на основе физико-химических полей мы использовали те же самые наборы соединений и схемы выравнивания, как и в случае непрерывных индикаторных полей (см. раздел III.1.2). Для описания молекулярных объектов мы использовали 2 типа физико-химических непрерывных полей - электростатическое и стерическое. При построении моделей поля гидрофобности и водородных связей мы не учитывали, поскольку в образовании комплексов соответствующие им типы взаимодействий участия не принимают.

Прогнозирующая способность 3D QSPR моделей, полученных с использованием НИП, существенно превышает прогнозирующую способность моделей, полученных на основе непрерывных молекулярных полей физико-химической природы. Так, в случае НИП (KRR2 с оптимизацией и без оптимизации коэффициента смешения полей) значение R²_pred равно 0.76, а в случае полей физико-химической природы (KRR2 с оптимизацией коэффициентов смешения h_f) - 0.50. Подобную разницу можно объяснить универсальным характером непрерывных индикаторных полей, которые способны в неявном виде аппроксимировать такие типы взаимодействия, которые недостаточно хорошо могут быть описаны с помощью стандартного набора физико-химических полей, оптимизированного для описания связывания органических лигандов с фармакологическими мишенями.

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

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

3.1 Построение регрессионных 3D QSAR/QSPR моделей методом CMF с использованием разработанных в дипломной работе наборов скриптов

Для построения регрессионных 3D QSAR/QSPR моделей методом СMF был использован набор скриптов, работающих в рамках среды статистического анализа данных R, разработанный ранее на химическом факультете МГУ в лаборатории физико-химических методов анализа строения вещества в группе, руководимой д.ф.-м.н., в.н.с. И.И. Баскиным. В ходе выполнения дипломной работы первоначальный набор скриптов был в существенной мере модифицирован и дополнен возможностью работать с индикаторными полями, а также проводить анализ моделей путём рассмотрения перекрывания полей регрессионных коэффициентов с молекулярными полями при помощи средств интерактивной компьютерной графики. Наборы скриптов были созданы в ходе дипломной работы с использованием программой среды R http://www.r-project.org/.

Для построения 3D QSAR/QSPR моделей при помощи этих скриптов необходимо, чтобы:

1. для каждой из молекул из базы данных была предварительно построена пространственная структура (т.е. для каждого атома должны быть заданы 3 декартовые координаты);

2. пространственные структуры всех молекул были выровнены;

3. структуры всех молекул из выборки были записаны в общий файл формата mol2 с расставленными на атомах значениями частичного заряда;

4. значения моделируемых активностей (свойств) для всех структур содержались в отдельном текстовом файле.

Ниже приведено описание методики построения регрессионных 3D QSAR/QSPR моделей при помощи этого набора скриптов.

3.1.1 Пространственное выравнивание молекул базы данных

В тех случаях, когда используемые при построении моделей пространственные структуры молекул не были изначально выровнены, то необходимое выравнивание осуществлялось при помощи имеющихся в данном пакете скриптов. Для этого сначала проводилось формирование подструктуры-шаблона при помощи скрипта mcmf-do-make-template.R, содержащего следующие параметры:

· mdb-fname - имя файла, содержащего базу данных молекул;

· imol - номер молекулы из базы данных, из которой должна быть извлечена подструктура-шаблон;

· atom_list - список номеров атомов, используемый при извлечении подструктуры-шаблона;

· template_fname - имя файла, в который будет записана подструктура-шаблон.

Далее сформированная подструктура-шаблон использовалась для выравнивания базы данных. Осуществлялось это при помощи скрипта mcmf-do-alignment.R, содержащего следующие параметры:

· mdb_fname - имя файла, содержащего базу данных молекул;

· templ_fname - имя файла, содержащего подструктуру-шаблон;

· mdb_a_fname - имя файла, в который будет записана выровненная база данных молекул;

· algorithm - используемый алгоритм выравнивания (в данном случае «arun» [113]).

3.1.2 Расчёт матрицы ядер для соединений из обучающей выборки