生物大分子模拟中的深度学习

深度学习的生物物理化学原理 - Notes Project Overview

引言

分子模拟广泛用于化学、凝聚态物理、材料、生物等学科。在生命过程的研究中，它不但能提供静态的最优结构，还能从原子尺度模拟蛋白和其他生物大分子如何通过构象变化、折叠、配体结合等行使其功能，获得热力学和动力学信息，是蛋白设计、药物设计、机理研究等流程中不可或缺的环节。

受限于计算能力，分子模拟的时长一般很难达到实验可观测的时间尺度，难以对构象变化等稀有事件直接进行模拟。为了解决这个问题，理论化学家们发展了多尺度计算方法（如获2013年诺贝尔化学奖的QM/MM）、粗粒化方法和增强抽样方法。近年这些方法都受益于深度学习从而大大拓展了适用范围，同时我们也能从中看出

笔记的这一章节架构如下：

Part I - 分子模拟的插画简介(Illustrated Molecular Dynamics)

• 体系简介 - scheme 1
• 势能面(力场)和运动方程 - scheme 2
• 构象空间探索、玻尔兹曼分布、遍历性原理 - scheme 3
• 宏观热力学、动力学性质统计 - scheme 4
• 时间尺度分离、多尺度建模 - scheme 5
• 反应坐标、自由能 - scheme 6
• 粗粒化 - scheme 7
• 增强抽样 - scheme 8

Part II - 监督学习用于粗粒化力场构建

Part III - 生成模型和强化学习用于求解统计力学

Part IV - 生成模型和强化学习用于增强抽样

---

7.1 分子模拟的插画简介(Illustrated Molecular Dynamics)

7.1.1 体系简介

我们用一个蛋白分子 1HD0 作为生物体系的样例，为了拟人化一些，给它一个简称叫 PHD（显然是精心选择过的哈哈哈）。

$N$ 个原子的蛋白分子可以用 $3N$ 个笛卡尔坐标 $\mathbf{R}=(\boldsymbol{R_1, R_2, \cdots, R_N})$，正如 PHD 们在生活中有相当多的行动可以选择。因为原子间有化学键等相互作用，实际上任意一个 $\mathbf{R}$ 就对应一个能量值 $U(\mathbf{R})$，实际上给足够时间的话可以通过量子力学原理将它算出来，称为势能面。将所有原子坐标的 $3N$ 个自由度 $\mathbf{R}$ 假想地画在一个二维平面上，纵坐标为 $U(\mathbf{R})$，就是势能面的一个示意图，分子结构的变化就相当于在势能面上移动。图上有山谷（能量局部极小值，对应一个亚稳态）、山峰（能量局部极大值）、连接山谷的鞍点（对应过渡态）等。

7.1.2 势能面(力场)和运动方程

然而实际不可能用计算机求解这么大体系的量子力学问题，得到精确的势能面。科学家们就假设体系运行遵循经典力学，通过一个有解析表达式的函数近似表达蛋白分子的能量，我们称它为力场，为简便起见同样记为 $U(\mathbf{R})$，其负梯度为原子的受力，再由牛顿第二定律得运动方程：

$$\boldsymbol{F}_{i}(t)=-\nabla_{\boldsymbol{R}_i} U(\mathbf{R})$$ $$\frac{\partial^2 \boldsymbol{R}_i}{\partial t^2}=-\frac{\nabla_{\boldsymbol{R}_i} U(\mathbf{R})}{m_i}$$

7.1.3 构象空间探索、玻尔兹曼分布、遍历性原理

分子模拟可以在不同控制条件(系综)下进行。如控制温度和体积不变时，化学平衡的原理告诉我们，经过无限长时间，探索过每处 $\mathbf{R}$，概率分布达到平衡后，两点间的概率之比为平衡常数 $K=\exp[-\Delta U/k_\mathrm{B}T]$，即在某一点 $\mathbf{R}$ 处停留的概率正比于 $\exp[-U(\mathbf{R})/k_\mathrm{B}T]$。这样的概率分布在物理化学中称为玻尔兹曼分布：

$$p(\mathbf{R})=\frac{e^{-\beta U(\mathbf{R})}}{Z},\quad \quad Z=\int e^{-\beta U(\mathbf{R})}\mathrm{d} \mathbf{R},\quad \quad \beta=\frac{1}{k_\mathrm{B}T}$$

