三维定量构效关系

2D-QSAR方法采用的参数，通常和分子的3D结构没有很大关系，所以也并不需要药物分子的活性构象，但是这也限制了2D-QSAR的准确性。通过分析药物分子甚至受体的3D结构，3D-QSAR方法可以提供更加合理的模型。目前世界上最主流的3D-QSAR方法是Co MFA（比较分子场分析）法和在其基础上发展的Co MSIA。

6.3.1 距离几何法

距离几何法（distance geometry approach，DGA），或称距离映射法（distance map，DM）是由加利福尼亚大学的Gordon M.Crippen在1979年提出的一种3D-QSAR方法，也是最早的成熟3D-QSAR方法。距离几何法通过计算受体的结合位点和药物分子的作用位点之间的能量参数预测化合物的结合能和生物活性。

距离几何法首先要找到配体分子中可能的作用位点。图6-3所示的是Crip-pen在1979年的文章中对一种糜蛋白酶抑制剂，m，m-二甲基苯氧基丙酮的作用位点的分析。深色的圆点表示配体作用位点。位点1位于羰基C原子上，位点2位于醚氧上，位点3和4位于甲基上，位点5位于苯环的中心。

图6-3 糜蛋白酶抑制剂m，m-二甲基苯氧基丙酮配体作用位点分解（引自Crippen）

随后计算配体分子各个作用位点的距离矩阵，如表6-1所示。不同的下标表示不同的配体作用位点类型。在本例中，由于位点4和5属同一类型，因此下标数字为1～4。

需要注意的是，表6-1并不是完全对称的，上三角阵和下三角阵中的数字大小并不完全相同。这是因为位点1在分子中的相对位置有较大的灵活性，它与位点3～5之间的距离并不固定。上三角阵中的数字是两个位点间距离的上限，而下三角阵中的数字表示两个位点间距离的下限。

表6-1 m，m-二甲基苯氧基丙酮配体作用位点的距离矩阵（引自Crippen）（单位：Å）

接下来要定义受体的结合位点，并计算出受体结合位点间的距离矩阵。与配体的结合位点矩阵不同，受体结合位点的矩阵是一个对称阵，因为在距离几何法的计算过程中，受体基本是保持固定的。受体结合位点的状态分为空（empty，配体作用位点可进入）和满（filled，配体作用位点不可进入）两种。例如，一个疏水的口袋当中就包含一个空位点，当配体和受体结合时，配体上的苯环作为一个作用位点，可以占据这一空位点。

下一步需要定义不同类型的配体分子作用位点和受体分子结合位点的结合能。表6-2与表6-3是交互生成的，通过拟合的方式得到。

表6-2 受体作用位点的距离矩阵（引自Crippen）（单位：Å）

表6-3 配体作用位点和受体结合位点的结合能（引自Crippen）（单位：kca I）

对于每一种结合方式，按表6-3将每一对作用位点和结合位点的结合能加和起来，便得到整个配体与受体的结合能。每一对作用位点和结合位点之间的结合为“全或无”的方式，只要结合，便是完全结合，不存在部分结合的情况。

如果对于已有的所有分子，通过距离几何模型预测的结合能都与实验值相仿（见表6-4），那么表6-2和表6-3中的各个数字就确定了，距离几何法的模型也就随之确定了。需要预测新分子的活性时，只要通过表6-3中的各种结合能参数找出最好的结合模式，算出结合能，就能得到新分子的大致活性。

表6-4 糜蛋白酶的预测值与实验值比较（引自Crippen）

注：p表示对位取代，m表示间位取代。

图6-4 二氢叶酸还原酶与其某一抑制剂结合的立体结构编号1～11的大球表示受体的结合位点（引自Crippen）

距离几何法将配体作用位点和受体结合位点都定义为距离矩阵有一个明显的好处，当分子平移或旋转时，它们的距离矩阵不会发生任何变化，这样就不必费心计算其平移和旋转过程了。距离几何法最后只是简单地给出一张表，列出哪些作用位点结合了哪些空结合位点。

Crippen最初将距离几何法用于酶蛋白酶抑制剂和DHFR（二氢叶酸还原酶）抑制剂的构效关系研究，效果良好。但是相比比较分子场等主流方法，距离几何法在配体作用位点的选取方面主观性较强，且计算操作繁复，因此近年来采用距离几何法进行3D-QSAR分析的文章数量很少，且绝大部分出自Crippen本人。

6.3.2 分子形状分析

分子形状分析（molecular shap eanalysis，MSA）方法是由A.J.Hopfinger在1980年研究DHFR的抑制剂的定量构效关系时提出的。

分子形状分析方法的第一步首先要确定配体分子的活性构象。假如受体的晶体结构已知，则其配体分子的活性构象可以用分子对接的方法得到。

接下来要选择一个分子（如药物活性最好的分子）作为参考结构，随后其他分子分别与该参考构象叠合，求出各个分子的分子形状参数。最早使用的分子形状参数是公共重叠体积（common overlap stericvalume）V0（i，j），即各个分子与作为参考结构的分子的重叠部分的总体积：

