阅读笔记-A Graph to Graphs Framework for Retrosynthesis Prediction

一、概要

逆向合成致力于寻找一组反应物以合成靶分子，对合成计划和药物开发至关重要。
本文的两个关键点：

• 反应中心识别模块，其将靶分子上的合成子分离，并将一对多图的翻译问题简化为多个一对一的翻译过程；
• 一个将合成子转换为最终的反应物的图结构的变分图翻译模型。
  这篇文章利用图的生成来进行逆向合成预测，其中必须为这个特定的任务设计新的组件。例如，反应物可以在不同的反应条件下进行反应。为了解决这种多模态问题，在这片文章中，作者还引入了低维的潜在变量来调节序列图的转换过程，以捕捉不同的输出分布。

二、Graph to Graphs 架构

这篇文章中，作者将逆合成任务表述为一对多的图到图转换的问题。具体而言，首先采用图神经网络来估计产物图中所有原子对的反应性得分，然后将具有高于我们所设定的阈值并且反应性得分最高的源自对作为反应中心。然后通过断开反应中心键的连接，将产物图拆分为合成子。最后基于所获得的合成子，通过一系列图变换生成反应物，其中使用一个潜在向量赖不活变换的不确定性从而产生不同的预测。整个框架如下图所示：

图中红色的反应的中心点由$G2Gs$标记出来，先将合成的图通过断开反应中心点来变成合成子。接下里根据得到的合成子，反应物通过一系列的图变换得到。生成的分子支架以蓝色边框为界。

2.1 预备知识

（1）Notation：

• $\mathbf{G}=(\mathbf{A},\mathbf{X})$；$\mathbf{A}$是邻接矩阵，$\mathbf{X}$是节点的特征矩阵；$n$为节点的数量；$b$为边的类型；则有 A ∈ { 0 , 1 } n × n × b \mathbf A\in\left\{0,1\right\}^{n\times n\times b} A∈{0,1}n×n×b， X ∈ { 0 , 1 } n × d \mathbf X\in\left\{0,1\right\}^{n\times d} X∈{0,1}n×d；则${\mathbf{A}}_{ijk}=1$表示第$i$个节点和第$j$个节点通过标签为$k$的边相连。

（2）Retrosynthesis Prediction：一个化学反应可以看作是两个分子集合 ( { G i } i = 1 N 1 , { G j } j = 1 N 2 ) ({\left\{G_{i}\right\}}_{i=1}^{N_{1}},{\left\{G_{j}\right\}}_{j=1}^{N_{2}}) ({Gi​}i=1N1​​,{Gj​}j=1N2​​)，并且$G_i$是反应物的图结构，$G_j$是生成物的图结构。本文中我们将$N_2$设定为1，即 ( { G i } i = 1 N 1 , G p ) ({\left\{G_{i}\right\}}_{i=1}^{N_{1}},G_{p}) ({Gi​}i=1N1​​,Gp​)
（3）Reaction Center and Synthon
反应中心是用来表示满足以下两个条件的原子对$(i,j)$：      \;\; 1）生成物图结构的第$i$个节点和第$j$个节点之间有边相连；2）生成物的图结构中对应的这两个原子没有边相连。另外，合成子是通过破坏反应中心的键从产物中提取的出来的子图，这些合成子稍后可以转化为反应物，并且合成子并不是有效的分子。
（4）Molecular Graph Representation Learning
作者使用了$\mathbf{L}$，层的$\mathbf{R}-\mathbf{G}\mathbf{C}\mathbf{N}$来获取节点表示矩阵${\mathbf{H}}^L$。之后用一个$Readout(\cdot )$函数（将${\mathbf{H}}^L$作为输入来获得图级别的embedding$h_G$）。

2.2 Reaction Center Identification

给定一个化学反应 ( { G i } i = 1 N 1 , G p ) ({\left\{G_{i}\right\}}_{i=1}^{N_{1}},G_{p}) ({Gi​}i=1N1​​,Gp​)，我们先针对生成物$G_p$中的每对原子定义一个二进制标记矩阵 Y ∈ { 0 , 1 } n × n \mathbf Y\in\left\{0,1\right\}^{n\times n} Y∈{0,1}n×n，并且$Y_{ij}=1$处表明生成物的图结构$G_p$第$i$个节点和第$j$个节点之间的化学键是反应中心。并且反应集（生成物和反应物）中的每个原子都是被单独标记的，这种性质使得我们能够简单地比较反应产物$G_p$，和反应物$G_i$中的每一对原子来识别反应中心，从而形成二进制标记矩阵$\mathbf{Y}$。有了这样的标记后$G2Gs$中的反应中心预测就可以转化为一个如下所示的二进制链接预测任务。
我们先用一个$L$层的$\mathbf{R}-\mathbf{G}\mathbf{C}\mathbf{N}$来计算节点的embeddings矩阵${\mathbf{H}}^L$以及整张生成物的图结构$G_p$的embedding$h_{G_P}$，即：

