基于神经网络的分子性质预测算法研究进展

朱洪翔, 傅钰江, 李雪, 陈博

(中石化(大连)石油化工研究院有限公司, 大连 116045)

分子性质预测对新型材料和药物设计具有重要意义,人们通常利用从头算方法,在此基础上逐一计算物理化学性质,这样的方法计算成本高,计算效率低,不利于相关研究的快速开展。随着计算机科学技术的发展,对于原子个数较少且重要的分子体系,计算机的算力已经可以足够支撑高精度的从头算能量点。2004年Brown等[1]提出交换对称多项式来拟合多原子分子势能面的方法,该方法利用从头算获得的能量点和能量梯度,进行交换对称多项式拟合,取得了良好的效果。2006年Qiu等[2]在多参考组态相互作用上,使用augcc-pV5Z基组并结合Davidson纠正,在使用分子对称的基础上,计算了15 000个从头算能量点,构造了描述氟原子和氢气反应的全域势能面。在石油化工领域,中国研究者[3-4]利用分子模型对页岩开采中的甲烷吸附过程和有机胺页岩抑制剂进行的分析,为油气开采提供理论依据。

近年来,利用机器学习和神经网络(neural network,NN)相关方法预测分子性质,引起了众多的关注[5-6],对于拟合精度、势能面维度以及相同原子的交换对称性问题都得到了一定程度的发展。机器学习在分子科学领域的发展主要从药物分子设计和性质计算开始[7-9],后来将应用领域扩展到分子计算领域,Rupp等[10]引入一种机器学习模型来预测大量有机分子原子化能。Faber等[11]针对118 000个有机小分子的13种分子性质展开研究,对应的机器学习模型可以有效地将预测误差从10 kcal/mol减小到3 kcal/mol。

神经网络则从势能面预测问题切入,2013年Chen等[12]提出通过系统增加样本点的方法来进行NN拟合势能面精度的方法。在H2+OH-反应体系中计算了17 000个从头算能量样本点,使用NN方法构造该反应体系的高精度全域势能面。在与Brown等[1]提出的交换对称多项式方法能量点数量相同的情况下,利用了更少的计算资源,取得了相当的效果。近几年,在生命科学领域,利用深度学习方法预测药物和蛋白质性质受到了广泛的关注。袁露[13]通过图神经网络构建出药物虚拟筛选模型;徐大杰[14]、丁小雨[15]通过生成模型和主动学习模型参与设计药物分子;顾耀文等[16]将注意力机制融入图神经网络中,对药物的毒性和代谢进行预测。曹晓勇[17]利用局部自由能对蛋白质结构进行优化。

通过文献调研,梳理近5年基于神经网络的分子性质预测算法的发展时间轴,如图1所示,分析分子性质预测领域相关关键词,绘制关键词关系矩阵图如图2所示。分子模拟(molecular dyna-mics,MD)权重最大,与分子模拟相关联的关键词中,第一性原理计算(ab-initio)权重最大。在分子模拟领域中,神经网络和深度学习是应用最广泛的智能算法,主要应用于力场和能量场的预测。本文结构如下:开篇介绍分子性质预测算法及其发展,第1节阐述基于深度学习的分子性质预测算法的两大类:多层感知机(multi-layer perceptron,MLP)和图神经网络(graph neural network,GNN),及文献中涉及的数据集和评价指标,第2、3节概述了基于MLP和GNN的分子性质预测算法,第4节通过实验复现不同框架下的算法并进行对比,最后讨论和展望分子性质预测领域未来的研究方向。

图1 基于神经网络的分子性质预测算法发展时间轴Fig.1 Development timeline of molecular property prediction algorithm based on neural network

图2 分子性质预测领域关键词检索关系矩阵图Fig.2 Keyword search relation matrix of molecular property prediction domain

1 概述

随着人工智能快速发展,在分子性质预测领域,人们提出了各种不同的算法,从浅层神经网络[18-21],到基于深度学习的最新方法[22-26]。有基于特征描述的[27],将原子所处的环境信息通过人工构建的描述符进行特征编码,该特征多用于基于前馈神经网络的输入;有基于消息传递的[28],将深度学习网络作用于单个原子之间的信息交换,从网络中学习其中的物理化学特征。

1.1 算法分类