由于只能通过平衡知道任意两点间概率的比值，求解统计力学的核心任务就是对这些比值进行归一化，表达为概率 $p(\mathbf{R})$。归一化常数 $Z$ 称为配分函数(partition function)，需要对所有可能的 $\mathbf{R}$ 积分。对对数配分函数 $\ln Z$ 求导我们就能获得一个体系所有的平衡态热力学性质。

“我太难了。”

然而这样的高维积分实际上是一个无法计算(intractable)的任务。分子模拟的策略是改积分计算为采样计算(对样本取平均)：有限时间的模拟，实际上是在探索构象空间，对玻尔兹曼分布进行采样。然而在有限时间的模拟中，不可能遍历整个构象空间，分子会陷在局部最小值附近的势阱中，短时间难以越过高势垒探索新的构象区域，采样就无法覆盖到一些重要变化，称为遍历性破缺(broken ergodicity)。

7.1.4 宏观热力学、动力学性质统计

始终应该明确一点，观测到的宏观结果都是许多微观结构、微观反应路径统计平均的结果。

7.1.5 时间尺度分离

PHD 们平时玩游戏、刷朋友圈是快自由度，决定他们未来方向的文献阅读、JC、idea、答辩、论文则是慢自由度。

正如刚刚 7.1.3 图中所示，生物大分子在势能面上移动时，可以发现它同时具有几种运动模式：在一些方向上进行快速的往返运动，如键的振动、部分溶剂分子的运动；在另一些方向上运动则非常缓慢，预示着一些反应、配体结合、构象变化等稀有事件的发生。这被称为时间尺度分离(separation of timescales)。生物大分子典型变化的时间尺度如图3。

7.1.6 反应坐标&自由能

分子模拟需要较好地表现运动最快的自由度，因此需要设定足够小的时间步长，导致模拟的总时长被快自由度限制。一个自然的想法就是，既然我们关注的稀有事件发生时，快自由度中的事件已发生了相当多次，可否将它们平均掉，只作为其他运动时的噪声出现？布朗运动和朗之万方程就是极好的例子：将溶剂分子的快速运动平均为随机力。生物大分子的各种动力学过程，同样可以选取一些核心的描述整个变化过程的量，称为集团变量(collective variables, CV) 或反应坐标(reaction coordinates, RC)：

$$\mathbf{s}(\mathbf{R})=\left(s_{1}(\mathbf{R}), s_{2}(\mathbf{R}), \ldots, s_{d}(\mathbf{R})\right),\quad\quad d \ll 3N$$

一个简单例子是选为 $\boldsymbol{R}_i$ 的线性组合。可以看出实际上是一个对高维运动降维的过程。对快自由度的能量平均后加入熵（多个微观状态 $\mathbf{R}$ 对应同一个反应坐标 $\mathbf{s}$，就有了微观状态数）就得到了自由能 $F(\mathbf{s})$：

$$p(\mathbf{s})=\int_{\mathbf{s}=\mathbf{s}(\mathbf{R})} p(\mathbf{R})\mathrm{d} \mathbf{R},\quad \quad p(\mathbf{R})=\frac{e^{-\beta U(\mathbf{R})}}{Z},\quad \quad$$ $$Z=\int e^{-\beta U(\mathbf{R})}\mathrm{d} \mathbf{R},\quad \quad \beta=\frac{1}{k_\mathrm{B}T}$$ $$\begin{aligned} F(\mathbf{s})&=-\frac{1}{\beta} \ln [p(\mathbf{s})\cdot Z]=-k_\mathrm{B}T \ln [p(\mathbf{s})\cdot Z]\\ &=-\frac{1}{\beta} \ln \int_{\mathbf{s}=\mathbf{s}(\mathbf{R})} e^{-\beta U(\mathbf{R})}\mathrm{d} \mathbf{R} \\ &=-k_\mathrm{B}T \int_{\mathbf{s}=\mathbf{s}(\mathbf{R})} p(\mathbf{R})\ln e^{-\beta U(\mathbf{R})}\mathrm{d} \mathbf{R}\\ &\quad\quad-T\cdot \left(-k_\mathrm{B} \int_{\mathbf{s}=\mathbf{s}(\mathbf{R})} p(\mathbf{R})\ln p(\mathbf{R}) \mathrm{d} \mathbf{R}\right)\\ &=\mathbb{E}_{\mathbf{R},\mathbf{s}(\mathbf{R})=\mathbf{s}} [U(\mathbf{R})]-T\cdot -k_\mathrm{B} \mathbb{E}_{\mathbf{R},\mathbf{s}(\mathbf{R})=\mathbf{s}} [\ln p(\mathbf{R})]\\ &=\mathbb{E}_{\mathbf{R},\mathbf{s}(\mathbf{R})=\mathbf{s}} [U(\mathbf{R})]-T\cdot S(\mathbf{s}) \end{aligned}$$

