Теория возмущений свободной энергии

Метод возмущений свободной энергии FEP (free energy perturbation) и термодинамическое интегрование TI (thermodynamic integration) — два важных подхода вычислительной химии, используемые для расчёта разностей свободных энергий, например, при оценке энергии связывания лиганда с белком, термодинамики растворимости, конформационной стабильности и т.д.

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

1 Возмущение свободной энергии

Информацию по разделу см. в [1, с. 312]

1.

Tuckerman M.E. Statistical mechanics : theory and molecular simulation [Electronic resource]. Oxford: Oxford University Press, 2010.

Пусть имеем термодинамическую систему, состоящую из молекулы лиганда \mathrm{L} и белка \mathrm{E}: в состоянии \mathcal{A} лиганд \mathrm{L} и белок \mathrm{E} и не связаны, в состоянии \mathcal{B} — имеет место комплекс \rm{ES}, соответственно введём потенциальные энергии U_{\mathcal{A}}(\mathbf{r}) и U_{\mathcal{B}}(\mathbf{r}) для несвязанного и связанного состояний. Свободная энергия Гельмгольца перехода в связанное состояние [1, с. 138] \Delta F_{\mathcal{A}\mathcal{B}} = F_{\mathcal{B}} - F_{\mathcal{A}} \tag{1}

может быть выражена через суммы по состояниям:

\Delta F_{\mathcal{A}\mathcal{B}} = F_{\mathcal{B}} - F_{\mathcal{A}} = - k T \ln \dfrac{Q_{\mathcal{B}}}{Q_{\mathcal{A}}}, \tag{2}

где, например,

\begin{aligned} Q_{\mathcal{A}} &:= \dfrac{1}{N! h^{3N}} \int_\text{ФП} e^{- \beta H(\mathbf{p}, \mathbf{q})} \mathrm{d} \mathbf{p} \mathrm{d} \mathbf{q} = \\ &= \dfrac{1}{N! h^{3N}} \int_\text{ФП} \exp \left\{ - \beta \left[\sum_{i=1}^N \dfrac{\mathbf{p}_i^2}{2m_i}\right] + U_{\mathcal{A}} (\mathbf{q}_1, \cdots, \mathbf{q}_N)\right\} \mathrm{d} \mathbf{p} \mathrm{d} \mathbf{q}, \end{aligned} \tag{3}

где H — гамильтониан системы \mathcal{A}, N — число частиц в ней, \beta = 1 / kT, \mathrm{d} \mathbf{p} \mathrm{d} \mathbf{q} — точка ФП; интегрирование ведётся по всему ФП, что эквивалентно интегрированию по области, соответствующей состоянию \mathcal{A}, т.е. там, где вероятность пронаблюдать микросостояние состояния \mathcal{A} не нулевая.

Сократив в 2 числитель и знаменатель на импульсозависимую часть \exp \{ \sum_i \frac{\mathbf{p}_i^2}{2 m_i} \}, перейдём к конфигурационным интегралам:

\Delta F_{\mathcal{A}\mathcal{B}} = F_{\mathcal{B}} - F_{\mathcal{A}} = -kT \ln \dfrac{Z_{\mathcal{B}}}{Z_{\mathcal{A}}}, \tag{4}

где

Z_{\mathcal{A}} := \int e^{ -\beta U_{\mathcal{A}} (\mathbf{q}_1, \cdots, \mathbf{q}_2) } \mathrm{d} \mathbf{q}. \tag{5}

Оценить \Delta F непосредственно по уравнениям 2 и 4 сложно, поскольку методы МД и Монте-Карло не позволяют непосредственно рассчитывать статистические суммы и интегралы, однако позволяют получить хорошие средние оценки ТД функций.