当前基于神经网络的分子性质预测算法主要呈现两种模式(图3):一种是基于前馈神经网络的多层感知机(MLP)架构算法,将分子势能看作每个原子势能之和,每个原子的势能通过一个独立的神经网络进行预测,多个原子的神经网络组成神经网络势,最后将所有原子的能量求和得到最终的分子能量;另一种是基于图神经网络(GNN)的分子性质预测算法,将分子看成不同原子之间相互连接的无向图,通过无向图中的节点与节点、节点与边、边与边的消息传递得到对应的原子和分子的势能。本文中将基于深度学习的分子性质预测算法分为基于MLP和GNN的分子性质预测算法,根据其内部的特征提取机制,总结上述两种框架的特点和功能,见表1。

表1 基于MLP和GNN两种框架的分子性质预测算法

图3 基于神经网络的分子性质预测算法分类示意图Fig.3 Classification diagram of molecular property prediction algorithm based on neural network

1.2 数据集

数据集是神经网络模型开发的基础,神经网络模型的成功多取决于数据集的质量和大小。研究者们设计并建立众多的数据集来适应不同的任务,如以数量著称的GDB[29]系列数据集,其建立的初衷是发现与合成药物,它包含了海量的稳定化学分子结构;以溶解度预测为研究目标的ESOL[30]数据集,包含1 128种化合物的水溶性数据;研究分子在水中特性的数据集FreeSolv[31],主要提供小分子在水中水化自由能的相关数据。

QM系列数据集,以GDB系列数据集为基础,选取其中某些分子利用密度泛函理论(density functional theory, DFT)计算得到相应的物化性质,其中具有代表性的是QM9[32]数据集。QM9数据集从GDB-17中1 660亿个分子中提取138 850个分子作为样本,包含C、H、O、N、F 5种元素,其中每个分子都有完整的空间结构,并在笛卡尔坐标系中进行描述。除空间结构信息外。QM9数据集还提供了每个分子的15种理化性质,每条性质均以密度泛函计算得到,即每个内层电子轨道由6个高斯型函数线性组合而成,每个价层电子轨道则会被劈裂成两个基函数,分别由3个和1个高斯型函数线性组合而成。上述15种理化性质分别为三个旋转常数、偶极矩、各向同性极化率、最高占据分子轨道的能量、最高与最低占据分子轨道的能量之差、电子空间范围、零点能、0 K下内能、298.15 K下内能、298.15 K下焓值、298.15 K下吉布斯自由能以及298.15 K下热容。更多的数据集[33-36]及其描述见表2。

表2 常见的分子性质特征数据集

1.3 评价指标

为了评价算法性能,设置合理的评价指标尤为关键。目前基于深度学习的分子模拟算法主要包含两种评价方法:一种是基于误差型的评价指标,一种是基于相关性的评价指标。其中基于误差型的评价指标分为平均绝对误差(mean absolute error,MAE)、均方误差(mean square error,MSE)和均方根误差(root mean square error,RMSE),基于相关性的评价指标为相关系数,又称皮尔逊相关系数(R)。

MAE误差又称L1 loss,计算公式为

(1)

MSE误差计算公式为

(2)

MSE越小代表预测值与标签值越接近,MSE的最小值为0,表示预测值与真实值相等。

RMSE误差又称L2 loss,其计算公式为

(3)

RMSE越小代表预测值与标签值越接近,RMSE的最小值为0,表示预测值与真实值相等。

三种评价指标对于预测值的敏感度不同,其中最敏感的是MSE误差,较为敏感的是RMSE,最不敏感的是MAE误差。

相关系数常用来描述预测值与原数据的拟合程度。常用的相关系数为R,其计算公式为

(4)

2 基于MLP的分子性质预测算法

基于MLP的分子性质预测算法一般流程如图4所示。网络输入通常为原子电荷数和位置信息,不同层级的特征进入不同的神经网络或模块中,模块中的特征通过全连接层进行特征交互,通过激活层进行非线性处理,最终的特征通过求和得到。

图4 基于MLP的分子性质预测算法流程图Fig.4 Flow chart of molecular property prediction algorithm based on MLP

2.1 基于RBF的分子性质预测算法

基于RBF(radial basis function)的分子性质预测算法主要针对径向基函数进行改进。由Hohenberg-Kohn第一定理可知,分子的一切物理化学性质均由该分子的核骨架决定。Behler等[20]将分子骨架信息融入神经网络中,通过神经网络预测分子的势能。该方法首先将分子的结构信息通过分段函数进行编码,经过编码后的信息通过径向对称函数处理后作为神经网络的输入,神经网络预测每个原子所具有的势能,最后通过求和得到分子的总势能。