自由能面(Free Energy Surface, FES)比势能面光滑许多。

7.1.7 粗粒化

两耳不闻窗外事，一心只读圣贤书。

“我变秃了，也变强了。”

粗粒化可认为是实现上述过程的一个方法：将数个原子视作一个粗粒，从而完成了对局部化学键振动这一快自由度的平均。由此实现降维，自由度大大减少，可采取的时间步长增大。粗粒化方法大致分为 bottom-up coarse-graining(拟合全原子模拟的受力或热力学参数) 和 top-down coarse-graning(拟合实验数据)。

7.1.8 增强抽样

月下柳梢星沉湾，伫马听浪入梦难。学海无涯苦作船。
春蚕吐丝千斤缆，蜡炬点灯万丈光。征帆不落桨声欢。

“老师真是世界上最好的导师，我一个人可能100年都发不出的文章，在Ta指点下1个月就完成了。”

与粗粒化不同，虽然一大类增强抽样方法也依赖于反应坐标，但它仍然采用全原子体系，而在反应坐标的维度上修改势能面使稀有事件更易发生，计算宏观性质的统计平均时只需根据重要性采样(importance sampling) 原理乘上两个势能面导致的概率分布之差。退火(annealing, tempering) 类方法可视作将 $U(\mathbf{R})$ 作为反应坐标。

$$\mathbb{E}_{\mathbf{s}\sim p(\mathbf{s})}[f(\mathbf{s})]=\int f(\mathbf{s}) p(\mathbf{s}) \mathrm{d}\mathbf{s} = \int \frac{f(\mathbf{s}) p(\mathbf{s})}{q(\mathbf{s})}\cdot q(\mathbf{s})\mathrm{d}\mathbf{s} = \mathbb{E}_{\mathbf{s}\sim q(\mathbf{s})} \left[\frac{f(\mathbf{s}) p(\mathbf{s})}{q(\mathbf{s})} \right]$$ $$p(\mathbf{s})=\frac{e^{-\beta F(\mathbf{s})}}{Z_p},\quad \quad q(\mathbf{s})=\frac{e^{-\beta [F(\mathbf{s})+V(\mathbf{s})]}}{Z_q}$$

增强抽样的一个代表性方法是 Metadynamics，它通过不断在当前位置加偏置势(bias potential) 跃出势阱。

7.2 监督学习用于粗粒化力场构建

Highlights

• 特征选取为粗粒化模型的内坐标，保证三维旋转、平移对称性
• 三种模型架构考虑的不同相互作用方式

7.3 生成模型和强化学习用于求解统计力学(Boltzmann Generator)

Highlights

• $\mathbf{x}$ 和 $\mathbf{z}$ 互相作为对方的“反应坐标”
• 神经网络（具体为深度生成模型中的流模型(flow model)）实现相互作用体系到无相互作用体系(如理想气体)的可逆变换
• 平行采样，计算自由能和自由能变
• 目标 AlphaFold Zero，只用物理知识预测大分子的结构系综和动力学
• 统计力学启发了 Boltzmann Machine 一类基于能量的生成模型，演化为如今的 Boltzmann Generator，反哺统计力学

7.4 生成模型和强化学习用于增强抽样

Discussion Highlights

• 核心是反应坐标的选取，这不是简单的对高维数据的降维，而是对高维运动（有了时间影响）。
• 神经网络虽然是通用函数拟合器，但需要人为设计架构用以保证泛化能力，不然就是人工智障。保证泛化能力的关键是抓住体系/数据的共性和结构，称为诱导偏置(inductive bias)，如对称性：
  • 三维平移、旋转对称性（通过内坐标等方法实现）
  • 同类原子交换对称性（极难满足）
  • Hamilton 力学的辛结构（较难满足）
• 统计力学（生物物理化学）与深度学习的 Cross-fertilization

---

• Previous
  引言 | 复杂无序系统的崎岖势能面和多尺度现象