Уравнение 4 можно переформулировать через средние ТД величины. Конфигурационный интеграл представим в виде интеграла от двух множителей-экспонент: Z_{\mathcal{B}} = \int e^{-\beta U_{\mathcal{B}}} \mathrm{d} \mathbf{q} = \int e^{-\beta U_{\mathcal{B}}} e^{\beta U_{\mathcal{A}}} e^{-\beta U_{\mathcal{A}} } \mathrm{d} \mathbf{q} = \int e^{- \beta (U_{\mathcal{B}}-U_{\mathcal{A}})} \underbrace{ e^{-\beta U_{\mathcal{A}}} }_{\rho_{\mathcal{A}}} \mathrm{d} \mathbf{q}, \tag{6} где \rho_{\mathcal{A}} — функция распределения, соответствующая состоянию \mathcal{A}. Среднее значение произвольной функции y по ансамблю \mathcal{A} можно вычислить через функцию распределения \rho_{\mathcal{A}}: \langle y \rangle_{\mathcal{A}} := \dfrac 1 {Z_{\mathcal{A}}} \int_\text{ФП} \rho_{\mathcal{A}} (\mathbf{p}, \mathbf{q}) \cdot y(\mathbf{p}, \mathbf{q}) \mathrm{d} \mathbf{p} \mathrm{d} \mathbf{q}, \tag{7}

учитывая 6 и 7, дробь в уравнении 4 примет вид:

\dfrac{Z_{\mathcal{B}}}{Z_{\mathcal{A}}} = \dfrac 1 {Z_{\mathcal{A}}} \int e^{-\beta U_{\mathcal{A}}} e^{-\beta (U_{\mathcal{B}} - U_{\mathcal{A})}} = \langle e^{-\beta (U_{\mathcal{B}} - U_{\mathcal{A}})} \rangle_{\mathcal{A}},

т.е. среднее берётся из распределения, соответствующего состоянию \mathcal{A}.

Наконец,

\Delta F_{\mathcal A \mathcal B} = -kT \ln \langle e^{-\beta (U_{\mathcal B} - U_{\mathcal A})} \rangle_{\mathcal A}. \tag{8}

Уравнение 8 известно как возмущение свободной энергии. Опишем, как следует использовать это уравнение. Состоянию \mathcal A отвечает функция распределения \rho_{\mathcal A} и потенциал U_{\mathcal A}, в ходе численного эксперимента мы получаем точки ФП, соответствующие этому распределению, и затем используем их для моделирования состояния \mathcal B, заменив потенциал на U_{\mathcal B}. Однако, распределению \mathcal B соответствует другая функция распределения \rho_{\mathcal B}, поэтому точки ФП следует перевзвесить¹, поделив на исходное e^{-\beta U_{\mathcal A}} и умножив на e^{-\beta U_{\mathcal B}}.

¹ англ. reweighting

Таким образом, при расчёте \Delta F по 8 мы собираем статистику по ФП, соответствующую начальному состоянию, и затем используем эти наблюдения для оценки \Delta F, как если бы они были взяты из распределения состояния \mathcal B. Ограничением этого подхода является то, что микросостояния из распределения \mathcal A могут не быть микросостояниями с высокой вероятностью в распределении \mathcal B, в этом случае разность U_{\mathcal B} - U_{\mathcal A} будет большой и экспоненциальный фактор e^{-\beta (U_{\mathcal B} - U_{\mathcal A})} \to 0, соответствующие состояния будут иметь малый вес в среднем по ансамблю \langle \cdots \rangle_{\mathcal A} и, как следствие, сходимость будет медленной. Другими словами, требуется, чтобы состояние \mathcal B было небольшим возмущением состояния \mathcal A и их области в ФП существенно перекрывались.

В постановке задачи, состоянию \mathcal{A} соответствует несвязанные белок \text{E} и лиганд \text{L}, состоянию \mathcal{B} — фермент-субстратный комплекс \text{EL}, на масштабах системы «лиганд—белок—гидратная оболочка» переход из несвязанного в связанное состояние в принципе можно считать небольшим возмущением.

2 Термодинамическое интегрирование

