【黄子按曰】大约15年前，我根据http://www.netsci.org上的一个网页，翻译形成了以下文字。最近，准备给背景是化学工程却因缘聚会到了我这里攻读生物信息学博士学位的留学生启蒙一下QSAR，却发现该网站及相应页面都已经不存在了。一声叹息！这里就只能给出我的中文译文了，不当之处敬请指正。

【黄子再曰】2015年4月11日，南开大学药学院药物化学CADD方向的殷璨同学云游至此，小试牛刀，找到了该文英文原文的PDF文件并邮件惠赠于我（An Introduction to QSAR Methodology.pdf）。在此谢过了！

定量构效关系方法学入门

 李艾伦
 29451 美国南卡罗来纳州棕榈岛卡罗来纳大街412号 网络科学公司
E-mail: netsci@netsci.org
杨坦利
27709 美国北卡罗来纳州研三角区五穆路葛兰素威康研究所
1997年10月发布

（黄健翻译） 

绪言

药物设计是一个反复的过程，始于一个具有某种生物学特性的化合物，止于该分子的化学合成与生物学活性的优化。分子或系列分子的化学特性与其生物学活性相关，当化学家形成这样一种假说时，药物设计的过程就开始了。我们无需详尽了解使药物具有活性的生化过程，只要比较活性和非活性分子结构的异同，通常就会使上述假说更加精确。然后选择具有最多活性相关功能基团或特性的化合物进行合成。

即使对简单的体系，组合策略的可能性也是呈爆炸性增多的。例如，为在不对称双取代苯环体系的4个开放位置放入10个取代基团，需要合成的化合物数目约为10,000个。替代这一实验密集型化合物优化的办法，是发展一种可以将生物学活性的差异与每一化合物容易获得的分子特性描述符(descriptors)的变化定量关联的理论。一个定量构效关系（Quantitative Structure Activity Relationship, QSAR）可帮助指导化合物合成。本文将逐步阐明用于推导QSAR的概念并回顾QSAR技术在医学研究中的应用。 

统计学概念

计算化学将分子结构描述成数字模型并用量子物理学和经典物理学方程模拟其行为。现有程序使科学家能够轻松生成和显示分子数据，包括几何结构、能量和相关属性（电属性、光谱属性和体积）。显示与操纵这些数据的通常范例为表格，其中化合物被定义为行而分子属性（描述符）被定义为相关的列。QSAR就是试图找出一系列分子的属性值与生物学活性变化的一致性，从而可以将这些“规则”用于评估新的化学实体。

QSAR通常为线性方程形式 Biological Activity = Const + (C1.P1)+(C2.P2)+(C3.P3)+…§ 

其中，参数P1到Pn可从该系列的每一分子计算，系数C1到Cn可通过协配参数与生物学活性的变化算出。因为上述方程通常通过统计学方法得出，所以下面将简要介绍推导QSAR时应遵循的原则。

一个简单QSAR例子是山德士(Sandoz)医学研究所在开发全新镇痛药物方面所报道的工作。该研究制备了与辣椒素(capsaicin)有关的香荚兰胺(vanillylamide)和香草基硫脲(vanillylthiourea)，并通过体外方法测定流入背根节神经元(dorsal root ganglia neuron)的45Ca2+来确定它们的活性。以EC50(μM)形式报道的数据见表1（注意，化合物6f在该系列中活性最高）。

表1 辣椒素类似物活性数据

Cmpd Number	Cmpd Name	X	EC50(μM)
1	6a	H	11.80 ± 1.90
2	6b	Cl	1.24 ± 0.11
3	6d	NO2	4.58 ± 0.29
4	6e	CN	26.50 ± 5.87
5	6f	C6H5	0.24 ± 0.30
6	6g	N(CH3)2	4.39 ± 0.67
7	6h	I	0.35 ± 0.05
8	6i	NHCHO	???

 

在没有更多信息的情况下，估计化合物6i活性的最好办法与唯一办法就是计算该系列化合物已知活性的均值。对化合物8活性的一个猜测是均值7.24，但这个猜测好吗？图1为数据的散点图表示。

数据的标准差 s 显示活性值分散在均值周围有多远。通过显示数据固有变异性的总和，s 值成为“猜测”质量的一种指标。标准差的计算如下：

与其依赖于上述有限的分析，不如找出影响该系列化合物活性的因素，并用以预测新化合物的活性。为达到上述目标，我们需要：

l 能够区分不同化合物的足够精确的测定的结合数据；

l 容易得出的一套参数并可能与受体的亲和力有关；

l 推导参数与结合数据关系的方法（QSAR）；

l 证实QSAR的方法；

QSAR方程是一种线性模型，它将一系列分子中计算出来的或实验测定的属性值变化与其生物活性的变化关联起来。要QSAR方法有效，实验选用的化合物（训练集）在化学上应具有多样性。在许多合成工作中，制备的化合物在结构上与先导结构相似。毫不奇怪，该系列化合物的活性值范围也常常有限。在这些情况下，必须额外添加化合物并在测试中填满训练集。

任何QSAR的质量将只与用于推导这个模型的数据的质量一致。剂量反应曲线应该平滑，包含足够点数以确保准确性，并应涵盖2个或更多数量级。对一特定观察实验的多次测试应可重现并且误差相对较小。我们所强调的问题就是信号-噪声比。同一化合物重复测试结果的变异应比该系列的变异度小得多。在生物学实验数据不符合上述原则时，从这些数据推导的QSAR模型会不可靠，所以应利用其它方法进行数据分析。

获得生物学数据后经常会发现，这些数据的表示方式不能用于QSAR分析。因为QSAR基于自由能与平衡常数的关系，所以用于QSAR研究的数据一定要按照发生在生物学应答过程中自由能变化的方式表式。当检测一种药物的药效（产生生物学效应所需剂量）时，自由能的变化可相应地由化合物浓度的反对数计算。