但是由于$\mathbf{L}$层的$R-GCN$只能聚合一个节点$\mathbf{L}$跳内的邻居信息，但反应中心的反应活性也可能受到远处原子的影响，为了确保边的$embedding$能够包含到远处原子的作用，作者将图级别的$embedding$也纳入了$edgeembedding$中，有：

• $||$：代表着串联；
• $H_i^L\in {\mathbb{R}}^k$代表着$H^L$的第$i$行；
• A i j = A [ i , j , : ] ∈ { 0 , 1 } b A_{ij}=A_{[i,j,:]}\in\left\{0,1\right\}^{b} Aij​=A[i,j,:]​∈{0,1}b代表着$G_p$中原子对$(i,j)$之间边的类别。

原子对之间最终的反应分数的计算方式为：${\mathbf{s}}_{\mathbf{i}\mathbf{j}}=\sigma ({\mathbf{m}}_{\mathbf{r}}({\mathbf{e}}_{\mathbf{i}\mathbf{j}}))$。$m_r$代表着一个将$e_{ij}$映射为一个标量的前向传递网络，$\sigma (\cdot )$代表着$Sigmoid$函数。在逆合成分析中已知某个反应类型的情况下，我们也可以通过将其嵌入到前馈网络的输入中来包含这些信息。接下来可以通过最大化二值标记矩阵Y的交叉熵来优化反应中心识别模块，如下所示：

实际上，在所有的原子对中，通常有几个反应中心。因此，前馈网络的输出对数值通常非常小。通过加权合成方法，可以缓解由于类分布不平衡而造成的问题。在推理过程中先计算所有原子对之间的反应性分数，然后选择超过阈值的原子对中分数最高的那一个作为反应中心。另外一种方法是选择得分最高的$k$个原子对作为反应中心。这种方法将会产生更多样化的反应路线但相应的需要花费更多的时间。接下来断掉反应中心之间的化学键，将反应产物$G_p$之间的每个子图看做一个合成单体。需要注意的是如果未识别出反应中心，整个生成物$G_p$将被看作为合成单体，然后制定一对一的图转换过程，将每个合成子转换成最终反应物的图结构。

2.3 Reactants Generation via Variational Graph Translation

给定一个化学反应 ( { G i } i = 1 N 1 , G p ) ({\left\{G_{i}\right\}}_{i=1}^{N_{1}},G_{p}) ({Gi​}i=1N1​​,Gp​)，我们把从$G_p$中提取出来的合成子表示为 { S i } i = 1 N 1 \left\{S_{i}\right\}^{N_{1}}_{i=1} {Si​}i=1N1​​。我们将其简化为$(S,G)$，接下来我们的目标就是学习到一个条件生成模型$p(G|S)$。值得注意的是，一个内在的问题是相同的合成子在不同的反应条件下可以转化为不同的反应物，这被称为多峰问题（感觉有点像cycle-gan中的问题，即一张图片只能翻译到另外一张图片，但其实一张图片可以翻译成很多内容风格都想相似但却有不同的图片）。为了减轻这种多模态问题并激发模块获得对反应物不同分布的建模，作者引入了一个地位的潜在变量$z$来捕获图转换过程中的不确定性。

2.3.1 生成模型

图$G$的生成是基于$S$和$latentvectorz$的条件分布。具体来说，我们首先自回归地生成一系列图变换动作$({\alpha }_1,...,{\alpha }_T)$，然后将它们应用到初始合成子图$S$上。这里${\alpha }_t$是用作图的修改的单步图变换（即动作）。$\mathcal{T}$来代表所有图变换动作的集合即：$a_1,...,a_T$，$\mathcal{T}$为合成子$S$转换为目标反应物$G$的所有轨迹的集合$(a_1,...,a_T)$。然后，我们将模型$p(G|z,S)$的建模转化为对图变换动作序列的联合分布$p(t|z,S)$建模。$p(G|z,S)$与$p(t|z,S)$之间的联系如下图所示：使用这种生成模型，可以通过从分布中采样行动序列，将合成子转化为反应物。
生成过程如下 ：

