Investigation of the behavior of dioxadet molecules in water by molecular dynamics simulation

Full Text

1. Введение

В конце 1980-х годов в НИИ онкологии им. Н.Н. Петрова был предложен препарат, показавший высокую противораковую активность [1]. К настоящему времени он прошел доклиническое и клиническое исследование и разрешен для медицинского применения под названием “Диоксадэт”. Однако изучение препарата и его гомологов продолжается. Как и все противораковые вещества, он имеет побочные эффекты, ограничивающие дозу применения, хотя в сравнении с препаратами аналогичного спектра действия, такими как цисплатин, побочные эффекты диоксадэта меньше [2, 3].

Несмотря на многолетнюю работу с ДХТ, совершенно не изучены его свойства на молекулярном уровне. Чтобы молекулы препарата попали внутрь клетки, они должны проникнуть через ее мембрану. Однако механизм прохождения ДХТ сквозь мембрану неизвестен. Для других молекул этот вопрос подробно исследуется методами молекулярной динамики [4–8]. Разработаны методы для расчета пассивной диффузии лекарственных молекул через липидные бислои [5, 6], с использованием которых рассчитана проницаемость празиквантела [7] и нифедипина [8].

С другой стороны, до сих пор неизвестно поведение молекул ДХТ в воде. Это касается доставки лекарственного препарата к мембране. С этой целью исследуются вопросы, связанные с растворимостью и ассоциацией, для чего также используется метод молекулярной динамики [9, 10].

Для молекулярно-динамического моделирования требуется задание межмолекулярных взаимодействий (поля сил) для моделируемого вещества. Для этого необходимо сконструировать трехмерную структуру моделируемой молекулы согласно типам атомов, установить их связи, задать ванн-дер-ваальсовские потенциалы, рассчитать частичные заряды на атомах для воссоздания правильного электростатического взаимодействия. Для многих молекул параметры взаимодействий найдены и апробированы в предыдущих исследованиях, и их можно найти в специализированных базах данных (библиотеках полей сил). Например, для моделирования органических жидкостей и растворов популярными являются поля сил OPLS, GROMOS, CHARMM, Amber. Они достаточно универсальны, но охватить всего многообразия молекул невозможно. Поэтому разрабатываются специальные подходы, которые помогают параметризовать новые молекулы. Они базируются на использовании какого-либо известного поля сил, однако приходится проводить специальные расчеты для уточнения строения новой молекулы и нахождения распределенных на ее поверхности зарядов на молекуле.

Насколько нам известно, молекулярно-динамическое моделирование ДХТ до сих пор не проводилось. В данной работе представлена первая попытка такого моделирования. Для этого, прежде всего, мы проводим параметризацию молекулы ДХТ с помощью трех разных доступных сервисов: ATB (Automatic Topology Builder, поле сил GROMOS), SwissParam (поле сил CHARMM) и AmberTools (поле сил Amber). В последнем мы используем также свои квантовохимические расчеты.

2. Модели молекулы диоксадэта

Трехмерная структура молекулы DXT, представленная на рис. 1 взята из базы данных PubChem [11]. Характерной чертой этой молекулы является наличие триазинового кольца, содержащего три атома азота. На нем присутствуют заряды, которые могут быть заметными по величине. По-видимому, эти заряды в значительной мере определяют основные черты поведения DXT в воде, его гидратную оболочку и ассоциацию. Заметим, что в молекуле имеются как гидрофобные участки, так и гидрофильные, поэтому, глядя на молекулу, трудно сразу сказать, какой тип взаимодействия у нее с водой.

Рис. 1. Структура молекулы диоксадэта. Цифрами отмечены атомы, на которых концентрируются наибольшие частичные заряды.

Для параметризации молекулы DXT мы использовали три разных подхода, которые широко применяются в настоящее время.

Подход ATB (Automated force field Topology Builder) [12]. На вход сервису подается файл со структурой молекулы. С использованием метода HF/STO-3G устанавливается оптимальная геометрия и рассчитываются распределенные заряды на атомах с помощью метода Малликена. Для подготовки файла топологии молекулы берется поле сил GROMOS. Используемые здесь квантовохимические расчеты весьма приближенные, поэтому рекомендуется всегда проводить проверку, насколько полученная молекулярно-динамическая модель подходит для решаемой задачи. Соответствующую модель DXT мы называем ATB.