其次，因为生物学数据通常有统计学倾斜，所以常采用对数转换将数据变成接近于正态分布。因此，当在平衡条件下测定反应时，最常用的转换是将浓度值（如IC50、EC50等）表示为log[C]或log1/[C]。辣椒素激动药转换后的数据见表2。
表2 辣椒素激动药转换后的数据

Cmpd Number	Cmpd Name	X	EC50	Log EC50	Log 1/EC50
1	6a	H	11.80 ± 1.90	1.07	-1.07
2	6b	Cl	1.24 ± 0.11	0.09	-0.09
3	6d	NO2	4.58 ± 0.29	0.66	-0.66
4	6e	CN	26.50 ± 5.87	1.42	-1.42
5	6f	C6H5	0.24 ± 0.30	- 0.62	0.62
6	6g	N(CH3)2	4.39 ± 0.67	0.64	-0.64
7	6h	I	0.35 ± 0.05	- 0.46	0.46
8	6i	NHCHO	?? ± ??	??	??

转换后数据相应于其均值的分布效果见图2。注意，投射在Y轴上的数据点分布已变得更加均匀了。

 

就转换后的数据而言，我们对6i活性的最佳猜测仍是该数据集的均值（0.40）。和前面一样，这种猜测的误差可由偏差平方均数的平方根计算。

 

 上面只是一个用于显示常用方法的例子；实际的数据集会有更多化合物和更多描述符。由于QSAR的目的就是突出结构特性与活性的关系，所以我们要找到一个或更多能将与这些分子与其活性相关联的结构特性。此外，我们还要找出能与系列中所有分子都能协调一致的参数。  

在QSAR研究中，有几类可能的参数。取代常数和其它物理-化学参数（例如Hammettσ常数）反映一个基团对分子的电效应。片段计数用于列举出现的特异性亚结构。其他参数还包括拓扑描述符及量子化学的计算值。 

在任何QSAR研究中，参数的选择都是重要的第一步。如果所选参数与活性相关性强，那么就有可能用于活性预测。如果相关性弱，就是知道参数的值也无助于活性的预测。因此，对一特定研究，应选择与所研究系列分子活性有关的参数，而且应以同样的方式得出这些参数的值。

 

山德士研究组将他们分析的辣椒素类似物分为3个区域：A区为一个芳香环；B区为一个酰胺键（amide bond）；C区为一个疏水侧链（见表1的图示）。对C区的假设是认为小而疏水的取代基团可能增加活性。根据这一假设，选择能最好地定义该特征的参数，它们是摩尔折射率(molar refractivity, MR)和疏水取代常数π(hydrophobic substituent constant)。表3中给出了它们的值。 

表3 辣椒素类似物参数值

Cmpd Number	Cmpd Name	X	Log EC50	π	MR
1	6a	H	1.07	0.00	1.03
2	6b	Cl	0.09	0.71	6.03
3	6d	NO2	0.66	- 0.28	7.36
4	6e	CN	1.42	- 0.57	6.33
5	6f	C6H5	- 0.62	1.96	25.36
6	6g	N(CH3)2	0.64	0.18	15.55
7	6h	I	- 0.46	1.12	13.94
8	6i	NHCHO	??	??	??

 对上述数据有两种分析方法：图形方式与统计方式。对于变量有限的问题，最常用的可视化解决方法是图形。对于本例，一幅活性对摩尔折射率或疏水性的图可增加我们对参数与活性关系的认识。山德士研究组推演出的图形见图3。

图3 辣椒素类似物参数值（黄子注曰：图3中2个X轴均为Log EC50，大谬！但原文如此。俺不得其解！）*

 上图对化合物6i的活性有所预示吗？知道6i的疏水性或摩尔折射率能够较好地估计出其活性吗？

这是一个仅涉及两个值的好例子，对上述问题的回答是有条件的肯定。在更复杂的情况下，多参数与活性相关。这就需要运用统计学方法推导出活性与参数集关系的方程。描述上述数据最好的线性方程为

Log EC50=0.764-0.817*π

我们对这一模型有多大信心？回答这一问题的第一步就是确定该方程预测系列中已知化合物活性的吻合程度。根据π值，上述方程可计算出EC50。由于方法不同，实测值与回归估计值会有所不同。实测值与计算值的差异叫做模型残差(residual)。活性的计算值及其残差（或估计误差）见表4。

表4 辣椒素类似物的计算值

Cmpd Number	Cmpd Name	X	Log EC50	π	Calculated Log EC50	Residual
1	6a	H	1.07	0.00	0.79	0.28
2	6b	Cl	0.09	0.71	0.21	- 0.12
3	6d	NO2	0.66	- 0.28	1.02	- 0.36
4	6e	CN	1.42	- 0.57	1.26	0.16
5	6f	C6H5	- 0.62	1.96	- 0.81	0.19
6	6g	N(CH3)2	0.64	0.18	0.65	- 0.01
7	6h	I	- 0.46	1.12	- 0.12	- 0.34
8	6i	NHCHO	??	- 0.98	1.60	??

残差是定量描述通过数据集回归方程估计个例值时误差的一种方法。残差标准误的计算就是求残差的均方根（该计算中分母减去2以反映对两个参数的估计）。

对一个改良的模型，从模型计算出的残差的标准差应小于原始数据的标准差。此前计算的均值的标准误为0.76，而QSAR模型的标准误为0.28。显然，运用线性回归可以提高我们分析的准确性。图4显示，测量值与计算值的关系呈一45度直线。

从一系列化合物中推导QSAR模型时，存在几个固有的假定。第一个假定是，计算（在某些情况下为测定）参数比测定活性更准确更经济。第二个假定是，离最佳拟合直线的偏差呈正态（高斯）分布。最后一个假定是，QSAR方程所描述直线的任何变化与活性及参数的数量级均无关。在这些假定下，我们可以运用多种方法评估模型的质量。

相关系数r估量观测数据与拟合回归直线的吻合程度，可以定量地描述数据的变异。模型或观测数据的误差都会导致点线不合。回归直线的吻合程度的指标如下计算：

其中回归方差为原始方差减去线周方差（黄子注曰：即残差的平方和）。原始方差为原始数据与均值离差的平方和。这在图5中有形象的表示。

上述计算如下：

原始方差 =（1.07-0.40）2 +（0.09-0.40）2 + …= 3.49

线周方差 = 0.282  +（-0.12）2  +（-0.36）2  + … = 0.40

回归方差 = 原始方差 — 线周方差 = 3.49 – 0.40 = 3.09

r2 = 回归方差 / 原始方差 = 3.09/3.49 = 0.89

图5：辣椒素类似物r2值的推导

可能的r2值范围在0到1之间。r2值为0，说明活性与研究所选参数无关。r2值为1，说明完全相关。对上面r2值为0.89的诠释是，辣椒素类似物LogEC50值89%的变化可通过疏水参数π的变化解释。

尽管数据点与回归直线非常吻合，但怎么才能确定这种相关性是否是偶然形成的呢？r2值越大，机会造成相关性的可能性越小。如果回归方程中有较多说明变量，那么，因拟合过程中的弹性而有可能与数据较好吻合。一条直线可以完全吻合两点，一条二次曲线可以吻合3点。如果有足够多的解释性变量，那么多元线性回归可以拟合观测数据2。假定数据呈高斯分布，下面的统计量F可以评估回归方程的统计学显著性。

通过r2和数据集中数据点的数目（或自由度）可以计算统计量F。辣椒素类似物的比值F计算如下：

该值常作为统计软件的标准输出，或可通过查阅统计表而确定回归方程的显著性。本例中，活性与π值无关的可能性不到1%（P=0.01）。

我们已经发现，生物学活性与疏水性有较好的相关性。加入一个反映分子大小的参数MR能否改进我们现有的模型呢？为了分析可能受几个变量（或属性）影响的关系，就有必要评估每一变量的贡献。在该数据集中，π和MR在某种程度上相关，所以拟合的顺序可能影响第二个变量对第一个的帮助程度。多元线性回归可用于决定多个变量在数据整体拟合中的相对重要性。

多元线性回归通过上调或下调现有的每一参数（最大化r2值）以使数据与回归方程最佳拟合（最小化线周方差）。回归程序通常以逐步逼近方式求解。即通过对参数的加减逐步回归使r2和s值优化。这种方式求出的系数反映了相应参数对生物学活性的相对贡献。

应用多元线性回归分析时，有两点需要说明。第一，如有足够多的参数，任何数据集都可以拟合一条回归直线。这样一来的结果就是，回归分析所需要的化合物通常要显著地比参数多；一种有用的经验方法是化合物数目为参数数目的3-6倍。困难的是，回归分析对内插值最为有效；但在合成中最有用的是外推法（即，回归分析所描述的实验区域已经是清楚的了，而将其用于一个新的尚未分析区域则可能问题多多）。

将多元线性回归用于辣椒素类似物，我们可以推导出生物学活性与疏水性及摩尔折射率的相关性方程如下。

LogEC50 = 0.762 - (0.819)π+ (0.011)MR

s = 0.313, r2=0.888

判断一个回归项的重要性需要考虑三点。

1. 回归系数的统计学显著性。

2. 典型的bixi效应（本例中为0.011、25.36*）的大小（*黄子曰：搞求不懂。）。

3. 任何与其它项的交叉相关性。

随着多元线性回归项的增多，r2值总是变大。我们按3位有效数字验算先前的计算（r2=0.89）以免舍入引起混淆。

上述分析的结果显示，在该系列内，原子的空间体积(steric bulk)对于活性而言并非一个重要的影响因子。疏水性常数的影响力支持疏水性结合位点的存在。由于该分析中，取代基团数目有限，不大可能通过进一步分析得出更多的推论。

本节中，我们逐步讨论了QSAR研究的数学基础。几位学者从不同侧面探讨了QSAR方法并已经发表了相关综述3-8。接下来的每个例子使用了这些技术来推导对活性很重要的化学因素的相关信息。 

QSAR诸方法

药物发挥其生物学功效涉及一系列事件，包括转运、与受体结合及代谢失活等。由于在多数情况下药物分子与推定受体间的相互作用机制尚不清楚（即，没有结合的晶体结构），这迫使我们只得从容易得到的属性（分子属性和描述符）进行推论，以此解释已知分子的这些相互作用。一旦关系确定，就可用于帮助预测新的或未知的分子。

推导定量构效关系的第一种方法将活性描述为基于热力学定律的化学结构的函数。从给定分子推导出的一系列参数，可用于表征自由能术语自由能变（ΔE）、焓变（ΔH）、熵变（ΔS）*。（*状态函数：在热力学中，把那些仅取决于状态而与变化途径无关的物理量，称状态函数。如温度、压力、体积、内能、焓、自由焓、熵等都是体系的状态函数。焓（enthalpy），表征在等压过程中物体吸收的热量，为一个体系的内能与体系体积和外界施加于体系的压强的乘积之和，用符号H表示；焓变（ΔH）为焓的变化。熵（entropy），指热能除以温度所得的商，标志热量转化为功的程度。用符号S表示）

电性效应如供失电子倾向、部分原子电荷、静电场密度等可用Hammettσ值、共振参数（resonance parameter, R值）、场参数（inductive parameter, F值）、Taft取代值（ρ*、σ*、ES）等表征。立体效应如分子体积和表面积可用摩尔折射率和Taft立体参数表征。焓效应可用分配系数（LogP）或由从分配系数导出的疏水性参数π计算。此外，各种结构指标可用于描述在分子内某些位点是否存在特异性功能基团。描述活性与上述参数集之间关系的线性方程就是Hansch方程

log1/[C] = A(logP)-B(logP)2+C(ES)+D(ρσ)+E+…

应用多元线性回归分析求出系数。通常，Hansch型研究针对这样的化合物，它们具有一个共同摸板（通常是刚性的如一个芳香环），结构上的差异仅限于特定位点上功能基团的变化。

Hansch使用上述方法分析了256个 4,6-二氨基-1,2-二氢-2,3-二甲基-1-(X-苯基)-s-三嗪类化合物，它们具有抗肿瘤二氢叶酸还原酶活性9。已发现，其中244个化合物的活性与N-苯环3和4位上存在疏水基团有关。推导这种关系时使用的参数包括：疏水性常数（π）和摩尔折射率（MR）针对N-苯环的间位和对位取代基，6个指示变量I1-I6用于表示有（值为1）或无（值为0）某种特异性结构特征。使用最小二乘法从这些数据中的出方程如下。

图6：Baker三嗪分析

log 1/[C] = 0.680(π3) - 0.118(π3)2 + 0.230(MR4) - 0.024(MR4)2 + 0.238(I1) - 2.530(I2) –

1.991(I3) + 0.877(I4) + 0.686(I5) + 0.704(I6) + 6.489

n = 244, r = 0.923, s = 0.377

求该方程的偏导数可得出MR4 的最佳值4.7和π3的最佳值2.9。注意，该数据集中化合物数目已减为244个。去除上述10个或更多变量不能正确预测的12个化合物，Hansch和Silipo报道的r和s值更佳。

尽管Hansch方法有许多限制，但它能用简单参数成功建立复杂生物学体系的模型。该方法已成功用于预测很多生物学检测的取代效应。该方法存在的主要问题是需要充分研究各种结构组合成的大量化合物。而且，这种分析方法没有考虑构象效应。有些学者已发表文章，提供了有关Hansch方法的更多背景10,11。

为避免Hansch型分析的组合问题，产生了另一种化合物设计方法。Free和Wilson使用了一系列取代常数，其相应特异性功能基团出现在母体化合物特定位点与生物学活性相关。

Activity = A + ΣiΣjGijXij

其中A为该系列化合物平均生物学活性，Gij为功能基团i在第j位时对活性的贡献，Xij为在第j位有（1）或无（0）功能基团i时对对活性的贡献。

使用上述方程的程序是建立系列矩阵并用以表示系列方程。然后推导出每个功能基团在每个位置的取代常数。接着使用统计学测试检验这些常数的重要性。如果显示模型有效，那么就可用于预测尚未合成化合物的活性值。总的说来，虽然需要大量化合物以研究多种替代类型的效应，但Free-Wilson方法实际上降低了所需类似化合物的数目。然而，该方法要求取代基的效应是相加性的。

1972年，John Topliss发表了一篇自动化Hansch途径详细方法的文章2。其方法假定所研究的先导化合物包含至少一个苯环，而苯环可以作为功能基团修饰的模板。对模板的第一种修饰是制备对氯衍生物以测定亲脂性。然后进行其它各种取代以探究并优化活性与分子疏水和电性特征的关系。尽管Topliss方法易于进行，但却有几大缺点。主要的问题在于其方法不能用于其它类型的研究，而且使用该方法具有很高的风险（它基本上不考虑取代基团间的相互作用，因为每次它仅改变一个取代基团）。

1960年代，作为将测定活性与化学属性相关联的工具，经典QSAR的应用巨增。但在许多领域，这些技术还不能运用或不能提供有用的关联。这些情况包括：活性取决于三维几何学、化合物的训练集不够、化合物数量太少、多样性不足、生物学活性不能很好定量等。通过扩展Hansch方法和发展新的QSAR方法，上述问题中许多已经得到处理。

有这样一些情况，由于各种原因，无法准确测得生物学活性值，如缺乏敏感的检测系统。在这些情况下，可以使用其它统计学方法；问题可以简化为一张分类表，表中化合物分为活性、部分活性、无活性等。然后在作为结果的数据集中搜寻可用于预测上述分类的模式。这种模式分析方法包括SIMCA(Soft, Independent Modeling of Class Analogy)13、ADAPT(Automated Data Analysis by Pattern recognition Techniques)14、CASE(Computer Automated Strucutre Evaluation)15和CSA(Cluster Significance Analysis)16。

模式识别方法17试图定义这样的参数集，这些参数可以使具有相似活性的化合物在n维空间的某些区域聚类。实现上述目标的方法可以是参数的或非参数的。参数法根据计算出的化合物属性在在n-维空间搜索聚类。这些方法不使用衍生值（例如，均值向量和协方差矩阵），而是在原始数据中找出聚类定义并应用迭代方法找出能最佳地确定分类表的线性参数集。

上述方法推导出的是判别式函数。而SIMCA方法使用主成分分析（Principal Component Analysis, PCA）来描述数据集。PCA的目的是减少描述生物学活性的变量数或将化学属性减到独立的较少的那么几个。通过分析生物学或化学属性的相关系数矩阵可以完成这一工作。

通过分解相关系数矩阵，主成分分析可用来为每个分类（如，活性与非活性）创建衍生变量；这可用于指出冗余或变量之间的相互关系。PCA通过将原始参数转化为一套称为“主成分”的互不相关的变量，试图从数据中找出简化的关系。通过本征值将对称的相关系数矩阵分解。用最大的本征值及其本征向量构建具有最大方差的原始变量的线性组合。然后逐步用稍小的本征值及其本征向量构建方差渐小的原始变量的线性组合。后续的本征向量之间无关。通过略去小本征值相应的本征向量而简化。总之，这些步骤将为数据找出正交座标轴集，而这可用于分辩数据中的差异。

Jurs研究组发展了检验化学结构影响活性的另一种方法。他们不是依赖多元统计分析来阐明构效关系，而是联合运用聚类分析和模式识别技术来找出这些关系。ADAPT程序产生一套分子描述（拓扑的、几何的、理化的）数据，这些数据源自三维模建。将这些数据点投射到n-维平面并应用模式识别方法进行分析。该分析的目的是从一系列化合物中区分活性与非活性化合物。

Jurs报道了ADAPT中方法学的几种应用。在一项对化学致癌原的研究中18，从一套含28个结构特征（包括片段描述符、亚结构描述符、分子连接性描述符、几何描述符）中推导出了一个线性的判别式函数。该研究选择了12个结构类的209个（其中致癌原130个、非致癌原79个）化合物。程序确定了一个含192个化合物的训练集，用以产生最佳描述符并分析整个数据集。在随机测试中，预测致癌性化合物的成功率为90%，预测非致癌性化合物的成功率为78%。

CASE程序使用拓扑学方法确定了活性必需亚片段，扩展了ADAPT中的技术。CASE能够区分位点性同分异构体。在分析结构相似的化合物的数据集时，CASE和ADAPT的应用均受限。

我们讲到现在为止，所有的方法在分析分子活性时都没有明确包含分子三维形状的贡献。虽然使用化学图形指数18、分子间结合距离19、分子表面积20、静电势21时包含了一些分子3D形状的信息，但率先详尽分析这些效应的是Hopfinger22研究组和Marshall23研究组。

在1979年，Marshall扩展2D-QSAR方法，明确地考虑了反映3D形状的构象柔韧性23。活性类似物方法（Active Analog Approach, AAA）的第一步是搜寻在某一生物学测试中有高度活性的化合物的所有构象。搜索的结果是一张原子间距离图，这可用于过滤该系列后续的构象搜索。该方法内在的假设是，具有相似活性谱的所有化合物具有相似的构象。一旦确定了“活性构象”，就可以计算每个分子的体积并重叠。用回归分析建立体积与生物学活性的关系Marshall及其同事将活性类似物方法商业化并与其它药物设计技术一起置入SYBYL分子模建软件包。

Hopfinger及其同事也在QSAR中使用了3-D形状。在对Baker三嗪的分子形状研究中24，计算出了该系列所有分子的共享空间及其静电场的差异。但这些计算与一套分子叠放规则结合时，就得到了不同分子的相对形状指数。将这些形状参数纳入标准的Hansch分析，可以改进计算参数与生物学活性的相关性，甚至可以不需从原始数据集中排除任何化合物。Hopfinger及其同事的该项技术已纳入CAMSEQ、CAMSEQ-II、CHEMLAB、CAMSEQ-M等程序中。

在1988年，Richard Cramer提出，可通过分析分子中取决于形状的立体参数和静电场研究分子生物学活性26。此外，比较分子立场分析（Comparative Molecular Field Analysis, CoMFA）采用偏最小二乘法（Partial Least Squares, PLS）和交叉验证(cross-validation)等新的数据分析方法而不是仅限与从数据拟合一条回归直线来建立活性预测的模型。

CoMFA要求科学家定义与推定药效基团叠放的系列分子的比对规则；活性构象与比对规则必须是特异的。比对时，程序把每个分子固定在一个三维网格里，计算三维网格里横切面网格点上探针原子（如一个SP3碳原子）与每个分子间静电与分子静力学力场的立体组分。得出的方程形式如下：

Act1 = Const1 + a1(stericxyz) + b1(stericxyz) + ... + a'1(estaticxyz) + b'1(estaticxyz) + ...

Act2 = Const2 + a2(stericxyz) + b2(stericxyz) + ... + a'2(estaticxyz) + b'2(estaticxyz) + ...

Actn = Constn + an(stericxyz) + bn(stericxyz) + ... + a'n(estaticxyz) + b'n(estaticxyz) + ...

 

传统回归方法要求数据集中的参数数目必须显著少于化合物数目（或数据的自由度）。而CoMFA分析得出的数据表中，参数数目远远多于化合物数目。PLS解除了传统方法的限制，用来推导所有立体和静电项的系数。PLS主要基于以下事实：一个分子相邻部分间的相互关系相似，故真实的维数小于网点的数目。因为这些系数取决于位置，所以该系列分子中与活性增减有关的立体体积与静电荷区域的取代模式就清楚了。作为一种预测工具，模型的大小（最佳模型所需组分数目）和有效性可通过交叉验证检验。

与传统回归方法不同，交叉验证通过模型对数据的预测程度而不是模型对数据拟合程度来评估其有效性。分析中使用余一法（leave-one-out）；用N-1个化合物建立模型，并用该模型对剩下的这一个（即第N个）化合物进行预测。依次将每个化合物排除在模型推导之外并预测。模型的性能可用交叉验证（或预测）的r2表示

r2 (cross-validated) = (SD - Press)/SD

SD为活性离均方差和。Press（Predictive Sum of Squares，预测方差和）为化合物除外时预测活性与实际活性差的平方和。

我们在前面已经讨论过，传统r2的取值范围为0~1。但交叉验证方法的r2的取值范围为-1~1。负值表示用活性的均值来估计比用模型估计好（即，与其用模型预测不如没有模型）。但得到交叉验证r2最高的模型时，就用它来推导传统的QSAR方程、r2和s值。模型的最终结果以系数的等高线图显示。

第一例CoMFA研究分析了21个甾体结构与人皮质类固醇结合球蛋白及人睾酮结合球蛋白的结合亲和力。该类化合物是刚性的，故研究中选用以消除构象依赖性效应。每一个甾体的模型使用来自剑桥晶体数据库的坐标构建并用Tripos力场最小化。侧链位置由系统构象搜寻完成。使用场配算法（Field Fit algorithm）在固定的网格（3-D网格用于计算CoMFA力场效应）比对每一结构。比较预测值和实际值，皮质类固醇交叉验证方法的r2为0.65（传统r2=0.897，s=0.397）；对睾酮结合甾体，交叉验证方法的r2为0.555（传统r2=0.873，s=0.453）。

如上所述，CoMFA从定义药效基团和叠合规则开始，推导可用于预测新化学实体活性的3-D模型。Accelrys（前分子模拟公司，Molecular Simulation Incorporated, MSI）的Apex和Catalyst程序可从化学结构与生物学活性的数据库中确定药效基团。这些模型而后可用于预测新化合物的活性。

Apex-3D是一个确定药效基团的自动化系统。它采用统计学技术与3-D模式匹配算法从一套生物活性分子中确定出可能的药效基团。该程序使用三种方法进行分子结构分类：一致归纳法从具有相似活性的化合物中确定共同的结构模式；差异归纳法确定造成化合物活性与非活性差异的结构模式；相伴变异归纳法阐明引起系列化合物生物学活性变化的结构变化特征。

