结构生物化学/分子建模/分子动力学

理论概述

一个简单的分子动力学笼子的例子，它展示了一个小的氨基酸链（色氨酸笼）在隐式溶剂中的折叠。

经典分子动力学 (MD) 是计算化学的一个分支，它专注于根据牛顿运动定律模拟粒子的运动。 原子被近似为“弹簧上的球”，这使得可以应用以下定律

$f=ma$

$v=v_{0}+at$

$x=vt+{\frac {1}{2}}at^{2}$

其中 $x$ 是位置，$v$ 是速度，$a$ 是加速度，以及 $t$ 是时间。

经典 MD 近似将键和角势视为谐波，并将二面角视为周期性的。这使得可以用类似于 AMBER 分子动力学程序使用的以下方程来表达系统的能量：

$E_{\rm {total}}=\sum _{\rm {bonds}}K_{r}(r-r_{eq})^{2}+\sum _{\rm {angles}}K_{\theta }(\theta -\theta _{eq})^{2}+\sum _{\rm {dihedrals}}{V_{n} \over 2}[1+{\rm {cos}}(n\phi -\gamma )]+\sum _{i<j}\left[{A_{ij} \over R_{ij}^{12}}-{B_{ij} \over R_{ij}^{6}}+{q_{i}q_{j} \over \epsilon R_{ij}}\right]+\sum _{\rm {H-bonds}}\left[{C_{ij} \over R_{ij}^{12}}-{D_{ij} \over R_{ij}^{10}}\right]$

模拟在一系列“时间步长”中进行。在每一步中，根据能量方程的导数计算每个原子的作用力。原子会移动一定的距离，这个距离与作用力的方向和大小以及时间步长的大小成正比。这个过程会重复很多时间步长，迭代地计算作用力和求解运动方程（基于作用力产生的加速度），并生成模拟复杂行为的结果，从药物化合物的对接到简单蛋白质的折叠。

这是一个简化的分子动力学模拟步骤。生产级 MD 代码使用更复杂版本的算法，包括温度和压力控制。

限制

量子效应

尽管已知原子是根据薛定谔方程而不是牛顿定律来运作的，但经典方法仍能用较高的精度进行模拟。可以通过检查参与粒子的德布罗意波长来了解这一点

$\Lambda ={\sqrt {\frac {2\pi \hbar ^{2}}{Mk_{B}T}}}$

其中 $M$ 是原子质量，$T$ 是温度。当 $\Lambda \ll$ 粒子之间的平均最近邻间距时，经典近似是合理的。因此，对于较重的系统，如有机分子，经典近似是合理的，但对于轻型系统，如 H₂、He 或 Ne，该方法的精度会大大降低。

由于原子被近似为离散粒子，因此无法用经典 MD 正确模拟较小粒子的相互作用。 化学键的形成和断裂、非共价反应中间体以及质子或电子隧穿都无法用这种方法模拟。在非常低的温度下进行的模拟也存在问题，因为在接近绝对零度时，经典力学越来越不受青睐，而量子解释则越来越受青睐。量子分子动力学方法可以用于模拟，这些模拟可以更电子地逼真，但计算成本要高得多。

系统大小和时间尺度

可以模拟的原子数量和时间量与可用的计算能力直接相关。虽然之前已经模拟过一个小型病毒的全部结构，但大多数精细的、全原子模拟都是针对蛋白质的规模，而不是细胞器或复杂的组装体。这是因为观察结果所需的模拟时间与系统大小直接相关。

表格 - 分子过程的时间尺度

秒	名称	过程
10−1		表面重组
10−3	毫秒	蛋白质折叠
10−8	10 纳秒	核磁共振
10−9	至 10+3	布朗动力学
10−9	至 100	格点蒙特卡罗
10−12	皮秒	快速振动
10−15	飞秒	电子激发

因此，系统越大，模拟运行的时间就需要越长才能观察到相关的行为。如今，专门为分子动力学模拟设计的计算机硬件 可以用一天的时间模拟大约 17 微秒的 23,000 原子系统，时间步长为 0.5 飞秒。在同一系统上，通用并行超级计算机可以获得大约几百纳秒的模拟时间。

恒温器

使用基于牛顿运动方程的方法，在“微正则系综”中获得了模拟，这是一个恒定能量系统 (nVE)。然而，这并不现实，为了与实验结果进行准确比较，模拟最常在正则系综中进行，该系综具有恒定温度 (nVT)。这个温度必须用算法来维持，这种算法被称为恒温器。常用的恒温器有很多，它们在可靠性、精度和计算成本方面各有优劣。

典型的模拟包括一个 MD 平衡阶段，在这个阶段，结构在 0K 恒定温度下模拟一定时间，一个加热阶段，在这个阶段，系统逐渐加热到所需的温度（通常为 300K），以及一个在该恒定温度下的最终生产模拟。最初的平衡不是为了使结构达到化学平衡，而是为了达到能量最小值（通常接近化学平衡状态）。如果省略平衡阶段，而在高能量区域运行模拟，模拟将继续对高能量相空间进行采样，从而产生与实验数据不一致的结果。

朗之万恒温器

这个恒温器遵循朗之万运动方程而不是牛顿运动方程，其中一个与速度成正比的摩擦力被添加到保守力中，它会调整粒子的动能，从而使温度正确。

然后系统由以下方程描述：

$ma=-v\xi +f(r)+f'$

其中 $m$、$v$ 和 $a$ 分别是粒子的质量、速度和加速度。 $\xi$ 是摩擦常数，而 $f'$ 是一个随机力，它会向粒子添加动能。它是从高斯分布中采样的，其中方差是设定温度和时间步长的函数。这平衡了摩擦力，以使温度保持在设定值。

安德森恒温器

这个恒温器将系统耦合到一个设定为所需温度的热浴。热浴由随机选择的系统粒子上的随机粒子碰撞来表示。

据已证实，将随机碰撞添加到牛顿 MD 系统会使模拟成为一个不可约且非周期的马尔可夫链。这意味着使用安德森算法生成的分布不是正则分布。

这个恒温器会随机影响粒子的速度，这会导致模拟系统中出现非物理动力学。

高斯恒温器

这个恒温器将系统与设定为所需温度的热浴紧密耦合。这里，每个粒子的动量变化表示为：

${\frac {dp_{i}}{dt}}=\displaystyle \sum _{j}^{N}(F_{ij}|q_{i}-q_{j}|)-p_{i}\alpha$

这种重新缩放会很快破坏动态相关性，因此不会保持正则系综。

贝伦森恒温器

这个恒温器首先于 1984 年提出，它涉及以控制温度的方式重新缩放每个粒子的速度。

与高斯恒温器类似，系统与一个设定为所需温度的热浴耦合，尽管这种算法使用的是弱耦合而不是强耦合，其中系统温度的校正使得温度偏差以某个时间常数 $\mathrm {T}$ 指数衰减，该时间常数与时间步长成正比。

${\frac {dT}{dt}}={\frac {T_{0}-T}{\mathrm {T} }}$

Nose-Hoover 热浴

该热浴在牛顿运动方程中添加了一个额外的项，$\eta$，产生了以下二阶微分

${\frac {dp_{n}}{dt}}=\displaystyle \sum _{i}^{N}{\frac {|p_{i}|^{2}}{2m_{i}}}-{\frac {3}{2}}nk_{B}T$

其中 $m$ 是质量，$T$ 是目标温度，$k_{B}$ 是玻尔兹曼常数。热浴被赋予了一定的质量，$Q$，对于 Nose 建议的 $6nk_{B}T$ 值，产生了上述每个粒子的动量变化方程。

使用该热浴，模拟系统的能量会发生波动，但系统和热浴的总能量是守恒的。如果系统是遍历的（所有可访问的微观状态的可能性相等），那么该系统将给出正则系综。

力场

分子动力学模拟需要一个“力场”，它列出了能量方程中每种相互作用类型涉及的参数。例如，AMBER 力场需要每个可能的两种原子类型的组合的键平衡长度、键力常数，每个可能的三个原子类型的组合的角度平衡值和角度力，每个可能的四种原子类型的组合的二面角平衡相位、二面角力常数和二面角周期性。额外的参数还描述了每个原子的非键合（静电）力。

这些参数的开发和分析，统称为力场，是分子动力学研究的重点。

力场开发

经典分子动力学的参数理想情况下会导致行为与原子水平上实际分子的行为非常接近，而实际分子在现实中受量子现象控制。经典参数的开发使得这种情况成为可能——一旦能量和力方程（例如 AMBER 方程）被用来建立系统的范式，参数就会被推导出来，使得使用它们计算的值与使用量子力学方法或实验计算的值相匹配。由于缺乏用于开发力场所需的各种分子和分子构象的实验数据，因此使用量子计算得出的数据来代替。

力场是为特定类型的分子开发的。例如，AMBER 为蛋白质提供了 Amber ff99SB 力场，为碳水化合物提供了 Glycam 力场，以及新开发的脂质力场 GAFFLipid。此外，还开发了“通用 AMBER 力场”或 GAFF，用于对不适合广泛力场类别的小配体和其他分子进行参数化。

力场的准确性对于模拟的准确性至关重要。诸如二面角力常数之类的简单项，可以对整体模拟产生极大的影响。事实上，蛋白质主链 $\phi$ 和 $\psi$ 角的不准确二面角参数会导致模拟过度倾向于α螺旋形成，这一事实影响了许多结果。

从历史上看，力场开发一直受到缺乏可用于拟合的能量和力数据的困扰。力场参数是通过拟合参数来生成的，使得使用这些参数预测的结构的能量或每个原子的力与从头算量子计算生成的能量或力匹配。这些计算非常昂贵，因此通常从一组参数推断到一个通用参数，以避免重复进行计算。例如，使用为四丙氨酸蛋白质主链生成的参数来描述任何蛋白质的主链。

计算能力的不断提高、更有效的算法以及对过去方法的了解为未来的参数开发带来了希望。例如，可以分别计算各种氨基酸的侧链参数，然后将它们添加到一个通用参数集中，以便分别描述每个氨基酸。

力场准确性

经典力场是通过拟合经典参数而开发的，这些参数使其产生的结果与实验和量子数据一致。它们按生物分子类型进行分类——蛋白质、碳水化合物和脂类力场已针对所有主要的 MD 程序单独开发。然而，这种分类仍然相当广泛，并且可能忽略系统特定的相互作用。

使用量子和经典计算方法比较了 N-甲基乙酰胺的两个主链二面角的势能面。当二面角接近平衡角时，经典近似将产生准确的结果，但在平衡角较大距离处，两种方法预测的能量偏差会显著增加。