CN114038516B - 一种基于变分自编码器的分子生成与优化方法 - Google Patents

一种基于变分自编码器的分子生成与优化方法 Download PDF

Info

Publication number
CN114038516B
CN114038516B CN202111414061.9A CN202111414061A CN114038516B CN 114038516 B CN114038516 B CN 114038516B CN 202111414061 A CN202111414061 A CN 202111414061A CN 114038516 B CN114038516 B CN 114038516B
Authority
CN
China
Prior art keywords
node
molecular
model
tree
substructure
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202111414061.9A
Other languages
English (en)
Other versions
CN114038516A (zh
Inventor
王爽
宋弢
庞善臣
李臻
魏志强
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Petroleum East China
Original Assignee
China University of Petroleum East China
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China University of Petroleum East China filed Critical China University of Petroleum East China
Priority to CN202111414061.9A priority Critical patent/CN114038516B/zh
Publication of CN114038516A publication Critical patent/CN114038516A/zh
Application granted granted Critical
Publication of CN114038516B publication Critical patent/CN114038516B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C20/00—Chemoinformatics, i.e. ICT specially adapted for the handling of physicochemical or structural data of chemical particles, elements, compounds or mixtures
    • G16C20/50—Molecular design, e.g. of drugs
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/045—Combinations of networks
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/08—Learning methods
    • G06N3/088—Non-supervised learning, e.g. competitive learning
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C20/00—Chemoinformatics, i.e. ICT specially adapted for the handling of physicochemical or structural data of chemical particles, elements, compounds or mixtures
    • G16C20/30—Prediction of properties of chemical compounds, compositions or mixtures
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C20/00—Chemoinformatics, i.e. ICT specially adapted for the handling of physicochemical or structural data of chemical particles, elements, compounds or mixtures
    • G16C20/70—Machine learning, data mining or chemometrics

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Computing Systems (AREA)
  • Chemical & Material Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • General Health & Medical Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Data Mining & Analysis (AREA)
  • Software Systems (AREA)
  • Evolutionary Computation (AREA)
  • Crystallography & Structural Chemistry (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Artificial Intelligence (AREA)
  • Biomedical Technology (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Molecular Biology (AREA)
  • Computational Linguistics (AREA)
  • Biophysics (AREA)
  • Medicinal Chemistry (AREA)
  • Pharmacology & Pharmacy (AREA)
  • Medical Informatics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Databases & Information Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明提出了一种基于变分自编码器的分子生成与优化，其特征在于，包括如下步骤：步骤1、将分子分解为子结构；步骤2、在变分自编码器和门控神经网络基础上构建生成模型，并通过该生成模型将子结构转变为分子子结构树，通过逐步添加子结构的方式生成化合物分子；步骤3、利用优化模型对化合物分子的隐向量进行优化，最终形成优化后的化合物分子，借此，本发明将分子转变为分子子结构树，在解码器过程中提出一种自回归模型，该模型通过逐步添加子结构的方式生成化合物分子，模型较好地学习到了分子子结构之间的连接关系，规避了节点添加过程中的化合价检测，解决了现有技术中的问题。

Description

一种基于变分自编码器的分子生成与优化方法
技术领域
本发明属于分子生成技术领域，特别涉及一种基于变分自编码器的分子生成与优化方法。
背景技术
分子生成与优化是药物设计的重要环节，基于深度学习的分子生成与优化一方面基于深度学习的药物生成与优化，从源头上实现了药物设计的智能化，药物分子从已有化学分子库筛选转变为从无到有地生成具有期望性质的药物分子，缩短了新药发现的周期。另一方面扩充虚拟分子库，为新药发现提供更多合理的、新颖的药物分子。
现有的基于深度学习的分子生成与优化的研究已取得一定的进展，但存在部分问题亟待解决。现有的分子生成模型中，多以原子为节点进行分子生成。为使得生成的分子符合化学价规则，生成过程的每一步都伴随原子化合价检测，这说明深度学习模型并未完全学习到分子的化合价规则，模型的学习能力和普适能力有待提升。
发明内容
本发明提出一种基于变分自编码器的分子生成与优化方法，解决了上述问题。
本发明的技术方案是这样实现的：一种基于变分自编码器的分子生成与优化方法，包括如下步骤：
步骤1、将分子分解为子结构；
步骤2、在变分自编码器和门控神经网络基础上构建生成模型，并通过该生成模型将子结构转变为分子子结构树，通过逐步添加子结构的方式生成化合物分子；
步骤3、利用优化模型对化合物分子的隐向量进行优化，最终形成优化后的化合物分子。
作为一种优选的实施方式，变分自编码器包括编码器、解码器和损失函数；
编码器，采用MSGG模型，通过三个通道提取分子特征，对于每个通道，MSGG模型采用门控神经网络和注意力机制结合的方式得到三个通道的输出，三个输出结合得到最终的分子隐向量表示；
解码器，采用自回归生成模型，该自回归生成模型由初始节点生成、分子子结构树动态表示、拓扑连接预测、新边生成、新节点生成和终止预测组成；
损失函数，给定一个训练样本，生成模型输出其对应的预测值，在训练的过程中计算预测值和真实值之间的差距数值。
作为一种优选的实施方式，初始节点生成的方法为，编码器通过MSGG模型首先将分子子结构树映射到d维的隐向量空间z，并为每个输入构建均值为μ，方差为σ的正态分布，再从d维的独立多元正态分布中采样初始向量zr，并将zr输入到一个多层感知机网络fr，为所有子结构的概率分布建模，即得到初始节点，初始节点的表达式为：
作为一种优选的实施方式，分子子结构树动态表示的方法为通过节点和边作为分子子结构树的基本单元，节点和边逐个添加的过程即为分子子结构树的生成过程，生成过程中的每一步对应一个GRU单元，一系列的GRU单元构成了GRU模型，GRU单元的数量随新节点的增加而动态增加，在t步生成的分子子结构树Tt经过GRU模型计算被记为GRU(Tt)，其表达式为：
其中，Ht是第t步分子子结构树的隐状态向量矩阵，由生成序列内的所有节点在t步的隐状态向量构成，Yt是GRU对截至t步时已生成的分子子结构树提取的整体特征向量表示，是在t步时输入到GRU模型中的第i个节点序列特征，/>是第i个节点在第t步对应的隐状态，/>是在t步第i个GRU单元的输出。
作为一种优选的实施方式，拓扑连接预测的方法为通过拓扑连接计算确定新生成节点，即父节点的位置，其拓扑连接计算的表达式为：
其中，Ct是节点Vt的拓扑连接序列，m是连接距离；
给定一个经过t步生成的分子子结构树，下一步t+1需生成的新节点的拓扑连接位置经过神经网络fl根据此时分子子结构树的状态计算得出，其计算的表达式为：
Ct+1＝Softmax(fl([Yt,zr]))
其中Ct+1是第t+1步的拓扑连接预测，Yt是上一时刻子结构树的状态表示，zr是隐空间向量。
作为一种优选的实施方式，新边生成的方法如表达式所示；
其中，是父节点特征，Nt(Vp)是节点Vp在步数t时的邻居顶点集合，Et+1,p是存在于新节点Vt+1和父节点Vp之间的边，fe是多层感知机网络，/>代表边Et+1,p的特征。
作为一种优选的实施方式，新节点生成的方法如表达式所示：
其中Vt+1是新节点，是步数t时父节点Vp已经生成的邻居节点的特征，/>是父节点特征，Yt是经过GRU模块处理后的生成子结构树在步数t时的状态描述，zr是隐向量特征，Nt(Vi)是节点Vi在步数t时的邻居顶点集合,fn是多层感知机网络。
作为一种优选的实施方式，终止预测的方法为，通过多层感知机网络fs进行计算是否终止分子子结构树的生成过程，计算表达式如下：
St+1＝sigmoid(fs([Yt+1,zr]))
其中St+1是是否终止的二值化标签，Yt+1是GRU模块处理的子结构树在步数t+1时的输出，zr是隐向量特征。
作为一种优选的实施方式，损失函数包括两部分，第一部分的表达式为：
其中表示子结构树T中的节点在节点排序π下变分自编码器的损失函数，DKL代表衡量隐向量分布与多元正态分布距离的KL散度，q(z|Tπ)代表子结构树T在节点排序π下经过编码器得到隐向量z的分布，Ν(0,Ι)代表多元正态分布；
第二部分的表达式为：
其中为解码过程的损失函数，/>表示初始节点的预测值与真实值之间的交叉熵，/>表示生成过程中第i步的拓扑连接的预测值与真实值的交叉熵，表示第i步的边的预测值与真实值的交叉熵，/>表示第i步的生成节点的预测值与真实值的交叉熵，/>表示第i步的终止预测值与真实值的二元交叉熵，n是子结构树中边的数量；
生成模型总的损失函数由第一部分和第二部分组成，其表达式如下：
作为一种优选的实施方式，优化模型优化的方法为，采用CycleGAN模型对分子隐向量空间进行优化，基于CycleGAN模型的分子优化损失由两部分构成，对抗损失函数和前向循环一致性损失函数，整体优化的损失函数为对抗损失函数与前向循环一致性损失函数之和，其表达式为：
其中，λ是对抗损失函数与前向循环一致性损失函数两者的平衡常数。
采用了上述技术方案后，本发明的有益效果是：
本发明将分子转变为分子子结构树，在解码器过程中提出一种自回归模型，该模型通过逐步添加子结构的方式生成化合物分子，模型较好地学习到了分子子结构之间的连接关系，规避了节点添加过程中的化合价检测，解决了现有技术中的问题。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图做简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为基于变分自编码器的分子生成过程示意图；
图2为BFS访问下图节点排序示例图；
图3为分子子结构树的动态特征表示图；
图4为生成分子子结构树与其拓扑连接关系图；
图5为新边生成模块示意图；
图6为新节点预测模块示意图；
图7为CycleGAN分子优化模型示意图；
图8为NeVAE模型和MSTAG模型对test分子集开展Penalized logP性质优化后数值最高的三个分子示意图；
图9为MSTAG模型从多元正态分布取样生成的新分子示意图；
图10为ZINC训练集分子示意图；
图11为针对Penalized logP指标的分子优化对比示意图；
图12为针对logP性质的分子优化对比示意图；
图13为针对Penalized logP指标开展的分子优化前后形态对比示意图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
如图1所示，一种基于变分自编码器的分子生成与优化方法，包括如下步骤：
步骤1、将分子分解为子结构；
步骤2、在变分自编码器和门控神经网络基础上构建生成模型，并通过该生成模型将子结构转变为分子子结构树，通过逐步添加子结构的方式生成化合物分子；
步骤3、利用优化模型对化合物分子的隐向量进行优化，最终形成优化后的化合物分子。
由于分子子结构具有特定的化学性质，本发明中以分子子结构作为基本单元，将分子表示为分子子结构树。分子生成过程转变为分子子结构树的生成过程，作为基本单元的子结构节点在生成过程中被逐步添加，并基提出了基于分子子结构树的自回归生成模型(Molecular Substructure Tree Autoregressive Generative Model,MSTAG)。
本发明在变分自编码器(VAE)和门控神经网络(GRU)基础上构建生成模型。变分自编码器由编码器、解码器和损失函数构成，在生成模型中，编码器将分子子结构树映射到隐向量空间，解码器从隐向量空间采样，生成新的子结构树。通过损失函数拉近隐向量分布与标准的多元正态分布的距离。变分自编码器训练完成后，解码器直接从多元正态分布中采样，生成新分子。在生成过程中，随着新节点的加入，已生成节点的邻居是不断变化的，生成树的状态也在不断改变。模型利用门控神经网络计算此过程中生成树的状态变化。为了更全面地利用分子结构信息，每一个新节点的添加都包含三个网络：拓扑预测网络，新边生成网络，新节点生成网络。后面网络的输入特征建立在前面网络的输出结果上，输入特征随着生成过程中生成树的状态作动态变化，待分子子结构树生成后，再将其转化为化合物小分子。优化模型对分子的隐向量进行优化，优化后的隐向量经解码器解码为优化后的分子。
将分子分解为子结构后，原始分子图将转换为由子结构构成的分子子结构树。在图论中，树是无向连接图，其中任何两个顶点通过一条路径连接。由树的结构可推论出，具有n个顶点的树包含(n-1)条边。在子结构树中，根顶点对应于包含原始分子图中原子序号0的子结构，该子结构决定了树顶点的唯一性。原始分子的原子序号借助于RDKit工具得到，通常序号为0的原子对应于SMILES字符串的首字符对应的原子。
子结构树是一种特殊的图，图通常被表示为G＝(V,E)，其中V是图内所有顶点的集合，V＝{v1,…,vn}，E是图内所有边的集合，每条边连接两个顶点E＝{(vi,vj)|vi,vj∈V}。拓扑连接关系是图特有的属性，在大多数有关图结构的研究中，图的拓扑连接关系用n行n列的邻接方阵M(n,n)来表示。如果在图结构中，节点vi和节点vj间有边连接，则设置邻接方阵中对应的位置为1(M(i,j)＝1,M(j,i)＝1)，其余位置设为0。这种用邻接方阵描述分子拓扑连接关系的方法常用在分子预测中作为分子图结构的描述符，后续结合频谱图卷积的方法提取分子特征，进而预测分子性质。对于分子性质预测这类问题来说，分子图的大小和节点间拓扑连接是固定的，因此对应的邻接矩阵的大小和取值也是固定的。在分子图生成问题中，若用邻接矩阵表示图并一次性生成邻接矩阵，由于原子节点之间的连接受化合价的约束，会导致生成的分子不能满足化合价规则，降低生成分子的有效性。同时采用一次性全部生成邻接矩阵的方法也会使得生成分子中原子数量固定不变，降低了生成分子的多样性。相比于一次全部生成所有节点，逐个增加节点的生成方式由于每一步增加了化合价约束，生成的分子有效性得到有效保障。在本发明的生成模型中，新的子结构树生成按照子结构逐步添加的方式，将树生成转变为子结构序列生成，且子结构序列的大小是动态可变的。
在图结构转化为节点序列的过程中，一个图结构对应n！种随机节点排序方式。对于分子子结构树的任意一个节点来说，中心节点和一阶邻居节点的连接受到化合价的约束，导致邻居节点比高阶子孙节点的关系更亲密，因此，本发明采用广度优先搜索(BFS)的方法遍历树的所有顶点。在传统图结构中，BFS首先随机选取一个节点作为根节点，优先遍历与根节点距离最近且相等的邻居节点，之后按照距离递增的顺序对图结构所有节点进行遍历。队列中先访问节点的邻居在下一个层次的子节点中被优先遍历。
起始节点选定后，在每次访问节点的过程中，相比于图的随机节点访问，BFS由于距离和邻居节点的限制，减少了下一个候选访问节点的数量。假设树结构中节点的最大度数是k(k<n)，则对于树结构中的任一节点，其访问一阶邻居最多有k！种排序方式，一个树结构对应不超过n×(k！)c(c<n)种BFS节点访问方式，其中n代表树结构的节点总数，c表示图结构中非端节点的节点数量。对于一个无向不循环图，树中的任一通路两端的节点有且只有1个邻居节点，由于BFS按照距离递增的顺序访问节点，逐级访问一阶邻居，二阶邻居，三阶邻居等，下一阶邻居节点的访问顺序的种类由上一阶的邻居节点数量决定。端节点的访问排列组合由上一阶节点决定。
如图2所示，展示了初始节点选定后无向图的BFS访问顺序。起始节点就是中心标号为0的节点，此图中，节点最大度数k＝3，节点数n＝14，非端节点c＝6。初始节点的一阶邻居、二阶邻居和三阶邻居分别用不同颜度进行标注。BFS根据初始节点0访问一阶邻居节点时，有k！(3！)种排序方式，按照BFS顺序访问二阶邻居节点时，由于二阶节点都是一阶节点的邻居，对于任意一个一阶邻居，其连接的二阶节点不超过k(3)个，每个一阶邻居节点连接的未访问的二阶邻居节点排序种类不超过k！(3！)。BFS遍历完二阶邻居节点后节点排序种类不超过k！×k！×k！(3！×3！×3！)，k！的数量与一阶邻居节点数量一致。三阶邻居节点排序种类的计算以此类推。在本实施例中，三阶邻居节点为端节点，其排序方式根据与其连接的二阶邻居节点的数量计算。因此，对于一个给定初始节点的图结构，按照BFS遍历节点，排序种类不超过(k！)c。对于未指定初始节点的无向图，初始节点的选择有n种可能。由此得出一个无向图对应的BFS节点访问排序不超过n×(k！)c(c<n)种。
本发明中，分子转化为初始节点固定的子结构树，由于其初始节点不需要选择，隐去第一项n，分子子结构树有(k！)c(c<n)种BFS节点遍历排序。初始节点固定的子结构树与BFS两者结合有效降低了节点排序组合数量。初始节点固定的子结构树与其节点的BFS排序并非一一对应关系。BFS遍历时由于每个节点的一阶邻居节点间的排序不固定，一个初始节点固定的子结构树对应多个BFS节点排序，因此，生成的不同节点排序可能对应相同的分子。在本发明中，在固定的BFS排序π下，包含n个节点的初始节点固定的子结构树的生成过程由三个序列表示：节点序列拓扑连接序列边序列/>其中/>是BFS排序π下第i个生成的子结构树顶点。/>表示在排序π下存在于第i个生成节点和前m个生成节点之间的连接关系。/>是在排序π下第i个生成的边。
类似于GraphRNN，在固定的节点排序π下，对初始节点固定的子结构树的概率分布建模任务分解为每一步生成时的节点、拓扑连接和边的联合概率分布，这个关联概率分布建立在此前已经生成的所有节点、拓扑连接和边的基础上。
其中Π是子结构树T的所有BFS节点排序，p(T)是子结构树T的概率分布。p(Cπ,Eπ,Vπ)是某种BFS节点排序π下的拓扑连接、边和节点的联合概率分布。是节点排序π下在第i个节点之前访问到的所有节点。/>指节点排序π下在第i次访问之前访问的所有拓扑连接、边。
变分自编码器(Variational Auto-Encoder,VAE)包含编码器和解码器。构成编码器的近似概率推理模型将输入/>编码为连续隐变量表示/>构成解码器的概率生成模型pθ(x|z)基于隐变量的前验概率分布p(z)将z解码为x。其中，/>和pθ是由两个深度神经网络参数化的分布。
在编码过程中，分子被编码器编码为隐向量空间的分布，本发明中的编码器采用上一章介绍的Multichannel Substructure-Graph GRU(MSGG)模型，通过三个通道提取分子特征。由上一章的实验结果可知，此模型能更全面地学习分子特征，从而更好地预测分子性质。这也说明基于此模型编码得到的隐向量能够获取分子的有效信息，在此对编码器作简略阐述：
首先，分子被分解为以化学键为基础的分子子结构，所有分子子结构构成了子结构库U,每个子结构拥有特定的序号。子结构作为节点，相邻子结构间共享的原子作为子结构树的边。利用神经网络聚合原子和化学键的信息得到节点特征Xo(s)和边特征Xe(e)。同时根据子结构在词汇库U中的序号，计算子结构的嵌入特征向量Xg(s)。
第二，在分子层面构建了节点通道、邻居节点通道和边通道。三个通道选取不同节点特征和边特征的组合作为输入。在节点通道中，输入的基本单位是子结构的节点特征Xo(s)，邻居节点通道聚集了节点和邻居节点的特征作为输入的基本单位，如公式所示。边通道不止包含了子结构树中的边特征，也融入了边连接的节点的嵌入特征，如公式所示：
其中s′是节点s的邻居节点，Xe(es,s′)是在节点s和节点s′之间的边特征。H(s)是邻居节点通道输入的基本单位，[,]是连接运算符。
其中和/>是边e连接的两个节点，/>是节点/>的嵌入特征，E(e)是边通道输入的基本单位。
第三，对于每个通道，MSGG模型采用双向的门控神经网络和注意力机制结合的方式得到三个通道的输出，三个输出结合得到最终的分子隐向量表示。
生长模型的解码器采用自回归生成模型，由六部分组成：初始节点生成、分子子结构树动态表示、拓扑连接预测、新边生成、新节点生成和终止预测。在生成过程中，生成的分子子结构树状态是不断变化的。本发明采用GRU模型计算动态变化的生成子结构树状态，为拓扑连接预测、新边生成、新节点生成和终止预测提供依据。
初始节点生成；
在每个分子子结构树中，初始节点都是第一个生成的，初始节点的选择影响到后续节点的生成。因此，初始节点对于整个分子子结构树是比较重要的。考虑到每个有机分子至少包含一个碳原子，在文献中碳原子被默认设置为生成过程的初始原子，但这在一定程度上限制了分子的多样性。在本发明的模型中，所有类型的子结构都有可能成为初始节点。
编码器MSGG首先将分子子结构树映射到d维的隐向量空间z，并为每个输入构建均值为μ，方差为σ的正态分布。从d维的独立多元正态分布中采样初始向量zr，并将zr输入到一个多层感知机网络fr，为所有子结构的概率分布建模，得到初始节点。
其中fr是一个多层感知机神经网络。
分子子结构树动态表示；
在生成过程中，节点和边是逐个添加的，已生成的子结构树的状态也随之更新。在动态生成过程中新添加的节点和边的预测，建立在当前已生成的子结构树状态的概率计算基础上。对生成树的状态进行合理准确的表达，是生成过程的重要基础和依据。在子结构树中，节点和边是子结构树的基本单元。生成的分子子结构树的状态不仅包含对于基本单元的特征描述，且包含基本单元之间的拓扑连接描述。针对节点和边的特征描述，本发明沿用第二章提到的节点特征和边特征描述方法。分子子结构树中第i个节点用字符Vi表示，其对应的节点特征表示为边用来连接相邻节点，在分子分解过程中，边通常对应于相邻两节点的公共原子。由于分子子结构的多样性，节点种类有接近四千种。相比之下，原子的种类有限，因此对应的边的类型有限，接近三十种。连接节点Vi和节点Vj的边被记为Ei,j,其特征描述为/>
在子结构树的生成过程中，节点是逐个生成并添加到已有树结构中的，这些生成的节点按照生成顺序构成一个节点序列。针对序列特性，模型引入在序列信息处理方面较为成熟的门控神经网络模型(Gated Recurrent Unit,GRU)提取子结构树的特征。在以文本、语音、音乐等为代表的普通序列里，序列信息是逐步增加的，但已经生成的序列信息不随信息的增加发生改变。与普通序列不同的是，本实施例提到的节点是树结构中的节点，拓扑连接关系是树结构的重要特征。在序列添加新节点的同时，要兼顾考虑新生成的节点与之前节点的拓扑连接关系，因为在已生成的分子子结构树中，新预测的节点必定与已生成的节点有连接。因此，已经生成的节点中还有新节点的一阶邻居，基于生成过程中拓扑连接的这一特点，将当前节点特征与直连邻居节点特征合并作为新的特征输入到GRU中。新生成的节点通过连接到已经生成的节点，将拓扑连接信息传递到所连接的节点，并更新整个序列的特征输入到GRU模型中。子结构树每增加一个新的节点，输入到GRU模型中的节点特征序列更新一次，节点特征序列在分子生成过程中迭代更新直至生成终止。
拓扑连接关系最直接反映在一个节点的邻居节点上，对于当前节点来说，所有的一阶邻居节点都与其有拓扑连接关系。生成过程中节点的邻居数量随着节点增加不断变化，在第t步，节点Vi的所有一阶邻居节点表示为Nt(Vi)。为了更好地描述子结构树的动态结构变化，输入到GRU的节点特征序列不仅要包含已生成的节点本身的特征，也要囊括在t步时已生成节点的所有一阶邻居节点特征。
在t步时对应的节点序列中第i个节点的输入特征计算如下：
其中是序列中第i个节点Vi的节点特征，Nt(Vi)是在t步时节点Vi的一阶邻居节点集合，/>对应于一阶邻居节点的特征，/>是在t步时输入到GRU模型中的第i个节点序列特征。
如图3所示，生成过程中的每一步对应一个GRU单元，第t步的第i个GRU单元的输入是 包含了t步第i个节点的特征和其邻居节点的特征。输入特征经过GRU单元的处理，得到该节点的隐向量状态和输出。
其中是t步节点特征序列中第i个GRU单元的输入，/>是在t步第(i-1)个GRU单元的输出。/>是第i个节点在第t步对应的隐状态。/>是在t步第i个GRU单元的输出。
一系列的GRU单元构成了GRU模型，GRU单元的数量随新节点的增加而动态增加。在t步生成的分子子结构树Tt经过GRU模型计算被记为GRU(Tt)，其表达式简化为：
其中Ht是第t步分子子结构树的隐状态向量矩阵，由生成序列内的所有节点在t步的隐状态向量构成，Yt是GRU对截至t步时已生成的分子子结构树提取的整体特征向量表示。
拓扑连接预测；
(1)拓扑连接计算
考虑到分子子结构树的节点间是互相连接的，当新的节点添加时，首先需要预测其连接到当前的子结构树中哪个节点。在GraphRnn中生成过程中每一步有且仅有1个新节点产生，在第t步子结构树应当有t个节点。新节点生成时，对应的拓扑连接序列为Ct＝{0,1}t，每个已经存在的节点都可能与新节点连接，且新节点所连接的节点数目最少为1。与图结构上的其他循环神经网络模型不同，本发明的生成模型目标是生成一个分子子结构树，每个节点(除初始节点外)有且只有一个父节点。在训练过程中，模型将分子子结构树延展为在BFS排序下的节点序列，先访问的节点在序列中的位置在后访问节点的前面，这意味着父节点往往出现在子节点的前面。按照BFS访问规则，先访问节点的子节点也被优先访问，序列中后访问节点的父节点的位置不可能排在先访问节点的父节点之前，具体证明过程如图4所示。因此新生成节点的父节点与前一个节点的父节点的位置有两种情况：两者父节点相同或者新生成节点的父节点在前一个节点的父节点之后。同理，子结构树展开的BFS节点序列中每一个节点与之前节点的连接距离都局限在前一个节点的父节点位置和它本身位置之间。为了减轻神经网络的计算负担，连接距离用m表示，在训练过程中，此值限定为12。其拓扑连接序列由以下公式表示，新节点的父节点仅存在于与之相邻最近的m个节点中，且在这m个节点中有且只有一个父节点。
其中Ct是节点Vt的拓扑连接序列。
给定一个经过t步生成的子结构树，下一步t+1需生成的新节点的拓扑连接位置经过神经网络fl根据此时刻子结构树的状态计算得出。
Ct+1＝Softmax(fl([Yt,zr]))
其中Ct+1是第t+1步的拓扑连接预测，Yt是上一时刻子结构树的状态表示，zr是隐空间向量。
(2)证明
命题：假设存在有n个节点的子结构树，按照BFS遍历子结构树的节点，得到节点序列v1,…,vn，对于
如果(vi,vj)∈E，那么并且/>
如果(vi,vj)∈E，那么
如果(vi,vj)∈E，但是那么/>并且存在(vi′,vk)∈E。
如果(vi,vj)∈E，但是那么/>并且存在(vk,vj′)∈E或者(vj,vj′)∈E。
观察1：在子结构树中，除了初始节点v1，所有节点有且仅有一个父节点，按照BFS对子结构树进行遍历得到节点序列，对于任一非初始节点vi(i>1)，其父节点在BFS节点序列中的位置总在节点vi之前。对于节点vi和vj(1<i<j)，节点vi的父节点在节点vj的父节点前面。
观察2：在子结构树中，除了叶节点之外的所有节点都有其子节点。若一个节点有不少于一个子节点，那么子节点在BFS序列中的序号是连续的。
证明1:若存在(vi,vj)∈E且i<j，可以推断出节点vi是节点vj的父节点，从观察1可得，节点vj的父节点有且只有一个，且位置在vj之前，已知(vi,vj)∈E，因此，不存在vj之前的其他节点与vj有拓扑边连接，即对于序号i′<j和序号k<j，满足并且
证明2：对于节点vj，vj′满足条件j<j′，从观察1可得vj的父节点在vj′的父节点之前，若存在(vi,vj)∈E，那么节点vj′的父节点的序号应该大于i，因此
证明3：若(vi,vj)∈E但是由观察2可推断，节点vj是节点vi的第一个子节点，对于任意的i<k<j，一定不存在(vi,vk)。若节点vi的子节点数量大于1，那么存在(vi,vj′)∈E。如果节点vi只有一个子节点，那么节点vj′的父节点一定在节点vi的后面，存在(vk,vj′)∈E或者(vj,vj′)∈E。
证明4：若存在(vi,vj)∈E但是由观察2可推断，节点vj是节点vi的最后一个子节点，对于序号j<j′，/>由于节点vj′的父节点在节点vi的后面，存在(vk,vj′)∈E或者(vj,vj′)∈E。
根据以上证明，可以划定新生成节点的连接范围。假设存在一个已生成序列，且生成序列的最后一个节点为vj，新预测的节点用vj′表示。从假设1可得，在已生成的序列中有且仅有一个作为新生成节点vj′的父节点，即生成序列中仅有一个节点与vj′连接。从假设2可推断，新预测节点的连接位置不能位于最后一个节点vj的父节点之前。由假设3和假设4可得出新节点vj′可能的连接位置位于最后一个节点vj的父节点和其后所有的生成节点之间。由此证明，新生成节点并非与所有已生成节点间都有连接的可能，其连接位置限定在一定范围内，用长度为m的拓扑连接序列预测其连接位置是合理可行的。
新边生成；
如图5所示，经过拓扑连接预测后，新生成节点的连接位置得以确定，称作新生成节点的父节点Vp。得到其父节点后，在子结构树中，边是由相邻两个子结构共有原子构成的，因此，根据父节点和现有子结构树的状态，可确定其连接的边类型。
在子结构树中，节点由几个原子和连接原子的化学键构成，相邻子结构之间共有的原子被描述为子结构树中的边。只有两个节点包含相同的原子，两个节点间才可能连接。因此，与边连接的节点必须包含边，即所预测的边必须包含在父节点Vp中，边预测的目标就是学习一个条件概率p(e|Vp)。选定父节点Vp后，边类型的概率p(e|Vp)分布在父节点Vp所包含的边的集合上。例如，若Vp包含碳原子C和氧原子O，那么边e一定包含C或O，且是节点Vp的子结构。因此，父节点Vp的特征是预测边类型的重要依据之一。
除此之外，边是存在于相邻节点间的公共子结构，由此可以推出，边暗示了节点内部的可能的连接位置。如果边的连接位置被占用了，那么新的边不能连接到相同的位置，因此，父节点Vp已经存在的边是预测即将生成的边类型的另一重要依据。
在本发明模型中，将父节点Vp连接的所有现存边的信息进行汇集，父节点Vp在步数t+1所预测的新边的概率分布用p(Et+1,p│Vp)表示，此概率分布的计算建立在步数t时父节点特征和已存在边特征的集合的基础上，结合两者的特征并送入多层感知机神经网络以计算新边类型的概率分布。针对父节点Vp在步数t时边类型预测如下：
其中Nt(Vp)是节点Vp在步数t时的邻居顶点集合，Et+1,p是存在于新节点Vt+1和父节点Vp之间的边，fe是多层感知机网络。
新节点生成
在生成过程中，当拓扑连接位置和连接边类型确定后，新节点的类型也随之确定。与边生成类似，考虑到边是两个相邻节点间的共有子结构，新顶点的预测也受连接边的限定，新生成的节点应当包含连接边。除此之外，两个节点的连接由于包含两者的连接位置，因此，其结合需满足化合价的约束，由此推理可得，新节点的生成受其父节点的限制。由于分子中节点和其邻居之间的所有连接都要受到化合价约束，在新节点预测中也要考虑其邻居节点的信息。由于分子中官能团的力场的存在，不同子结构之间互相影响，所以新节点的预测也受到已生成的子结构树的影响。新节点预测的目标是学习包含边结构Ep,t+1的子结构集合上的条件概率分布，t+1步时对应的条件概率分布描述为p(Vt+1|Vp,Ep,t+1,Tt)。条件概率分布的计算建立在父节点特征、其连接的边特征和父节点已有的邻居节点特征、生成子结构树的特征和初始的隐向量特征基础上，同样通过多层感知机网络fn为条件概率分布建模。新节点预测如下：
其中是步数t时父节点Vp已经生成的邻居节点的特征，Yt是经过GRU模块处理后的生成子结构树在步数t时的状态描述，zr是隐向量特征。
终止预测；
伴随新的节点和新边的生成，子结构树的状态也随之更新，新的子结构树状态通过GRU模块重新计算，用GRU(Tt+1)表示。每进行一次更新，都要预测一下是否终止生成过程，同样通过多层感知机网络fs进行计算。
St+1＝sigmoid(fs([Yt+1,zr]))
其中St+1是是否终止的二值化标签，Yt+1是GRU模块处理的子结构树在步数t+1时的输出。
损失函数；
在训练的过程中，给定一个训练样本，模型输出其对应的预测值，期望在训练过程中预测值与真实值得差距越来越小，这就需要给模型指明一个前进的方向，而损失函数承担了这一个角色，其用来衡量预测与真实值的差异。随着训练数据的不断输入，模型向着损失函数越小的方向学习参数权重。因此损失函数定义的好坏，与模型的性能表现息息相关。
本发明中生成模型的最终目标为，给定一个多元正态分布的采样值，通过生成模型学到的化学编码知识和分子的结构信息，输出满足化学价规则的新颖的化学小分子。从中可分析出，模型表现性能的好坏与两方面有关，一是输入与多元正态分布的偏差程度，二是生成过程中结构预测的准确性。鉴于本章提出的生成模型整体架构为变分自编码器，通过编码器将输入的小分子映射到隐向量空间，在隐向量空间中，用Kullback-LeiblerDivergence(KL散度)衡量隐向量与多元正态分布的差异程度，将隐向量拉近多元正态分布，从而为后续的随机采样生成奠定基础。在解码的过程中，本章的模型包含四个预测，根节点预测、拓扑连接预测、边预测、节点预测，将每个阶段的预测值与真实值的差异加入到损失函数的构成中，作为损失函数的一部分，预测值与真实值的差异用交叉熵损失进行评估。
其中，变分自编码器部分的KL损失函数表示如下：
其中表示子结构树T中的节点在节点排序π下变分自编码器的损失函数，DKL代表衡量隐向量分布与多元正态分布距离的KL散度，q(z|Tπ)代表子结构树T在节点排序π下经过编码器得到隐向量z的分布，Ν(0,Ι)代表多元正态分布。
解码器由自回归生成模型构成，其中，初始节点的生成只与隐向量有关，而拓扑连接预测、边预测和节点预测的预测网络输入都随着生成过程中生成树状态的变化而更新。一个子结构树对应一个初始节点，但拓扑连接、边和节点的预测都与节点总数量有关系，这三部分的损失函数是生成过程每一步的损失函数加和。
其中为解码过程的损失函数，/>表示初始节点的预测值与真实值之间的交叉熵，/>表示生成过程中第i步的拓扑连接的预测值与真实值的交叉熵。/>表示第i步的边的预测值与真实值的交叉熵，/>表示第i步的生成节点的预测值与真实值的交叉熵。/>表示第i步的终止预测值与真实值的二元交叉熵。n是子结构树中边的数量。拓扑连接预测、边预测和节点预测都是在节点排序π下开展的。
生成模型总的损失函数由上述两部分构成，其综合考虑了分子子结构树的多种排序方式：
训练的过程结合了Teacher forcing训练策略，在生成过程中的每一步，都伴随预测结果的产生，新的拓扑连接、新的边、新的节点，但预测的结果并不是一定准确的，尤其在模型训练的初始阶段。本章提出的生成模型的解码部分由自回归模型构成，在自回归模型中，先生成的输出结果将作为后续生成步骤的输入。若将不准确的预测结果传入后续生成过程，将影响后续拓扑连接、边和节点的类型预测。利用Teacher forcing训练策略用训练集的期望输出(真实值)替代模型的预测输出，传入到后续的生成过程，以确保后续生成过程的计算建立在正确的输入上。
当新生成过程结束后，生成的子结构树将转变为分子。子结构节点之间的连接由其内部标定的连接位置决定，连接处要满足化合价的规则。对于最终的分子，若原子的化合价还没被填满，则用氢原子进行化合价填充。
基于CycleGAN的分子优化模型；
上述提到的生成模型的目标是学习分子的化合价规则，最终达到模型在正态分布中随机采样即可生成有效化学分子的目的。生成符合化学价规则的新的化合物分子可以扩充化合物分子库，但为了使得生成的分子更具有应用价值，我们希望生成的分子具有期望的化学性质。因此，需要对模型增加新的规则以生成具有特定性质的化学分子。在本实施例中，采用了CycleGAN模型对分子隐向量空间进行优化，由于不同性质的分子映射得到的隐向量也不同，通过对抗学习的方法先生成具有特定性质的分子对应的隐向量，再利用解码器的自回归生成模型将隐向量解码为化学分子。
CycleGAN的提出是为了解决图像的风格迁移问题，与熟知的Pix2Pix对比可以发现，Pix2Pix在风格转移的过程中两个域之间有相同的训练样本，即要求两个域中的数据需要两两配对。但很多现实场景中两个域并非有配对的数据，CycleGAN的提出解决了在源域和目标域之间没有一对一映射情况下的风格转换，例如将一类图片转换成另一类图片。在分子优化的问题上，训练集的分子之间不存在配对的关系，因此，采用CycleGAN作为优化模型，模型架构如图7所示，将不具备期望性质的分子对应的隐向量分布定义为源域X，具备期望性质的分子对应的隐向量分布定义为目标域Y。将目标域上的分子“风格”迁移到源域上去，达到优化分子性质的目的。
具体算法过程如下：针对源域X(不具备期望性质的分子隐向量)和目标域Y(具备期望性质的分子隐向量)，其中的分子隐向量数据分别记为x和y。设立两个生成器GXY和GYX，从源域X的数据中选择数据x输入到生成器GXY中，并生成隐向量GXY(x)，而生成器GYX从目标域Y中选择数据y输入到生成器GYX中，并生成隐向量GYX(y)。此时设立两个判别器DX和DY，DX负责分辨源域的真实数据x和生成器GYX从目标域转换而来的数据GYX(y)，DY则负责分辨目标域的真实数据y和生成器GXY从源域转换而来的数据GXY(x)。
基于CycleGAN的分子优化损失由两部分构成：对抗损失函数和前向循环一致性损失函数，其中对抗损失函数分为两部分：
其中，GXY从源域X选取数据生成隐向量GXY(x)，并使GXY(x)尽量接近于目标域Y中的数据。GXY致力于最小化损失函数，而DY尝试最大化。
同理，对于目标域Y到源域X的风格转变，其损失函数表示如下：
其中，GYX从目标域Y选取数据生成隐向量GYX(y)，并使GYX(y)尽量接近于源域X中的数据。GYX致力于最小化损失函数，而DX尝试最大化。
除了上述两个目标函数外，CycleGAN还设置了前向循环一致性的损失函数，评估输入x经过生成器GXY和生成器GYX后是否能够转变回原来的x，即x→GXY(x)→GYX(GXY(x))≈x，以保障两个生成器在目标域和源域上有普适性，其损失函数描述如下：
最终整个的优化损失函数表示为：
整体优化的损失函数为对抗损失函数与前向循环一致性损失函数之和，λ作为两者的平衡常数。
实验结果及分析
实验设置
数据集
本实施例采用ZINC数据集作为训练数据集，MSTAG模型在ZINC数据集上学习小分子结构特性和化合价规则。ZINC数据集由25万个类药化合物分子构成，数据集中的分子的原子数量不固定，平均每个分子约由23个非氢原子构成，ZINC数据集是在生成模型中常见的数据集，同样较常出现的还有QM9数据集，此数据集中的分子大小都是固定的，常用在对分子大小有要求的模型算法中，ZINC数据集相比于QM9数据集，分子多样性更大，结构更为复杂。对于ZINC数据集来说，将其中的分子转变为分子子结构树后，得到的子结构数量为3794，对应的边集合中边的数量为27。
评价指标
考虑到分子的复杂性和特殊性，生成模型的主要目标是学习分子的语法和语义特性，并生成满足化合价规则的新分子。这里的语法指模型在生成过程中要遵守分子的化学规则和化合价约束，模型只有学习到这些规则才能生成有实际应用价值的新分子。分子的语义特性在本实施例中指分子的基本属性，例如分子量、原子数量、原子类型等构成属性和logP、SA、QED等性质属性。分子的性质与分子的结构息息相关，不同的子结构组成的分子表现出的性质也有所不同。一个好的生成模型不仅可以学习分子“语法”特性，也应学习分子内部的这种“语义”特性。在本实施例中，对生成模型这两方面的表现通过两个衡量标准进行评估，第一个衡量标准由分子的可用性、新颖性和独特性组成，评价生成模型是否学到了分子的语法特性。第二个衡量标准是比较训练集中分子的属性分布和生成集合中分子的属性分布，通过性质分布的比较评价生成模型是否学到了分子的语义特性。
新颖性：在生成模型中，模型经过对训练集的分子数据的学习，已较为熟悉训练集中分子的特性，有较大可能生成与训练集中分子相同的化学分子。但生成模型的主要作用之一是扩充现有的虚拟分子库，为虚拟筛选等提供充足的备选分子。若生成分子与原有分子重合较多，则无法起到扩充分子库的作用。因此，设立了新颖性的衡量指标，用以评价模型生成新分子的能力。训练集中的分子集合用Vt表示，生成的分子集合用Vg表示，新颖性的计算方式如下：
其中set(Vg)表示对生成分子集Vg去重，set(Vg)∩Vt表示生成分子与训练集中分子重复的数量。如公式所示，新颖性表示生成分子集合中未在训练集出现的生成分子所占的比例。
独特性：新颖性用来衡量生成分子与训练集中分子的差异性，我们希望生成模型生成的分子不仅不在训练集中出现，而且在新生成的分子集合中，重复的越少越好，独特性用来衡量生成分子集合内部的差异性。其计算方式如下所示：
有效性：分子的有效性衡量生成模型学习分子化合价规则的能力，其计算方式为生成分子的集合中满足化合价规则的分子所占的比例。
无化合价检验的有效性(Validity w/o check)：在大多数基于序列的分子生成模型中，所生成分子都达到了100％的有效性。这是因为在序列生成的每一步生成过程中，新节点的添加都要做一次化合价检验，只有满足化合价规则的节点才能被添加到生成队列中。在以原子为基本实施例点的分子图结构中，化合价检测运用更为便捷，即检查所连接的化学键是否超出了原子化合价即可。在本章生成模型中，基本实施例点为分子子结构，由于一个子结构包含不止一个原子，化合价检测的步骤不便于在每个原子上展开，为此根据片段之间的连接特点设计了一个初步的粗筛查。
在MSTAG模型的生成过程中，新的节点添加时首先要预测在已生成的分子子结构树中新节点的添加位置，然后根据当前连接的已生成节点的状态预测连接边的类型，由于分子子结构树中的边是两个相邻子结构的公共结构，因此若边结构被包含在所连接的节点结构内，那么此边可以被添加。当边被添加后，新节点类型确定时，同理，若新节点包含所连接边结构，才可能被添加到分子子结构生成树中。这种粗略的过滤可以进行初步的化合价筛选。而Validity w/o check指在边生成和节点生成的步骤中均未实施这种初步的化合价筛选。
实验性能
生成模型的化学规则学习性能比较
为了与其他方法进行对比，模型从多元正态分布中随机取样，生成了1万个化合物分子，同时，从训练集中随机选取了相同数量的分子进行对比。对于生成的新分子和从训练集中选取的1万个分子分别计算了新颖性、独特性、有效性、Validity w/o check四个指标。
参与对比的模型有GraphVAE,GraphNVP,MRNN,GraphAF和GCPN。GrahNVP模型以分子图为操作对象，基于流的生成模型(flow based model)一次生成分子图的特征矩阵和邻居矩阵。MRNN利用循环神经网络以原子逐渐添加的方式生成分子。GraphAF利用基于流的自回归模型以原子逐个添加的方式生成分子。GCPN将强化学习的方法应用到分子生成模型中，同样采取了原子逐个添加的方式。这几个模型都是以分子图为对象，对比结果如下表所示。
从上表中可以看出在分子的有效性比较中，按照序列逐个生成原子的方式进行分子生成的模型表现优于一次性生成所有原子和连接关系的模型，例如GraphVAE和GraphNVP。这是因为在原子按照序列生成的过程中，每添加一个新的原子，都对它和与其连接的原子进行化合价检测，这就保证了分子生成过程中不满足化合价约束的原子不会被添加到序列中，确保最终得到的分子是满足化合价约束的有效的分子。至于在序列生成的模型的比较中，评价模型是否学到了相应的化学规则的标准是validity w/o check。在所有序列生成的模型中取消化合价检测的步骤，根据最终生成的分子中满足化合价规则的分子所占的比例评估模型的性能。一个充分学习到分子化学规则的模型即使去掉化合价检测也能生成满足化合价约束的分子。从上表中可以看到，在没有化合价约束的情况下，模型MSTAG生成的分子达到了100％的有效性，而在其他序列生成的模型中，表现最好的模型是GraphAF，性能仅达到了68％，模型表现超出其47％。充分说明本章提出的模型已学习到了分子内部的化学价规则。
生成模型的属性学习比较
关于生成模型的性能度量，模型除了学习分子的化学规则外，我们还希望模型能学习到分子的内部“语义”信息，即分子的性质和与性质关联的结构特征。在本小节中，列出了分子的5个属性值：分子量(Molecule Weight,MW)、辛醇-水分配系数(water-octanolpartition coefficient,logP)、可合成性得分(Synthetic accessibility score,SA)、类药性定量估计(Quantitative Estimation of Drug-likeness,QED)和每个分子的平均原子数量(Number of atoms,NA)。统计了训练集中随机选取的1万个小分子的五个属性的平均值，同时，也对生成的1万个小分子做了5个属性的数据统计，结果呈现于下表中：
从表中观察可得，分子的形状信息类似于分子量(MW)和原子数(NA)，训练集和生成分子对应的属性均值在较小的差异范围内。分子的性质信息类似于辛醇-水分配系数(logP)、可合成性(SA)和类药性定量估计(QED)都在合理的分子限定范围内。
基于Penalized logP的分子优化性能比较
本实施例以分子性质Penalized logP为优化目标，从ZINC数据集中选择了800个Penalized logP值最低的小分子，用基于CycleGAN的分子优化模型进行优化。
Penlized logP是一个综合性的性质指标，其计算方式如下：
penalizedlogP(m)＝logP(m)-SAS(m)-ringpenalty(m)
其中m表示分子，logP(m)代表辛醇-水分配系数，是衡量药物分子类药性的一个重要指标。SAS(m)衡量了分子合成的可行性。ring_penalty(m)的设置是对分子中碳环加以约束，避免分子产生不切实际的大碳环。
在基于Penalized logP指标对分子优化中，我们与其他方法进行了比较，包含ChemVAE，GrammarVAE和NeVAE。ChemVAE是第一个基于变分自编码器实现分子从无到有的设计的生成模型，其以SMILES字符串作为分子的描述符，以CNN作为编码器，RNN作为生成器，基于高斯过程(Gaussian process,GP)模型和贝叶斯优化在分子隐向量空间进行优化。GrammarVAE模型以SMILES字符串作为分子描述符，构建语法规则解析树，提取SMILES编码规则，对分子的规则树通过CNN编码，并利用RNN解码，训练变分自编码器学习SMILES的语法规则。同样基于高斯过程和贝叶斯优化对隐向量空间进行优化。NeVAE与前两者不同，是建立在变分自编码器之上的针对无向分子图的逐步生成模型。在编码过程利用原子的特征和边权重，每个原子的特征汇集了邻近的多个邻居原子的信息，得到嵌入向量表示。在解码过程中，根据隐向量，先确定分子图中原子节点的数量和特征，之后逐个添加边，随着边和边权重的动态更新解码出完整分子图。利用贝叶斯优化方法对隐向量进行优化，再解码出优化后的分子。
MSTAG模型与三个模型的优化对比结果展示在下表中，经过模型优化后，对产生的最优Penalized logP结果进行比较，表中前三个模型的数据来自NeVAE。观察可得，本文MSTAG模型取得了最佳的Penalized logP值。
此外，为了避免优化结果的偶然性，更客观地评价模型优化能力，与NeVAE类似，将Penalized logP得分最高的三个分子进行展示，结果如图8所示，经过MSTAG模型优化后得分最高的三个分子的Penalized logP值均高于3，高于NeVAE最佳Penalized logP值。
本发明提出了一个基于变分自编码器的分子生成和优化模型。首先，将分子重构为分子子结构树，利用MSGG编码器将分子子结构树编码到隐向量空间；其次在解码器中构建了MSTAG分子子结构树生成模型，对隐向量空间进行采样，利用MSTAG自回归生成模型生成分子子结构树；最后再将分子子结构树转变为化合物分子。在优化模型中，利用CycleGAN模型对分子的隐向量空间进行优化，利用训练好的MSTAG生成模型基于优化向量生成优化后的化合物分子。结果表明，MSTAG模型对分子的化合价规则做了充分的学习，在没有化合价检验的情况下生成分子的有效性达到100％。同时优化模型明显提升了分子性质。
以上仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种基于变分自编码器的分子生成与优化方法，其特征在于，包括如下步骤：
步骤1、将分子分解为子结构；
步骤2、在变分自编码器和门控神经网络基础上构建生成模型，并通过该生成模型将子结构转变为分子子结构树，通过逐步添加子结构的方式生成化合物分子；
步骤3、利用优化模型对化合物分子的隐向量进行优化，最终形成优化后的化合物分子；
编码器，采用MSGG模型，通过三个通道提取分子特征，对于每个通道，MSGG模型采用门控神经网络和注意力机制结合的方式得到三个通道的输出，三个输出结合得到最终的分子隐向量表示；
解码器，采用自回归生成模型，该自回归生成模型由初始节点生成、分子子结构树动态表示、拓扑连接预测、新边生成、新节点生成和终止预测组成；
损失函数，给定一个训练样本，生成模型输出其对应的预测值，在训练的过程中计算预测值和真实值之间的差距数值；
所述终止预测的方法为，通过多层感知机网络fs进行计算是否终止分子子结构树的生成过程，计算表达式如下：
St+1＝sigmoid(fs([Yt+1,zr]))
其中St+1是是否终止的二值化标签，Yt+1是GRU模块处理的子结构树在步数t+1时的输出，zr是隐向量特征；
所述损失函数包括两部分，第一部分的表达式为：
其中表示子结构树T中的节点在节点排序π下变分自编码器的损失函数，DKL代表衡量隐向量分布与多元正态分布距离的KL散度，q(z|Tπ)代表子结构树T在节点排序π下经过编码器得到隐向量z的分布，Ν(0,Ι)代表多元正态分布；
第二部分的表达式为：
其中为解码过程的损失函数，/>表示初始节点的预测值与真实值之间的交叉熵，/>表示生成过程中第i步的拓扑连接的预测值与真实值的交叉熵，表示第i步的边的预测值与真实值的交叉熵，/>表示第i步的生成节点的预测值与真实值的交叉熵，/>表示第i步的终止预测值与真实值的二元交叉熵，n是子结构树中边的数量；
生成模型总的损失函数由第一部分和第二部分组成，其表达式如下：
2.根据权利要求1所述的一种基于变分自编码器的分子生成与优化方法，其特征在于，所述初始节点生成的方法为，编码器通过MSGG模型首先将分子子结构树映射到d维的隐向量空间z，并为每个输入构建均值为μ，方差为σ的正态分布，再从d维的独立多元正态分布中采样初始向量zr，并将zr输入到一个多层感知机网络fr，为所有子结构的概率分布建模，即得到初始节点，初始节点的表达式为：
3.根据权利要求1所述的一种基于变分自编码器的分子生成与优化方法，其特征在于，所述分子子结构树动态表示的方法为通过节点和边作为分子子结构树的基本单元，节点和边逐个添加的过程即为分子子结构树的生成过程，生成过程中的每一步对应一个GRU单元，一系列的GRU单元构成了GRU模型，GRU单元的数量随新节点的增加而动态增加，在t步生成的分子子结构树Tt经过GRU模型计算被记为GRU(Tt)，其表达式为：
其中，Ht是第t步分子子结构树的隐状态向量矩阵，由生成序列内的所有节点在t步的隐状态向量构成，Yt是GRU对截至t步时已生成的分子子结构树提取的整体特征向量表示，是在t步时输入到GRU模型中的第i个节点序列特征，/>是第i个节点在第t步对应的隐状态，是在t步第i个GRU单元的输出。
4.根据权利要求1所述的一种基于变分自编码器的分子生成与优化方法，其特征在于，所述拓扑连接预测的方法为通过拓扑连接计算确定新生成节点，即父节点的位置，其拓扑连接计算的表达式为：
其中，Ct是节点Vt的拓扑连接序列，m是连接距离；
给定一个经过t步生成的分子子结构树，下一步t+1需生成的新节点的拓扑连接位置经过神经网络fl根据此时分子子结构树的状态计算得出，其计算的表达式为：
Ct+1＝Softmax(fl([Yt，zr]))
其中Ct+1是第t+1步的拓扑连接预测，Yt是上一时刻子结构树的状态表示，zr是隐空间向量。
5.根据权利要求1所述的一种基于变分自编码器的分子生成与优化方法，其特征在于，所述新边生成的方法如表达式所示；
其中，是父节点特征，Nt(Vp)是节点Vp在步数t时的邻居顶点集合，Et+1,p是存在于新节点Vt+1和父节点Vp之间的边，fe是多层感知机网络，/>代表边Et+1,p的特征。
6.根据权利要求5所述的一种基于变分自编码器的分子生成与优化方法，其特征在于，所述新节点生成的方法如表达式所示：
其中Vt+1是新节点，是步数t时父节点Vp已经生成的邻居节点的特征，/>是父节点特征，Yt是经过GRU模块处理后的生成子结构树在步数t时的状态描述，zr是隐向量特征，Nt(Vi)是节点Vi在步数t时的邻居顶点集合,fn是多层感知机网络。
7.根据权利要求1所述的一种基于变分自编码器的分子生成与优化方法，其特征在于，所述优化模型优化的方法为，采用CycleGAN模型对分子隐向量空间进行优化，基于CycleGAN模型的分子优化损失由两部分构成，对抗损失函数和前向循环一致性损失函数，整体优化的损失函数为对抗损失函数与前向循环一致性损失函数之和，其表达式为：
其中，λ是对抗损失函数与前向循环一致性损失函数两者的平衡常数。
CN202111414061.9A 2021-11-25 2021-11-25 一种基于变分自编码器的分子生成与优化方法 Active CN114038516B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111414061.9A CN114038516B (zh) 2021-11-25 2021-11-25 一种基于变分自编码器的分子生成与优化方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111414061.9A CN114038516B (zh) 2021-11-25 2021-11-25 一种基于变分自编码器的分子生成与优化方法

Publications (2)

Publication Number Publication Date
CN114038516A CN114038516A (zh) 2022-02-11
CN114038516B true CN114038516B (zh) 2024-04-19

Family

ID=80145538

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111414061.9A Active CN114038516B (zh) 2021-11-25 2021-11-25 一种基于变分自编码器的分子生成与优化方法

Country Status (1)

Country Link
CN (1) CN114038516B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116130036B (zh) * 2023-01-09 2024-03-01 四川大学 一种基于图表示的金属有机框架的逆向设计方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110970099A (zh) * 2019-12-10 2020-04-07 北京大学 一种基于正则化变分自动编码器的药物分子生成方法
CN112071373A (zh) * 2020-09-02 2020-12-11 深圳晶泰科技有限公司 药物分子筛选方法及系统
CN112397157A (zh) * 2020-10-28 2021-02-23 星药科技(北京)有限公司 基于子图-变分自编码结构的分子生成方法
CN113327651A (zh) * 2021-05-31 2021-08-31 东南大学 一种基于变分自编码器和消息传递神经网络的分子图生成方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110970099A (zh) * 2019-12-10 2020-04-07 北京大学 一种基于正则化变分自动编码器的药物分子生成方法
CN112071373A (zh) * 2020-09-02 2020-12-11 深圳晶泰科技有限公司 药物分子筛选方法及系统
CN112397157A (zh) * 2020-10-28 2021-02-23 星药科技(北京)有限公司 基于子图-变分自编码结构的分子生成方法
CN113327651A (zh) * 2021-05-31 2021-08-31 东南大学 一种基于变分自编码器和消息传递神经网络的分子图生成方法

Also Published As

Publication number Publication date
CN114038516A (zh) 2022-02-11

Similar Documents

Publication Publication Date Title
Popova et al. MolecularRNN: Generating realistic molecular graphs with optimized properties
Fu et al. Core: Automatic molecule optimization using copy & refine strategy
CN110147892B (zh) 基于变分轨迹上下文感知的人类移动模式推测模型、训练方法及推测方法
Fu et al. Mimosa: Multi-constraint molecule sampling for molecule optimization
Faez et al. Deep graph generators: A survey
CN107992976B (zh) 热点话题早期发展趋势预测系统及预测方法
CN108647226B (zh) 一种基于变分自动编码器的混合推荐方法
CN112364975A (zh) 基于图神经网络的终端运行状态预测方法及系统
CN111428848B (zh) 基于自编码器和3阶图卷积的分子智能设计方法
CN109214503B (zh) 基于kpca-la-rbm的输变电工程造价预测方法
CN109583588B (zh) 一种短期风速预测方法及系统
Xu et al. Graph partitioning and graph neural network based hierarchical graph matching for graph similarity computation
CN114676687A (zh) 基于增强语义句法信息的方面级情感分类方法
CN114038516B (zh) 一种基于变分自编码器的分子生成与优化方法
CN108733921B (zh) 基于模糊信息粒化的变压器绕组热点温度波动范围预测方法
CN115345222A (zh) 一种基于TimeGAN模型的故障分类方法
CN114528971A (zh) 一种基于异质图神经网络的图谱频繁关系模式挖掘方法
CN114385910A (zh) 基于知识追踪的在线学习内容推荐方法及系统
Wang et al. Reinforcement-enhanced autoregressive feature transformation: Gradient-steered search in continuous space for postfix expressions
CN112183721B (zh) 一种基于自适应差分进化的组合水文预测模型的构建方法
CN113076545A (zh) 一种基于深度学习的内核模糊测试序列生成方法
Petkov et al. Dag-wgan: Causal structure learning with wasserstein generative adversarial networks
CN105512249A (zh) 一种基于紧凑进化算法的本体匹配方法
Singh et al. ChemoVerse: Manifold traversal of latent spaces for novel molecule discovery
CN117133116B (zh) 一种基于时空关联网络的交通流预测方法及系统

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant