Теория возмущений свободной энергии
Пусть имеем термодинамическую систему, состоящую из молекулы лиганда \(\rm L\) и белка \(\rm E\): в состоянии \(\cl A\) лиганд \(\rm L\) и \(\rm E\) и белок не связаны, в состоянии \(\cl B\) — имеет место комплекс \(\rm{ES}\), соответственно введём потенциальные энергии \(U_{\cl A}(\bf r)\) и \(U_{\cl B}(\bf r)\) для несвязанного и связанного состояний. Свободная энергия Гельмгольца перехода в связанное состояние [1, с. 138] \[\Delta F_{\cl A \cl B} = F_{\cl B} - F_{\cl A} \tag{8.1}\] может быть выражена через суммы по состояниям: \[ \Delta F_{\cl A \cl B} = F_{\cl B} - F_{\cl A} = - k T \ln \dfrac{Q_{\cl B}}{Q_{\cl A}}, \tag{8.2}\] где, например, \[ \begin{aligned} Q_{\cl A} &:= \dfrac{1}{N! h^{3N}} \int_\t{ФП} e^{- \beta H(\bf p, \bf q)} \dd \bf p \dd \bf q = \\ &= \dfrac{1}{N! h^{3N}} \int_\t{ФП} \exp \left\{ - \beta \left[\sum_{i=1}^N \dfrac{\bf p_i^2}{2m_i}\right] + U_{\cl A} (\bf q_1, \cdots, \bf q_N)\right\} \dd \bf p \dd \bf q, \end{aligned} \tag{8.3}\] где \(H\) — гамильтониан системы \(\cl A\), \(N\) — число частиц в ней, \(\beta = 1 / kT\), \(\dd \bf p \dd \bf q\) — точка ФП; интегрирование ведётся по всему ФП, что эквивалентно интегрированию по области, соответствующей состоянию \(\cl A\), т.е. там, где вероятность пронаблюдать микросостояние состояния \(\cl A\) не нулевая.
Сократив в Ур. 8.2 числитель и знаменатель на импульсозависимую часть \(\exp \{ \sum_i \frac{\bf p_i^2}{2 m_i} \}\), перейдём к конфигурационным интегралам: \[\Delta F_{\cl A \cl B} = F_{\cl B} - F_{\cl A} = -kT \ln \dfrac{Z_{\cl B}}{Z_{\cl A}}, \tag{8.4}\] где \[ Z_{\cl A} := \int e^{ -\beta U_{\cl A} (\bf q_1, \cdots, \bf q_2) } \dd \bf q. \tag{8.5}\]
Оценить \(\Delta F\) непосредственно по уравнениям Ур. 8.2 и Ур. 8.4 сложно, поскольку методы МД и Монте-Карло не позволяют непосредственно рассчитывать статистические суммы и интегралы, однако позволяют получить хорошие средние оценки ТД функций.
Уравнение Ур. 8.4 можно переформулировать через средние ТД величины. Конфигурационный интеграл представим в виде интеграла от двух множителей-экспонент: \[Z_{\cl B} = \int e^{-\beta U_{\cl B}} \dd \bf q = \int e^{-\beta U_{\cl B}} e^{\beta U_{\cl A}} e^{-\beta U_{\cl A} } \dd \bf q = \int e^{- \beta (U_{\cl B}-U_{\cl A})} \ub{ e^{-\beta U_{\cl A}} }_{\rho_{\cl A}} \dd \bf q, \tag{8.6}\] где \(\rho_{\cl A}\) — функция распределения, соответствующая состоянию \(\cl A.\) Среднее значение произвольной функции \(y\) по ансамблю \(\cl A\) можно вычислить через функцию распределения \(\rho_{\cl A}\): \[ \langle y \rangle_{\cl A} := \dfrac 1 {Z_{\cl A}} \int_\t{ФП} \rho_{\cl A} (\bf p, \bf q) \cdot y(\bf p, \bf q) \dd \bf p \dd \bf q, \tag{8.7}\] учитывая Ур. 8.6 и Ур. 8.7, дробь в уравнении Ур. 8.4 примет вид: \[ \dfrac{Z_{\cl B}}{Z_{\cl A}} = \dfrac 1 {Z_{\cl A}} \int e^{-\beta U_{\cl A}} e^{-\beta (U_{\cl B} - U_{\cl A})} = \langle e^{-\beta (U_{\cl B} - U_{\cl A})} \rangle_{\cl A}, \] т.е. среднее берётся из распределения, соответствующего состоянию \(\cl A\). Наконец, \[ \Delta F_{\cl A \cl B} = -kT \ln \langle e^{-\beta (U_{\cl B} - U_{\cl A})} \rangle_{\cl A}. \tag{8.8}\]
Уравнение Ур. 8.8 известно как возмущение свободной энергии. Опишем, как следует использовать это уравнение. Состоянию \(\cl A\) отвечает функция распределения \(\rho_{\cl A}\) и потенциал \(U_{\cl A}\), в ходе численного эксперимента мы получаем точки ФП, соответствующие этому распределению, и затем используем их для моделирования состояния \(\cl B,\) заменив потенциал на \(U_{\cl B}\). Однако, распределению \(\cl B\) соответствует другая функция распределения \(\rho_{\cl B}\), поэтому точки ФП следует перевзвесить1, поделив на исходное \(e^{-\beta U_{\cl A}}\) и умножив на \(e^{-\beta U_{\cl B}}\).
1 от англ. reweighting
Таким образом, при расчёте \(\Delta F\) по Ур. 8.8 мы собираем статистику по ФП, соответствующую начальному состоянию, и затем используем эти наблюдения для оценки \(\Delta F\), как если бы они были взяты из распределения состояния \(\cl B\). Ограничением этого подхода является то, что микросостояния из распределения \(\cl A\) могут не быть микросостояниями с высокой вероятностью в распределении \(\cl B\), в этом случае разность \(U_{\cl B} - U_{\cl A}\) будет большой и экспоненциальный фактор \(e^{-\beta (U_{\cl B} - U_{\cl A})} \to 0\), соответствующие состояния будут иметь малый вес в среднем по ансамблю \(\langle \cdots \rangle_{\cl A}\) и, как следствие, сходимость будет медленной. Другими словами, требуется, чтобы состояние \(\cl B\) было небольшим возмущением состояния \(\cl A\) и их области в ФП существенно перекрывались.
В постановке задачи, состоянию \(\cl A\) соответствует несвязанные белок E и лиганд L, состоянию \(\cl B\) — фермент-субстратный комплекс EL, на масштабах системы «лиганд—белок—гидратная оболочка» переход из несвязанного в связанное состояние в принципе можно считать небольшим возмущением.
Термодинамическое интегрирование
Выбор начального и конечного состояний не ограничивается связыванием в комплекс, состояния могут соответствовать любым областям в ФП. В случае, когда \(\cl A\) и \(\cl B\) достаточно далеки, вводят промежуточные состояния, переходы между которыми можно считаться малым возмущением. Уравнение Ур. 8.8 в этом обобщении имеет вид: \[ \Delta F_{\cl A \cl 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_{\cl A} + g(\lambda) U_{\cl B}, \\ f(\lambda) \t{ монотонно убывает от $1$ до $0$}, \\ g(\lambda) \t{ монотонно возрастает от $0$ до $1$}. \end{cases} \]
Из выражения свободной энергии, соответствующей состоянию \(U(\lambda)\): \[ F(\lambda) = - k T \ln Q(\lambda), \] можно получить соотношение для производной свободной энергии \(F\) по параметру пути \(\lambda\): \[ \begin{aligned} \pdv{F}{\lambda} &= -\dfrac{kT}{Q} \pdv{Q}{\lambda} = - \dfrac{kT}{Z} \pdv{Z}{\lambda} = - \dfrac{kT}{Z} \pdv{}{\lambda} \int_\t{ФП} e^{- \beta U(\bf q; \lambda)} \dd \bf q = \\ &= - \dfrac{kT}{Z} \int_\t{ФП} \left(-\beta \pdv{U(\bf q; \lambda)}{\lambda}\right) e^{- \beta U(\bf q; \lambda)} \dd \bf q = \left\langle \pdv{U(\lambda)}{\lambda} \right\rangle_\lambda, \end{aligned} \tag{9.1}\] т.е. \(\partial F / \partial \lambda\) может быть рассчитана через среднее по ансамблю с распределением \(e^{-\beta U(\lambda)}\) (\(\lambda = \const\)).
Изменения свободной энергии при переходе из \(\cl A\) в \(\cl B\) затем может быть рассчитано через интеграл по пути \(\lambda\): \[ \Delta F_{\cl A \cl B} = \int_0^1 \left\langle \pdv{U}{\lambda} \right\rangle_\lambda \dd \lambda. \tag{9.2}\]
Уравнение Ур. 9.2 называется формулой термодинамического интегрирования [2], этот интеграл не зависит от выбора \(f(\lambda)\) и \(g(\lambda)\), однако, при выборе удачного вида (т.н. выпуклая линейная комбинация) \[ U(\lambda) = \lambda U_{\cl A} + (1 - \lambda) U_{\cl B}, \] уравнение Ур. 9.2 упростится до \[ \Delta F_{\cl A \cl B} = \int_0^1 \langle U_{\cl B} - U_{\cl A} \rangle \dd \lambda. \]
На практике, однако, от интегрирования по \(\lambda\) вновь переходят к суммированию по некоторым промежуточным состояниям \(\lambda_1, \cdots, \lambda_n\): \[ \Delta F_{\cl A \cl B} = \sum_{0 \le \lambda_k \le 1} \left\langle \pdv{U}{\lambda} \right\rangle_{\lambda_k}, \tag{9.3}\] точки ФП для вычисления средних \(\langle \partial U / \lambda \rangle_{\lambda_k}\) берутся из распределения \(e^{-\beta U(\lambda = \lambda_k)}\).