随着深度学习的不断发展并在其他领域取得巨大成功,基于神经网络势(neural network potential,NNP)逐渐成为分子模拟的主流算法之一。欧式距离作为典型特征虽然可以衡量原子间距,但将角度和空间信息掩盖。Smith等[37]尝试将分子骨架中的二面角引入,将势能的影响因素从原子间的距离扩展为距离和二面角。而Zhang等[24]则将每个原子的位置信息进行多维度编码,包含距离信息、x方向分量、y方向分量、z方向分量组成特征向量,经归一化后输入神经网络中预测分子的势能。但前馈式神经网络所能容纳的特征有限,且随着网络层次的加深,特征丢失严重,针对RBF的分子性质预测算法性能提升有限。

2.2 基于Residual Block的分子性质预测算法

前馈神经网络特征丢失严重,从数学本质上看,是特征在多层网络的前向传递时梯度发生衰减,残差网络则将特征进行跃层连接如图5所示,将上一层特征的梯度传递到下一层,在保证模型表达能力的同时,尽量减少梯度消失。研究者考虑将残差机制引入,同时将普通残差结构扩展为多特征融合残差结构。

图5 SchNet[23]残差机制示意图Fig.5 Schematic diagram of residual block in SchNet

其中具有代表性的是Schütt等[25]针对原子能量的预测任务提出的SchNet网络,该网络构建具有多尺度感知的残差原子交互模块,增加了原子之间的局部相关性。Unke等[38]受SchNet的启发提出PhysNet,该网络在引入残差机制的基础上进一步改进了内部交互残差模块,与SchNet中只在特征采样前后进行融合不同,PhysNet在网络中采取多级残差设计。Lu等[36]建立了一个新的数据集Frag20,并对PhysNet进行了简化得到sPhysNet,利用sPhysNet在新数据集进行测试。结果表明复杂的残差结果在一定程度上可以提高模型的预测精度,但残差网络架构对原子间的内部交互表达不兼容,特征大多仍为单原子的性质融合,对复杂分子性质预测能力不足。

2.3 基于Interaction的分子性质预测算法

不同于将分子性质仅看做原子性质的加权和,研究者认为一些预测效果较差的性质可能受原子间的交互影响,故将原子间的交互机制引入[39]。若将分子看成图结构,则构成分子的每个原子则组成图中的一个结点,原子和原子之间的联系可以通过关系矩阵进行关联,如图6所示,关联度越大该原子的能量在空间受到的影响越大。但在实际训练时很难将所有性质在有效时间步中训练到合适精度,所以有些研究者增加内部交互模块的训练时长。

图6 相邻时间步原子间交互矩阵示意图Fig.6 Schematic diagram of interaction matrix between atoms in adjacent time steps

在此基础上,Schütt等[33]还将原子扩张影响引入,具体是将原子矢量特征与原子间距离进行非线性耦合。经过若干时间步修正后的能量输入至神经网络预测该原子的能量。Yao等[40]将这种原子间的影响直接独立出来,利用独立的网络进行训练,一个电荷网络势负责预测库伦能和范德华能,一个能量网络势负责由原子结构影响的总势能,三者之和即为分子的总势能。Li等[41]引入主要官能团的影响,对于不同类型的体系利用不同的算子进行特征描述,在分子能量的预测中具有良好表现。

3 基于图神经网络(GNN)的分子性质预测算法

图神经网络(graph neural network,GNN)是近年来出现的一种利用深度学习在图结构中进行学习的框架,其优异的性能引起了学者高度的关注和深入的探索。通过在图中的节点和边上制定一定的策略,GNN 将图结构数据转化为相应的特征表示,在多种任务中取得优良的效果。基于图神经网络的分子性质预测算法一般流程为模型通过构建好的输入特征,通过多层图卷积进行消息传递,如图7所示,基于图神经网络的分子模拟算法其主要特征是将分子结构看成图结构,并利用图神经网络的结构不变性进行预测。

图7 GNN算法基本流程图Fig.7 Basic flow chart of algorithm

3.1 基于MP的分子性质预测算法

消息传递本来指代计算机间用于实现同步的通信机制,在分子性质预测中,研究者将不同层级的图中包含的信息更新机制称为消息传递。组成分子的不同原子构成图或者超图结构,图的节点包含描述结点本身特征信息,图的边则包含两个节点的连接关系特征,消息传递机制规定了下一层图中的节点和边的更新规则。