式中：（vi，k∩vj，l*）表示分子i中的k原子与分子j中的l*原子的重叠体积，如图6-5所示。

图6-5 分子i的k原子与分子j的l*原子重叠体积的2D示意图（引自A.J.Hopfinger）

除了V0（i，j）之外，Hopfinger还设计了一些其他的分子形状参数，例如：

S0（i，j）=[V0（i，j）]2/3

L0（i，j）=[V0（i，j）]1/3

这些参数称为第一代分子形状参数。

1983年，Hopfinger进一步改进了分子形状法，提出了第二代分子形状参数。首先是分子势能场公式：

式中：u是分子序号；R，θ，φ是空间坐标（见图6-6）；n是分子u中的原子个数；ri是探针原子或分子和药物分子中第i个原子的距离；a（T）i是探针类型相关的参数；b（T）i是药物分子中原子类型的参数；Qi是第i个原子的电荷；Q（T）是探针电荷；ε（ri）是介电常数（在Hopfinger的原文中，ε（ri）=3.5）。

图6-6 探针分子相对于分子u的坐标R，θ，φ（自A.J.Hopfinger）

根据分子势能场，Hopfinger构造了几个分子参数：

在最后建立线性模型时，采取这些分子参数的对数形式，例如Hopfinger的DHFR抑制剂的QSAR模型：

lg（1/C）=-2.34F+0.29F2+0.37（π3+π4）+9.39 　（6-15）

上面的模型中，F=lg[ΔPα，β（12，2）1/2]，探针为氢离子；C是半数抑制摩尔浓度；π3+π4表示DHFR抑制剂5-苯甲基-2，4-二氨基嘧啶上苯环的3，4位取代疏水常数之和。Hopfinger认为这一模型（r=0.961，s=0.105）比用第一代参数V0建立的模型（r=0.931，s=0.137）更好。

6.3.3 虚拟受体法

虚拟受体法是3D-QSAR中较晚出现的分支，是一系列3D-QSAR方法的统称，其共同点是采用某些探针原子或基团在分子表面建立虚拟受体环境。包括Compass法、GERM法、RSM法等。

6.3.3.1 Compass

Compass方法是Jain在1994年提出的3D-QSAR方法，其流程如图6-7所示。

图6-7 Compass方法的过程

整个过程分成3个步骤：①生成每个分子的多种构象，将每种构象的每一种取向定义为分子的一种姿势（pose），通过药效团模型初步叠合分子；②从分子叠合结果的质心出发，向三维空间投射162条射线，每条射线各分子平均表面外一定距离处作为一个样本点（sampling point）[见图6-8（a）]，对于每个分子，其范德华表面到各个样本点的距离[见图6-8（b）]组成一个矢量，用于表征该分子的形状，或者，也可以通过各个分子上离样本点最近的氢键给体或受体基团到各个样本点的距离来定义分子的形状。得到这些参数之后，通过三层BP网络（神经网络方面的内容稍后介绍）进行模型训练，其输入是该分子的形状表征矢量，即162维的距离矢量或氢键给体受体矢量，输出是某种值处于0～1之间的生物活性数据。神经网络输入层也就是162个神经元，中间层可取3个神经元。神经网络训练若干步之后，重新调整分子姿势，姿势调整与神经网络训练交替进行，直至网络收敛；③神经网络收敛之后，将新分子代入模型，即可预测其生物活性及可能的活性构象。

图6-8 （a）样本点的确定；（b）距离矢量计算（引自Jain）

6.3.3.2 GERM

1994年，Walters提出了GERM方法（genetically evolvedreceptor models），该方法因其采用遗传算法得到受体活性位点的原子水平模型而得名。

首先，在配体分子表面的三维空间放置若干虚拟受体原子，这些原子可以模拟配体表面并与配体相互作用。这些虚拟受体原子依据其原子类型写成一个线性字符串，如图6-9所示。原子类型与字符串数字的对应关系如表6-5所示。由于配体结构预先并不知道，这些虚拟受体原子的初始数量和种类都是任意的，要通过计算配体结合能来选取最佳组合。而原子的不同种类的组合几乎是无限的（一个由60个虚拟受体原子，8种原子类型组成的模型有超过1054种组合），因此不可能采取枚举的方法，而要使用遗传算法（稍后介绍）完成筛选。遗传算法当中的适应性函数是计算出虚拟原子和所有配体分子的相互作用能的负对数相对于生物活性数据对数的回归系数（见图6-10）。

图6-9 虚拟受体原子的生成（引自Walters）

表6-5 虚拟原子类型及其编码（引自Wa Iters）

6.3.3.3 受体表面模型方法

受体表面模型方法（receptor surface model，RSM）是由MSI公司（Molecular Simulations Incorporated）的Mathew Hahn提出的。RSM方法首先通过一个或多个经叠合的结构（通常是样本中活性较高的那些分子），作出虚拟受体表面，这个虚拟表面代表了这些分子形状的叠加。叠合的正确与否，对于模型的正确性及其泛化能力都有很大的影响。