Подход SwissParam [13]. На вход сервису также подается файл со структурой молекулы, однако принципиальная особенность данного метода в том, что здесь не проводятся квантовохимические расчеты. Вместо этого оптимальная геометрия и заряды подбираются из имеющейся базы данных, содержащей типичные фрагменты молекулы. Сконструированная таким образом молекула оптимизируется под заданные критерии, например, чтобы молекула была нейтральной. Преимущество этого сервиса заключается в малых временны́х затратах. Однако надежность его нельзя считать высокой; в частности, следует с осторожностью относиться, например, к распределению частичных зарядов по молекуле. Модель DXT, созданную с его помощью, мы называем SwissParam.

Подход AmberTools [14]. В данном сервисе атомы молекулы типизируются в поле сил Amber. Преимущество этого сервиса состоит в том, что методы и базисы для квантовохимического расчета могут выбираться самостоятельно. По умолчанию, для расчета электронной плотности и оптимальной геометрии молекулы используется метод Хартри–Фока (HF). Однако он обычно оказывается недостаточно точным, поэтому мы также использовали метод B3LYP — один из наиболее популярных для органических молекул. Таким образом, в данном сервисе мы создали две модели DXT, которые назвали AmberHF и AmberB3LYP соответственно. В обоих случаях использован один и тот же базис–6-31G*.

Итак, мы получили четыре модели DXT. Файлы топологии для них доступны на сайте: https://doi.org/10.5281/zenodo.7251232. Заряды на атомах триазинового кольца для всех моделей приведены в табл. 1. Прежде всего отметим их большие значения. Они больше, чем на соответствующих атомах в других частях молекулы. Это может говорить об особой роли триазинового кольца во взаимодействии DXT с окружающей водой и друг другом. Видно, что разные модели дают различные значения заряда для тех же атомов.

Таблица 1. Заряды на атомах триазинового кольца DXT (нумерацию атомов см. на рис. 1)

Отметим необычно большие заряды для модели AmberHF. Мы склонны считать, что они завышены из-за невысокой точности методаHF. Для модели SwissParam наблюдаются наименьшие заряды, равные между собой для всех атомов азота (атомы N1, N2, N3) и углерода (С4, С5, С6). Равенство зарядов кажется странным, так как данное кольцо не является полностью симметричным. Полученный результат, очевидно, связан с тем, что здесь не проводится квантовохимических расчетов, а используются готовые данные от фрагментов молекул. Имеющаяся симметрия молекулы проявляется в зарядах обеих моделей Amber. Действительно, атомы азота N2 и N3 можно считать эквивалентными, так же как атомы углерода С4 и С6 (см. рис. 1 и табл. 1). Однако эта симметрия молекулы не нашла отражения в модели ATB, в которой все атомы кольца заряжены по-разному. Это мы считаем важным на этапе расчета частичных зарядов на молекуле.

3. Молекулярно-динамическое моделирование

Для всех созданных моделей молекул DXT проведено их молекулярно-динамическое моделирование в SPC-модели воды [15]. Использован пакет Gromacs 2021.6 [16]. После помещения молекулы в раствор сначала проводилась минимизация энергии системы, затем –релаксация в NVT-ансамбле с термостатом v-rescale [17] и в NPT-ансамбле с термостатом Нозе–Хувера [18] в течение 2 нс. После этого рассчитывалась равновесная траектория, соответствующая 100 нс. Температура и давление для всех моделей были взяты равными 310 К и 1 бар соответственно. Количество молекул воды составляло 2150 для растворов с одной молекулой DXT (ребро модельного куба равно 4 нм), 2950 и 2920 — для двух и трех молекул (ребро куба — 4.5 нм) и 3996 молекул для четырех молекул (ребро — 5 нм). При моделировании нескольких молекул они помещались на максимально далекое расстояние друг от друга.

4. Результаты

Энергия растворения

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

$\Delta H_{DXT}=⟨H_{solution}-⟨H_{water}⟩⟩.$

На каждом шаге траектории рассчитывается разность между текущей энергией раствора, H_solution, и средним значением энергии того же количества воды без растворенных молекул, ⟨H_water⟩, и проводится усреднение по всей траектории. Такой расчет мы сделали для одной, двух и трех молекул в растворе. В табл. 2 показаны соответствующие значения для наших моделей DXT.

Таблица 2. Энергия растворения DH_DXT разного числа (одна, две, три) молекул DXT в воде (столбцы с номерами 1, 2, 3) для разных моделей DXT в единицах кДж/моль. Указанный разброс значений отражает стандартное отклонение при усреднении по траектории

Заметим, что эта величина не является, строго говоря, энтальпией растворения. В рассчитанную энергию дает вклад также взаимодействие между молекулами DXT. Кроме того, каждая молекула DXT влияет на окружающую воду, которая может несколько иначе взаимодействовать с другой молекулой DXT.

Мы видим, что приведенные значения возрастают (с учетом большого разброса) пропорционально числу молекул для всех моделей. Однако мы хотим обратить внимание на различие между моделями. Видно, что для модели SwissParam значения положительны, т.е. выпадают из общего ряда. Можно сказать, что эта модель “проявляет гидрофобные свойства”, в отличие от других моделей, что, по-видимому, есть следствие относительно малых зарядов на триазиновом кольце (табл. 1). Однако этот вопрос требует дополнительных исследований. Обратим также внимание, что модель AmberHF дает наибольшие по модулю отрицательные значения , что, очевидно, связано с наибольшими значениями зарядов на ее триазиновом кольце.

Водородные связи DXT

Вода может завязывать водородные связи с атомами азота и кислорода DXT, а также атомами водорода гидроксильной группы и амина (–NH–). Водородные связи определялись по стандартному критерию, используемому в пакете GROMACS (угол водород–донор–акцептор больше 30° и максимальная длина связи меньше 0.35 нм). Количество водородных связей на одну молекулу DXT, найденное усреднением по всей траектории, приведено ниже:

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

Мы также изучили водородные связи между молекулами DXT. Несмотря на то, что молекула содержит атомы, способные образовать межмолекулярные водородные связи, оказалось, что они практически отсутствуют. За время, составляющее 20% от всего времени моделирования, удается наблюдать одну водородную связь между молекулами. В остальное время, несмотря на то, что молекулы могли находится рядом, водородные связи отсутствовали. Это означает, что преимущественная взаимная ориентация молекул DXT в растворе не обеспечивает контактов между атомами, способными к образованию водородных связей.

Гидратная вода вокруг триазинового кольца

Анализ проводился для одиночной молекулы DXT, растворенной в воде. Рассчитывалась функция радиального распределения (ФРР) центров молекул воды относительно поверхности атомов триазинового кольца DXT (рис. 2) и изучена их ориентация с использованием подхода, описанного в работе [19].

Рис. 2. Функция радиального распределения (ФРР) молекул воды относительно поверхностей ближайших атомов триазинового кольца DXT; ФРР всех моделей стремятся к единице на больших расстояниях, как для модели AmberHF. Кривые для AmberB3LYP, ATB и SwissParam смещены вниз для удобства их рассмотрения.

Обратим внимание на ближайшие пики на рис. 2: первый при 0.1 нм и второй вблизи 0.2 нм, четко проявляющиеся для моделей AmberHF, AmberB3LYP и АТВ. Дополнительный анализ показал, что самые близкие молекулы воды связаны водородной связью с атомами азота кольца, а второй пик образован молекулами воды, связанными водородной связью с молекулами первого пика подобно тому, как это показано на рис. 6 из работы [19]. Заметный пик при 4.5 нм соответствует молекулам воды, часть которых связана водородными связями с молекулами, образующими второй максимум. С увеличением расстояния, ФРР для этих моделей неплохо совпадают. Это понятно, поскольку влияние особенностей разных моделей молекул DХТ на далекую воду уменьшается. Модели AmberB3LYP и AmberHF имеют практически совпадающие функции радиального распределения. Некоторые различия для модели АТВ проявляются в меньшей интенсивности первых пиков и небольшом сдвиге второго максимума в сторону больших расстояний, что можно связать с тем фактом, что все заряды на кольце в этой модели, в отличие от первых двух, различны (табл. 1). Для модели SwissParam вместо двух первых пиков наблюдается один общий пик при 0.15 нм. Это обусловлено, по-видимому, меньшими зарядами на кольце (табл. 1), что приводит к ослаблению водородных связей воды с DXT.

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

Ассоциаты