• 令$S^i$来表示对$S$施加$a_{1:i}$的结果。并且有$S^0=S$以及$p(S^i|S^{i-1},z)=p(a_i|S^{i-1},z)$。
• 作者用$MDP(MarkovDecisionProces)$来建模这个过程即$p(S^i|S^{i-1},z)=p(S^i|S^{i-1},...,S^0,z)$。用 马尔科夫决策过程意味着每一次行动仅取决于目前为止已修改的图结构。因此图变换模型可以被以下公式表示：$\mathbf{p}(\mathbf{t}|\mathbf{z},\mathbf{S})=\mathbf{p}({\mathbf{a}}_{1:\mathbf{T}}|\mathbf{z},\mathbf{S})={\prod }_{\mathbf{i}=1}^{\mathbf{T}}\mathbf{p}({\mathbf{a}}_{\mathbf{i}}|\mathbf{z},{\mathbf{S}}^{\mathbf{i}-1})$
• 接下来介绍$a_i$的含义：首先我们有$a_i=(a_i^1,a_i^2,a_i^3,a_i^4).$假设有$m$种原子，并且有：1） a i 1 ∈ { 0 , 1 } 2 a^{1}_{i}\in\left\{0,1\right\}^{2} ai1​∈{0,1}2预测了图转换过程的终止；2） a i 2 ∈ { 0 , 1 } n a_{i}^{2}\in\left\{0,1\right\}^{n} ai2​∈{0,1}n表示第一个需要注意的节点；3） a i 2 ∈ { 0 , 1 } n a_{i}^{2}\in\left\{0,1\right\}^{n} ai2​∈{0,1}n表示第二个需要注意的节点；4） a i 2 ∈ { 0 , 1 } b a_{i}^{2}\in\left\{0,1\right\}^{b} ai2​∈{0,1}b表示两个节点之间边的类型。则这个分布$p(a_i|z,S^{i-1})$可以被分解为三个部分：1）终点预测$p(a_i^1|z,S^{i-1})$；2）节点选择$p(a_i^{2:3}|z,S^{i-1},A_i^1)$；3）边标记$p(a_i^4|z,S^{i-1},a_i^1)$
• 终点预测：作者按照下式来参数化这个分布：
  
  $\tau$代表着$softmax函数，m_t(\cdot )代表着前馈网络$。在每次转换时，我们先采样$a_i\sim p(a_i^1|z,S^{i-1})$，如果$a_i^1$代表着终结，那么整个过程将停止，现有的图结构$S^{i-1}$将会是最终被此模型生成的反应物$\mathbf{G}$。
• 节点选择：首先给定一个所有可能在转换中添加的原子的集合(比如碳，氧等原子) { v 1 , . . . , v m } \left\{v_{1},...,v_{m}\right\} {v1​,...,vm​}，并且我们把这个集合表示为$V={\bigcup }_{i=1}^m{v}_i$。我们先将现存的合成子图结构$S^{i-1}$通过添加新的孤立的节点（原子）来变为${\tilde{S}}^{i-1}$，即$S^{I-1}={\tilde{S}}_{i-1}\bigcup V$。第一个原子先从$S^{i-1}$中选择，第二个节点（原子）根据选择的第一个原子从${\tilde{S}}^{i-1}$中选择。计算方式如下：
  
  ${\beta }_1,{\beta }_2$是掩码它们是掩码，用来消除某些原子被选中的概率。并且${\beta }_1$强迫模型去从$S^{i-1}$中选择节点，${\beta }_2$防止第一个节点被再次选择。所以只有第二个节点可以从$V$中选择，并且它意味着为当前的合成子图结构$S^{i-1}$添加新的原子。所以，节点选择分布$\mathbf{p}({\mathbf{a}}_{\mathbf{i}}^{2:3}|\mathbf{z},{\mathbf{S}}^{\mathbf{i}-1},{\mathbf{a}}_{\mathbf{i}}^1)$是综合以上两个条件分布而来的。
• 边标记：边标记分布的参数化建模如下：
  

整个过程如下所示：

最后，可以通过列举将合成子$S$转化为$G$的所有可能的图转换序列来计算将合成子$S$转化为最终反应物$G$的概率: