药物分子设计简介

传统药物设计从总体上来讲，缺乏成熟完善的发现途径，具有很大的盲目性，一般平均要筛选10000种化合物以上才能得到一种新药，因此开发效率很低，很难迅速得到合适的新药来治疗越来越多的疑难杂症。
  随着计算机技术及计算化学、分子生物学和药物化学的发展，药物设计进入了理性阶段，其中药物分子设计是目前新药发现的主要方向。它是依据生物化学、酶学、分子生物学以及遗传学等生命科学的研究成果，针对这些基础研究中所揭示的包括酶、受体、离子通道及核酸等潜在的药物设计靶点，并参考其它类源性配体或天然产物的化学结构特征，设计出合理的药物分子。计算机辅助药物设计方法（CADD）是药物分子设计的基础。从20世纪60年代构效关系方法（QSAR）提出以后，经过40多年的努力和探索，尤其是20世纪90年代以后，随着多种新的方法的出现，CADD方法已经发展成为一门完善和新兴的研究领域，它大大提高了药物开发的效率，为人们攻克一些顽症提供了崭新的思路和成功的希望。
计算机辅助药物设计方法（CADD）大体可以分为三类：1.基于小分子的药物分子设计方法，这类方法主要是针对受体结构未知的药物分子，主要包括定量构效方法和药效团模型方法；2.基于受体结构的药物分子设计方法。随着分子生物学和结构生物学的发展，越来越多的生物大分子结构被解析。因此对于一些未知三维结构的受体大分子，它们的结构常常也可以通过同源蛋白质建模的方法得到。在这种情况下，就可以采用基于受体结构的药物分子设计方法来寻找新的先导化合物。基于受体结构的药物分子设计方法的思路是通过研究受体结构的特征以及受体和药物分子之间的相互作用方式来进行药物设计。常用的方法是分子对接方法和从头设计方法；3.计算组合方法。主要包括两部分的内容，一方面是采用计算机技术设计合成组合库的构造块，通过计算机生成包含足够分子多样性的虚拟组合库；另一方面则是把得到的虚拟组合库和其它分子设计方法结合起来进行药物分子设计。当然上面介绍的这些方法之间也并不是完全独立的。对于某个体系常常会采用多种药物设计方法进行先导化合物的设计。如知道了蛋白质和药物分子相互结合的复合物结构时，可以从蛋白质结构出发用基于分子对接的数据库搜索方法进行先导化合物的寻找；也可以从复合物的相互作用信息，得到药效团模型，然后从药效团模型出发进行药物设计；同时还可以从蛋白质结构出发，用片断生长的方法进行全新药物分子的设计[也就是从头设计方法]。
  计算机辅助药物设计方法（CADD）中目前还有几个难题：1.受体和配体之间自由能的评估，因为分子活性的大小常常是由受体和药物分子之间的结合自由能来决定的，所以体系的自由能得计算非常重要。虽然目前已有很多方法，如线性相互作用能方法等，但总的来看结合自由能的评估还需要做大量的工作；2.大分子的构象问题。确定构象的方法有两种：实验方法和理论计算方法。用实验方法解析的大分子结构的数目有限，而常用的构象分析方法，如系统搜索方法、距离几何方法、蒙特卡洛模拟方法、遗传算法和分子动力学方法等，对于大分子还很难确定出其优势构象，所以发展新的蛋白质结果预测方法是计算化学的一个重要课题；3.溶剂效应的考察。如何处理溶剂问题是计算机辅助药物分子设计中一个难题。我们所关心的大量生命过程都是在溶剂条件下进行的，在溶剂和真空中，分子的构象可能会存在较大的区别，所以发展快速而准确的溶剂效应计算方法是计算机辅助药物分子设计的一个重大挑战。
  随着计算机技术、大分子计算方法和生物信息学的发展，有理由相信计算机辅助药物设计方法（CADD）会进入一个新的发展时期，人类古老的梦想—理论设计好所需物质然后再进行合成的时代正在成为现实！

 

药物设计常用方法QSAR方法介绍
一QSAR方法定义和发展过程简介
定量构效方法（quantitative structure-activity relationship, QSAR）是应用最为广泛的药物设计方法。所谓定量构效方法就是通过一些数理统计方法建立其一系列化合物的生理活性或某种性质与其物理化学性质之间的定量关系，通过这些定量关系。可以预测化合物的生理活性或某些性质，指导我们设计出具有更高活性的化合物。
早在1867年，Crum-Brow 和Fraser就提出了构效关系的概念，1900年前后，Overton 和Meyer等提出了麻醉作用的类脂学说，即化学结构各异的麻醉剂其活性随着脂-水分配系数增加而增加的现象，这可能是最早提出的化合物生理活性和物理化学性质之间的定量分配关系模型。但只是到近几十年，尤其是Hansch法提出后，随着计算机技术的发展和多变量解析技术的引入，定量构效关系方法才逐渐发展和应用起来，现在它已经成为药物设计和药物开发中不可缺少的工具。
二：常用QSAR方法介绍
1 二维定量构效关系方法（2D-QSAR）
 传统的二维定量构效关系方法很多，有Hansch法、模式识别Free-Wilson法和电子拓扑法。 其中最为著名应用最为广泛的就是Hansch和Fujita提出的Hansch法。它假设同系列化合物某些生物活性的变化是和它们某些可测量物理化学性质的变化相联系的。这些可测量的特性包括疏水性、电性质和空间立体性质等，都有可能影响化合物的生物活性。Hansch法假定这些因子是彼此孤立的，故采用多重自由能相关法，借助多重线性回归等统计方法就可以得到定量构效关系模型。Hansch法最初可以表达为下面的公式