虚拟受体表面可以是“开放”（open）的或“闭合”（closed）的。开放的虚拟表面有一些孔洞，用于代表溶剂可及区域或未知区域；闭合的虚拟表面则完全覆盖表面所有区域。实际上，受体表面模型未必是连续的，也可以由已知区域的许多小片段组成。

图6-10 适应性系数的求算（引自Walters）

可以用范德华场函数或Wyvill场函数定义受体表面形状。范德华场函数：

V（r）=r-r VDW 　（6-16）

式中：r为某个网格点到某个原子的距离；r VDW是这个原子的范德华半径。对于位于范德华表面上的点，V（r）=0；范德华体积内部的点，V（r）＜0；范德华体积外的点，V（r）＞0。对所有原子周围的网格点均计算V（r），若某一网格点同时与多个原子相邻，则在其V（r）中选取最小值。范德华场函数无边界，因此要定义一个截断距离，对于每个原子，只计算一定半径范围之内的场函数值。

Wyill场函数：

与范德华场函数不同，在有限的距离R内，Wyvill函数有边界并完全衰减。显然，V（0）=1，V（R）=1，V（R/2）=1/2，这样表面得以平缓连续地展开（见图6-11）。场函数值是所有原子场值贡献的加和。R是场值衰减为0时的距离，对于距离超过R的点，其场值一律不予计算。在Hahn的原作中，R定义为原子范德华半径的2倍，即不同类型原子有不同的R。

图6-11 范德华场函数（a）和Wyvill场函数（b）生成的表面模型（引自Hahn）

在受体表面上的每一个格点上，储存着部分电荷、静电式、形成氢键倾向和疏水性的信息。将这些信息对应的标量值算出受体与配体之间的相互作用能Einter，再用Einter与生物活性数据进行回归。

6.3.4 比较分子场分析法和比较分子相似指数分析法

1988年，Cramer提出了比较分子场法（comparative molecular field anlysis， Co MFA）。Co MFA方法是目前最为通用的3D-QSAR方法。

Cramer方法首先需要药物分子的活性构象。如果受体晶体结构已知，则可以通过分别作分子对接得到这些分子的活性构象；如果受体晶体结构未知，那么一种办法是进行药效团叠合，或者也可以通过同源建模模拟受体的结构再进行对接。有时对接后药物分子的构象和位置差别较大，不适合作叠合，可以将生物活性最好的分子的最低能量构象作为模版，其他分子的活性构象以该模版为标准选取。

活性构象得到以后，用这些构象作药物分子的叠合，叠合通常需要挑选一个分子作为模版，可以选择生物活性数据最好的分子作模板。

分子叠合完成之后，将其所在的空间按一定步长划分成若干个格点，在每个格点中放置一个探针原子。探针原子的种类可自选，如碳正离子、水分子等，但是也受到场势能函数的限制（如静电场和立体场通常用碳正离子）。随后分别计算每个探针离子和各个分子间的相互作用能，用这个能量作为PLS的参数，计算出最佳主成分。这个过程如图6-12所示。

图6-12 Co MFA的基本流程（引自Cramer）

图6-12所示的是立体场和静电场的相互作用能，这两者经常被综合考虑，不过如有需要也可以分别建立PLS模型，此外一般还有氢键场和疏水场可选择。

Co MFA方法可以通过比较直观的方法指导新的药物分子的改造工作。包含Co MFA模块的软件如SYBYL等，能够绘制出等势面图，并且用不同颜色的区域表示此处引入相应基团对活性的可能影响，如图6-13所示。

例如，黄色区域可能表示在此处引入疏水性基团可以提高分子的活性，红色区域可能表示此处引入大体积基团会降低分子活性等。如今，借助SYBYL等软件，可以作出更清晰的等势面图。

图6-14中，绿色（G）区域表示大体积基团提高分子活性；黄色（Y）表示大体积基团降低分子活性；蓝色（B）区域表示正电基团提高分子活性；红色（R）区域表示负电基团提高分子活性。

图6-13 Co MFA等势面图（引自Cramer）

图6-14 醌类似物的Co MFA立体场和静电场等势图（引自Haifeng Chen等）

1994年，Klebe等在Co MFA的基础上提出了比较分子相似指数分析（com-parativemolecularsimilarityindicesanalysis，Co MSIA）方法。Co MSIA方法中的分子场（立体场、静电场、疏水场、氢键给体场、氢键受体场）由下式算出：

式中：j为分子序号；i为分子j中的原子序号；wik为原子i中的第k类理化参数值；wprobe，k为电荷+1、半径1Å、疏水性+1的探针原子的第k类理化参数值；α为衰减因子；riq则是格点q中的探针原子和原子i间的距离。

Co MSIA中用Gauss函数代替了Co MFA当中的Coulomb和Lennard-Jones函数。Co MFA的势函数需要定义截断值并可能在某些区域急剧变化，而Co MSIA则可以避免这些问题。格点的大小、位置发生改变后Co MFA模型可能产生很大变化，而Co MSIA则稳定得多（见图6-15）。

图6-15 Co MFA和Co MSIA对格点位置变化的敏感性（引自Klebe等）

由于最常用的Co MFA工具SYBYL同时也包含Co MSIA模块，因此这两种方法经常会被同时采用，比较分析。