Вычислительная алхимия
В классическом молекулярном моделировании решается задача наиболее полного описания исследуемой системы: хорошая модель должна быть столь физична, сколь позволительно этого добиться при имеющихся вычислительных ресурсах.
КМ можно достичь субатомного описания относительно небольших молекул; применение КЗ полей позволяет приближенно рассчитывать взаимодействия целых ансамблей макромолекул; насыщенные мицеллярные растворы, можно описывать уравнениями фазового состояния. В описанных ситуациях реальным физическим явлениям ставится в соответствие некоторый математический формализм, объясняющий их. Принципиально иной подход положен в основание методов вычислительной алхимии, суть которых заключается в моделировании нефизичного процесса превращения атомов одной природы в другую, при этом по крайней мере нарушается закон сохранения материи (массы).
Алхимический формализм одинарной и двойной топологии [1, с. 50]
От примера с оценкой энергии связывания фермента E и субстрата L мы перешли к интегрированию по пути \(\lambda\) между двумя произвольными состояниями в ФП. По определению: точки ФП — это \(6N\)-мерные векторы, состоящие из компонент импульсов \(\bf p\) и координат \(\bf q\) всех частиц в системе, ФП описывает состояние физичной системы в смысле классической механики, однако мы можем обогатить ФП новыми степенями свободы, включив в него массу, заряд и другие свойства частиц. В расширенном ФП можно построить путь \(\lambda\), соединяющий состояния \(\cl A\) и \(\cl B\), отвечающие разным химическим молекулам, фактически — смоделировать алхимическую трансформацию.
На практике расчёт энергии перехода из одного состояния в другое проводится по Ур. 9.2, однако этот расчёт можно реализовать по-разному. В дальнейших рассуждениях без ограничения общности будем считать, что алхимической трансформации \(\cl A \to \cl B\) подвергается молекула лиганда.
В МД потенциальная энергия \(U(\bf r_1, \cdots, \bf r_N)\) вычисляется по Ур. 1.3. Вводимый путь алхимической трансформации \(\lambda\), параметризующий превращение лиганда, может включаться в выражение СП, и тогда топология (числовые параметры атомов СП и порядок их связей) трансформируемой молекулы плавно изменяется при движении в ФП, т.е. свойства атомов (масса \(m\), заряд \(q\) и др.) находятся в функциональной зависимости от пути \(\lambda\): \[ \begin{cases} q^i(\lambda) = \lambda q^i_{\cl A} + (1 - \lambda) q^i_{\cl B}, \\ \quad \vdots \\ m^i(\lambda) = \lambda m^i_{\cl A} + (1 - \lambda) m^i_{\cl B}. \\ \end{cases} \]
В другой реализации, расчёт проводится непосредственно по Ур. 9.2, тогда в процессе трансформации молекула лиганда будет представлять собой суперпозицию начального \(\cl A\) и конечного \(\cl B\) состояний, каждая из которых будет описываться своей топологией. В этой постановке всё ФП или его область являются функцией от \(\lambda\).
В обоих подходах в начале (\(\lambda = 0\)) и в конце (\(\lambda = 1\)) пути существую физичные молекулы, однако в формализме одинарной топологии в процессе трансформации параметры атомов, будь то заряд, масса или др., будут принимать промежуточные значения, т.е. не будут соответствовать ни одному настоящему атому. В формализме двойной топологии очевидно возникает проблема: потенциальная энергия молекулы в процессе трансформации согласно Ур. 9.2 описывается в виде взвешенной суммы. При попытке поместить две молекулы в одну область пространства неизбежно произойдёт расталкивание электронных облаков1, поэтому также необходимо указывать список атомов, общих для обеих структур, чтобы «исключить» их из взаимодействий между собой (и позволить занять одну область пространства), но оставить взаимодействие с окружением, масштабируя его согласно пути \(\lambda\).
1 в МД аппроксимирующихся зарядами и потенциалами взаимодействия атомов
Интегрирование по пути
Для проведения алхимической трансформации необходимо просемплировать ФП в разных точках алхимического пути \(\lambda = 0..1\) и затем оценить значение интеграла по формуле Ур. 9.3. Поскольку интегрирование проводится численно, узлы интегрирования (значения \(\lambda\)) могут выбираться из разных предпосылок. При интегрировании методом параллелограммов значения \(\lambda\) могут браться с определённым шагом \(\Delta \lambda\) или вовсе произвольно. Более продвинутые методы интегрирования могут накладывать ограничения на выбор \(\lambda\); часто используемый метод Гаусса2 уже предоставляет набор значений \(\lambda\), который необходимо использовать в расчётах, итоговое значение интеграла вычисляется по простой формуле \[ \int_a^b f(\lambda) \dd \lambda \approx \sum_{i=1}^n w_i f(\lambda_i), \] где \(w_i\) — известные из метода Гаусса значения весов. Например, для числа узлов \(n = 9\), значения \(\lambda_i\) и \(w_i\) приведены в Табл 9.1. Значения \(\lambda_i\) лежат симметрично относительно середины пути \(\lambda = 0.5\); веса \(w_i\) крайних точек в итоговом интеграле минимальны.
2 англ. Gaussian quadrature
| \(i\) | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 |
|---|---|---|---|---|---|---|---|---|---|
| \(\lambda_i\) | 0.0159 | 0.082 | 0.1933 | 0.3379 | 0.5 | 0.6621 | 0.8067 | 0.918 | 0.9841 |
| \(w_i\) | 0.0406 | 0.0903 | 0.1303 | 0.1562 | 0.1651 | 0.1562 | 0.1303 | 0.0903 | 0.0406 |
Расчёт энергии связывания
Метод алхимической трансформации позволяет оценивать энергию превращения одной молекулы в другую путём изменения параметров атомов или преобразованиями функции внутринней энергии \(U\). Получаемые в результате алхимии значения энергии \(\Delta_\t{alc} F\) не могут быть проверены в эксперименте непосредственно, но могут использоваться для получения ТД параметров связывания.
Опишем ТД цикл (Рис 9.1) для определения разности энергий связывания \(\Delta \Delta_b G\) пары лигандов A и B, в активном центре фермента E. Для трансформаций лигандов A \(\to\) B и комплексов EA \(\to\) EB соответствующие значения алхимических энергий \(\Delta_\t{alc} F_L\) и \(\Delta_\t{alc} F_C\) могут быть рассчитаны, а значения энергий связывания \(\Delta_b G_{\rm A}\) и \(\Delta_b G_{\rm B}\) определены в эксперименте. Исходя из ТД цикла и с учётом знаков
\[ \Delta_\t{alc} F_C = -\Delta_b G_{\rm A} + \Delta_\t{alc} F_L + \Delta_b G_{\rm B}, \] откуда \[ \ub{\Delta_\t{alc} F_C - \Delta_\t{alc} F_L}_{\Delta_\t{calc} \Delta_b F} = \ub{\Delta_b G_{\rm B} - \Delta_b G_{\rm A}}_{\Delta_\t{exp} \Delta_b G}, \] т.е. разница энергий связывания \(\Delta \Delta_b G\) может быть определена экспериментально или рассчитана.
Моделирование белка (в данном случае в комплексе с лигандом) (Рис 9.1 б) требует соответствующего размера ячейки, заполненной водой. Для корректного сопоставления энергии, размер систем (количество частиц3) при превращении лигандов должен соответствовать системе белка (комплексов), однако, при превращении свободных лигандов фермент должен быть удалён на бесконечное расстояние ({Рис 9.1}а), а поскольку нас интересует только разность энергий \(\Delta_\t{alc} F_L\) (энергия алхимического превращения), из системы можно исключить фермент и всю необходимую его воду, имеющие одинаковый энергетический вклад до и после превращения; проще говоря, моделирование лигандов можно проводить в ячейке малого размера, что значительно экономит время расчёта.
3 уже нарушается в процессе алхимической трансформации
Рассчитанные значения \(\Delta \Delta F_b\) могут быть использованы как поправки: энергия связывания с лигандом B может быть рассчитана через экспериментальное значение энергии связывания \(\Delta_b F_A\) и вычисленную поправку — энергию алхимической трансформации \(\Delta \Delta_b F\):
\[ \Delta_b G_{\rm B} = \Delta_b G_A + \Delta_{\rm A \to \rm B} \Delta_b F, \] последнее верно в силу того, что и \(F\), и \(G\) являются функциями состояния, более того, в конденсированных фазах \(\Delta F \approx \Delta G\).
Резюмируем: имея в распоряжении экспериментально измеренную энергию связывания с помощью метода алхимической трансформации можно рассчитать поправки и определить энергию связывания для целого набора химически близких лигандов, что делает описанный метод мощным инструментом для поиска ингибиторов и лекарственных молекул.