上述方法遵循的逻辑与医学化学家工作时的逻辑相似。什么样的药效基团模式存在于活性分子中而不存在于非活性分子中？药效基团由不同化学中心（以原子为中心的功能基团）及这些中心间的距离来描述。这些中心可以是芳香环、供电子能力、氢键形成位点、亲脂性区域、局部原子电荷。每个分子的信息以规则形式的知识库存贮，知识库可用于预测新结构的活性。

Apex-3D含有一个专家系统，可以根据已确定的药效基团自动选择最佳构象并匹配结构。当有定量的生物学数据时，可以生成任何可能确定的药效基团的3D-QSAR模型。根据生物学活性数据的类型，该程序可确定不同结合取向、不同受体亚型、激动剂或拮抗剂的药效基团。

运用该方法要求科学家将化合物分配到一个或多个活性类别的训练集。如果需要，也可以自定义特异的化学中心。程序中自动确定药效基团的部分按以下步骤建立知识库：

l 确定数据集中每个化合物所有可能的相互作用中心；

l 根据变元集产生拓扑学（2D）或拓扑图形（3D）的距离矩阵；

l 采用聚类算法从所有分子对中确定可能的药效基团；

l 根据药效基团在每一活性分类的化合物中的采用贝氏统计的出现率与非机会出现率对药效基团进行分类；

l 设定与药效基团相关的统计学概率与置信区间，使药效基团规则能对训练集所有分子进行适当分类；

l 将极可能含有药效基团的化合物与药效基团匹配。

建立知识库后，科学家就可以用它来预测训练集外化合物的生物学活性。

把计算的“母核”（biophore）位点、次级位点或整个分子属性的指标关联起来，上面所定义的药效基团可用于建立3D-QSAR模型。Apex中生成的3D-QSAR模型采用选择变量改良的多元线性回归分析进行筛选。特殊的随机检验评估偶然相关。涉及步骤有：

l 交互式地选择药效基团用于分析；

l 从程序可能选择的参数库中交互式地选择参数

l 采用逐步多元回归计算每一药效基团的最佳3D-QSAR模型，并用统计学方法评估其有效性与预测能力。

至此，3D-QSAR模型的知识库已可用于计算新化合物的活性。对于一个新结构的一套构象，根据每一构象所含药效基团与3D-QSAR模型的匹配情况可以计算出其活性范围。

迄今为止，Apex的应用见于文献的还很少。下面讨论一个具有启发价值的将其运用于非肽血管紧张素II拮抗剂研究的例子。文献中31-44报道几种结构不同的该活性类别的化合物。在所有高活性分子是否具有相同结合位点存有争议的情况下30假定出药效基团模型。然后使用药效基团自动化确定过程分析这些化合物并评估它们作用于同一位点的可能性。研究中采用了特异性结合活性(IC50)变化范围达6个数量级的一套共55个化合物。使用几何结构经AMPAC优化含多个构象的3D结构。根据活性，化合物分为四类：活性最高(<100nM)的有27个化合物。Apex-3D生成没有假阳性和假阴性的可以正确分类全部最高活性类化合物的规则。存在几个“母核”的事实与存在多个结合位点是一致的。

对其中某些药效基团，还可以推导出定量的3D-QSAR模型。有一个含48个化合物的模型，它的统计学参数如下：预测R2=0.83，预测RMSE=0.86。训练集外化合物的预测活性在统计学范围内。

与Apex一样，Catalyst也从活性不同的一套分子中产生结构-活性假说。当所有分子的连接性与活性值确定了，Catalyst在相对特异的位置推导出由一套通用化化学函数（表面疏水区域、氢键向量、电荷中心或其它用户自定义特征）构成的假说。Catalyst可以生成多达10个假说并将其按估计的统计学显著性排序。这些假说可用图形方式检验，或与新的分子匹配、或直接用于柔性3D数据库搜寻。

其过程的第一步是找出每一个分子覆盖低能构象空间的代表性构象。选择代表性结构以最大化构象空间的采样率。第二步是列出在活性化合物中常见而在非活性化合物中罕见的候选特征（由于训练集中所有成员不大可能都具有全部的重要结合特征，所以活性类无需纳入全部确定的特征）。假说的“代价”（cost）由描述假说所需“笔”数（number fo bits）和用假说预测活性时的错误定义。最小复杂度理论表明，一个预测性假说应使其代价最小化。上面确定的所有候选假说都要进行最小化。通过非参数随机检验对各种医学训练集结果的统计学显著性（找出机会相关的概率很低）进行验证。

构建QSAR模型的一项新技术是遗传函数近似算法（Genetic Function Approximation, GFA）45。该方法尤其适用于数据集所含变量远多于样本数，或数据集中变量与活性间有非线性关系的情况。

GFA从一组随机构建的QSAR模型总体开始；通过评估每个模型的相对预测力作为一个衡量误差的指标，并以此对上述随机模型进行排序。选择排在前两位的较好模型作为母模型，沿用每个母模型的参数项创建新的子模型，如此反复进行使总体逐步进化。总体中最差的模型被新的模型取代。随着进化，总体中质量越来越高的模型就丰富起来。

不同的模型实际上是对数据的多种拟合。然后，科学家能够利用他们的科学知识和直觉从最终的模型中进行选择。比较模型间的异同可能促使科学家考虑其它能够解释数据的机制或设计新的实验以决定可能的各种模型。由于模型通常仅使用数据集全部变量中的一部分，GFA表现得象变量简化算法类中的一员；数据中关键信息集中于少数变量时，该方法最为适用。

实验结果与发表数据的比较证实，GFA推导的模型与逐步回归、线性回归、偏最小二乘法回归等标准方法推导的模型相当，在有些例子中甚至更优。随着计算源和仪器源实验数据的增加，GFA简化变量的能力、GFA发现非线性关系的能力、GFA找出代表数据集多种解释的多种模型的能力，将在数据分析中变得越来越有用。

除结合之外，还有许多与分子生物学活性有关的属性，包括运输、分布、代谢和清除。由于上面所讨论的技术强调结合，通常没有考虑其它任何过程，所以也就常常失败。但是，我们应该注意到，它们仍是一种有效而且可能的应用。它们通常比专门或直觉的方法更有效也更成功。由于没有一种方法适用于所有情况，所以很多研究组仍在开发能构建含电子变量和拓扑学指标SAR方程的方法。

Peter Goodford报道了一种基于能量的化合物设计方法——网格法。GRID程序46是检测已知结构分子在能量上有助于结合的位点。用特异性探针基团与目标结构间的静电、氢键及Lennard-Jones相互作用计算能量。GRID也用于区分不同探针的选择性结合位点。

总之，拓扑学方法始于与分子结构不同的理论代表图；该图用顶点表示原子，用边表示化学键。接着就是定义结构的数字指标，这些指标包含有亚基团变量与原子变量（原子类型、原子重量、原子数、化合价与核心电子比等）的概括信息。此外，还会生成一系列指标（χ指标），它用一套不同大小和复杂度的片段描述分子。回归分析应用中，χ指标将化学结构与色谱停留时间、摩尔折射率、电离能、原子化热（heats of atomization）等理化特性关联。虽然将拓扑学指标应用于QSAR已广为报道47-49，但这些方法通常局限于预测多聚物和烃类。 

小结

推导定量构效关系并不容易。分子通常是柔性的，可能有很多有用的属性与活性相关。在一个研究项目的早期，用于建摸的化合物通常很少。因此，在一个极为多元的变量空间，我们只有很少的样本（化合物）。哪些变量重要，我们怎么优化它们？必须清楚的是，需要许多训练化合物覆盖整个活性空间，而且所需拟合技术强调的不仅是推导拟合，而且也强调拟合的预测质量。尽管这些方法尚未发现一个新的化合物，但它们已经帮助科学家检验了研究项目中产生的大量数据。随着方法的发展，它们在诸如组合化学等领域中的应用将会越来越广阔。

参考文献

1. Christopher S.J. Walpole, Roger Wrigglerworth, Stuart Bevan, Elizabeth A. Campbell, Andy Dray, Iain F. James, Kay J. Masdin, Martin N. Perkins and Janet Winter, J. Med. Chem., 36, 2381 (1993).

2. John G. Topliss, Utilization of Operational Schemes for Analog Synthesis in Drug Design, J. Med. Chem., 15, 1006 (1972).

3. Yvonne C. Martin, A Practitioner's Perspective of the Role of Quantitative Structure Activity Analysis in Medicinal Chemistry, J. Med. Chem., 24, 229 (1981).

4. John G. Topliss, Quantitative Structure-Activity Relationships of Drugs, Academic Press, New York, (1983).

5. Franke, R., Theoretical Drug Design Methods, Elsevier, Amsterdam, (1984).

6. Seydel, J.K., QSAR and Strategies in the Design of Bioactive Compounds, VCH, Weinheim, (1985).

7. Yvonne C. Martin, Accounts of Chem. Res., 19, 392 (1986).

8. Corwin Hansch, Accounts of Chem. Res., 2, 232 (1969).

9. Corwin Hansch and Carlo Silipo, J. Amer. Chem. Soc., 97, 6849 (1975).

10. Robert F. Gould (ed.), Biological Correlations -- The Hansch Approach, Advances in Chemistry Series, No. 114, American Chemical Society, Washington, D.C., 1972.

11. Yvonne C. Martin, Quantitative Drug Design, Marcel Dekker, New York, (1978).

12. Spencer M. Free and James W. Wilson, A Mathematical Contribution to Structure-Activity Studies, J. Med. Chem., 7, 395 (1964).

13. Svante Wold, Pattern Recognition by Means of Disjoint Principal Components Models, Pattern Recognition, 8, 127 (1976).

14. A.J. Stuper, W.E. Brugger and P.C. Jurs, Chemometrics: Theory and Application, B.R. Kowalski (Ed.), American Chemical Society, Washington, D.C., 1977.

15. G. Klopman and M.I. Dimayuga, J. Comput.-Aided Mol. Design, 4, 117 (1990).

16. J. McFarland, J. Med. Chem.,

17. Peter C. Jurs, Chemometrics and Multivariate Analysis in Analytical Chemistry, in Reviews in Computational Chemistry, Volume 1, K.B. Lipkowitz and B.B. Boyd (Eds.), VCH Publishers, Inc., New York, 1990.

18. Peter C. Jurs, J. T. Chou and M. Yuan, J. Med. Chem., 22, 476 (1979).

19. Lemont B. Kier, Molecular Orbital Theory in Drug Research, Academic Press, New York, NY, 1971.

20. Gordon M. Crippen, J. Med. Chem., 22, 988 (1979).

21. M. Mabilia, R.A. Pearlstein and A.J. Hopfinger, Computer Graphics in Molecular Shape Analysis, in Molecular Graphics and Drug Design, A.S.V. Burgen, G.C.K. Roberts and M.S. Tute (Eds.), Elsevier Science Publishers, Amsterdam, 1986.

22. H. Weinstein, in Chemical Applications of Molecular Electrostatic Potentials, Peter Politzer and Donald G. Truhlar (Eds.), Plenum Press, New York, NY, 1981.

23. R. Potenzone, E. Cavicchi, H.J.R. Weintraub and A.J. Hopfinger, Comput. Chem., 1, 187 (1977).

24. G.R. Marshall, C.D. Barry, H.E. Bosshard, R.A.. Dammkoehler and D.A. Dunn, The Conformational Parameter in Drug Design: The Active Analog Approach, in Computer Assisted Drug Design, ACS Symposia, 112, E.C. Olson and R.E. Christofferson (Eds.), American Chemical Society, Washington D.C., 1979

25. Anton J. Hopfinger, J. Amer. Chem. Soc., 102, 7196 (1980).

26. Richard D. Cramer III, David E. Patterson and Jeffrey D. Bunce, J. Amer. Chem. Soc., 110, 5959 (1988).

27. Components are defined as the linear combination of all parameters (the independent data) which are used to fit the shape of the parameters to the shape of the activity data (the dependent data).

28. V. E. Golender and A. B. Rozenblit, Logical Structural Approach to Computer Assisted Drug Design, in Drug Design, Vol IX, Academic Press (1980).

29. V. E. Golender and A. B. Rozenblit, Logical and Combinatorial Algorithms for Drug Design, Research Studies Press, UK (1983).

30. V. E. Golender and E. R. Vorpagel, Computer Assisted Pharmacophore Identification, in 3D QSAR in Drug Design: Theory, Methods and Applications, ESCOM Science Publishers, Netherlands (1993).

31. A. P. Thomas, C. P. Allott, et. al., "New Nonpeptide Angiotensin II Receptor Antagonists. 1. Synthesis, Biological Properties, and Structure-Activity Relationships of 2-Alkyl Benzimidazole Derivatives", J. Med. Chem., 35, 877-885 (1992).

32. R. H. Bradbury, C. P. Allott, et. al., "New Nonpeptide Angiotensin II Receptor Antagonists. 2. Synthesis, Biological Properties, and Structure-Activity Relationships of 2-Alkyl-4-(biphenylylmethoxy)quinoline Derivatives", J. Med. Chem., 35, 4027-4038 (1992).

33. R. H. Bradbury, C. P. Allott, et. al., "New Nonpeptide Angiotensin II Receptor Antagonists. 3. Synthesis, Biological Properties, and Structure-Activity Relationships of 2-Alkyl-4-(biphenylylmethoxy)pyridine Derivatives", J. Med. Chem., 36, 1245-1254 (1993).

34. P. R. Bovy, D. B. Reitz, et. al., "Nonpeptide Angiotensin II Antagonists: N-Phenyl-1H-pyrrole Derivatives Are Angiotensin II Receptor Antagonists", J. Med. Chem., 36, 101-110 (1993).

35. N. B. Mantlo, P. K. Chakravarty, et. al., "Potent, Orally Active Imidazo[4,5-b] pyridine-Based Angiotensin II Receptor Antagonists", J. Med. Chem., 34, 2919-2922 (1991).

36. W. T. Ashton, C. L. Cantone, et. al., "Nonpeptide Angiotensin II Antagonists Derived from 4H-1,2,4-Triazoles and 3H-Imidazo[1,2-b][1,2,4]triazoles", J. Med. Chem., 36, 591-609 (1993).

37. K. S. Atwal, S. Z. Ahmed, et. al., "Dihydro pyrimidine Angiotensin II Receptor Antagonists", J. Med. Chem., 35, 4751-4763 (1992).

38. J. V. Duncia, A. T. Chiu, et. al., "The Discovery of Potent Nonpeptide Angiotensin II Receptor Antagonists: A New Class of Potent Antihypertensives", J. Med. Chem., 33, 1312-1329 (1990).

39. B. De, M. Winn, et. al., "Discovery of a Novel Class of Orally Active, Non-Peptide Angiotensin II Antagonists", J. Med. Chem., 35, 3714-3717 (1992).

40. K. Kybo, Y. Kohara, et. al., "Nonpeptide Angiotensin II Receptor Antagonists. Synthesis and Biological Activity of Potential Prodrugs of Benzimidazole-7-carboxylic Acids", J. Med. Chem., 36, 2343-2349 (1993).

41. I. Sircar, J. C. Hodges, et. al., "Nonpeptide Angiotensin II Receptor Antagonists. 2. Design, Synthesis, and Structure-Activity Relationships of 2-Alkyl-4-(1H-pyrrol-1-yl)-1H-imidazole Derivatives: Profile of 2-Propyl-1-[[2'-(1H-tetrazol-5-yl)-[1,1'-biphenyl]-4-yl]-methyl]-4-[2-(trifluoroacetyl)-1H-pyrrol-1-yl]-1H-imidazole-5-carboxylic Acid (CI-996)", J. Med. Chem., 36, 2253-2265 (1993).

42. J. Weinstock, R. M. Keenan, et. al., "1-(Carboxybenzyl_imidazole-5-acrylic Acids: Potent and Selective Angiotensin II Receptor Antagonists", J. Med. Chem., 34, 1514-1517 (1991).

43. R. M. Keenan, J. Weinstock, et. al., "Imidazole-5-acrylic Acids: Potent Nonpeptide Angiotensin II Receptor Antagonists Designed Using a Novel Peptide Pharmacophore Model", J. Med. Chem., 35, 3858-3872 (1992).

44. R. M. Keenan, J. Weinstock, et. al., "Potent Nonpeptide Angiotensin II Receptor Antagonists. 2. 1-(Carboxybenzyl)imidazole-5-acrylic Acids", J. Med. Chem., 36, 1880-1892 (1993).

45. D. Rogers and A. J. Hopfinger, "Applications of Genetic Function Approximation (GFA) to Quantitative Structure-Activity Relationships (QSAR) and Quantitative Structure Property Relationships (QSPR)", J. Chem. Inf. Comp. Sci., 34, 854-866 (1994).

46. Peter J. Goodford, J. Med. Chem., 28, 849 (1985).

47. A.T. Balaban, Applications of Graph Theory in Chemistry, J. Chem. Inf. Comput. Sci., 25, 334 (1985).

48. Lowell H. Hall and Lemont B. Kier, The Molecular Connectivity Chi Indexes and Kappa Shape Indexes in Structure-Property Modeling, in Reviews in Computational Chemistry, Volume 2, K.B. Lipkowitz and D.B. Boyd (Eds.), VCH Publishers, Inc., New York, 1991.

49. I.B. Bersuker and A.S. Dimoglo, The Electron-Topological Approach to the QSAR Problem, in Reviews in Computational Chemistry, Volume 2, K.B. Lipkowitz and D.B. Boyd (Eds.), VCH Publishers, Inc., New York, 1991.