Выбор начального и конечного состояний не ограничивается связыванием в комплекс, состояния могут соответствовать любым областям в ФП. В случае, когда \mathcal{A} и \mathcal{B} достаточно далеки, вводят промежуточные состояния, переходы между которыми можно считаться малым возмущением. Уравнение 8 в этом обобщении имеет вид: \Delta F_{\mathcal{A}\mathcal{B}} = -kT \sum_{i} \ln \langle e^{-\beta (U_{i+1} - U_i)} \rangle_i, где i — номер промежуточного состояния, U_i — соответствующий ему потенциал, усреднение ведётся по функции распределения e^{-\beta U_i} этого состояния.

Однако, можно пойти дальше и от суммирования дискретных состояний перейти к непрерывному интегрированию. Введём путь трансформации \lambda = 0..1 — величину по сути аналогичную координате реакции между двумя состояниями ТД системы. Пользуясь \lambda мы можем определить потенциальную энергию U(\lambda) промежуточных состояний:

\begin{cases} U(\lambda) = f(\lambda) U_{\mathcal{A}} + g(\lambda) U_{\mathcal{B}}, \\ f(\lambda) \text{ монотонно убывает от $1$ до $0$}, \\ g(\lambda) \text{ монотонно возрастает от $0$ до $1$}. \end{cases}

Из выражения свободной энергии, соответствующей состоянию U(\lambda): F(\lambda) = - k T \ln Q(\lambda), можно получить соотношение для производной свободной энергии F по параметру пути \lambda: \begin{aligned} \frac{\partial F}{\partial \lambda} &= -\frac{kT}{Q} \frac{\partial Q}{\partial \lambda} = -\frac{kT}{Z} \frac{\partial Z}{\partial \lambda} \\ &= -\frac{kT}{Z} \frac{\partial}{\partial \lambda} \int_{\mathrm{ФП}} e^{-\beta U(\mathbf{q}; \lambda)} \, \mathrm{d}\mathbf{q} \\ &= -\frac{kT}{Z} \int_{\mathrm{ФП}} \left(-\beta \frac{\partial U(\mathbf{q}; \lambda)}{\partial \lambda}\right) e^{-\beta U(\mathbf{q}; \lambda)} \, \mathrm{d}\mathbf{q} \\ &= \left\langle \frac{\partial U(\lambda)}{\partial \lambda} \right\rangle_\lambda \end{aligned} \tag{9}

т.е. \partial F / \partial \lambda может быть рассчитана через среднее по ансамблю с распределением e^{-\beta U(\lambda)} при \lambda = \mathrm{const}.

Изменения свободной энергии при переходе из \mathcal{A} в \mathcal{B} затем может быть рассчитано через интеграл по пути \lambda:

\Delta F_{\mathcal{A}\mathcal{B}} = \int_0^1 \left\langle \frac{\partial U}{\partial \lambda} \right\rangle_\lambda \, \mathrm{d}\lambda \tag{10}

Уравнение 10 называется формулой термодинамического интегрирования [2], этот интеграл не зависит от выбора f(\lambda) и g(\lambda), однако, при выборе удачного вида (т.н. выпуклая линейная комбинация) U(\lambda) = \lambda U_{\mathcal{A}} + (1 - \lambda) U_{\mathcal{B}}, уравнение 10 упростится до \Delta F_{\mathcal{A}\mathcal{B}} = \int_0^1 \langle U_{\mathcal{B}} - U_{\mathcal{A}} \rangle \mathrm{d} \lambda.

2.

Zwanzig R.W. High‐Temperature equation of state by a perturbation method. I. Nonpolar gases // The Journal of Chemical Physics. 1954. Vol. 22, № 8. P. 1420–1426.

На практике, однако, от интегрирования по \lambda вновь переходят к суммированию по некоторым промежуточным состояниям \lambda_1, \cdots, \lambda_n: \Delta F_{\mathcal{A}\mathcal{B}} = \sum_{0 \le \lambda_k \le 1} \left\langle \frac{\partial U}{\partial \lambda} \right\rangle_{\lambda_k}, \tag{11} точки ФП для вычисления средних \langle \partial U / \lambda \rangle_{\lambda_k} берутся из распределения e^{-\beta U(\lambda = \lambda_k)}.