lg1/c=algp+bσ+cE+……+constant    
既活性和疏水性参数π或lgp、电负性参数σ以及立体参数E有关。后来Hansch发现药物要交替穿过水相和类脂构成的体系，其移动难易程度和lgp呈现出函数关系。如果经过一定时间后药物在最末一相中为浓度lgc，则以lgc对lgp作图，可以发现它们之间呈抛物线关系，因此上式又可以写成下面的形式
lg1/c=a（lgp）2+b lgp +cσ+cE+……+constant   第一个式子适用于体外活性数据，而第二个式子适用于体内活性数据。
 Hansch和Fujita等人最初所采用的构效关系模型中，仅采用了一些简单的分子参数。但对于一个分子来说，可以用很多分子参数来表示分子的不同特征，比如各种拓扑参数、热力学参数、量化计算得到的参数以及分子形状参数等。研究结果表明用这些参数往往能得到更好的结果。因此在实际应用过程中，我们总是尽量选择最佳参数来得到最有效的模型而不必局限于Hansch和Fujita所提出的参数。
  此外除了传统的线性回归方法之外，一些新的数理统计方法非数值算法也被用于构效研究中，如偏最小二乘法、人工神经网络及遗传算法等，这些新方法的应用大大推动2D- QSAR方法的发展。
2 三维定量构效关系方法（3D-QSAR）
近些年来，随着构效关系理论和统计方法的进一步发展，又出现了一些三维定量构效关系（3D-QSAR）方法。譬如分子形状分析（molecularshape analysis, MSA）、距离几何方法(distance geometry, DG)和以及比较分子场分析（comparative molecular field analysis, CoMFA）方法。
与2D-QSAR比较，3D-QSAR方法在物理化学上的意义更为明确，能间接反映药物分子和靶点之间的非键相互作用特征。因此近十多年来3D-QSAR方法得到了迅速的发展和广泛的应用。在3D-QSAR方法中，比较分子场分析（comparative molecular field analysis, CoMFA）方法是目前最为成熟且应用最为广泛的方法。
2.1 比较分子场分析（comparative molecular field analysis, CoMFA）方法简介
  CoMFA的基本原理是：如果一组相似化合物以同样的方式作用于同一靶点，那么它们的生物活性就取决于每个化合物周围分子场的差别，这种分子场可以反映药物分子和靶点之间按的非键相互作用特性。其计算可以简单的分为三个步骤（1）首先确定药物分子的活性构象，再按一定的规则（一般为骨架叠加或场叠加）进行药物分子的叠合；（2）然后，在叠合好的分子周围定义一定的步长均匀划分产生格点，在每个格点上用一个探针离子来评价格点上的分子场特征（一般为静电场和立体场，有时也包括疏水场和氢键场）；（3）最后通过偏最小二乘方法建立化合物活性和分子场特征之间的关系并给出各种分子面的等势能面。
  近年来，研究人员对传统的CoMFA进行了大量的改进，其中涉及到活性构象的确定，分子叠加规则、分子场势函数的定义以及分子场变量的选取等等，在很大程度上提高了CoMFA计算的成功率。其中最具有代表性的可能就是比较分子相似因子分析（comparative molecular similarity indices analysis, CoMSIA）方法。
2.2 比较分子相似因子分析（comparative molecular similarity indices analysis, CoMSIA）方法简介
  与CoMFA方法相比，最大的不同就是分子场的能量函数采用了与距离相关的高斯函数的形式，而不是传统的Coulomb 和Lennard-Jones 6-12势函数的形式。CoMSIA方法中共定义五种分子场的特征，包括立体场、静电场、疏水场以及氢键场（包括氢键给体场和氢键受体场）。这五种分子场可以通过公式计算得到。在CoMSIA方法中，由于采用了与距离相关的高斯函数形式，可以有效地避免在传统CoMFA方法中由静电场和立体场的函数形式所引起的缺陷。由于分子场能量在格点上的迅速衰退，不需要定义能量的截断（cutoff）值. 对一些实际体系进行了这两种方法的比较结果分析，在计算中采用了不同的格点数，且对体系均采用全空间搜索策略。结果表明CoMFA计算对不同的格点大小值以及叠合分子不同的空间取向非常敏感，采用不同的空间取向时，回归系数的差值最大可以达到0.3以上。而CoMSIA方法在计算不同格点大小取值以及分子空间取向下得到的结果则稳定的多，在一般情况下，CoMSIA计算会得到更加满意的3D-QSAR模型。
3. 4D-QSAR方法简介
 1997年，hopfinger等提出了4D-QSAR的概念。作者首次采用遗传算法选择分子动力学产生的构象来产生最佳的构效关系模型。在这个方法中，作者用每个格点对用的原子占有率来作为PLS的变量，作者根据原子的不同特征定义了七种不同种类的原子模型。在4D-QSAR方法中，作者考虑了药物分子的整个构象空间，而不是一个分子，而且考察了多种原子叠合方式，因此在概念上比传统的CoMFA方法有一定的进步

三 建立定量构效关系的方法
3.1 线性回归方法
在传统二维构效研究中，多重线性回归（multiple linear regression, MLR）方法是最为常见的统计方法。一个分子可以用很多分子参数来表达，但在建立多重线性回归模型的时候，为避免过拟合（overfitting），我们只能从这些物理化学参数中选择一部分参数来建立回归模型。一般来讲，同系物数目和所选取参数数目的比应大于3～5，也有人提出应大于2的n次方（n表示选取的参数个数），怎样选取合适的参数一直是定量构效关系研究中的一个难题，而且对于线性回归来说，当体系噪声较强或干扰严重时，有可能导致所得的模型失真，为了克服多重线性回归的不足，在数学上可采用主成分回归方法。
所谓主成分回归就是采用主成分分析方法(principle component analysis, PCA)对活性影响最大的几个主要成分建立定量构效关系模型。所谓主成分是一组新的变量，它是原来变量Xij的线性组合，第一个主成分所能解释原量的方差最大，第二个次之，第三个再次之，往下依此类推。也就是说，主成分是一种线性组合，用它来表示原来变量所产生的平方误差最小。运用主成分分析，原变量矩阵X可以表达为得分（即主成分）矩阵T, 该矩阵由本征矢量上的投影所得。主矩阵与矩阵的本征矢量一一对应，即XP=T。主成分回归可以有效解决共线问题，同时由于去掉了不太重要的主成分，因而可以削弱噪声（随机误差）所产生的影响。主成分回归可以分为两步：1 测定主成分数，并通过主成分分析将X矩阵降维；2 对于降维的X矩阵再做线性回归分析。但是如果在第一步消去的是有用的主成分，而保留的是噪声，则第二步多元线性回归所得的结果就将偏离真实的数学模型。
在主成分回归法的第一步中，我们处理的仅是X矩阵，对于矩阵Y中的信息并未考虑。事实上，Y中可能包含有用的信息，所以一种很自然的想法是测试矩阵X因子时同时考虑矩阵Y的作用。偏最小二乘法在考虑自变量的同时也考虑因变量的作用，同时通过折衷各自空间内的因子使模型较好的同时描述自变量和因变量，从而有效地减少相关因素的影响。
在构效关系研究过程中，还用到了许多其他的数理统计方法，比如说逐步回归方法、非线性最小二乘方法以及一些模式识别方法等，在此就不一一介绍。
3.2 遗传算法
3.2.1遗传算法简介
遗传算法（genetic algorithm, GA）基本思想源自于达尔文的自然选择学说。它表明遗传和变异是决定生物进化的内在因素，生物的遗传特性，使生物界的物种保持相对稳定，而变异特性则使生物个体产生新的性状，以至于形成新的物种，推动了生物的进化和发展。
遗传算法是模拟达尔文遗传选择和自然淘汰生物进化过程的计算模型。它的思想源自于生物遗传学和适者生存的自然规律，是具有“生存+检测”迭代过程的搜索算法。遗传算法以一种群体中的所有个体为对象，并利用随机化技术指导对一个被编码的参数空间进行高效搜索。其中，选择、交叉和变异构成了遗传算法的遗传操作。参数编码、初始群体的设定、适应度函数的设定、遗传操作设计、控制参数设定五个要素组成了遗传算法的核心内容。
遗传算法的主要特点是1直接对结构对象进行操作，不存在求导和函数连续性的限定；2具有内在的隐并行性和更好的全局寻优能力；3 采用概率化的寻优方法，能自动获取和指导优化的搜索空间，自适应的调整搜索方向，不需要确定的规则。遗传算法的这些性质已被人们广泛的用于组合优化、机器学习、信号处理、自适应控制和人工生命等领域。作为一种新的全局优化搜获方法，遗传算法以其简单通用、适用于并行处理以及高效、实用等显著特点在各个领域得到了广泛的应用，取得了良好的效果，并逐渐成为最为重要的智能算法之一。
遗传算法的兴起是在80年代末至90年代初，最初只是用于多肽的构象分析，近些年来，遗传算法在工业分析、光谱分析、蛋白质三级结构预测和数据分析等方面也取得了广泛应用。
3.2.2 遗传算法基本原理
遗传算法的原理比较简单，传统遗传算法常常用一个二进制线性数字串表示一个个体（individual）,若干个体构成一个种群（population）.一般我们用随机的办法产生初始种群。初始种群产生以后，就用定义的“适应性函数”（fitness function）去评价种群中的每个个体。然后通过个体的得分选择部分个体并通过交叉和变异去繁衍它们的后代。当然得分高的被选中的可能性会大些，同时那些未选中的个体就自然消亡。留下来的个体和它们新产生的子代就形成了新的种群。经过一些列交叉和变异过程后，就可以得到总体得分比较高的种群，这就是我们所要寻找的具有特殊性能的构型。
我们习惯上把Holland 1975年提出的GA称为传统的GA。其主要步骤如下：1编码. GA在进行搜索之前先将解空间的解数据表示成遗传空间的基因型串结构数据，这些串结构数据的不同组合便构成了不同的点；2 产生初始种群。随机产生初始串结构数据。每个串结构数据称为一个个体，多个个体构成一个种群。GA以这个串结构数据作为初始点开始迭代优化过程；3 进行适应性指评估检测。适应性函数表明个体或解的优劣性。不同的问题，适应性函数的定义方式也不同；4 选择。选择的目的是为了从当前群体中选出优良的个体，使它们有机会作为父代为下一代繁殖子孙。选择实现了达尔文的适者生存原则；5 交叉。交叉操作是遗传算法中最主要的遗传操作。通过交叉操作可以得到新一代个体，新个体组合了其父辈个体的特性，交叉体现了信息交换的意思；6 变异。变异首先在群体中随机选择一个个体，对于选中的个体以一定的概率随机地改变串结构数据中的某个串的值。同生物界一样，变异为新个体的产生提供机会。
3.2.3 遗传算法基本特点
遗传算法作为一种快捷、简便、容错性强的算法，在各类结构对象的优化过程中显示出明显的优势。与传统搜索方法相比，遗传算法具有以下特点：1 搜索过程不直接作用在变量上，而是作用在参数集进行了编码的个体上。此编码操作使得遗传算法可直接对结构对象（集合、序列、矩阵、树、图、链和表）进行操作；2 搜索过程是从一组迭代到另一组解，采用同时处理群体中多个个体的方法，降低了陷入局部优化的可能性，并易于优化；3 采用概率变迁的规则来指导搜索方向，而不采用确定性搜索规则；4对搜索空间没有任何特殊的要求（如连通性、凸性等），只利用适应性信息，不需要导数等其他辅助信息，适用范围更广。
3.2.4 遗传算法的工作原理
遗传算法是一个以适应度函数（或目标函数）为依据，通过对群体中个体施加遗传操作实现群体内个体结构重组的迭代优化过程。它的理论主要有模式定理、积木块假设、最小欺骗问题和隐并行性问题等。我们主要介绍模式理论和最小欺骗问题。
模式理论是研究较广、影响较大的GA理论。它将GA的运算过程理解为模式操作过程并从模式运算的角度解释GA的性能特点。GA的计算过程是从一个个体组成的初始群体出发，循环的执行选择、交叉、变异过程。表面上看运算过程直接作用于群体，但其中隐含了模式集合中的操作过程。模式定理从模式操作的角度论证了在选择、交换和变异因子作用下某一模式量的变化。模式定理揭示出具有低阶、短定义矩以及平均适应度高于群体平均适应度的模式在子代中将以指数级增长，由于重组、变异的作用，并非低阶优于平均的模式都是有效模式，但在规模为N的群体中，有多个模式是有效的，即每一代GA处理个体的同时，获得了对多个模式的并行处理，并且无须额外的存储量，这一性质称为隐并行性。也就是GA在群体中搜索的过程可看作是隐含的对模式的抽样过程，然后通过遗传算子的作用将短的低阶模式的信息组合起来，形成高阶模式，直至搜索到最优解，这说明同传统的优化算法相比，GA有很强的处理能力。
研究欺骗是为了预测给定问题用GA求解的难易程度，这是我们所关注的问题。由模式理论，一个问题能否用遗传算法有效地求解，取决于问题编码是否能满足积木块假设，满足者用遗传算法求解率高，不满足者求解率较低，甚至找不到满意解，这一现象称为欺骗。最小欺骗问题旨在尽可能使问题的编码违背积木块假设，使高阶积木块不能生成最优解。但实验结果表明，对包含欺骗的问题，用GA求解时得到最优解和得不到最优解的情况都有可能发生。这说明欺骗不是衡量用GA解决问题难易程度的唯一标准。目前我们还无法判断一个问题包含欺骗的多少和问题相对于GA难易程度间的关系。

