Химия · Аналитическая и физическая химия · 2 курс

Компьютерная химия и докинг

Как новое лекарство «испытывают» на компьютере до синтеза в колбе?

Ситуация

В Институте биоорганической химии АН РУз в Ташкенте изучают вещества, выделенные из растений, и их производные как лекарства. Прежде чем синтезировать новое соединение, химик «строит» его на компьютере: рисует молекулу, получает её 3D-модель, минимизирует энергию, рассчитывает колебательный спектр и «вставляет» в карман фермента. Эта работа выполняется в программах вроде Avogadro, HyperChem, Gaussian. Насколько точно компьютер предсказывает эксперимент, какой метод подходит для какой задачи и чему в результатах можно доверять?

Тема

Компьютерное моделирование в химии. Компьютерный класс Института биоорганической химии АН РУз: кандидат в новое лекарство «строится» на компьютере ещё до того, как попадёт в лабораторию — рисуется молекула, переводится в 3D…

Ключевые понятия

Силовое поле и оптимизация
Квантово-химический расчёт
Докинг и QSAR

Формулы

E = \sum k_{b} {(r - r_{0})}^{2} + \sum k_{θ} {(θ - θ_{0})}^{2} + \sum \frac{V_{n}}{2} (1 \pm \cos n φ) + \sum ε [{(\frac{r^{*}}{r})}^{12} - 2 {(\frac{r^{*}}{r})}^{6}] + \sum \frac{q_{i} q_{j}}{4 π ε_{0} ε_{r} r}

Силовое поле (молекулярная механика): слагаемые связей, углов, торсионные, ван-дер-ваальсовы и электростатические.

𝐫_{k + 1} = 𝐫_{k} - s \nabla E

Оптимизация геометрии: наискорейший спуск (в методе сопряжённых градиентов направление «помнит» и предыдущий шаг).

\hat{H} Ψ = E Ψ, 𝐅 𝐂 = 𝐒 𝐂 𝜺

Уравнение Шрёдингера и уравнения Хартри–Фока–Рутаана: орбитали как сумма базисных функций.

\tilde{ν} = \frac{1}{2 π c} \sqrt{\frac{k}{μ}}

Волновое число гармонического колебания: силовая постоянная $k$ и приведённая масса $μ$ (в многоатомной молекуле — собственные значения матрицы Гессе).

Δ G_{св} = R T \ln K_{d}

Докинг: энергия связывания и константа диссоциации ( $K_{d}$ мала — лиганд связывается прочно), $T$ = 298 К.

Теория

Вычислительная химия: что рассчитывается и какие программы используются

Компьютерная химия рассчитывает по теории величины, измеряемые на опыте: длины связей, валентные и торсионные углы, энергии конформеров, энергию связи и порядок связи, теплоту образования, потенциал ионизации (теорема Купманса: $I \approx - ε_{HOMO}$ ), сродство к электрону, энергетическую щель (HOMO–LUMO), дипольный момент, спектры. Результат расчёта всегда сравнивают с экспериментом — отсюда столбцы «расчёт» и «эксперимент» в таблицах на сцене. Программы бывают двух видов: графические (ISIS Draw, ChemWin, ChemDraw) — рисование формул, названия, 2D-изображение; расчётные (Avogadro, HyperChem, Gaussian, ChemOffice/Chem3D) — 3D-геометрия, энергия, спектр. Расчёт начинается с начальной геометрии: 2D-рисунок «поднимается» в 3D (шаблоны VSEPR, z-матрица — длина связи, угол и торсион), а затем оптимизируется.

Молекулярная механика (эмпирический метод) и оптимизация

Молекула рассматривается как шарики, соединённые пружинами: силовое поле — сумма слагаемых растяжения связей, изгиба углов, торсионных, ван-дер-ваальсовых и электростатических. Его параметры ( $k$ , $r_{0}$ , $θ_{0}$ , $V_{n}$ ) подбираются по опыту (силовые постоянные из ИК-спектров, геометрии из электронографии, энергии конформеров) — это параметризация. $r_{0}$ — идеальная (ненапряжённая) длина; реальная длина в оптимизированной молекуле немного отличается от неё из-за влияния соседних атомов. Оптимизация геометрии — поиск точки, где градиент энергии равен нулю: метод наискорейшего спуска прост, но медлен («зигзаг»), сопряжённые градиенты и квазиньютоновский метод (BFGS) быстрее. Оптимизация находит только локальный минимум; для глобального минимума нужен поиск конформеров (систематический скан или много случайных стартов). Силовое поле в этом эксперименте подобрано так, что барьер в этане равен 12,0 кДж/моль, а в бутане гош − анти = 3,8 кДж/моль. В программе Avogadro то же самое делает «Auto Optimize» (силовые поля MMFF94, UFF).

Квантово-химические методы: ab initio, базисы, выбор метода

В методах ab initio («с нуля») экспериментальных параметров нет: в методе Хартри–Фока каждый электрон движется в усреднённом поле остальных, а орбитали ищут как сумму базисных функций (гауссовых функций $e^{- α r^{2}}$ ). С ростом базиса энергия по вариационному принципу уменьшается, подходя сверху к пределу (на сцене: атом H, n = 1…8). Но у гелия даже предел Хартри–Фока на 1,14 эВ выше точной энергии — это энергия корреляции электронов; её учитывают методы MP2, CCSD(T) и теории функционала плотности (DFT, например B3LYP). Правила выбора базиса: минимальный (STO-3G) — быстро, качество низкое; 6-31G* — с поляризующими d-функциями, для большинства органических задач; для анионов и водородных связей — диффузные (+) функции; для точных энергий — cc-pVTZ. Стоимость расчёта быстро растёт: HF ~ $N^{4}$ , MP2 ~ $N^{5}$ , CCSD(T) ~ $N^{7}$ ( $N$ — число базисных функций). Выбор метода: тысячи атомов — молекулярная механика; сотни — полуэмпирические методы или DFT; точная термохимия — CCSD(T). Gaussian — самый распространённый программный пакет для расчётов ab initio и DFT (геометрия, частоты, ЯМР, УФ, термохимия).

Полуэмпирические методы, заряды и УФ-спектр

Полуэмпирические методы упрощают квантовую механику и заменяют часть интегралов параметрами, взятыми из опыта: простейший — π-электронный метод Хюккеля (параметры $α$ и $β$ ), затем расширенный метод Хюккеля, CNDO, MNDO, AM1, PM3, PM7 (HyperChem, MOPAC). Они быстры и позволяют получить энергетические параметры (теплоту образования, потенциал ионизации), распределение заряда, диполь, графическое изображение орбиталей. Парциальные заряды на атомах зависят от метода: Гастайгера (выравнивание электроотрицательности), Малликена, ESP-подгонка; для воды заряды Гастайгера дают диполь 1,2 Д, а силовые поля жидкой воды (TIP3P, SPC/E) — 2,3–2,4 Д (опыт, газ — 1,85 Д): если сменить силовое поле, изменятся заряды и результаты. УФ-спектр: по Хюккелю энергия перехода HOMO→LUMO $Δ E = (x_{HOMO} - x_{LUMO}) | β |$ ; если $| β |$ подогнан по бутадиену (217 нм), то для длинных полиенов длина волны завышается — метод не учитывает чередование связей (точнее — ZINDO, TD-DFT). π-орбитали Хюккеля и ароматичность — в эксперименте «Молекула: орбитали, симметрия».

Расчёт и истолкование спектров

ИК-спектр: матрица вторых производных энергии (матрица Гессе) делится на массы и диагонализуется — из собственных значений получаются частоты нормальных колебаний, из собственных векторов — движение атомов; если при колебании дипольный момент меняется, оно активно в ИК (интенсивность $\propto {(\partial μ / \partial Q)}^{2}$ ), если не меняется — проявляется только в спектре комбинационного рассеяния (Рамана). Гармоническое приближение завышает частоты, поэтому рассчитанные частоты масштабируют (≈ 0,9 для HF, ≈ 0,96 для B3LYP) или подгоняют силовое поле прямо под экспериментальные частоты (в этом эксперименте так сделано для воды). УФ — электронные переходы (TD-DFT, ZINDO), ЯМР — химические сдвиги (метод GIAO; ChemOffice оценивает сдвиги ¹H и ¹³C по таблицам инкрементов — аддитивности), масс-спектр — пути фрагментации. Рассчитанный спектр помогает истолковать экспериментальный: какая линия какому колебанию принадлежит (подробнее — в эксперименте «Спектры: определение строения»).

Биологическая активность: QSAR, PASS, докинг; границы модели

Зависимость структура–активность (SAR/QSAR): активность предсказывают регрессией по дескрипторам молекулы (lgP, молекулярная масса, доноры и акцепторы H-связей, электронные эффекты) (уравнение Ганша); для лекарств, принимаемых внутрь, действует «правило пяти» Липински: M ≤ 500, lgP ≤ 5, доноров ≤ 5, акцепторов ≤ 10. PASS (В. Поройков с коллегами, Москва) по структурным фрагментам молекулы байесовским методом даёт вероятности сотен видов активности: если Pa (активно) > Pi (неактивно), активность вероятна, а при Pa > 0,7 — высоко вероятна. Молекулярный докинг — размещение лиганда в кармане фермента: поиск (Монте-Карло, генетический алгоритм) + оценочная функция (стерическое соответствие, H-связи, десольватация, энтропия); программы вроде AutoDock Vina. Виртуальный скрининг так отсеивает тысячи соединений, выбирая для синтеза лучшие. Границы модели: силовое поле в этом эксперименте упрощено в учебных целях и подогнано под несколько молекул; расчёт ab initio — только для атомов с s-гауссовыми функциями; карман фермента — не настоящий белок, а схематичная модель; энергии докинга и $K_{d}$ — модельные, не экспериментальные. Результат расчёта — гипотеза, её подтверждает эксперимент.

Знаете ли вы?

В 2013 году Нобелевскую премию по химии получили М. Карплюс, М. Левитт и А. Уоршел «за создание многомасштабных моделей сложных химических систем» — в их работах объединены молекулярная механика и квантовая химия.
Нобелевская премия по химии 1998 года была разделена между Вальтером Коном (теория функционала плотности) и Джоном Поплом (вычислительные методы квантовой химии; один из основателей программы Gaussian).

Что внутри эксперимента

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

Другие эксперименты раздела

Спектры: определение строения — 2 курс
Сплавы: диаграммы состояния — 2 курс
Химическая кинетика и уравнение Аррениуса — 2 курс
Катализ и ферменты — 2 курс
Коллоиды и адсорбция — 2 курс