Две и более молекул DXT в воде склонны подходить друг к другу и находиться рядом в течение достаточно долгого времени. Для количественного описания такого поведения мы следили за изменением расстояния между центрами масс пары молекул со временем. На рис. 3 показана динамика относительного движения двух молекул DXT. Изначально мы помещали молекулы на большом расстоянии друг от друга. Однако уже за время релаксации они обычно сближаются, а сблизившись, остаются вместе, о чем говорит расстояние между ними, равное примерно 0.5 нм. Из рис. 3а видно, что молекулы из модели АТВ в течение всего времени моделирования (100 нс), находились рядом. Однако на рис. 3б–3г видны резкие выбросы, указывающие на то, что молекулы отходили друг от друга на большие расстояния. Однако основную часть времени молекулы находятся рядом, в первую очередь это относится к модели AmberB3LYP (рис. 3б).

Рис. 3. Зависимости расстояния между центрами масс двух молекул DXT от времени для разных моделей: а — ATB, б — Amber B3LYP, в — Amber HF, г —SwissParam.

Как было отмечено выше, водородные связи между молекулами DXT практически не образуются, т.е. ассоциация в данном случае полностью определятся электростатическим взаимодействием. Мы заметили, что пара близких молекул DXT часто располагается так, чтобы их триазиновые кольца были параллельны. В случае ассоциатов из трех и четырех молекул обычно присутствует хотя бы одна такая пара. Более того, эта тенденция иногда приводит к образованию стопок из трех и четырех молекул. На рис. 4 показаны кластеры из четырех молекул DXT, взятые из молекулярно-динамической траектории. Как привило, возникает бесформенный ассоциат (рис. 4а), однако иногда можно видеть хорошо сформированные стопки (рис. 4б). Приведенный структурный анализ показан на примере модели ATB. Однако для других моделей наблюдается похожий вид ассоциатов.

Рис. 4. Скриншоты, демонстрирующие агрегаты из четырех молекул DXT в воде (молекулы воды не показаны), возникающие в модели ATB.

5. Заключение

В работе проведена параметризация молекулы диоксадэта для полноатомного классического молекулярно-динамического моделирования. Для этого использованы различные доступные сервисы: ATB, SwissParam, а также AmberTools, в стандартном виде и с использованием дополнительных квантовохимических расчетов. Проведено моделирование молекулы диоксадэта в воде и сравнение результатов для разных моделей между собой. Несмотря на то, что полученные частичные заряды на триазиновом кольце молекулы и энергии растворения заметно различаются для разных моделей, структурные свойства качественно совпадают, что видно по схожему поведению ФРР гидратной воды вблизи триазинового кольца. Дополнительно проведено моделирование двух, трех и четырех молекул в воде. Показано, что молекулы диоксадэта в воде склонны образовывать лабильные ассоциаты, среди которых иногда встречаются стопки молекул. Все модели демонстрируют сходные свойства. Для выбора лучшей модели молекулы диоксадэта необходимы дальнейшие исследования. Файлы с топологией молекулы доступны для использования на сайте https://doi.org/10.5281/zenodo.7251232.

About the authors

E. A. Yakush

Email: nikmed@kinetics.nsc.ru
Russian Federation, Novosibirsk; Novosibirsk

A. V. Kim

Email: nikmed@kinetics.nsc.ru
Russian Federation, Novosibirsk; Novosibirsk

N. N. Medvedev

Author for correspondence.
Email: nikmed@kinetics.nsc.ru
Russian Federation, Novosibirsk

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Fig. 1. Structure of the dioxade molecule. The numbers indicate the atoms on which the largest partial charges are concentrated.

Download (50KB)

Indexing metadata

3. Fig. 2. Radial distribution function (RDF) of water molecules relative to the surfaces of the nearest atoms of the triazine ring DXT; The RDFs of all models tend to unity at large distances, as for the AmberHF model. The curves for AmberB3LYP, ATB and SwissParam are shifted down for ease of viewing.

Download (107KB)

Indexing metadata

4. Fig. 3. Time dependence of the distance between the centers of mass of two DXT molecules for different models: a - ATB, b - Amber B3LYP, c - Amber HF, d - SwissParam.

Download (191KB)

Indexing metadata

5. Fig. 4. Screenshots showing aggregates of four DXT molecules in water (water molecules not shown) arising in the ATB model.

Download (133KB)

Indexing metadata

Х Международная конференция им. В.В. Воеводского “Физика и химия элементарных химических про­цессов” (сентябрь 2022, Новосибирск, Россия).