3.2.5 遗传算法的一个计算实例
遗传算法在程序上实现并不复杂，我们可以通过考虑如下简单的非线性优化问题来具体讨论遗传算法的实现过程。
函数定义如下：f(x)=x*sin(10Л* x)+1.0     (-1≤x≤2)
我们要用遗传算法得到上面这个函数的最大值。已知方程的极值为f(1.85)=2.85。整个实现过程可分为以下几个步骤：
1 编码：我们使用二进制的字符串表示方程的自变量x。字符串的长度取决于所需要的精度。定义域的范围是3。精度要求小数点后取6位有效数字，也即是要把这个定义域范围分为3×1000000个相等的区域。这样就需要22位的二进制字符串来表示每个个体，因为2的21次方小于3×1000000，所以要取到3×1000000格点[也就是前面所说的区域]，至少要取到2的22次方。
 把一个二进制字符串转化为x的实数解需要两个步骤：首先把二进制字符串转化为十进制的整数；然后把得到的整数化为实数解。
2 产生初始种群：这一步比较简单，所有的个体构成初始种群[个体就是我们所说的方程的解]，为每个个体产生一个22位二进制字符串，需要强调每个字符串对应着方程的一个解。
3 适应性函数：在这里每个个体的适应性函数就是方程f(x)的值，由于我们在这里求的是极值，所以每个个体所对应f(x)的值越大，就表明个体的适应性就越好，则被选出的几率就越大。适应性函数相当于自然进化中的环境，它会直接控制种群的演化过程。
4 遗传操作：在遗传算法中有两个主要的遗传操作：突变和交叉。突变是在群体中随机选择一个个体，然后以一定的概率随机地改变串结构中某个串的值[这相当于改变了x的值，也就改变了f(x)的值，最终导致适应性函数值的变化。交叉操作的结果类似也是改变了个体的适应性函数值，这个步骤可以加快寻找最大值的速度。
5 参数的选择：在整个迭代优化过程中种群大小、交叉几率和突变几率对计算过程会产生较大的影响。
6 计算结果：经过上面的编码操作产生初始种群，然后带入每个个体值算出f(x)的值根据适应性函数[也就是我们的要求，在这里我们的要求是取极大值，所以值越大表明个体值的适应性就越好]判断个体值的好坏，这一过程辅之以遗传操作也确保更快更好的找到所需的值。循环这一过程，直到找到使f(x)最大的值，任务完成。
 当然对不同的问题在实际操作上会有一定的差别，比如可能会对种群中的个体采用不同的数据结构，会采用其他的遗传操作。但总体上来讲，其基本流程是类似的，最大的差别就是根据不同的问题采用不同的适应性函数。

 

药物设计常用方法分子对接方法

一：概述
  分子对接是指两个或多个分子通过几何匹配和能量匹配相互识别的过程，在药物设计中有十分重要的意义。药物分子在产生药效的过程中，需要与靶酶相互结合，这就要求两个分子要充分接近并采取合适的取向以使二者在必要的部位相互契合，发生相互作用，继而通过适当的构象调整，得到一个稳定的复合物构象。通过分子对接确定复合物中两个分子正确的相对位置和取向，研究两个分子的构象特别是底物构象在形成复合物过程的变化是确定药物作用机制，设计新药的基础。
  分子对接计算把配体分子放在受体活性位点的位置，然后按照几何互补、能量互补以及化学环境互补的原则来评价药物和受体相互作用的好坏，并找出两个分子之间最佳的结合模式。由于分子对接考虑了受体结构的信息以及受体和药物分子之间的相互作用信息，因此从原理上讲，它比仅仅从配体结构出发的药物设计方法更加合理。同时，分子对接筛选的化合物库往往采用的是商用数据库，比如可用化合物数据库（ACD）、剑桥晶体结构数据库（CSD）、世界药物索引（WDL）、药用化合物数据库（CMC）以及可用化合物搜索数据库（ACDSC）等等，因此筛选出来的化合物都为已知化合物，而且相当大一部份可以通过购买得到，这为科研提供了很大的方便，近年来，随着计算机技术的发展、靶酶晶体结构的快速增长以及商用小分子数据库的不断更新，分子对接在药物设计中取得了巨大成功，已经成为基于结构药物分子设计中最为重要的方法。
  分子对接的最初思想源自于“锁和钥匙”的模型，即“一把钥匙开一把锁”。不过分子对接，也就是药物分子和靶酶分子间的识别要比“钥匙和锁”的模型要复杂的多，首先表现在药物分子和靶酶分子是柔性的，这样就要求在对接过程中要相互适应以达到最佳匹配；再者，分子对接不仅要满足空间形状的匹配，还要满足能量的匹配，底物分子与靶酶分子能否结合以及结合的强度最终是由形成此复合物过程的结合自由能的变化值决定。互补性和预组织是决定分子对接过程的两个重要原则，前者决定识别过程的选择性，而后者决定识别过程的键和能力。互补性包括空间结构的互补性和电学性质的互补性。受体和底物分子在识别之前将受体中容纳底物的环境组织的愈好，其溶剂化能力就越低，则它们的识别效果愈佳，形成的复合物越稳定，这就是分子识别的预组织原则。

 二：分子对接的种类
 1.刚体对接：指在对接过程中，研究体系的构象不发生变化。适合考察比较大的体系，如蛋白质和蛋白质间以及蛋白质和核酸之间的对接。
 2.半柔性对接：指在对接过程中，研究体系尤其是配体的构象允许在一定的范围内变化。适合处理大分子和小分子间的对接，对接过程中，小分子的构象一般是可以变化的，但大分子是刚性的。
 3.柔性对接：指在对接过程中，研究体系的构象基本上可以自由变化的。一般用于精确考虑分子间的识别情况。由于计算过程中体系的构象可以变化，所以计算耗费最大。
   三：分子对接中的重要问题
  分子对接的目的是找到底物分子和受体分子间的最佳结合位置，所以要面对的重要问题是如何找到最佳的结合位置和如何确定对接分子间的结合强度？
  如何找到最佳的结合位置牵涉到优化的问题。底物分子和受体分子都是可以自由转动和平动的，同时两个分子自身的构象也存在变化，因此它们之间可能的结合方式是非常复杂的，所以简单的系统搜索方法是不够的，要引入其他高效的优化方法，常用的有遗传算法、模拟退火以及禁忌搜索等。
  如何确定对接分子间的结合强度涉及到底物分子和受体分子间结合能力的预测，牵涉到结合自由能的计算。结合自由能包括以下几个方面的贡献:1 受体分子和底物分子气态下分子对接过程的自由能变化，约为对接过程中的函变;2 受体分子、底物分子以及复合物分子的溶剂化自由能；3 对接过程中的熵变。在这几项中，气态下分子对接过程的函变可以通过分子力学的方法简单求算，但去溶剂化能准确而快速的求算还存在一定的问题，不过熵变的计算可能是最大的问题，因为它的计算需要耗费大量时间，而在实际的药物设计过程中，研究人员总是希望能快速筛选成千上万的分子。所以目前采用的是较为简单的自由能评价方法。

四：几种有代表性的分子对接方法
 4.1 DOCK
  DOCK是Kuntz研究小组发展的分子对接程序，可能是目前应用最为广泛的分子对接程序之一.它能自动地模拟配体分子在受体活性位点的作用情况，并把理论预测最佳的方式记录下来。而且该方法能够对配体的三维结构数据库进行自动搜索，因此被广泛应用于基于受体结构的数据库搜索的药物设计中，并取得了巨大的成功。用DOCK进行药物设计以及数据库的搜索基本上可以分为下面几个步骤:配体和受体相互作用位点的确定，评分系统的生成，DOCK计算及DOCK结果的处理与分析。活性位点的确定和表达是DOCK最重要的特点之一。活性位点特征的确定对于DOCK研究是非常重要的，因为配体分子和受体相互作用过程的模拟主要就是参考几何位点的几何特征进行的。在DOCK中，活性位点的确定通过sphgen程序来完成。DOCK软件包中sphgen程序生成受体表面所有的凹陷的负像，并对这些负像进行聚类分析。在DOCK程序中，表面点采用了Richards提出的模型。在这些表面点的基础上，采用sphgen程序生成了负像，它实际上由一些与分子表面点相切的圆球叠加而成。
  在生成负像的基础上，就可以进行配体分子和活性口袋之间的匹配。在这里，配体也
 采用一组球集来表示，和负像不同的是，配体所用的球集表示配体所占的空间区域。如果
 配体分子能和活性口袋形成比较好的匹配，那么配体的球集一定能和活性口袋中的负像形
 成好的叠合。配体分子和负像之间的匹配原则是基于配体和受体之间球集的内坐标的比较。
 按照匹配原则得到了配体和受体之间的匹配情况之后.就要通过合理的得分函数来选择最优的结果。DOCK提供了多种得分函数来评价配体和受体之间的结合情况，包括原子接触得分以及能量得分。
  DOCK进行分子对接时，配体分子可以是柔性的。对于柔性的分子，其键长和键角保持
 不变，但可旋转二面角是可以发生变化的。在DOCK中，柔性分子的构象变化通过下面的
 操作实现:首先是刚性片断的确定，然后是构象搜索。构象搜索采用两种方法:一种是锚优
 先搜索(anchor-first search),第二种方法是同时搜索(simultaneous search).

 4.2 AUTODOCK
  AUTODOCK是Scripps的Olson科研小组开发的分子对接软件包，最新的版本为3.05，AUTODOCK采用模拟退火和遗传算法来寻找受体和配体最佳的结合位置，用半经验的自由能计算方法来评价受体和配体之间的匹配情况。在AUTODOCK中，配体和受体之间结合能力采用能量匹配来评价。在1.0和2.0版本中，能量匹配得分采用简单的基于AMBER力场的非键相互作用能。非键相互作用来自于三部分的贡献:范得华相互作用，氢键相互作用，以及静电相互作用。在3.0版中，AUTODOCK提供了半经验的自由能计算方法来评价配体和受体之间的能量匹配
  在最早的AUTODOCK版本中，作者采用了模拟退火来优化配体和受体之间的结合。在
 3.0版本中，Morris等发展了一种改良的遗传算法，即拉马克遗传算法(LGA)。测试结果表
 明，LGA比传统的遗传算法比模拟退火具有更高的效率。在LGA方法中，作者把遗传算法和局部搜索(local search)结合在一起，遗传算法用于全局搜索，而局部搜索用于能量优化。在AUTODOCK中，局部搜索方法是自适应的，它可以根据当前的能量调节步长大小。LGA算法引入了拉马克的遗传理论，LGA最大的特点就是通过进化映射(developmental mapping)把基因型转化为表现型而实现局部搜索和遗传算法的结合。基因型空间通过遗传算子突变和交叉来定义;而表现型则通过问题的解来定义，这里表示体系的能量得分。

 4.3 FIexX
  FIexX 是德国国家信息技术研究中心生物信息学算法和科学计算研究室的Matthias
 Rarey等发展的分子对接方法，现在己经作为SYBYL分子模拟软件包中的一个模块实现了商业化。FIexX中结合了多种药物设计的方法进行配体和受体之间的对接。在FIexX中，配体和受体之间结合情况的评价采用了类似Bbhm提出的基于半经验方程的自由能评价方法。在FlexX中，分子对接的流程主要分为下面的步骤:
 (1)核心结构确定
  对接的第一步为核心片断(base fragment)的选择。核心片断是指能对配体和受体之间
 相互作用起决定作用的基团，而且核心片断的构象要尽量少一些。核心基团的正确选择对分子对接的计算结果有非常重要的影响。因为如果核心基团和受体之间不存在明显优势的相互
 作用时，则很难正确预测正确的结合模式。随着核心基团的增加，则核心基团和受体之间的
 相互作用也会相应增强，那么结合模式准确预测的机会就会大大增加。因此在选择核心基团时，核心基团包含的基团要尽量多一些。而且核心基团的构象数要尽量少一些。当核心基团选定以后，就可以把配体分子划分为多个片断。
 (2)核心结构的放置
  当选择好了核心基团以后，就要把核心基团放置在活性位点的正确部位。在放置核心基
 团的时候，FlexX采用了一种形态聚类算法(pose clustering algorithm)算法)，在这个算法中，一个核心基团可以看作为一个具有明确相互作用点的刚性物体，而受体分子的活性口袋也可以看作为一个具有明确相互作用点的刚性物体。把核心基团放置在活性口袋中的过程就相当于把配体中的二个相互作用点叠合在活性位点的三个相互作用点上(假设这三个点不共线)。在匹配中，两个三角形三个顶点所具有的相互作用特征应该是符合的，同时三角形对应边长的差别应该在一定的范围内。放置核心结构的第一步就是找出所有相匹配的这些三角形，而且对配体的位置进行坐标转化。当所有的转化完成以后，检查配体是否和受体产生了碰撞，去掉一些不合理的核心结构取向。对于得到这些核心结构的可能位置，通过核心结构空间的位置均方根位移(r.m.s.D )进行聚类分析。对应那些r.m.s.D值小于一定阀值的空间位置进行归并，仅仅保留那些相差较大的空间位置。最后，检查这些核心结构和受体之间的相互作用情况，对结构进行简单的修饰.
 (3)配体分子的生长
  当核心结构在活性口袋中的位置确定以后，配体分子的其它部分可以分为小的片断，依
 次“生长”在核心结构上。片断生长采用树形搜索(tree search)的方法，算法和SYBYL中的系统搜索所采用的算法基本类似，在搜索的过程中要尽量删除那些无用的分支。搜索树的第一层是核心结构在活性口袋中的不同放置位置。在下一层中，片断要采用尽可能多的形式连接到核心结构上。然后，按照树形结构，片断依次连接。如果片断能和受体形成氢键或盐桥，则优先连接，因为相互作用形式越明确，片断的几何定位越容易。在对接过程中，一个完全无遗漏的树形结构在造作上是很困难的，因此对于树的每个节点，我们仅仅只考虑最佳的k种结果，这样会为下面的生长节约所要耗费的计算开支。k种最佳的结构采用上面介绍的得分函数来选择。
  增加新的基团以后，如果发现新的相互作用，或配体和受体之间存在重叠，则需要对配
 体的位置进行优化。FIexX采用点的加权叠合法(weighted superposition of points)来进行优化。
  当配体生长结束以后，我们可以得到k个最佳的配体和受体的结合形式，从中用户可以
 选择需要的结果。FIexX的分子对接过程中，配体的生长以及柔性的考察和DOCK中采用的锚优先方法基本类似。FIexX己经实现了商业化，而且程序可以自动地对数据库中的多个分子进行分子对接的计算，然后给出最佳的结果。同时FIexX提供了友好的图形界面，易于操作，因此可能在药物设计中具有较好的应用前景。

 4.4 Affinity
  Afinity是Accelrys (MSI)和杜邦联合开发的分子对接方法，也是最早实现商业化的分子对接方法。Afinity中提供了多种分子对接的策略，可以根据用户的需要提供多种方法的组合。在Afinity中，分子对接可以大致分为两个步骤:首先通过蒙特卡罗或模拟退火计算来确定配体分子在受体活性口袋中可能的结合位置:然后，在第一步的基础上，采用分子力学或分子动力学方法进行进一步细致的分子对接。和其它分子对接方法比较，Affinity具有自己的特色。首先，Affinity中提供了多种对接方法的结合，比如蒙特卡罗方法和分子力学、分子动力学以及模拟退火方法的结合;这些方法结合的灵活性为多种分子对接问题提供了解决方案。第二，在Affinity中，不仅仅配体是柔性的，受体的重要部位，比如活性位点中的某些残基也可以定义为柔性的区域。第三，Affinity提供了精确和快速计算配合和受体之间非键相互作用的两种有效方法，一种是基于格点的能量计算方法，而另一种则是单元多偶极(cell multipole method)方法。第四，Afinity采用了Stouten提出的溶剂化模型来考察配体和受体在堆积过程中溶剂化能的变化。
  在Affinity中，配体和受体之间匹配主要采用能量得分的评价方式。对于能量得分，Afinity采用两种力场，即CVFF力场和CFF力场。Afinity中提供了基于格点的能量计算方法
 以及不基于格点的能量计算方法。对于这两种不同的方法，其能量评分函数略有区别。在基于格点的能最计算中，可以考虑溶剂效应的影响。
  Afinity提供了一套多种方法相结合的分子对接流程，整个流程都是自动化的，不需要用户给予任何的千预。在这个流程中，采用蒙特卡罗的方法对每个优化后的结构进行取样，得到配体在尽量多样性的空间取向。在Afinity中，对接流程不是固定的，用户可以按照自己的需要设计相应的流程。比如，用户可以把蒙特卡罗取样和分子动力学模拟结合起来，用蒙特卡罗模拟先得到若干个配体在活性口袋可能的对接位置，然后对这些对接结构用分子动力学进行更加细致的采样，最后对从不同起点得到的分子动力学轨迹进行分析来得到正确的结果。
 Afinity方法采用了多种方法的结合为用户解决不同的问题提供了不同的解决方案。Afinity适合对配体和受体之间的相互作用模式进行精细地考察，但它不太适合对大量的配体分子进行基于分子对接的虚拟筛选，因为Affinity对配体和受体都采用了柔性的策略，需要消耗较大的计算量。

 

五：软对接的实现过程
  把受体分子和配体分子放在空间中的任意位置，想通过软对接找到可能的结合构象，需要采用下面的操作步骤:
 1. 把两个分子中的大分子作为目标分子,小分子作为探针分子，分别产生目标分子和探针分子的溶剂可极性表面。分子溶剂可及性表面表征方法把分子的表面描绘成一个很平滑的空间，它较好地表征了分子在溶液中和溶液的接触情况。在分子对接中，采用溶剂可及性表面是很合适的，平滑的分子表面使对接过程操作起来比较容易.可及性表面点矢量的匹配性可以很好地描述底物分子和靶酶分子的局部匹配情况。
 2. 计算目标分子和探针分子的坐标重心，把目标分子坐标重心作为坐标系的零点。把探针分的坐标重心平移到坐标系的零点位置作为对接的起点计算两个分子的几何尺寸，确定搜索空间。搜索空间用一个矩形来定义，它包含整个受体分子或配体分子。在对接的过程巾.探针分子的坐标重心到目标分子的坐标重心之间的距离不能超过二者半径之和。
 3. 上面的对接范围定义适用于进行全局搜索的情况。在对接过程中，也可以把探针分子
 的运动定义在特定的空间范围内，使得探针分子只能在特定空间区域(比如活性位点)和目标分子对接。这样做的目的是为了约束探针分子的运动范围，而只在特定的范围进行局部搜索
 4. 探针分子转动和平动六个自由度去对接目标分子，并对每种结合情况进行评估。评估采用上面介绍的表面匹配得分。
 5. 采用优化方法来优化分子对接的过程，得到最佳的分子对接模式。在SFDOCK程序中，可以采用不同的优化方法，包括单纯型法、蒙特卡罗模拟退火方法以及遗传算法。
 6. 对于优化得到的结果，保留表面匹配最佳的且在空间上具有较大差别的构象。在这里，
 我们采用坐标均方根位移(r.m.s.)来评价两个探针分子在空间取向上的差别。对于得到的
 基于表面匹配的最佳分子对接模式，可以采用能量匹配进行进一步的分析。

  所谓柔性对接是主要是指在分子对接过程中，底物和受体的构象是可以允许发生变化的。在SFDOCK方法中，目前我们仅仅允许配体小分子的可旋转二面角发生变化。因此，柔性对接和软对接比较，需要优化的变量个数是不同的。软对接只有六个变量: 三个平动自由度，三个转动自由度。而柔性对接除了这六个自由度以外，还包括了底物分子的部分二面角变量.至于柔性对接的优化和软对接的优化方法则基本相同。具体的计算步骤也基本相似。

 六：分子对接过程中的优化方法
  分子对接过程实际上是一个复杂的优化过程，它的复杂之处可以从两个方面来考虑。首先，对接时受体分子和底物分子的构象会发生一定程度上的变化，同时考虑底物和受体的构象变化是很困难的，软对接只是部分解决了这个问题，对于复合物前后构象变化较大的情况仅仅采用软对接是无法得到好的结果的，在这种情况下只能采用柔性对接。但考虑到计算效率的问题，因此在一般情况下，我们实际上仅仅只考虑小分子的柔性变化，而不考虑靶酶构象的变化(靶酶构象的变化相对底物分子而言变化不大)。其次，在分子对接时，构象空间中存在大量的局部极小，并且对接目标函数在某些区域常常是不连续的，即使在结合位点附近，采用常见的梯度优化方法，比如最陡下降法、共扼梯度法以及牛顿叠代法一般是很难得到最佳解的。解决此类问题的最佳方法是采用一些启发式的随机方法如Monte Carlo方法、模拟退火方法以及遗传算法等。
 下面我们简单介绍一下在SMOCK采用遗传算法时的优化过程，它分为下面的步骤:
 产生初始种群:首先随机地产生初始种群中的若干个个体，每个个体代表一种对接情况。在软对接中，个体由六个数的一维数组组成，这六个数代表探针分子的平动和转动值。在柔性对接中，个体中除了六个探针分子的平动和转动值以外，还包括用户定义的探针分子的可旋转单键的二面角数值。产生了初始种群后，就可用得分函数来来评价每个个体。在遗传算法的优化过程中，目标分子和探针分子之间的表面匹配得分或能量匹配得分为遗传优化的得分函数。
 选择操作:种群中的所有个体被评价后，就可以根据种群中个体的得分结合随机方法来选择被新种群保留的个体。
 交叉操作:选择操作结束后，就可以对种群中的个体进行交叉操作，用新个体替代被淘汰的个体。交叉操作是这样进行的:如果一个随机数小于交叉几率，则在种群中随机地选择两个个体作为母体，然后将这两个母体随机地分为两段，并将相应的部分进行交叉换位，得到两条新链。在交叉操作的过程中，可以根据需要进行单点交叉或多点交叉。
 突变操作:突变操作是这样进行的。先产生一个随机数，如果此数小于突变几率，则在种群中随机地选择一个个体进行突变操作，突变操作时在这个个体中随机产生一个突变位点，并将此位点的数值用一个随机数代替。通过交叉和突变而产生的新的种群中的所有个体均要用得分函数来予以评价。
 比较操作:在遗传算法操作过程中，为了将最好的个体标记下来，我们用了一个“精华”种群来保存它们，在进行每次遗传操作后，逐一比较新种群中的个体和“精华”种群中的个体，如果新种群中存在更好的个体，就把它们拷贝到“精华”种群中去。在比较的过程中，我们要求“精华”种群的个体空间取向上差别尽可能大。具体比较操作如下:如果种群中的某个个体比“精华”种群的所有个体都具有更好的得分，而且它和“精华”种群中的个体之间的坐标r.m.s.满足预先定义的条件，则用这个个体来替代“精华”种群中的最差个体。如果某个个体比“精华”种群中的所有个体都具有更好的得分，但它和某个个体之间的坐标r.m.s.较小，则用这个较好的个体来替代精华种群中的这个较差的个体。
 考虑到遗传算法在局部极值附近会振荡较长时间，为了帮助遗传算法脱离局部极值，每次交又和突变操作后。把种群中的若干个得分最差的个体用随机产生的个体替代。在经过足够长的优化过程后，当精华种群的个体不再发生变化时，计算收敛。