基于MP的分子模拟算法主要针对消息传递机制进行改进。Gilmer等[28]在消息传递阶段分别借鉴GG-NN[42]和Interaction Networks[43]的消息函数,对于消息传递距离的改进,消息传递神经网络(message passing neural networks,MPNNs)提出两种方法,一种是将未连接的节点增加虚拟边,另一种是增加与所有节点连接的全局节点,这样的设计类似卷积神经网络中的增加感受野。Lu等[44]在改进消息传递机制的基础上,考虑不同数量级原子间相互作用,即若干2个原子间的相互作用Pair-wise、若干3个原子间的相互作用Triple-wise等。仅针对分子结构进行编码从某种程度上忽略了原子在空间上的排布信息,一方面可以将空间方向信息融入图结构中[45],另一方面可以建立3D图结构[46],在一定程度上增加了模型的精度。但基于消息传递机制的图神经网络存在特征传递瓶颈问题,有些研究者通过增加注意力机制来克服[47],有些则通过等变机制来解决[48],但效果提升有限。

3.2 基于Equivariant的分子性质预测算法

分子中原子结构具有旋转不变性、平移不变性和镜像反转不变性等几何性质。其中,具有旋转不变性的称为SO(3),除SO(3)外具有平移不变性的称为SE(3),除SE(3)外具有镜像反转不变性的称为E(3),如图8所示。基于Equivariant的分子性质预测算法主要将分子间的位置关系看成旋转等变的特征描述。在特征不断传递更新时,保持特征旋转等变,这和卷积神经网络有所区别,卷积神经网络不能保持特征的旋转等变。

图8 三种不变性分子示意图Fig.8 Schematic diagram of three invariant molecules

大多旋转的等变网络考虑原子轨道间的电子相互作用,一方面引入物理运动算子[49],学习原子间的势和力。利用来自可训练的潜在力向量方向信息,以及受牛顿物理学启发的物理注入算子的优势,整个模型保持旋转等变,并通过更可解释的物理特征推断出多体相互作用。Qiao等[50]利用有效的紧束缚模拟和学习映射预测量子化学性质。Schütt等[51]提出了极化原子相互作用神经网络PaiNN。

信息传递神经网络已成为图形学习的首选方法,特别是在预测化学性质和加速分子动力学研究方面。虽然它们很容易扩展到大型训练数据集,但以前的方法已被证明不如内核方法有效。不变表示的局限性是一个主要原因,并将消息传递公式扩展到旋转等变表示。Qiao等[26]提出了一般多体等变神经网络——UNiTE,对于N-body的张量,一方面通过对角化进行简化,另一方面通过线性映射,将两部分的特征进行高阶卷积后,再进行消息传递,经过等变归一化使用对称消息池读出预测结果。Takamoto等[52]提出了通用神经网络势PFP(preferred potential),该网络面向更广阔的应用空间,为了这种普适性,PFP提供了复杂的数据集生成模块,并在LiFeSO4F中的锂扩散、金属有机框架中的分子吸附、Cu-Au合金的有序-无序转变,以及费托催化剂的材料发现等方面进行有效性验证。

3.3 基于李群的分子性质预测算法

李群指光滑可微的群,通常可以认为综合了群和光滑流形的概念。光滑、可微的流形,指领域和欧式空间同构的线性空间,也就是李群的每个元素存在线性空间或者向量空间作为切空间。由于单位元在群中的特殊地位,所以李群单位元的切空间李代数是非常重要的概念。由于李群的光滑性,每个元素的切空间结构相同,都可以通过线性变换变化至单位元的切空间,也就是李代数。

基于李群的分子性质预测算法将分子间的作用看成对称平移等变问题[53],通过过向旋转等变网络中增加角度信息,利用角度信息使模型预测准确度提升更明显。或者增加对称适应原子轨道(symmetry-adapted atomic orbitals,SAAO)的分子特征[50],通过对该矩阵进行对角化处理得到具有旋转不变性的自适应原子轨道基,如图9所示。将上述特征映射到对应的图结构中,经过多层的消息传递层和特征编码最终通过解码得到最终的势能预测值。

图9 三维空间上的等变特征示意图Fig.9 Isovariant feature diagram in 3D space

除了在三维空间中实现特征的旋转、平移、反射和置换不变性,更高维空间的等变特征也受到了广泛的关注[54],甚至重新定义网络中的乘法[55]和卷积[56]。一方面扩展特征维度[23],另一方面增加物理化学性质[57],使网络特征更具解释性。上述算法均未考虑特征的异向传递,虽然在预测精度上已经具有较大的优势,但忽略非对称结构在不同空间方向具有不同的理化性质。

其他新颖的改进方法则受其他深度学习任务启发,比如:Liu等[58]将注意力机制引入到消息传递过程中,提出DeepMoleNet网络,该网络通过加权不同原子的贡献,将以原子为中心的对称函数(atom-centered symmetry functions,ACSF)作为teacher描述符,而不是以传统方式使用ACSF作为输入,使化学可解释的见解能够融合到多任务学习中。Godwin等[59]提出了NoisyNodes策略,该策略通过简单的噪声正则化解决过度平滑问题。NoisyNodes用噪声破坏输入图,并添加一个噪声校正节点级损失。多样的节点级损失鼓励潜在的节点多样性,噪声节点可以作为GNN中的补充构建模块。

4 实验分析

近年来,基于深度学习的分子性质预测算法呈现爆发式增长,基于GNN的深度学习分析模拟算法受到更广泛的关注。为了对当前基于深度学习的分子模拟算法进行比较,选取0 K下的内能U0、298.15 K下的内能U、298.15 K下焓值H、298.15 K下吉布斯自由能G的预测结果进行对比。

根据目前公开的结果,如表3所示,在预测精度上, UNiTE[26]的预测精度最高,在算法大类对比中,如表4所示基于MLP的算法精度低于基于图神经网络的算法精度。从图1所示趋势图不难看出,基于MLP的分子性质预测算法在2018年后逐渐被GNN相关算法所取代,基于GNN分子性质预测算法已逐渐成为研究热点。

表3 几个算法的公开代码链接

表4 不同分子性质预测算法结果对比表

本文中对相关算法进行调研,将每个算法的代码公开情况,汇总形成表5。在几个表现较突出的算法中,基于深度学习相关算法预测时间要远小于密度泛函理论相关算法,如图10(a)所示,其中几个算法在准确度上可以媲美DFT相关算法。与机器学习相关算法相比,准确度上有较大的优势。

表5 大类对比表

图10 实验结果图Fig.10 Diagram of experimental result

本文中还对几个模型的模型特征表达能力进行对比,如图10(b)随着训练集的增加,模型的测试误差降低,不难看出模型的特征表达能力仍然有较大的扩展空间,随着新型数据集的开发,基于深度学习的分子性质预测算法将更具优势。

5 总结与展望

概述了近些年来基于神经网络的分子性质预测算法,分为MLP和GNN两大类进行分析,结合公开数据集对算法进行对比验证,发现在精度上基于GNN的分子性质预测算法要高于基于MLP的分子性质预测算法,并且基于GNN的分子性质预测算法在近三年受到更广泛的关注。基于MLP的分子性质预测算法将分子性质看作若干原子分子性质之和,每个原子的相关信息通过独立的神经网络开展预测任务,相关联的分子之间通过残差结构进行特征融合,得到最终的预测结果。基于GNN的分子性质预测算法将分子结构看成图结构,原子与原子之间的关系映射为无向图之间的边,原子映射为图的节点。原子的性质通过消息传递机制,传递给下一层的图结构,典型的GNN分子性质预测算法存在节点到节点的消息传递,边到节点的消息传递,全局节点到各节点的消息传递,每个节点与上一层的若干节点相关联,多层消息传递后,特征将像图像中的卷积中的特征张量具有一定的感受野,最终通过读出层读出最终的预测结果。

基于以上的调研,未来可以开展的工作包含以下几个方面。

(1)超大规模数据集的构建。目前的算法在特征表达上仍未到达瓶颈,大规模和超大规模的数据集将有利于当前的算法形成预测更为精确的模型,使算法在可接受的时间复杂度内达到可以与第一性原理计算相媲美的预测结果。此外,模型的损失函数与模型的评价指标高度契合,即模型训练用的损失函数大多由评价指标或评价指标的变体构成,这导致模型在一定程度上过度向评价指标方向倾斜从而忽略其他关键因素。

(2)图神经网络的各向异性特征传递。当前的图神经网络在预测精度上已经具有较大的优势,当前的图神经网络较为依赖等变的特征传递而忽略非对称结构在不同空间方向具有不同的理化性质,该问题将制约图神经网络在非对称性结构中的性质预测。此外,目前的模型针对过拟合问题缺乏系统的优化,结合图神经网络特征,一方面考虑增加池化、噪声等结构,另一方面设计更适合的损失函数。

(3)材料科学与生命科学中的实际应用。在材料科学中,基于势能面搜索的材料结构预测,目前的算法更关注搜索的时间和空间复杂度,而不是材料的性质,在基于势能面搜索的材料性质预测仍然是需要解决的问题。此外,目前的算法针对材料结构预测往往基于简单体系,当考虑材料的强相互作用、磁性时,仍不能起到实际指导作用。在生命科学中,目前相关研究只能针对小蛋白质结构进行预测,目前距离真正理解蛋白质结构与功能的关系并在此基础上开展具有特定结构或功能蛋白质的设计仍有较大差距。