欢迎光临！

此项目（包括此网站）由社区管理。Github地址点这里。

此项目创始人：Abdusemi

此项目中的资料来自中国药科大学孙慧涌、邹毅等老师的CADD课程。

内容：

计算机药物设计绪论(Introduction to Computer-Assisted Drug Design)

参考书目：

• 计算机辅助药物设计

• Comprehensive Medicinal Chemistry II Volume 4

• Computer-Assisted Drug Design

相关期刊：

• Journal of Chemical Information and Modeling (JCIM)

• Journal of Chemical Theory and Computation (JCTC)

• Journal of Computational Chemistry (JCC)

• Journal of Cheminformatics

• PLoS Computational Biology (PLoS CB)

• Briefings in Bioinformatics (BIB)

• Bioinformatics

• Nucleic Acids Research (NAR)

• Nature Machine Intelligence

• Journal of Medicinal Chemistry (JMC)

• Journal of American Chemical Society (JACS)

• Nature, Science, Nature Communications, etc

计算机辅助药物设计的发展史

CADD技术诞生于药物研发的发展过程中。

药物研究与开发的历史，是个由粗到精，由盲目到自觉，由经验性的试验到科学的合理设计的过程，大致可以分为3个阶段：发现阶段(Discovery)，发展阶段(Development)，设计阶段(Design)。

药物研发的早期发现阶段：

• 偶然性大，缺乏目标

• 对毒副作用不够重视

• 对疾病的发病机制特别是分子层面的机制缺乏足够的认识

药物研发的发展阶段：

• 药物化学的发展阶段大致是在20世纪30年代到60年代。合成药物的大量涌现，是药物发展的“黄金时期”

• 分子药理学形成和酶学的发展，对阐明药物的作用原理起了重要的作用

药物研发的设计阶段：

• 一方面，爆发式增长过后，药物研发陷入低潮，系统性疾病的复杂性导致研究人员难以通过传统的方法获得有效的药物。研发经费激增，而成功率骤降

• 另一方面，“反应停”事件的出现导致了人们对于药物安全性的日益重视，也致使药物的研发难度大幅上升

  反应停事件：1959年，西德各地出生过手脚异常的畸形婴儿。伦兹博士对这种怪胎进行了调查，于 1961年发表了“畸形的原因是催眠剂反应停”。反应停是治疗女性妊娠时早期孕吐的一种药物。

• 随着数据的积累,比如：靶标数目的增加、化合物实体库规模增长、复杂相互作用网络，常规实验手段不能满足指数增长的化合物的筛选和活性评价

• 药物设计的诞生得益于分子生物学、结构生物学、有机化学、物理化学、计算机科学等诸多学科的发展，其本身充分体现了学科交叉的特征

计算机辅助药物设计的概念

计算机辅助药物设计（CADD）是以计算机作为操作媒介，利用计算化学、计算生物学、分子图形学、数理统计以及数据库等技术，研究药物和受体的相互作用，以期发现、设计和优化创新药物分子的方法学集合，其技术核心在于分子模拟，研究目的在于合理药物设计.

分子模拟是CADD的技术核心与实现手段

分子模拟是利用计算机以原子（量子）尺度模型来模拟分子的结构与行为，进而模拟分子体系的各种物理、化学性质的方法。它是在实验数据基础上，通过基本物理、化学原理，构筑的一套数理模型和算法。《中国大百科全书》

简单的说：计算理论和计算机图形学的结合

分子模拟的历史

分子模拟的思想溯源于上世纪初量子力学的兴起：

• 1925年，Heisenberg发表了第一篇现代量子力学的文章

• 1926年，Schrodinger发表了著名的波动方程

• 1927年，采用量子力学的方法计算了氢分子的轨道，量子化学由此诞生

• 1930年，D. H. Andrews提出分子力学的基本思想

• 1933年，J. D. Bernal和R. H. Flowler提出了水的原子模型

• 1946年，Frank Westheiner完成了第一次分子力学计算

分子模拟的理论基础

1964年，加州大学的沃尔特·科恩（Walter Kohn）教授和霍恩伯格（Hohenberg）博士从薛定谔方程出发，严格证明了一个重要性质:分子在基态的能量及其它所有性质可以由基态的电子密度的分布得到。即分子结构决定分子性质

Walter Kohn获得1998年诺贝尔化学奖

分子模拟的目的就是要用理论方法（模拟分子结构） 去实现以往用实验才能证明的东西（推导分子性质）

分子模拟的发展及意义

• 1974年，Allinger等人发布了分子力场MM

• 1981年，Peter Kollman发布了AMBER力场

• 1983年，Martin Karplus发布了分子模拟程序CHARMM

• …

分子模拟曾经限制在部分科学家手中，现在随着理论和计算机技术的发展，任何有兴趣的人都可以应用。而且应用范围越来越广，尤其在生命科学、药学、高分子、材料科学等多个领域取得了重要成果。21世纪也正好是这些领域大展身手的世界，分子模拟的作用也将更加重要。

分子模拟的计算理论

根据基本原理的不同，分子模拟可分为量子力学模拟和分子力学模拟

	考虑电子运动状态		考虑核运动状态（电子运动作近似假设）
量子力学模拟	用波函数表示	分子力学模拟	用质点运动轨迹表示
	解薛定谔方程		解力学方程（经典牛顿力学方程）

量子力学模拟（QM、QM/MM）

分子力学模拟

• 全原子模拟（all-atom simulation）

• 粗粒化模拟（coarse-grained simulation）

介观模拟（mesoscopic simulation）

连续介质模拟（continuum simulation）

From the thousand-atom size scale, in the early 1990s, to a full protocell, on the billion-atom size scale,nowadays!!!

合理药物设计是CADD的研究目的

合理药物设计（rational drug design）是指在分子病理学的基础上，依据靶点（target）的结构信息，发现和设计与之契合的药物分子，并充分考虑药代动及选择性等性质的药物设计过程。

合理药物设计包括：

• 基于结构的药物设计（sturucture-based drug design，SBDD）

• 基于性质的药物设计（property-based drug design，PBDD）

• 基于机理的药物设计（mechanism-based drug design，MBDD）

靶点（target）

靶点主要指受体，广义上的受体包括酶、离子通道、膜蛋白受体、抗原-抗体、核酸等

基于结构的药物设计（SBDD）

基于结构的药物设计包括：

• 基于受体结构的药物设计(target structure-based drug design, TSBDD)

  直接药物设计方法

  借助X射线晶体衍射、核磁共振或同源模建获得受体或受体-配体复合物的三维结构数据，经计算机图形学再现并设计药物分子

  特点：目标清晰，设计准确，是首选的方法之一

• 基于配体结构的药物设计(ligand structure-based drug design, LSBDD)

  间接药物设计方法

  在受体三维结构不确定的情况下，从一系列已知活性分子出发（作用于同一靶标），分析其结构变化与生物活性强弱的关系，在计算机的辅助下，找出药效团，据此特征设计药物分子

  特点：效率低，选择性差，结构新颖性不佳

计算机辅助药物设计的应用 The application of CADD

CADD的应用范围

在现代药物研发过程中，CADD技术已被广泛应用于各个环节：

• 靶点识别（bioinformatics）

• 生物靶标的结构特性表征

• 靶标可靶性评估

• 苗头化合物发现（hit）

• 药物-靶标相互作用特征分析

• 苗头分子结构优化与改造（hit to lead）

• ADME/T性质的优化

• 临床分析（AI）

CADD的应用环境

software: Algorithm and Logic

hardware: Computation and Visualization

分子模拟的常用软件

常用的综合软件平台:

• Discovery Studio

• Schrodinger

• MOE

• SYBYL

常用的开源软件:

• Modeller，免费（https://salilab.org/modeller/）

• Autodock，免费（https://www.autodock.scripps.edu/）

• AMBER（Assisted Model Building with Energy Refinement）

• GROMACS，免费（https://www.gromacs.org）

• CHARMM（Chemistry at HARvard Macromolecular Mechanics）

• NAMD，免费（https://www.ks.uiuc.edu/Research/namd/）

CADD的技术组成:

• 分子力学（Molecular Mechanics，MM）

• 分子动力学（Molecular Dynamics，MD）

• 分子对接（Molecular Docking）

• 同源模建（Homology Modeling）

• 量化计算（Quantum Mechanical Calculation，QM）

• 定量构效关系（Quantitative Structure-Activity Relationships，QSAR）

• 药效团（Pharmacophore）

• 人工智能（Artificial Intelligence，AI）

• 数据库（Database）

分子力学（Molecular Mechanics）：分子力学是分子模拟的基础之一，是连接微观机制与宏观表象的桥梁。

分子动力学（Molecular Dynamics）：分子动力学是分子力学的应用拓展，其可进一步从时间尺度把握分子作用的微观机制。

分子对接（Molecular Docking）：分子对接是分子力学的又一重要应用（基于力场的分子对接），其可快速识别可与目标靶点相结合的配体分子，是基于结构药物设计中最为重要的方法之一。

同源模建（Homology Modeling）：同源模建通常以已知同源结构为模板，通过分子力学优化获得目标序列结构。

量化计算（QM Calculation）：量化计算是分子模拟的又一重要分支，通过对电子波函数的精确描述，可有效把握化学反应的分子机制，是分子模拟中最为精确的理论研究方法。

定量构效关系（QSAR）：定量构效关系（quantitative structure-activity relationships），是研究一组化合物的活性、毒性、药代性质与其结构（structural）、物理化学性质（physicochemical）、拓扑结构（topological）之间的相关关系，并用数理统计模型加以表征的研究方法。

药效团（Pharmacophore）：药效团是产生特定药理作用所必须的物理化学特征及其在空间的分布，其是一组离散的物理化学特征在空间特定位置的分布，亦属3D-QSAR的一种。

人工智能（AI）：基于人工智能的药物设计是近年最为热门的研究领域之一，其亦属广义QSAR模型。通过对已知数据建立复杂非线性映射关系来有效预测未知数据属性。

数据库（Database）：数据库很多时候作为分子模拟数据来源，通过数据共享，极大节约时间、空间、资源/劳力成本，在CADD中发挥着极为基础的作用。

CADD中几个基础概念 Basic Concepts in CADD

分子模拟中常用度量单位

• 距离单位：埃（Å，10^-10m）、纳米（nm，10^-9m）

• 角度单位：度（degree，°）、π（π = 180°）

• 时间单位：飞秒（fs，10^-15s）、皮秒（ps，10^-12s）、纳秒（ns，10^-9s）、微秒（µs，10^-6s ）、毫秒（ms，10^-3s）

• 能量单位：kcal/mol、kJ/mol、hartree（a.u.）、eV（1 hartree = 27.2 eV = 627.5029 kcal， 1 kcal = 4.184 kJ）

• 温度单位：K（通常分子模拟温度为298.15 K或310 K）

坐标体系（Coordinate System）

笛卡尔坐标（Cartesian Coordinates）即用原子的xyz值来表征分子构型。 如：甲烷分子的笛卡尔坐标表示法

内坐标（Internal Coordinates）用原子之间的距离、角度和二面角表征分子构型。通常用Z矩阵（Z-matrix）表示

两种坐标体系的比较：笛卡尔坐标一般用于定义含有大量原子的体系，如蛋白质、DNA等。对于一个含有n个原子的体系而言，其坐标数量为3N；笛卡尔坐标被广泛应用于分子力学模拟之中。分子内坐标一般用于定义含有较少原子的体系，如有机小分子化合物。对于一个含有n个原子的体系而言，其坐标数量为3N-6；分子内坐标一般应用于量子力学模拟之中。

常见的分子存储格式：

• MDL公司的SDF（MOL）格式（无子结构信息、无原子类型、无部分电荷）

• Tripos公司的MOL2格式（有子结构信息、有原子类型、有部分电荷）

• 蛋白晶体库的PDB（ENT）格式（有子结构信息、无原子类型、无部分电荷）

• Accelrys公司的SMI格式（SMILES，仅字符串形式，无坐标相关信息）

目前使用的分子格式中，pdb格式广泛应用于生物大分子文件的存储（如蛋白质或核酸）

mol2格式常用于存储小分子信息，亦可用于生物大分子文件的存储（不常见）

其他形式的文件（sdf、mol、smi等）均用于小分子文件的存储（缺乏子结构定义形式）

SDF（MOL）文件格式（苯分子）：

MOL2文件格式（ALA氨基酸）：

分子电荷形式：

• 形式电荷（formal charge）：分子整体带电量（-1）

• 部分电荷（partial charge）：又名原子电荷（atomic charge）

PDB（ENT）文件格式（ALA氨基酸）：

SMI文件格式：分子线性输入规范（Simplified Molecular Input Line Entry System，SMILES）

• 计算效率极高，便于存储

• 使用传统化学符号（B, C, N, O, P, S, F, Cl, Br, I）

• 分子表述忽略H原子

SMILES键型

CC	ethane	(CH3CH3)
C=O	formaldehyde	(CH2O)
C=C	ethene	(CH2=CH2)
O=C=O	carbon dioxide	(CO2)
COC	dimethyl ether	(CH3OCH3)
C#N	hydrogen cyanide	(HCN)
CCO	ethanol	(CH3CH2OH)
[H][H]	molecular hydrogen	(H2)

SMILES成环规则:

• 脂肪族或非芳香碳：C

• 芳香碳：c

• 成对数字表述闭合环结构 如c1ccccc1（或C1=CC=CC=C1）表示苯，而 C1CCCCC1表示环己烷

• SMILES金属 [Al] [As] [Au] [Be] [Bi] [Cd] [Ca] [Fe] [Hg] [K] [Li] [Mg] [Na] [Ni] [Pt] [Sb] [Sn] [Zn] [Zr]

异构体与手性:

• 异构体原子使用“/”和“”，如：反式1,2-二氟乙烯：F/C=C/F 顺式1,2-二氟乙烯：F/C=CF

• 手性原子使用“@”

SMILES在人工智能模型中的应用：由于SMILES格式的文本格式，其在基于人工智能的分子生成模型中扮演着十分重要的角色

分子表面：

• 溶剂可及表面（solvent accessible surface）：指探针分子在目标分子表面进行“假想滚动”时，其球心形成的轨迹（并非真实表面）。通常使用半径为1.4 Å的水分子作为探针分子。探针分子的中心可以放置于可及表面的任一点上，但不能穿入分子中任何原子的范德华球。

• 范德华表面（vdW surface）：以原子的范德华半径生成的分子表面，由于缺少探针半径弥补，范德华表面通常更加离散

分子表面比较：

溶剂可及表面（solvent accessible surface）

范德华表面（vdW surface）

常用视图软件:

• PYMOL

• VMD

分子力学基本原理(Molecular Mechanics)

分子力学的概念及发展(The concept and development of molecular mechanics)

分子力学的起源：1930，D. H. Andrews

• 在分子内部，化学键都有“自然”的键长值和键角值。分子要调整它的几何形状（构象），以使其键长值和键角值尽可能接近自然值，同时也使非键作用（van der Waals）处于最小的状态，从而给出原子核置的最佳排布，即分子的平衡构型。

参考: Andrews, D.H. The relation between the Raman spectra and the structure of organic molecules, Phys. Rev. 1930, 36, 544

分子的经典力学模型：1946，T. L. Hill

• Hill提出用van der Waals作用能和键长、键角的形变能来计算分子的能量，以优化分子的空间构型

• Hill指出：“分子内部的空间作用是众所周知的，（1）基团或原子之间靠近时则会相互排斥；（2）为了减少这种作用，基团或原子会趋于相互离开，但是这将使键长伸长或键角发生弯曲，又会引起相应的能量升高。最后的构型将是这两种力折衷的结果，并且是能量最低的构型”

现代分子力场：1976， Norman L. Allinger

• 1976年，Allinger发布了MM1力场，1977年发布了MM2力场

虽然分子力学的思想和方法在19世纪30-40年代就已经建立，但是直到60-70年代以后，随着电子计算机的发展，用分子力学来确定和理解分子结构和性质的研究才越来越多。随着并行计算和超级计算机的普及，分子力学方法已经广泛应用于多个研究领域，如有机化学、生物化学、药物设计、材料设计等

分子力学（Molecular Mechanics）是使用基于原子核坐标的势函数和一系列相关势参数描述体系能量的方法

由于势函数和参数都是由经验和实验值拟合获得，故亦称经验力场方法（empirical force field method）

分子力学方法的建立和正确使用在于它的背景和若干假设

在分子模拟中许多问题涉及较大体系（如生物大分子、高分子体系等），无法用量子力学方法求解。其原因在于量子力学计算涉及体系的电子，以至于计算量远超计算机承受范围

分子力学方法忽略电子的运动，只计算以原子核的位置为函数的体系能量，因此可用于计算包含大量原子的体系 由于模拟参数的精确拟合，分子力场有时甚至可以获得近乎量子力学的精确解，而所需的模拟时间仅为量化计算的几十分之一甚至几百分之一，极大的提高了计算效率

注意：分子力学不能提供分子中依赖电子分布的性质

分子力学基本假设：

Born-Oppenheimer（波恩—奥本海默）近似假设（BO近似）

• 要点：核的振动运动比电子的运动慢得多，因此可近似地假设核的运动不影响电子的运动，他们是独立的

• 具体讲，∵M（质量）核＞＞M电子， V（速度）核＜＜V电子，

  ∴ 在解电子运动方程时可以近似地认为核的构型（位置）保持不变，即可把电子运动与核运动分开处理，核在电子的平均场中运动

• 显然，现实中不存在BO近似，把能量完全写做原子核坐标的函数并不准确

分子力学的模型表示形式：

分子模型假设：

• 球表示原子

• 弹簧表示化学健

• 弹簧的形变表示化学键的伸缩、弯折和扭转

• 非键结合的原子通过范德华和静电相互作用表示

分子力学基本特点：

• 分子力场不仅描述结构，也描述性质（结构决定性质）

• 一个完整的力场包括一套势函数以及相关的力场参数

• 不同的力场可以有相同的函数形式，但不同的参数

• 不同力场的参数（即使描述相同对象）不能混用

• 注意：力场都是经验的（相对从头计算而言），即没有绝对的正确与错误之分；分子力场参数都是拟合特定分子的数据而生成，因此只能说某个力场更适用于某些体系

计算的能量代表什么？

一个分子力学模型能简单地通过增加在所有键上的张力和所有原子间的范德华和库仑相互作用给出一个“能量”值。这个量称作E_MM，即分子力学能（也可称作空间（位）能）

1. E_MM最接近反映真实事物性质的是一个分子的内能：U标准的热力学方程，使U与其它量有如下关系：

\[H = U + PV （H: 焓、P: 压力、V: 体积）\]

\[\Delta H = \Delta U + P \Delta V\]

2. 对一个简单分子力学模型没有外部压力（真空中），有：

\[\Delta H ≈ \Delta U ≈ \Delta E_MM\]

然而，自由能的变化（ΔG）与焓变（ΔH）和熵变（ΔS）有关，如下：

\[\Delta G = \Delta H - T \Delta S (T：绝对温度、 \Delta S：熵变)\]

3. 如果能使一个过程熵变足够小（ΔS ≈ 0），那么，

\[\Delta E_MM ≈ \Delta G\]

即：分子力学能量的变化ΔE_MM可以合理的近似等于自由能的变化ΔG

分子力场的组成形式 The composition of molecular force field

分子力学的能量表达形式包括:

1. 成键相互作用项

• 键伸缩能bond stretching/compression

• 键角弯折能angle bending

• 二面角扭转能torsion rotation

• 交叉作用项crossing terms

2. 非键相互作用项

• 静电相互作用能electrostatic interactions

• 范德华相互作用能van der Waals interactions

3. 附加项

以AMBER力场函数为例说明:

键伸缩能（Bond Stretching）：Harmonic势函数

最基本的描述方法就是胡克定律（Hook’s law）

参考键长：当所有其它项在力场中值为零时，键长所采用的值

平衡键长：当体系处于能量最小时，键长所采用的值，此时其它项目对体系存在贡献

谐振子势函数在势井底部比较符合真实情况，即基态位置上比较符合实际，远离则不准确

键伸缩能（Bond Stretching）：Morse势函数

AMBER 、CHARMM等采用谐振子函数形式

MM系列力场、MMFF94等采用莫斯函数

键角弯折能（Angle Bending）：Harmonic势函数

• 谐振子模型在偏离平衡位置不大的情况下（10°以内）可以取得很好的结果

• 采用谐振子的力场包括：AMBER 、CHARMM等

二面角扭转能（Torsion Rotation）

注意：

• 由于二面角的扭转对总能量的贡献小于键长和键角的贡献，一般情况下二面角的改变要比键长和键角的变化自由得多

• 因此在一些处理大分子的力场中常保持键长、键角不变，只考虑二面角及其他的作用来优化整个分子的构象和能量。比如构象搜索就指定扭转角个数来寻找能量最低点

交叉相互作用项（Crossing Terms）

• AMBER、CHARMM等力场中没有交叉相互作用项

• MM2和MMFF94支持键伸缩-键角弯折相互作用项

静电相互作用（Electrostatic Interactions）：由于不同元素吸引电子的能力不同，从而产生分子中电荷的不均衡（不规则）分布。描述这种电荷分布的方法有许多种。

电荷分布的表征和分布计算：

• 点电荷法

• 偶极矩法

点电荷法：

偶极矩法：根据一定规则计算出每个化学键的偶极矩，通过计算偶极-偶极相互作用来描述静电相互作用。

范德华相互作用（Van der Waals Interactions）：

静电相互作用不能完全概括一个体系中所有非键相互作用。如在惰性气体原子中，多极矩为零，不存在偶极-偶极或偶极-诱导偶极相互作用，但实际存在原子间相互作用。这种偏差即是由范德华相互作用引起的

范德华相互作用是存在分子间的一种相互作用，主要有三种来源，即诱导力、色散力和取向力

范德华相互作用的计算：原子及分子间的色散吸引及交换排斥相互作用可通过量子力学来计算，这种计算通常考虑电子相关及采用较大基组，计算量巨大。在分子力学中通常使用较简便的方法计算，其中最有名的范德华势函数是Lennard- Jones 6-12函数（兰纳-琼斯势）

氢键作用（Hydrogen Bond）

• H同时与两个原子相互作用

• 氢键能量在15~20 kJ/mol，而一般的共价键能量在400 kJ/mol

• 在生物体系中H一般与O、S和N形成氢键

由不同方法计算得到的能量绝对值无实际意义。只有当它与同体系的其它构象能量相比较时才有意义

• 比较不同程序计算得到的能量值无意义

• 用同一种程序时，比较不同分子的能量值无意义

力场的参数化

分子力场由两部分组：势函数形式，配套的力场参数

力场参数拟合的主要方法：实验拟合，量化计算拟合，此外，还有基于经验规则产生的力场参数（如Dreiding力场）

目标参数来源：

• 各类键长、键角的“本征值”一般取自晶体学、电子衍射或其他谱学数据

• 键伸缩和角变力常数主要由振动光谱数据确定

• 二面角扭转力常数通常由分子内旋转位垒（NMR谱带和弛豫时间、从头计算）来推算

• 非键相互作用参数则主要由晶格参数和液体的物理性质数据拟合得到

分子力场的种类 The type of molecular force field

• 全原子力场（all-atom force field）：考虑分子体系中的每一个原子，并对其分子参数进行精确定义。模拟最为精确，但计算量大

• 联合原子力场（combined force field）：相对全原子力场进行了简化，忽略了体系中的非极性氢原子，将其参数整合到与它们成键的相邻原子上

• 粗粒化力场（coarse grained force field）：进一步精简分子结构的力场参数，种类较多。例如，将氨基酸侧链看作一个颗粒而进行力场拟合。这一类力场通常针对特定体系

MMX力场

Allinger group, 1976

• 包括MM、MM2、MM3和MM4

• 主要适用于非极性有机小分子

• 函数形式比较复杂，包含交叉项

• 也可用于生物大分子体系，但是速度较慢

AMBER力场

Kollman group, 1981 * 最初仅为蛋白质和核酸体系提供相应的原子类型和力场参数 * 1990，发展了适用于多糖模拟的力场参数（Homan 1990） * 1995、2000、2004加入了适用于有机小分子的原子类型和参数 * 使用最广的分子力场之一（ff14SB+gaff） * 被SYBYL等综合商业软件所采纳

Wang, JM; Wolf, RM; Caldwell, JW; Kollman, PA; Case, DA, Development and testing of a general amber force field, J. Comput. Chem. 2004, 25, 1157. (citation>7500)

CHARMM力场

Karplus group, 1983

• 适用于各种分子性质的计算和模拟

• 对于从孤立的小分子到溶剂化的大生物体系都可以给出较好的结果

• 2009年加入CHARMM小分子力场使得应用范围更广泛（charmm36+CGenFF）

• 被Discovery Studio等综合商业软件所采纳

OPLS力场

Jorgensen group, 1985

• 将模拟液态小分子非键相互作用参数加入AMBER力场而产生的新力场

• 特别适用在液相系统中模拟体系物理性质

• 被不同分子模拟软件包广泛采纳，如GROMACS以及Schrödinger中的改良版

GROMOS力场

Gunsteren group, 1984

• 联合原子力场，力场参数主要以纯流体或混合流体体系在凝聚态下的热力学特征为实验数据拟合生成

• 适用于大分子体系的力场

• 计算速度较快

• GROMACS软件原生力场

MMFF力场

Halgren group, 1996 (Merck Molecular Force Field)

基于MM3力场发展而来，定义了非常完备的原子类型 既适用于有机小分子，也适用于大分子体系，如蛋白质 被多个著名分子模拟软件包采纳，如Discovery Studio、MOE、SYBYL、RDKit

注意事项：

• 不同力场定义原子类型不同。因此，在选择力场时，需充分考虑所选力场能否完整表征目标体系中的所有原子

• 不同力场势能函数组成不同。力场对于不同的势能函数分项各有侧重，因此，应该充分考虑目标研究体系的主要特征，选择适应的力场。如有的力场考虑氢键，有氢键函数

• 不同力场势能函数参数不同。力场参数通过不同的实验体系拟合而成，如amber蛋白力场主要从蛋白大分子体系的实验数据中拟合参数，而 MM力场则从有机小分子体系中拟合参数注

常用分子力学优化算法 Molecular mechanics minimization

分子力学为精确结构分析奠定了基础

• 众所周知，分子的物理、化学及生物学性质依赖于其三维结构（尤其是生物大分子中），即分子构象

• 虽然实验技术的发展使得分子结构测定变得普及，但是即使如此，仍然效率低、开销大。因此，分子稳定构象或优势构象的预测是分子模拟的重要任务之一

• 分子力学的发展为精确构象分析奠定了理论基础，其中能量优化算法为稳定分子构象的识别提供了算法保障分

能量优化的方法及特点

单纯形法（simplex algorithm）：无需设置力场，通过调整原子位置排除位阻（如根据原子范德华半径判断位置是否 合适），精度较低，一般用于调整体系初始结构（如模建蛋白侧链调整）

导数法（derivative algorithm）：需要设置力场，可提供势能面信息，导数的方向指向能量最小化的方向。可分为一阶导数法及二阶导数法

一阶导数方法

• 最陡下降法（Steepest Descent）：最陡下降法是最常用的优化方法之一。能量梯度是进行搜索的方向，每次搜索之后旧的方向被新点处的梯度所取代。该方法对能量梯度依赖性较大，因此当体系能量远离最低点时，优化效率很高；但在最低点附近时，由于能量梯度接近于零，收敛很慢。适用于优化的初级阶段

• 共轭梯度法（Conjugate Gradient）：通过当前的梯度和先前的最小化梯度比较确定能量最低点。在能量最低点附近时，收敛速度比最陡下降法快。

对于复杂体系（如蛋白体系），考虑到最陡下降法和共轭梯度法各有优势，因此实际应用中常采用二者联合的优化策略。如 在amber模拟中：首先以最陡下降法进行优化，以期快速降低体系总能量。随后采用共轭梯度法，加快在能量低点附近的收敛速率，获得能量最优构象。

二阶导数方法

• 牛顿-拉森法（Newton-Raphson）：牛顿-拉普森法是一种二级微商算法，需要求解二阶导数矩阵（Hessian矩阵）。一级微商可确定优化方向，而二级微商可确定沿梯度在何处改变方向，因此这种算法的效率优于一阶导数法。由于在计算中需储存大量的二级微商，该算法计算量较大，速度慢，不适用于大体系。可首先以其他算法优化到能量最低点附近，再以此方法进行优化。

各优化方法的特点：

• 最陡下降法：方向变化大，收敛慢，优化幅度大

• 共轭梯度法：收敛快，易陷入局部最小值，对初始结构偏离不大

• 牛顿-拉森法：计算量大，当微商小时收敛快优化策略：可以多方法联合使用

注意事项：

• 结构优化往往只是局部优化，因此获得稳定构象只是初始结构附近的“最优构象”，而往往不是全局最优构象，因此结构优化的结果与初始构象的选择密切相关

• 若想找到全局最优构象，则需要将所有可能的初始构象进行优化，然后通过能量的比较获得全局最优构象（全局构象分析）

• 对于复杂分子体系，由于体系初始构象可能随自由度的增加而指数增加，因此往往很难通过简单优化策略获得全局最优构象

构象搜索与分析 Conformation searching and analysis

构象分析简史

1950年前后，德里克·巴顿（Derek Barton）在取代环已烷反应性的研究中指出其反应性受取代分子的轴向性质影响，即受取代环已烷构象影响

1960年，亨得利克莱森（Hendricleson）用理论方法计算了甲基环己烷（methycyclohexane）扭船式与椅式构象转化过程的能量为6.7 kJ/mol，表明扭船式是可以存在的。直到1975年研究者才从实验中证实了环已烷扭船式的存在，表明基于计算的构象分析可以有效辅助实验观测。

事实上，甲基环己烷还有更多的局部最小构象（local minima）

构象分析方法（conformational search）：

定义：指在构象空间中寻找所有可能出现的比较稳定的构象，即在分子的势能面上寻找所有极小值位点

目的：为了确定分子的优势构象

特点：不同于能量最小化，能量最小化的一个特征是只能优化到最接近初始结构的极小值位点；构象搜索可以搜索全空间（理论上）的极小值位点

能量势垒（Energy Barrier）：当一个构象转变成另一个时，通常会遇到能量势垒，其可能来自于分子内部或分子间。势垒越高，结构转换越慢越难。

常用构象搜索算法 * 系统搜索法（仅用于小分子体系） * 片段连接法（常用于小分子体系） * 随机搜索法（常用于小分子体系） * 遗传算法（常用于小分子体系） * 模拟退火法（大分子、小分子体系均适用） * 分子动力学模拟方法（大分子、小分子体系均适用）

系统搜索法（格点搜索法）：即在构象空间中以小的间隔变量进行逐点搜索，如果变量足够小，则可能搜索到全空间

步骤：

• 固定键长、键角

• 确定分子中可旋转键

• 设定步长

• 变化二面角，产生新构象并计算能量

案例：丙氨酸二肽势能面搜索

• 确定两个旋转角φ、ψ

• 固定键长键角

• 设定步长θi=300

• 计算体系能量（利用AMBER力场）

由图可见，两个区域特别重要,对应于α-螺旋（helix）和β-折叠（strand）结构

优越性：对于自由度较少的分子，系统搜索法得到的构象可以离散地覆盖整个势能面，即可以覆盖整个构象空间 缺点：计算效率低，只能应用于小分子。随着可自由旋转键的增加，构象数目呈指数增加。

\[Number\_of\_Conformation = \prod_{i=1}^N \frac{360}{\theta_i}\]

• θ_i是对键所选的两面角步长，N为旋转键的数

• 如N=5，θ=30 ，构象数为125，即248,832，69 hr（设1s计算一个构象）

• 如N=7，θ=30 ，构象数为127 ，即35,831,808 ，415天

片断连接法

• 系统搜索方法不能用于包含很多二面角的体系，在这个方法基础上发展起来的“片段连接（build-up）”方法可以在一定程度上缓解这个限制

• 基本思路：每种分子片段有其优势构象，通过把几个三维的“结构片断”相互连接组成一个完整分子来实现构象分析的过程

• 方法优势：这种方法比系统搜索更加有效的原因在于，与分子中可旋转单键相比，其组成“结构片断”的数目会少一些。这种方法尤其对于一些环系片断可能更为有效。常用于小分子体系

步骤：

• 确定构造整个分子需要哪些“结构片断”

• 产生所有片断的构象模式（储存于片断数据库中）

• 把“结构片断”相互连接

片断连接方法依赖条件：

• 结构片断的构象要相对独立。它的构象不依赖于与之连接的其他“结构片断”构象的变化

• 储存在数据库中的“结构片断”的构象模式应该覆盖这些“结构片断”在不同分子中所有可能的构象模式

• 内在缺陷：对于复杂分子，分子中基团之间的构象往往相互影响，“结构片断”的典型构型未必能够很好替代其在分子中的构象模式

随机搜索法：

• 改变笛卡尔坐标：识别分子中二面角，在每个原子的坐标x、y、z上随机加上随机量

• 改变内坐标：随机改变分子中可旋转键二面角

• 两种方法效率基本相等

随机搜索法的优势与缺陷：

• 比使用分子动力学的方法效率提高一个数量级

• 随机性较大，在有限搜索时间内，很难判断是否找到分子最佳构象

• 一般用于小分子体系

遗传算法

是一种借鉴生物界自然选择和自然进化的概率选择搜寻算法，具有高度并行、随机、自适应等特点

基本思路：以构象搜寻为例，个体由一系列可旋转键角度组成（染色体），通过大量个体的复制、变异、交叉互换产生新种群，并测量新种群（个体）的适应值（构象能量）来保留优势群体继续迭代搜索，最终得到优势构象

方法优势：并行度高，可同时获得一组优势构象

其中：

编码并产生初始化种群（产生大量分子构象）

计算个体适应值（如通过分子力学计算各构象能量）

复制、互换、变异产生新群体（不同个体的角度互换或随机更新）

遗传算法的特点与不足：

• 随机性较大，如初始种群的大小、突变概率、交叉概率的选取均可能导致十分不同的结果，参数的选择亦处于经验水平

• 对于自由度较多体系，所需初始种群的数量亦需指数增加，因此对于较大体系可能显著降低收敛效率，且计算量较大

• 常用于小分子体系

模拟退火法

也称作蒙特卡罗退火法，统计降温法或随机弛豫法。模拟退火方法是蒙特卡罗方法（Metropolis算法）搜集构象空间，其不同于传统蒙特卡罗方法之处在于运用体系能量的同时还把温度也作为体系的一个变量

整个过程分为两个步骤：

1. 升温熔化体系

2. 逐渐降温

在任一个温度下，体系的初始构象1，相应能量为E(1)，构象发生微小的随机变化产生新的构象2，相应的能量为E(2)，能量变化为ΔE=E(2)–E(1)

当ΔE<0时，接受构象变化； 当ΔE>0时，则在（0，1）间选择一个随机数R（如0.35），将其与P(ΔE)=exp(-ΔE/k_BT)相比较（Metropolis准则）若P(ΔE)>R则接受变化，新的构象成为下一次随机变化的起始点；否则拒绝变化，旧的构象仍是下一次随机变化的起始点。不断循环往复，在每个温度下都进行这个过程，直至体系温度降到足够低，即体系“冻结”在某个固定的构象上。

模拟退火法特点:

• 模拟退火效率如何取决于一些参数的设置，例如初始温度、降温因子以及随机数种子等，这些参数设定了退火的进程，决定了退火的效率

• 模拟退火方法的优点在于它取舍构象时不仅接受能量下降的变化，同时也接受部分能量上升的变化，因而有可能跳出局部势阱，寻找到新的能量最低点。同时模拟退火方法不依靠于体系起始构象，这样就消除了人为因素带来的影响

• 适用于大分子构象搜索

分子动力学方法

基本思想：根据分子的势函数，得到作用在每个原子上的力，利用牛顿第二定律求解运动方程，得到原子在势能面上的运动轨迹，从而达到构象搜索的目的

• 基于分子动力学模拟的构象搜索方法

  1. 淬火模拟（Quenching Simulation）：高温模拟

  2. 退火模拟（Annealing Simulation）：低温模拟

• 淬火模拟可以找到更多分子构象（高温易于翻越势垒），而退火模拟可以搜索出能量更低的构象

• 适用于大分子构象搜索

分子动力学模拟 (Molecular Dynamics Simulation)

分子力学与分子动力学区别

• 分子力学方法：研究体系处于某一稳定状态的结构和性质，与“时间”无关

• 分子动力学模拟：研究原子和分子层面的动态行为，获得体系随时间变化的性质和规律，是真正的“计算机实验”

分子动力学基本原理 The principle of molecular dynamics simulation

分子动力学里程碑

1957年，Alder首次运用MD模拟，对刚性球体系进行研究（无势能函数，J. Chem. Phys. 1957, 27, 1208）

1974年，Rahman完成了第一个真正意义的MD模拟，即液态水模拟（J. Chem. Phys.1974, 60, 1545）

1977年，McCammon完成了第一个蛋白大分子的MD模拟（J. Andrew McCammon,Bruce R. Gelin & Martin Karplus, Nature, 1977, 267, 585–590）

分子动力学方法：利用牛顿力学基本原理，通过积分运动方程产生体系的连续构型，从而得到所有原子的运动轨迹，即体系中各粒子的位置和速度随时间变化的过程（坐标）

可通过对动力学轨迹的进一步分析获得体系的各种性质

牛顿运动定律：

• 物体总是连续地作匀速直线运动，除非一个外力作用于它

• 力等于动量与时间的比率，即 𝑓 = dp/dt = d(mv)/dt = m dv/dt = ma

• 对每一个作用都有一个大小相等方向相返的对抗作用

分子动力学一般步骤

1. 给定条件参数（温度、粒子数、模拟时间等）

2. 体系初始化（初始位置和速度）

3. 计算作用于所有粒子上的力

4. 解牛顿运动方程，计算短时间内（time step）粒子的新位置

5. 计算粒子新的速度和受力情况

6. 重复3-5直至体系达到设置模拟时间。体系平衡后，等间隔保存原子的坐标，这些信息称为轨迹（trajectory）

7. 分析轨迹数据，得到体系的统计性质

体系初始化:

• 初始化粒子的坐标（确立体系的初始构型）：依据实验数据（X-Ray、NMR、冷冻电镜）、同源模建等

• 初始化粒子的速度（对每个原子指定初始速度）：根据温度T时各粒子遵守Maxwell-Boltzmann速度分布随机生成

粒子受力求算：

• 按照经典力学理论，一个粒子受到的力等于体系势能函数对粒子坐标的一阶导数，其中势能由力场决定

• 在每个时间步长中，各原子的受力情况通过势能函数求算

• 两个原子之间的受力大小相等，方向相反

分子动力学中积分方法：

• 在一个连续势的作用下，所有质点运动被偶连在一起，引起一个多体问题，不能用解析法求解。在这种情况下，运动的方程只能利用有限差分方法来积分

• 有限差分方法的基本思想是把整体划分成许多小段，每一段通过固定的时间被分隔。在时刻t，体系中每个质点的总受力情况等于与其他质点相互作用的矢量加和

\[X_i = \frac{F_{xi}}{2m_i} (t - t_0)^2 + v_{i,0} (t - t_0) + X_{i,0}\]

• 经过𝛿𝑡=(𝑡 − 𝑡_0)变化后，力𝐹_xi也发生改变

• 所以新的力𝐹_xi又被重新计算，从而又导出新的位置和新的速度

采用有限差值方法的分子动力学计算可以按照以下步骤进行：

模拟时间步长选择

在一个分子动力学模拟中，步长选择没有硬性规则

• 步长太小，有限时间内覆盖像空间不足

• 步长太大，不稳定性在积分算法中可能升高，由于原子间高能量的重叠，这种不稳定性可能导致模拟体系的崩溃

当模拟可柔分子时，一个有效的控制方式是时间步长应该设置为运动最快周期的十分之一的时间

最高频率振动是由于键伸缩，特别是与氢原子键合的键，一个C-H键振动具有一个近似10 fs的重复周期

因此一般时间步长应设置为1fs

一个真实的模拟通常会固定C-H键的伸缩，因此可以将模拟时间步长设置为2fs

边界条件（显示溶剂模拟）

在实际均一的溶液体系中，每个溶质分子（如蛋白质）周围都可以看作为一个无限的周期性环境

基本假设：

1. 把溶质分子和周围的部分溶剂分子当作一个相对独立的单元（模拟体系）

2. 如果研究体系受周围溶质或溶液的影响，则应采用周期性边界条件

3. 当所关心的问题受周围环境影响较小时，模拟可以采用简单的非周期性边界条件

周期性边界条件：

周期性边界条件中基本单元形状的选择与模拟体系的形状有关：

• 最常用的是立方体和截顶八面体，采用这两种形状时，分子的考察最为简单。如果模拟的是一个核酸分子，则采用六方柱

• 如果模拟的是一个球形的蛋白质分子，则采用截顶八面体或菱形十二面体

非周期性边界条件：

• 当采用非周期性边界条件考虑溶剂效应时，在研究的重要部位外加入非周期的水分子环境，比如球形水环境

• 下图所示的水分子环境是以活性位点的重心为球心的半径为20Å的水分子球

• 在实际的分子力学或分子动力学模拟中，为了加快计算速度，可以把体系划分为柔性区（R1以内）、约束区（R1-R2之间）和完全刚性区（R2以外）

分子动力学轨迹分析

分子动力学轨迹可以给出任何分子性质随时间的演化特征（含时性质），如体系构象转变过程、药物结合/解离过程等，常用于分析体系动力学性质的方法包括均方根位移（RMSD）、均方根振动（RMSF）、分子间相互作用能等

从单个轨迹点分析出的热力学性质为静态性质（特殊案例：分子力学优化能）

从多个轨迹点分析得出的热力学性质称为含时性质

能够预测含时性质是分子动力学最大优势之一

均方根位移 Root Mean Square Deviation（RMSD）：计算体系在动力学模拟中随时间变化情况，可用于分析体系总体构象变化，判断其是否达到稳定状态

均方根振动 Root Mean Square Fluctuation（RMSF）：计算动力学轨迹中每个片段（如氨基酸残基）的平均波动情况，可用于分析体系局部稳定性等情况

分子间相互作用计算

计算动力学轨迹中两片段相互作用情况（如药物-靶标结合自由能），由于考虑了分子间相互作用动态变化过程（含时性质），其计算结果比单点计算（如分子对接打分）更为严谨

分子力学优化及动力学模拟实验

实验目的：

1. 掌握分子动力学模拟，观测体系构象变化。

2. 熟悉 Discovery Studio 的基本操作。

3. 熟悉 Protein Data Bank 数据库。

实验原理：

使用 Discovery Studio 软件进行分子力学优化及动力学模拟。

本实验所用软件环境：

DS Version：19.1.0.18287 PP Version：19.1.0.1963 DS Client Version：19.1.0.18287 OS Distribution：Windows OS Version：10.0.22000

分子动力学操作流程 The process of MD simulation

分子初始构象获取:

小分子构象获取:

• ChemDraw等软件画出结构

• 数据库获得

大分子结构获取:

• PDB数据库获取

• 同源模建

分子预处理:

小分子结构预处理

• 分子状态重生成（异构体、质子化等）

大分子结构预处理

• 蛋白结构检测（去除非必须元素，主、侧链修复等）

• 质子化状态确定

蛋白质结构的预处理：通过实验方法所测得的生物大分子三维结构或多或少存在一些缺陷，比如分辨率低导致模型质量不佳、晶体结构无H原子和电荷信息、在模型构建步骤中相对粗糙的力场带来的结构误差等，这就需要在做SBDD前对生物大分子（如蛋白）的三维结构进行准备和能量最小化。

蛋白结构准备的一般原则：

• 缺失氨基酸及其侧链的修复

• 正确设置键级、计算电荷和加氢

• 限制性优化

• 水的处理

• 辅因子的处理

• 其它小分子的处理

Standard MD Cascade：

• 分子力学优化（Minimization）：体系能量优化，结构优化的目的在于优化分子中因实验（低精度结构）或模建（加H或侧链修复等）产生的结构不合理问题，使用最陡下降法与共轭梯度法联合优化分子结构，以获得稳定初始状态。

• 升温期模拟（Heating）：为了保持体系稳定，模拟初期需需从低温开始模拟（如50 K），逐渐升温至加室温（300 K）或体温（310 K）。体系首先根据Maxwell-Boltzmann分布随机生成初速度。升温过程需要逐步慢速进行，以防止体系不稳定。一般需要对体系进行一定限制。

• 平衡期模拟（Equilibration）：升温过程结束后，即进入平衡期，平衡期模拟的主要作用是避免在加热过程中产生的局部或者全局的不稳定构象（T=300 K）。一般而言，经过平衡后，体系趋于稳定，可以开始后续的采样及分析。一般需要对体系进行一定限制。

• 产生期模拟（Production）：产生模拟又称为采样模拟（sampling），是一个MD模拟有效数据的来源，也是MD模拟整个过程耗时最长阶段。此阶段通常不加任何限制。

数据分析：

动力学轨迹观察

RMSD计算

RMSF计算

实验步骤：

1. 获得分子起始构象：从 PDB 数据库网站（https://www.rcsb.org/）获取 HP53 蛋白晶体结构(PDB code: 1YRF)。(下载)

2. 分子预处理：

   （1）点击 Discovery Studio 软件上的 Macromolecules → Prepare Protein → Clean Protein 进行蛋白结构修复，矫正。

   （2）点击 Discovery Studio 软件上的 Macromolecules → Prepare Protein → Prepare Protein 进行蛋白结构修复，矫正，质子化状态确定。设置参数如下：

3. Standard MD Cascade：点击 Discovery Studio 软件上的 Simulation → Run Simulation → Standard Dynamic Cascade 进行分子力学优化、升温期模拟、平衡期模拟、产生期模拟。设置参数如下：

4. 结果分析 — RMSD、RMSF 计算：点击 Discovery Studio 软件上的 Simulation → Analyze Trajectory → Analyze Trajectory 进行 RMSD、RMSF 计算.设置参数如下：

实验结果：

分子预处理实验结果 , Standard MD Cascade 实验结果, 结果分析实验结果

定量构效关系（Quantitative Structure-Activity Relationships）(QSAR)

构效关系（Structure-Activity Relationship）的弊端：

• 非定量

• 过于经验化，难以标准化

• 需要大量的化合物样本和实验数据

• 难以预测新化合物的活性

概念（concept）

定量构效关系（quantitative structure-activity relationships，QSAR），是研究一组化合物的活性、毒性、药代性质与其结构（structural）、物理化学性质（physicochemical）、拓扑结构（topological）等之间的相关关系，并用数理统计模型加以表征的研究方法。

QSAR的发展简史

1863年，法国斯特拉斯堡大学Cros博士的论文提到哺乳动物体内醇类毒性会随着其水溶性的降低而增加，并可以达到一个最大值。这是关于分子结构和生物活性之间关系的最早文献报道

1868年，英国药理学家Fraser和化学家Crum Brown在研究生物碱的碱性N原子甲基化前后的生物效应时，提出化合物的生理活性依赖于其组分的变化，即生物活性Φ是化合物组成C的函数： \(Φ＝f(C)\) 这就是著名的Crum-Brown方程。

二十世纪30年代Hammett研究了有机物毒性与分子电性（σ）之间的关系，提出芳香环上取代基的电性效应可表示为 \(σ=log(K_x/K_H)\)，其中 \(K_H\) 和 \(K_x\) 分别表示苯甲酸和相应的取代苯甲酸在25摄氏度水溶液中的电离常数，而且认为化合物的毒性或麻醉作用与溶解度相关。

二十世纪50年代Taft提出立体参数的概念，把Hammett的思想扩展到脂肪族化合物，指出 \(E_s=log(K_s/K_H)\)，其中 \(K_s\) 和 \(K_H\) 分别表示酰基上取代的乙酸甲脂的酸水解反应速率常数和乙酸甲脂的酸水解速率常数，即分子体系的空间效应可决定化合物的水解速率常数。

经典的QSAR模型：Hansch模型 1964年，Hansch和Fujita（藤田）指出，分子的生物活性主要是由其疏水效应（\(logP\)）、立体效应（ \(E_s\)）和静电效应（ \(σ\)）决定的，并且假设这三种效应彼此可以独立相加，据此提出Hansch模型： \(log(1/C) = a logP + b E_s + ρ σ + d\)。Hansch模型揭开了经典QSAR研究的篇章，成为QSAR发展历史中的里程碑。 Hansch模型的提出标志着药物定量构效关系研究的开始，也被认为是从盲目药物设计过渡到合理药物设计的重要标志。

QSAR的意义：在受体结构未知的情况下，揭示化合物的结构与活性的依赖关系，建立表征这种关系的数学模型，以预测未知化合物的活性，演绎受体与药物相互作用关系。

QSAR建模的基本要求：2002年在葡萄牙的Setubal召开的一次国际会议上，与会的科学工作者们提出了关于QSAR模型有效性的几条规则，被称为“Setubal Principles”，这些规则在 2004年11月得到了进一步详细的修正，并被正式命名为“OECD Principles”。规定一个QSAR模型要达到调控目的（regulatory purpose），必须满足以下5个条件：

1. a defined endpoint（明确目标）

2. an unambiguous algorithm（明确算法）

3. a defined domain of applicability（明确的使用范围）

4. appropriate measures of goodness-of-fit, robustness and predictivity（稳定）

5. a mechanistic interpretation, if possible（可解释）

实际上这五个条件已经全面概括了QSAR研究工作的基本要求。

QSAR的分类：

根据配体分子属性：

1. 二维定量构效关系（2D-QSAR）

2. 三维定量构效关系（3D-QSAR）

3. 混合型定量构效关系（XD-QSAR）

二维定量构效关系（2D-QSAR）

二维定量构效关系方法（2D-QSAR）是将分子整体的结构性质作为参数（不考虑分子三维结构），对分子生理活性（活性、毒性、药代性质等）进行回归分析，并量化其化学结构与生理活性相关关系

经典的二维定量构效关系方法包括Hansch法、Free-Wilson法等，以 Hansch法最为著名

根据用途不同可分为：QSAR、QSTR（toxicity）、QSPR（pharmacokinetics）等（QSXR）

特点：不考虑分子三维构象、无需优化、无需叠合、计算效率高

2D-QSAR的一般步骤

数据采集 –> 数据处理 –> 模型构建 –> 外部检测

Paper Database –> Descriptor calculation and selection –> MLR/PLS cross-validation –> External testing

QSAR的数据组成：化学结构特征数据和生物活性表征数据

化学结构特征数据（分子描述符）又称分子描述符（molecular descriptor），是描述分子结构属性的变量，其可以是物理化学参数、图形拓扑结构参数等。分子描述符在QSAR模型的构建中通常作为自变量出现，其可以由实验测得或计算获得。化学描述变量由组成性描述符，理化参数描述符、分子片断、连接描述符和其它参数构成。

组成性描述符：分子量，各类原子数目，芳香环数目，氢键供体、受体数，可旋转键数

理化参数描述符：疏水性参数（Lipophilicity Parameters），电性参数（Electronic Parameters），立体参数（Steric Parameters）

疏水性：非极性溶质与水混合时会形成互不相溶的两相，即非极性分子有离开水相进入非极性相的趋势，即所谓的疏水性（Hydrophobicity），所产生的效应称为疏水效应（Hydrophobic effect）

疏水性在药物的吸收、转运、药物-靶标相互作用以及药代动力学中都有十分重要的作用

由于疏水效应与溶质在溶剂中的分配比例有关，因此常用疏水常数或脂水分配系数（LogP）来表示，疏水常数在QSAR研究中是一个重要参数

疏水性参数（Lipophilicity Parameters）：脂水分配系数（ logP ），正辛醇/水分配系数 \(logP = log(\frac{C_o}{C_w})\)。P值越高，疏水性（hydrophobic）越强。P值越低，亲水性（hydrophilic）越强。

Lipinski Rule of Five（类药五原则）又称五倍率规则，其可解释约70%上市药物

• 化合物的分子量不大于500道尔顿

• 化合物中氢键供体（包括羟基、氨基等）数量不大于5个

• 化合物中氢键受体的数量不大于10个

• 化合物的脂水分配系数的对数值（logP）不大于5

• 化合物中可旋转键的数量不大于10个

五倍率成药规则揭示大部分成功药物不会有太大亲脂性（logP<5），其原因在于： * 脂溶性过大的化合物容易与血浆蛋白结合成大分子，不易透过生物膜，从而降低了其分配特性 * 脂溶性过大易造成药物在脂肪组织中的堆积，从而产生蓄积毒性 * 适当的调节LogP，可以有效改变药物半衰期，从而改变给药剂量和给药半衰期

LogP的理论预测方法：在QSAR研究中，尽管分配系数的实验值较计算值更有价值，但对于没有实验值或很难用实验法来测定的化合物，可靠的估算方法变得十分必要。

预测方法包括碎片加合法和基于分子性质的计算方法

碎片加合法：碎片加合法的基础是假定分子的疏水常数具有加合性。把分子划分为基本片段，每种特定的基本片段具有特定的贡献值，整个分子的logP值是其所含的所有片段贡献的总和。片段的贡献值为取代基疏水常数（substituent hydrophobic constant），用 π 表示，可经计算获得，亦可通过查表获得。碎片加合法的优势在于其概念清楚，计算快捷，结果精度较高，适用范围广，因此是目前应用最广泛的方法。

基于分子性质的计算方法：

分子表面积法：耶洛夫斯基（Yalkowsky）最早报道了分子的表面积与logP的关系。

分子极性表面（Polar Surface Area， PSA）：极性原子以及与之相连氢原子的表面积加和

PSA广泛用于吸收或扩散性质的预测

电性参数（Electronic Parameters）

Hammett电性常数（𝝈）：芳香环间位或对位上侧链取代基对分子反应性的影响（表示芳香取代基的诱导和共轭效应之和），用参数𝝈表示，正值表示为吸电子基，负值表示为推电子基。 \(σ=log(K_x/K_H)\) \(K_H\)， \(K_x\) ：苯甲酸和相应的取代苯甲酸在25摄氏度水溶液中的电离常数。值的正负可以反映基团吸电子或给电子能力。值的大小可以反映反应能力的强弱。

立体参数（Steric Parameters）

Taft立体参数（\(E_s\)）：取代基团的大小对酸性介质中脂肪族化合物水解速率的影响，其值均≤0（氢为最小基团）。 \(E_s=log(K_s/K_H)\) 其中: \(K_H\) =乙酸甲脂的酸水解速率常数。\(K_s\) =酰基上取代的乙酸甲脂的酸水解反应速率常数。取代基越大，水解速率越慢，\(E_s\) 越负。氢的 \(E_s\) 值最大，为零。

\(-CCl_3\) \(E_s\) =-2.06

\(-CF_3\) \(E_s\) =-1.16

分子片断描述符：分子片断描述符将分子中某一特征片断，如原子片断、环片断以及亚结构片断作为描述符代码，是一种拓扑学范畴的描述符。

由于分子片断描述符仅考虑彼此独立的分子片断，而可能丢失分子结构内部各基团的排列位置与相互联系的信息，因此产生了分子连接指数描述符

分子连接性指数：分子连接性指数反应了分子中各原子排列状况、分支大小，其与多种理化常数及生物活性相关

分子片段组成和分子片段连接关系进一步拓展构成分子片段与连接描述符，即分子指纹（molecular fingerprint）

其它参数：此外还有位置描述符、环境描述符、几何描述符等，有时为了尽可能地减少信息损失，可同时并用几种描述符。

活性数据又可称为应变量，由实验测定，可以是连续的如“y=pC”，也可以是离散的如“活性-非活性”、“弱-中-强”等。在QSAR中，应变量活性参数通常以产生标准生物效应时药物的物质的量剂量或物质的量浓度的负对数（log1/C）表示。

药物的生物活性定义为产生预定的生物效应时所需剂量或浓度：

• 半数有效量EC/D_50：产生最大生物效应一半时的浓度/剂量（effective）

• 半数致死率LD_50：一半死亡时的剂量（lethal）

• 抑制活性IC_50：活性被抑制50%时抑制剂的浓度（inhibitory）

• 全抑制浓度MIC：完全抑制所需最低浓度（最低抑制浓度，minimum）

对生物活性数据的要求：

• 准确

• 有代表性

• 同源

• 数量尽可能多

QSAR的建模（分析）方法

目前，几乎所有探索化合物结构-活性关系的分析方法都是以统计学为基础的。进行QSAR数据分析，最常采用的建模种类包括：根据用途分为：回归分析和分类分析（模式判别）。根据算法分为：线性方法和非线性方法。此类方法均属于化学计量学（Chemometrics）范畴。

建立模型的方法是影响 QSAR模型质量的关键因素。目前最常用的建模方法：线性方法包括多元线性回归，主成分回归，偏最小二乘法。非线性方法包括人工神经网络，支持向量机，朴素贝叶斯。线性方法通常适用于回归分析，但很少用于分类分析。非线性方法通常适用于回归与分类分析，一般认为非线性方法是人工智能的基础。

线性回归分析是指对一组数据进行最小二乘拟合并建立函数关系的过程。当有几种性质可能对活性有贡献时，可用多元线性回归来处理。事实上，因变量只受一个自变量影响的情况非常少见，通常由几个自变量共同影响一个因变量。

在QSAR建模中，经典的多元线性方法包括Hansch分析法和Free-Wilson分析法，其可用于同源先导化合物活性的优化和预测，分析药物作用机制，推测受体模型结构等。其最大优点是可获得物理意义明确的因果模型。

回归模型评判方式：\(R^2\) 和S

有很多内部检验的方法可以用来评估一个模型的拟合能力、稳定性和内部预测能力，如相关系数（决定系数）、交互检验以及各种残差分析（均方根误差、标准偏差等）。

相关系数(R) \(R = \sqrt{1-\frac{\sum (y_{pred}-y_{exp})^2}{\sum (y_{exp}-y_{mean})^2}}\) R越高，s越小，表明模型的拟合能力越强。

标准偏差(s) \(s = \sqrt{\frac{\sum (y_{pred}-y_{exp})^2}{n-k-1}}\) 其中 \(n ≥ 5k\), n为样本数，k为自变量数目（分子描述符数目）

R值注意问题：

• 虽然R是衡量总回归效果的重要标志，但是R值的大小与回归方程中因变量个数n（样本数量）以及自变量个数k（分子描述符数量）有关。

• 当n相对于k不大时，会获得较大的R值，即容易产生偶然相关（过拟合， overfitting），特别是当n=k+1时，即使k个自变量与因变量Y完全不相关，亦有R=1的结果。

• 因此进行多元线性回归时要注意n与k的比例。一般认为，参与回归分析的化合物数目n与所得到的关系式中参数项数目k（即分子描述符个数）之比应不小于5:1（至少4:1）

分子描述符取舍原则:若交叉相关系数>0.9，说明两参数高度相关，即回归方程中保留一个即可，对两个条件相似参数可删除与目标值相关性小的参数。

模型检验

内部验证（交互验证）（cross validation）：

• Leave-One-Out (LOO)

• Leave-Many-Out (LMO)

• Bootstrapping

最常用交互检验方法：留一法（Leave-One-Out, LOO）

交互检验：依次从N个样本中抽出1个样本，用剩下的N-1个样本来建立构效关系模型，然后用建立的模型预测抽取出来的1个样本的活性，重复这个操作，直到所有样本都被抽去和预测。然后计算内部测试集预测误差的平方和（PRESS）和交互检验相关系数（QLOO）

\[PRESS=\sum{(y_{pred}-y_{exp})^2}\]

\[Q_{LOO} = \sqrt{1-\frac{PRESS}{\sum{(y_{exp}-y_{mean})^2}}}\]

外部验证

外部测试集（test set）预测：从研究的化合物中挑选出足够多的样本组成预测集，预测集中的样本不参加模型的构建，然后通过模型对预测集中的分子的预测结果来检验模型真实的预测能力（Rtest）。

QSAR模型构建的注意事项：

• 模型应用范围：定量构效关系的研究只能应用于作用机制相同的同源化合物，一般认为，结构相近的同源物，其在体内作用机制相同

• 实验数据选择：当一组同源物的生物活性变化幅度若小于一个对数单位（即小于10倍）时，往往难以得到满意的相关结果，这是由相关系数（R）的计算方法决定的

• 模型可靠性判断条件：一般认为，当 \(R^2>0.8\) 或 \(Q^2>0.5\) 时，模型具有较好预测能力；模型外部测试结果（Rtest）应与内部交互验证结果（Q）相当

• 模型样本调整：在进行回归分析计算中，若有偏离较大的化合物需剔除时，剔除的数据点不能过多，以避免过多人为因素干扰。对于被去除的数据点应给予合理的解释

定量构效关系不能代替传统药物设计，其不能发现全新结构先导化合物。这一方面工作可以通过基于结构的药物设计进行（如分子对接）

在预测生物活性方面并不总是成功的，这是因为所得到的定量关系式不能完全解释化合物与受体间作用，即简单的理化参数不足以完全描述化合物生物活性的本质

QSAR的应用:

• 预测同源物的生物活性：由一组同源物所得到的QSAR模型可用来预测同系物生物活性。对于作用机制相同的同源物，定量构效关系的研究常可得到满意的结果。对于作用机制不同的化合物，其QSAR模型表达形式往往不同，因此不能随意外推到其他类型的分子

• 避免合成过多的化合物：在化合物的设计方面，比较常见的系列化合物为同系物，如甲基、乙基、丙基、丁基等。这样的同系物对于探索哪种理化参数对生物活性有显著的影响方面只能获得少量的信息，因此对指导进一步的化合物设计帮助不大

• 更有目的地提高化合物的选择性：对一系列同源物的两种不同活性（如：亚型选择性）分别进行定量构效关系的研究，往往可以得到不同的关系式，根据不同的关系式所提供的差异性信息可有目的地提高化合物的选择性

• 预测化合物的成药性质：化合物在体内的ADMET性质与其结构的关系不像配体-受体结合时关系那样紧密，因此，许多针对药物药效学而开发的方法难以适用。研究表明，药物结构的某些描述符与其ADMET性质存在定量关系，其可起到预测的作用。QSPR研究，有针对性地改造化合物的结构以期改进药动学性质。QSTR研究，改善分子的毒性。

• 协助理解药物的作用机制：化合物结构的基本骨架必须与受体相适应才能发挥生理作用。取代基团的改变可显著影响化合物与受体的结合。研究表明，取代基是通过疏水性、电性和立体效应等因素来影晌化合物与受体相互作用的。定量构效关系的研究，如Hansch法，可定量描述化合物理化参数与生物活性关系，因而起到协助了解药物作用机制的作用。

三维定量构效关系

3D-QSAR是引入分子三维结构信息进行定量构效关系研究的方法，这种方法可间接反映药物分子与生物大分子相互作用中两者间非键相互作用特征，相对于二维定量构效关系有更加明确的物理意义和更丰富的信息量。1980年以来，三维定量构效关系也逐渐成为基于机理的合理药物设计的主要方法之一。

3D-QSAR分析方法：3D-QSAR是QSAR与计算化学和分子图形学相结合的研究方法，是研究药物与受体间的相互作用、推测模拟受体结构、建立药效-结构活性关系，并进行药物分子优化的有力工具。常用的3D-QSAR包括：分子形状分析（Molecular Shape Analysis，MSA）。假想受点点阵（Hypothetical Active Site Lattice，HASL）。距离几何法（Distance Geometry Methods，DGM）。比较分子场分析（Comparative Molecular Field Analysis，CoMFA）。比较分子相似性因子分析（Comparative Molecular Similarity Indices Analysis，CoMSIA）。

比较分子场分析（CoMFA）

由Cramer于1988年提出，其指出：引起生物学效应的药物分子与受体间相互作用大多是可逆的非键相互作用，如范德华相互作用、静电相互作用、氢键相互作用等，并称这些相互作用为分子场。因此，在受体三维结构不明确的情况下，可以通过研究这些药物分子周围的作用势场的分布情况，并以此作为化合物的结构特征变量，建立其与化合物活性之间的关系，定量预测新设计化合物生物活性。CoMFA的提出是QSAR研究领域中的一项重大突破

按照CoMFA原理，如果一组结构类似的化合物以相同方式作用于同一靶点，那它们与受体分子之间的各种作用场应具有一定相似性，而其活性差异取决于化合物周围不同的分子场。

3D-QSAR建模的主要步骤

Processes for developing a CoMFA model:

• Preparation of the data（数据准备）

• Optimization and overlap of the 3D structures（分子优化与叠合）

• Calculation of molecular fields（分子场计算）

• Model training with cross validation（模型训练）

• Validating the prediction capability with the test set（外部检测）

Training set -> Geometry optimization -> Descriptor Calculations -> Model Development

Test set -> Geometry optimization -> Descriptor Calculations -> Model Validation

CoMFA操作基本过程

1. 分子结构处理，确定活性构象（药效构象）：当有受体的晶体结构或受体与小分子的作用位点信息清楚的情况下，可采用分子对接的方法确定分子的活性构象

2. 计算分子中的原子电荷：获取原子净电荷以便能计算分子的静电场

3. 分子叠合（关键步骤）：叠合一般分为骨架叠合和场叠合。对于结构差异较大的化合物，叠合规则的选取对研究结果影响显著

   骨架叠合规则：

   1. 在了解作用机理的前提下，可用已知活性构象为模板，优化、叠合其他分子

   2. 在活性构象未知的情况下，用活性最高分子的低能构象作为模板，优化、叠合其他分子

   3. 亦可依据活性类似法（Active Analogue Approach，AAA）对所有分子进行系统构象搜索，找出其共同构象，从而确立活性构象进行叠合

4. 分子场计算：在叠合的分子周围产生一个包容所有分子的矩形盒子，并划分成规则排 列的格点。用某种基团或小分子作为探针（ \(H^+\)、\(Csp^{3+}\)、\(H_2O\)、\(CH_3\) 等），在网格中以一定的步长移动（0.4~2 Å），计算分子在格点中的各种分子场，用以作为QSAR建模时的自变量

常用的探针有：

• \(Csp^{3+}\) — 计算立体能和静电能

• \(H_2O\) — 计算疏水场和氢键

• \(CH_3\) — 计算van der Waals场

• \(H^+\) — 计算静电场

5. 模型构建：由于分子场格点较多（往往采集到>2000个分子场值），其通常远超化合物数量，因此不能采用多元线性回归法构建模型，而需采用偏最小二乘法（Partial Least Squares，PLS），以克服自变量数目过多所引起的随机相关问题。模型精度评价亦采用交互验证系数Q，当Q^2>0.5时表示模型具有较好精度。

6. 外部检测并利用模型预测新设计化合物活性

7. CoMFA模型的可视化：CoMFA可用等高线图（contour map）的方式将各种场的分布用图形直观表示。从图上可以清楚地观察到各种场分布强弱对化合物活性的影响，从而可以根据模型进行现有化合物的结构修饰和改造，设计新的分子结构。

3D-QSAR的限制

• 在3D-QSAR中，分子的活性构象的确定以及叠合很大程度上会受到人为因素的影响，模型的有效性可能受到质疑

• 相较而言，2D-QSAR无需对分子构象进行处理，因此其所构建的模型更为客观，稳定性更好

3D-QSAR的优势 * 3D-QSAR可以更有效、直观的分析化合物与受体之间的潜在作用模式

胰蛋白酶抑制剂2D-QSAR实验

实验目的：

1. 掌握数据集中训练集与测试集的拆分方法。

2. 掌握分子描述符的选择和建模方法的确定。

3. 掌握模型结果的分析。

4. 掌握未知化合物活性预测。

实验原理：

使用 Discovery Studio 软件进行 2D-QSAR 模型的构建、外部数据集检测、未知活性化合物预测。

本实验所用软件环境：

DS Version：19.1.0.18287

PP Version：19.1.0.1963

DS Client Version：19.1.0.18287

OS Distribution：Windows

OS Version：10.0.22000

QSAR建模一般流程：

• 已知活性数据收集

• 数据集准备（训练集与测试集拆分等）

• 分子描述属性计算（传统分子描述符、分子指纹等）

• 模型构建（多元线性回归、偏最小二乘等）

• 外部数据集检测

• 未知活性化合物预测

实验步骤：

1. 已知活性数据收集：本实验使用指导老师提供的 dataset-qsar.sdf 数据集。 下载

2. 数据集准备（训练集与测试集拆分等）：点击 Discovery Studio 软件上的 SmallMolecules → Create QSAR Model → Generate Training and Test Data 进行训练集与测试集拆分。设置参数如下： 完成后，点击打开报告中result中的Test Set和Training Set。

3. 分子描述属性计算（传统分子描述符、分子指纹等）：Discovery Studio 会在模型的构建中自动计算。在构建模型时，只需在 Calculable Properties 中挑选要计算的描述符。

4. 模型的参数设置与构建：点击 Discovery Studio 软件上的 Small Molecules → Create QSAR Model → Create Multiple Linear Regression Model 进行多元线性回归模型的构建。设置参数如下：

5. 未知活性化合物预测：点击 Discovery Studio 软件上的 Small Molecules → Calculate Molecular Properties → Calculate Molecular Properties进行未知活性化合物预测。设置参数如下：

实验结果

训练集与测试集拆分结果, 模型的构建结果 ， 未知活性化合物预测的结果

基于人工智能的药物设计(Artificial Intelligence-based Drug Design)

问题：多样性抑制剂如何建模？ 解决方法：模式识别是以建立分类模型为目的算法统称，常以非线性方法做为算法基础，通过迭代优化获得稳定模型（机器学习过 程）。当生物活性可定性分为有效和无效时，判别分析可以根据化合物属性（分子描述符）进行分类，建立的判别模型可用来预测新化合物的属性类别。

机器学习的概念

机器学习（Machine Learning，ML）是一门多领域交叉学科，涉及概率论、统计学、算法复杂度理论等多门学科，其专门研究计算机如何模拟人类学习行为，以获取新知识或技能，并重新组织已有知识结构，不断改善自身性能。机器学习是人工智能的核心（Artificial Intelligence，AI ），其应用遍及人工智能的各个领域（ML≈AI）。机器学习主要使用归纳、总结而不是演绎、推理。

人工智能的发展简史

人工智能在20世纪五六十年代被正式提出，至今发展超过60年

1946年，第一台电子数字计算机ENIAC（埃尼阿克）问世

1950年马文·明斯基（Marvin Lee Minsky）建造了世界上第一台神经网络计算机，这也被看做是人工智能的起点。

1956年，在由达特茅斯学院举办的一次会议上（Dartmouth Conference），计算机专家约翰·麦卡锡（John McCarthy）提出了“人工智能”一词，后来被认为是人工智能正式诞生的标志。其他与会专家：马文·明斯基（Marvin Lee Minsky，人工智能之父）、克劳德·香农（Claude Shannon，信息论的创始人）、赫尔伯特·西蒙（Herbert Simon，诺贝尔经济学奖得主）等科学家。

1997年，IBM的DeepBlue（“深蓝”）战胜国际象棋世界冠军。它存有70万份大师对战的棋局数据，可搜寻并预测随后的12步棋

2016年3月，AlphaGo以4:1击败韩国围棋冠军李世石，成为近年来人工智能领域重要的里程碑事件（2017年击败国际冠军柯洁）

2017年，AlphaGo Zero（第四代AlphaGo）在无任何输入数据的情况下，自学围棋3天后便以100:0横扫第二版AlphaGo，学习40天后又战胜了在人类高手看来不可企及的第三版AlphaGo —“大师”（AlphaGo Master）

AlphaGo是由Google旗下DeepMind团队开发的人工智能围棋程序，具有自我学习能力，其能够搜集大量围棋对弈数据和名人棋谱，学习并模仿人类下棋

科研工作中的人工智能

Marwin H. S. Segler, et al., Nature, 2018, 555, 604-610.

Chemputer

Leroy Cronin*, Nature, 2018, 559, 377-381.

Leroy Cronin*, Science, 2020, 370, 101-108.

A mobile robotic chemist

Andrew I. Cooper*, Nature, 2020, 583, 237-241.

人工智能在药物发现中的应用

• 靶标的识别：药物靶标的鉴别通常也是药物设计中最难的一步，其通常以大数据为数据基础，通过整合系统生物学（生物网络分析）、生物信息学/计算生物学（数据挖掘与建模）、结构生物学（蛋白三维结构表征）、药理/毒理学（药物的作用机制及代谢途径等）等多学科的交叉来实现新靶点的发现。

  • 生物标记物的鉴别（biomarker identification）

  • 靶标结构建模（protein modeling）

    • 蛋白建模：AlphaFold2 在第14届蛋白结构预测大赛中（Critical Assessment of Structure Prediction，CASP）甩第二名于无形（Baker Group），一举解决了蛋白的结构预测问题，并对人类全蛋白组结构进行了预测（>98%）。

    

    • 基于深度神经网络的从头预测

      • 成对的氨基酸之间的距离

      • 连接这些氨基酸的化学键之间的角度

      • 进化信息

      • 在整个PDB数据库进行训练

  • 蛋白-蛋白相互作用（protein protein interaction，PPI）

• 苗头化合物的发现：

  • 虚拟筛选（Virtual Screening）

  • 老药新用（Drug Repurposing）

• 苗头化合物的优化（先导化合物的设计）：

  • 高活性先导化合物的发现：基于组合库，非线性回归（SVR）建模优化候选活性剂结构，使其成为先导化合物

• 先导化合物的优化（ADMET性质的预测）：对先导化合物进行代谢、吸收、毒性等方面的预测和进一步优化，使其获得更好的体内活性

基于机器学习的判别模型构建（Constructing Classification Model based on Machine Learning Approaches）

• 监督学习（supervised learning）：监督学习，即在机器学习过程中提供标识（目标已知）。监督学习可以从给定训练集中学习出一个函数，通过迭代优化（自我学习）减少模型误差，达到精确预测目的。常见的监督学习算法包括回归分析（regression）和判别分类（classification）

• 非监督学习（unsupervised learning）：非监督学习又称归纳学习，其一般在机器学习过程中没有给定明确标识（目标未知），需要通过迭代优化来确定相关属性。常见方法包括聚类分析（clustering）。

判别模型构建的一般流程

采集数据 -> 处理模型 -> 构建 -> 外部检测

Database Paper(ChemDraw) -> Descriptor calculation and selection -> Machine learningCross-validation -> Model testing

数据集构建:

正样本活性剂的获取：（BindingDB、ChEMBL、PubChem等）

负样本非活性剂的获取：

一般认为，筛选模型的构建需要扩大活性剂与非活性剂的数量比值，如1:100，其原因在于自然界中天然抑制剂的数量要远小于化合物分子的化学空间。

描述符计算

• 化学描述符的计算与传统2D、3D回归QSAR类似，可分为2D描述符（组成性、理化性质、分子片段、拓扑结构描述符等）和3D描述符（各种分子场等）

• 在QSAR研究中常遇到所选用分子描述符之间存在不同程度的相关性，因而使得提供的信息发生重叠，掩盖了要分析的问题本质。因此在数学上可以通过降维的方式来提取主要信息，常用的方法包括系统选择法、主成分分析方法（Principal Component Analysis，PCA）等

判别模型的构建

交互验证与外部检测

模型检验

• 内部验证（5倍或10倍交叉验证）

• 外部验证（预留的部分数据集或外部数据集）

判别模型的精度指标（定性指标）

• TP（true positive）：判断正确的正样本数

• FP（false positive）：判断错误的正样本数

• TN（true negative）：判断正确的负样本数

• FN（false negative）：判断错误的负样本数

• Positive = TP + FP

• Negative = TN + FN

• SE（sensitivity）：敏感度 \(SE = \frac{TP}{TP + FN}\) (0,1)

• SP（specificity）：特异性 \(SP = \frac{TN}{TN + FP}\) (0,1)

• GA（global accuracy）：全局准确率（仅用于平衡样本判别） \(GA = \frac{TP + TN}{TP + TN + FP + FN}\) (0,1)

• MCC（Matthews correlation coefficient）：马修斯系数（平衡样本和非平衡样本均适用） \(MCC = \frac{TP × TN - FN × FP}{\sqrt{(TP + FN)(TP + FP)(TN + FN)(TN + FP)}}\) (-1,1)

非平衡样本中MCC判定的重要性：

假设1:100的数据集中负样本全部判断正确，正样本全部判断错误，则 GA>99%，因而在非平衡数据集中应采用MCC等判断标准，如图，当正样本多数判断错误时，MCC趋近于0，与事实相符合。

定量预测指标—AUC

AUC（area under curve，AUC）：线下面积（0~1），由受试者工作特征曲线（receiver operating characteristic curve，ROC）得出。适用于两组不同类别数据对同一种打分形式的排序（如分类模型对每个对象的打分、概率等）

不限于机器学习，其适用于各种平衡、非平衡样本的二分类问题

判别模型可用标准

• 对于一个平衡数据集（1:1），当一个模型的GA值在内外测试中均达到0.8以上时表明模型的精度较为可靠，可用于推广研究（GA≥0.8）

• 当一个模型的MCC值在内外测试中均达到0.5以上时表明模型的精度较为可靠，可以用于实际的虚拟筛选或预测研究（MCC≥0.5）

• 若内测与外测MCC或GA值相差悬殊（如ΔMCC>0.3；ΔGA>0.2 ），则表明所构建模型过拟合比较严重，需要重新构建模型

判别模型的特点及注意问题

• 可用于多样性分子判别，不强调分子结构是否相似或作用机制是否相同

• 模型构建时，训练集的样本容量一般较大

• 分类模型的构建不考虑一个活性分子的具体活性值，只区分是否有活性，如活性剂（active）或非活性剂（inactive），一般认为IC50<10 μM的分子为活性剂

• 平衡样本（1:1）的分类模型可以使用全局准确率（global accuracy，GA）或马修斯系数（MCC）来判断模型的精度，非平衡样本（如1:100）只能使用MCC值等来判断模型精度

常用的人工智能算法 Commonly used Machine Learning Algorithms

• 朴素贝叶斯（Naive Bayesian，NB）

• 决策树与随机森林（Decision Tree and Random Forest，DT and RF）

• 支持向量机（Support Vector Machine，SVM）

• 人工神经网络（Artificial Neural Network，ANN）

• 深度学习（Deeping Learning）

  深度神经网络（Deep Neural Network，DNN）

  卷积神经网络（Convolutional Neural Network，CNN）

  循环神经网络（Recurrent Neural Network，RNN）

  对抗神经网络（GAN）、强化学习（RL）、变分自编码器（VAE）等

朴素贝叶斯（Naive Bayesian，NB）

• 朴素贝叶斯法是基于贝叶斯定理与特征条件独立假设的分类方法

• 朴素贝叶斯分类器（Naive Bayes Classifier，NBC）发源于古典数学理论，具有坚实的数学基础，以及稳定的分类效果

NBC工作原理

\[条件概率： P(A|B) = \frac{P(AB)}{P(B)}\]

\[全概率公式： P(A) = \sum_iP(A|B_i)P(B_i)\]

\[贝叶斯公式： P(B_i|A) = \frac{P(A|B_i)P(B_i)}{\sum_jP(A|B_j)P(B_j)} \Rightarrow p(C|F_1, F_2, ..., F_n) = \frac{p(F_1, ..., F_n|C)p(C)}{p(F_1, ..., F_n)}\]

贝叶斯公式在药物设计预测模型中的含义：当某分子含有某特征时（F1~Fn），它成为活性剂或非活性剂的概率，概率大的一方为该分子所属类别。

NBC算法举例

8个分子，5个活性剂，3个非活性剂。已知活性剂带芳香环的概率为0.8，非活性剂带芳香环的概率为0.3，现从8个分子中随机挑选一个，结果其带芳香环。求该分子为活性剂的概率是多少？ 0.8163

NBC特点

优点：NBC模型不需要额外的核函数映射输入数据，所需估计的参数较少，对缺失数据不敏感，算法清晰便于理解

缺点：NBC模型假设属性之间相互独立，这个假设在实际应用中通常较难成立，其会对NBC模型的分类精度产生一定影响

决策树与随机森林（Decision Tree and Random Forest，DT and RF）

• 决策树是一种树型结构，其中每个内部结点表示在一个属性上的测试，每个分支代表一个测试输出，每个叶结点代表一种类别

• 决策树学习是以实例为基础的归纳学习

• 决策树学习采用的是自顶向下的递归方法，其基本思想是以信息熵为度量构造一棵熵值下降最快的树，到叶子节点处的熵值为零，此时每个叶节点中的实例都属于同一类

DT特点

优点：决策树算法的最大优点是它可以自学习，其在学习的过程中不需要使用者了解过多背景知识，只需要对训练实例进行较好的标注

缺点：决策树对训练集通常有较好的分类能力，但对未知的测试数据未必有好的分类能力，即泛化能力弱，容易发生过拟合现象

解决方法：随机森林（RF）

RF工作原理

• 从样本集中随机采样选出n个样本

• 从所有属性中随机选择k个属性，选择最佳分割属性作为节点建立决策树

• 重复以上两步m次，即建立m棵决策树

• 这m棵决策树形成随机森林，通过投票表决结果决定数据属于哪一类

上述案例也可以使用NBC等其它分类器进行操作，习惯上，这些分类器组成的“总分类器”，仍然叫做随机森林

支持向量机（Support Vector Machine，SVM）

发展简史：

• 1963年，Vapnik在解决模式识别问题时提出了支持向量方法

• 1995年，Vapnik正式发布了支持向量机算法

理论基础：统计学习理论（VC维理论，Vapnik-Chervonenkis Dimension）

特点：小样本、非线性、高维模式

SVM工作原理

SVM是从线性可分情况下的最优分类面发展而来

• 图中方形点和圆形点代表两类样本，H为分类线，H1、H2分别为过各类中离分类线最近的样本且平行于分类线的直线，它们之间的距离叫做分类间隔（margin）

• 所谓最优分类线就是要求分类线不但能将两类正确分开（训练错误率为0），而且使分类间隔最大

• 推广到高维空间，最优分类线就变为最优分类面

• 一般情况下，分类样本情况复杂，很难用简单线性算法对目标体系进行区分，因此需要引入复杂数理统计方法（如核函数）对数据进行处理

• 核函数（kernel function）的基本功能为接受两个低维空间的向量，计算其在高维空间中经过某种变换的内积非线性映射到线性

• 核函数种类：线性核函数，多项式核函数 Sigmoid核函数，径向基核函数（Radial Basis Function，RBF）

人工神经网络（Artificial Neural Network，ANN）

• 人工神经网络（ANN）是20世纪80年代以来人工智能领域兴起的研究热点

• 人工神经网络从信息处理角度对人脑神经元网络进行抽象，建立简化模型，根据不同的连接方式设计不同网络架构

ANN工作原理

神经元模型:

• x1~xn为输入向量的各个分量

• w1~wn为神经元各个突触的权值

• b为偏置

• f为传递函数（激励函数），通常为非线性函数

• o为神经元输出

• 数学表示o=f(WX’+b)

常用激励函数：

• Sigmoid

• tanh

• ReLU（Rectified Linear Units）

• Softplus

神经网络是一种运算模型，由大量节点（或称神经元）相互连接构成，每个节点代表一种特定的输出函数，称为激励函数（activation function）

每两个节点间的连接都代表一个通过该连接信号的加权值，称之为权重，其相当于人工神经网络的记忆

网络的输出则依据网络的连接方式（拓扑结构）、权重值、激励函数的不同而变化

多个神经元相互作用形成神经网络

深度学习（Deep Learning，DL）

深度学习的概念由Hinton于2006年提出，其概念源于人工神经网络

深度学习可以认为是含多个隐含层的人工神经网络

深度学习通过组合低层特征形成更抽象的高层特征，以发现数据的分布式规律

深度学习特点

传统模型特点：固定特性+简单分类

深度学习模型特点：自学习给定输入数据特征+学习分类器

与传统人工智能相比，深度学习最大特点在于其能够自学习与任务相适应的特征

卷积神经网络（Convolutional Neural Network，CNN）

卷积神经网络（CNN）是一类包含卷积或相关计算且具有深度结构的前馈神经网络（Feedforward Neural Network），是深度学习的代表算法之一

卷积神经网络的研究始于二十世纪80至90年代，在二十一世纪后，随着数值计算设备的改进，卷积神经网络得到了快速发展，并被大量应用于计算机视觉、自然语言处理等领域

CNN工作原理

输入层

隐含层

（1）卷积层（convolutional layer）：卷积核（convolutional kernel），激励函数（activation function）等

（2）池化层（pooling layer）

（3）全连接层（fully-connected layer）

输出层

卷积层（convolutional layer）：卷积层的功能是对输入数据进行特征提取，其内部包含多个卷积核，组成卷积核的每个元素都对应一个权重值

卷积层内每个神经元都与前一层中位置接近的区域的多个神经元相连，区域的大小取决于卷积核的大小，被称为“感受野（receptive field）”，其含义可类比视觉皮层细胞的感受野

卷积层作用：特征提取与放大

池化层（pooling layer）：又称汇聚层，在卷积层进行特征提取后，输出的特征图会被传递至池化层进行特征选择和信息过滤

池化层作用：降低输出规模、增加可解释性、防止过拟合

全连接层（fully-connected layer）：卷积神经网络中的全连接层等价于传统前馈神经网络中的隐含层

全连接层通常搭建在卷积神经网络隐含层的最后部分，并只向其它全连接层传递信号

特征图在全连接层中会失去三维结构，被展开为向量并通过激励函数传递至下一层

LeNet-5：手写字体识别模型（LeNet-5）诞生于1994年，是最早的卷积神经网络之一

基于机器学习的判别模型构建实验

实验目的：

1. 了解数据集中正负样本处理方法。

2. 掌握模型构建与结果分析。

3. 掌握未知化合物活性预测。

实验原理：

使用 Discovery Studio 软件进行，以朴素贝叶斯为例对 FXR 活性剂与非活性剂进行机器学习判别模型构建。

本实验所用软件环境：

• DS Version：19.1.0.18287

• PP Version：19.1.0.1963

• DS Client Version：19.1.0.18287

• OS Distribution：Windows

• OS Version：10.0.22000

判别模型建模一般流程:

• 已知活性数据收集

• 数据集预处理（正样本/负样本、训练集/测试集准备等）

• 分子描述属性计算（传统分子描述符、分子指纹等）

• 模型构建（机器学习算法）

• 外部数据集检测

• 未知化合物类别预测

实验步骤：

1. 已知活性数据收集：本实验使用指导老师提供的 dataset-AI-2.sdf 数据集。 下载

2. 数据集预处理（正样本/负样本、训练集/测试集准备等）：本实验中，指导老师已经做好了正样本和负样本的分类。训练集/测试集的准备：点击 Discovery Studio 软件上的 Small Molecules→Create QSAR Model→Generate Training and Test Data 进行训练集与测试集拆分。完成后，点击报告中的 test set 和 training set 设置参数如下：

3. 分子描述属性计算（传统分子描述符、分子指纹等）：Discovery Studio 会在模型的构建中自动计算。在构建模型时，只需在 Calculable Properties 中挑选要计算的描述符。

4. 模型的构建与内外部验证：点击 Discovery Studio 软件上的 Small Molecules → Create QSAR Model → Create Bayesian Model 进行朴素贝叶斯模型的构建。设置参数如下：

5. 未知活性化合物预测：未知活性化合物数据集用的是已知活性数据收集，点击 Discovery Studio 软件上的 Small Molecules → Calculate Molecular Properties → CalculateMolecular Properties 进行未知活性化合物预测。设置参数如下：

实验结果：

数据集预处理结果, 模型的构建结果 ， 未知活性化合物预测的结果

讨论：

可以从在内外部验证中所得出的模型精度指标看出，内部验证MCC为0.991，外部验证MCC为0.931。因为内外部验证ΔMCC<0.3，所以此模型可以用于实际的虚拟筛选或预测研究（MCC≥0.5）

基于结构的药物设计 (Structure‐based drug design) (SBDD)

介绍

早期药物发现的途径

1803年德国药师萨顿从阿片中分离到黄色结晶(吗啡)

1971年提取、分离、鉴定紫杉醇 “卵巢癌和转移性乳腺癌治愈率达33%，总有效率达75%以上”

一直到20世纪70年代，随机筛选/偶然发现仍然是最主要的先导化合物发现方式，至今仍不可替代

计算机辅助药物设计方法：

基于结构的药物设计的基本概念

基于结构的药物设计（Structure‐based drug design, SBDD），广义上讲是基于配体结构和受体蛋白结构的药物设计的统称；狭义上讲就是基于受体结构的药物设计，即依据与药物作用的靶点（广义上的受体，如酶、受体、离子通道、抗原、核酸、多糖等）的三维结构，运用分子识别原理（互补性），设计对受体进行调控的先导物，或根据已有药物作用力大小和构效关系判断来推测新化合物的药效，达到发现活性分子的目的。

基于配体的药物设计也称为间接药物设计；基于受体结构的药物设计也称为直接药物设计。

蛋白质的结构与预测

蛋白质的三维结构

蛋白质是生物体中含量最丰富的生物大分子，参与了细胞中大部分的生命过程，是细胞最重要的组成物质。

当一个或多个蛋白的结构、功能或其参与相互作用的信号通路发生异常，将引起人类疾病的发生。经过严格的靶标验证实验，那么它就可能成为药物设计的靶标蛋白，通过设计药物分子与靶蛋白的结合便可达到治疗疾病的目的。

对靶蛋白三维结构的获取及分析则是SBDD的先决条件。

蛋白质的基本结构单元：氨基酸

组成蛋白质的种类繁多，结构各异，但元素组成基本相似，主要有碳（C）、氢（H）、氧（O）、氮（N）和硫（S）。除此以外，有些蛋白质发挥功能还需要少量金属离子的参与，如铁（Fe2+）、镁（Mg2+）、钙（Ca2+）和锌（Zn2+）等。

氨基酸的分类

疏水性氨基酸：

这些疏水侧链参与了范德华（van der Waals, vdW）相互作用，它在维持蛋白质构象中起着主要的作用，是蛋白质折叠的主要驱动力。

苯丙氨酸Phe的苯基侧链有时能参与较弱的极性相互作用。

Pro比其它任何氨基酸残基都具有更强的立体化学效应，它的侧链是和主链氮共价相连的，封闭的环常常形成一个转角，改变了主链的方向。

带电氨基酸:

带电荷的氨基酸都属于亲水性氨基酸，包括酸性氨基酸，如天冬氨酸和谷氨酸和碱性氨基酸，如精氨酸、赖氨酸和组氨酸。这些氨基酸主要参与了氢键相互作用和正负电荷的吸引作用（即盐桥）。

组氨酸在电中性状态下，其氮原子上的氢可在两个氮原子上互变（τ互变和π互变），pKa均在6左右，具体以哪个异构体为主取决于其周围的环境；当两个氮均被质子化时，氨基酸整体带正电荷。

在pdb中，组氨酸的残基名都用“HIS”表示，而在Amber力场中，τ互变、π互变及带正电的His则分别表示为“HIE、HID和HIP”。

极性中性氨基酸：

Ser、Thr、Asn和Gln主要参与氢键相互作用，属于亲水性氨基酸。

芳香环类侧链如Tyr和Trp虽然是疏水的，但它们都能形成氢键，可以位于蛋白质的内部或外部。

具有完全疏水性侧链的Phe是蛋白质内部疏水核心的主要成分。

成对的Cys残基能形成二硫键，从而增强蛋白质结构的稳定

Cys也与His一样，常见于酶活性位点

当侧链的巯基去质子化时具有较强的亲核性，常作为共价药物设计的靶向残基。

蛋白中辅因子的作用：辅助蛋白质折叠，物质运输，辅助酶催化

氨基酸通过肽键组成多肽

• 一个氨基酸的氨基与另一个氨基酸的羧基通过脱水形成的酰胺键称为肽键。

• 由两个氨基酸形成的肽为二肽，以此类推。不超过10个氨基酸形成的肽也被统称为寡肽，超过10个氨基酸的肽也称为多肽。

氨基酸通过肽键组成多肽

• 位于肽键平面外的另外两个单键（N-Cα和Cα-C）在没有侧链位阻的情况下是可以自由旋转的。

• 绕肽键一边的N-C_α键旋转的角称为phi（φ），绕着同一C_α原子的C_α-C键旋转的角被称为psi（ψ）。

• 其主链具有可旋转键和刚性平面相互交替的特征，自由度仅与两个构象角密切相关。

拉氏图（Ramachandran图）

• 将蛋白或多肽主链的φ和ψ值可视化作成的图形称为拉氏图，由 G. N. Ramachandran 等人提出的，以φ和ψ的角度为横纵坐标，规定φ和ψ角允许的构象区域的一个图形

• 主要的应用包括了辅助结构生物学中蛋白晶体结构的解析以及对同源建模后模型质量的检测。

在线服务：PROCHECK（https://www.ebi.ac.uk/thornton-srv/software/PROCHECK/）

蛋白质的高级结构

在生理温度下，水溶液中蛋白质的多肽链在大多数情况下折叠成呈球状。蛋白质中不同氨基酸的序列是它的一级结构（primary structure），直接由编码它的基因序列决定，序列又决定了蛋白折叠成更高级别结构的方式，包括二级结构、三级结构和四级结构。

蛋白质的高级结构：二级结构

α螺旋（α-helix）

上：β转角（β-turn）

下：β折叠（β-sheet）

卷曲或环形区域（coil or loop）

蛋白质结构的获取方法

• X‐射线晶体学（X‐ray crystallography）:通过对衍射的位置、强度计算，读出原子坐标值（电子密度），解析结构得到晶体空间结构。在SBDD中，同一靶点的晶体结构要尽可能选择分辨率（Resolution）高（数值低）的结构作为靶标结构进行后续设计。

• 冷冻电镜技术（Cryo‐EM）：冷冻电镜，是用于扫描电镜的超低温冷冻制样及传输技术(Cryo‐SEM)，可实现直接观察液体、半液体及对电子束敏感的样品，如生物、高分子材料等。

• 核磁共振技术（NMR）：核磁共振能够给出单个功能基团的信息,如离子化状态、pKa和氢键；提供溶液中主链和侧链的动态运动的位点特异信息。NMR的优势在于能在接近自然生理环境的条件下对生物分子进行研究。

• 蛋白质三维结构预测：

基于理论的预测方法：从头预测法（ab initio prediction）

• 原理：采用理论计算（分子力学、分子动力学、量子化学）方法，直接从分子和原子参数模拟肽段在三维空间中所有可能的姿态，并计算出自由能最低的一个。

• 特点：计算量极大，不常用。

• 代表工具：QUARK https://zhanglab.ccmb.med.umich.edu/QUARK/

基于知识的预测方法：穿针引线法，折叠识别（threading, fold recognition）

• 原理：通过比较目标序列与已知的折叠库，找到目标序列的折叠模式，并使用相关的方法将各个部分的折叠连接起来。

• 特点：折叠法是近年来发展起来的一种比较新的方法。它可以应用到没有同源结构的情况中，且不需要预测二级结构，即直接预测三级结构，从而可以绕过现阶段二级结构预测准确性不超过 65% 的限度，是一种有潜力的预测方法。

• 代表工具：I-TASSER (Iterative Threading ASSEmbly Refinement)https://zhanglab.ccmb.med.umich.edu/I-TASSER/

基于知识的预测方法：同源模建法（homology modeling）

• 原理：根据待建模序列（目标）与一个或多个已知蛋白结构的序列（模板）间的同源性（序列一致性），直接由目标序列的一级结构预测其三级结构

• 理论基础：蛋白质三级结构的保守性远超过一级序列的保守性

• 特点：模板蛋白和目标蛋白的序列一致性需要大于30%，越大建模准确性越有保障。

• 代表工具：(1) MODELLER (https://salilab.org/modeller/) (2) SWISS-MODEL (https://swissmodel.expasy.org/)

同源模建的步骤（以DS MODELER为例）

1. 搜索并识别模板：

   • 查找目标蛋白序列（https://www.uniprot.org/）

   • 搜索并识别模板

     • 载入序列

     • 使用BLAST搜寻模板：每行表示一条命中的氨基酸序列。命中的序列按照E值（序列无缝比对存在偶然性的可能性大小，表征了序列比对的可行度）进行排序。E值最低的序列，结果最可靠，排在第一行。

     • 挑选一个或多个合适的同源模板（templates）。一个理想的template 需要涵盖整个target 的长度，具有较高的序列一致性（Sequence identity），并且E 值要够小（小于 1×10^-5）。一般而言，若有多条模板与target 都具有相似的同源性，但模板之间相似性并不非常高，那么通常使用多模板来构建同源模型，使建模过程中模型的每个部分都采用最合适的模板。

2. 序列比对：通过多序列比对将目标序列直接比对至模板序列

3. 产生全原子模型（自动产生）

   • 模型的PDF Total Energy越低，表明该模型在同源约束条件下优化的越好；模型同限定的同源约束条件偏差越小，该模型的可信度越大

   • DOPE是一个基于原子统计势能的程序，主要用于模型评估。它的分数可以认为是衡量同一分子不同构象可信度的标准，能够帮助选择预测结构的最优模型。分数越低，模型质量越可靠。一般可以粗略地选取PDF Total Energy最低的模型作为最合理的初始模型；若其值相似，则可以利用DOPE score作为衡量模型质量的依据。

   • loop区建模：主要是目标蛋白和模板蛋白的比对结果中存在缺口的部分如何处理的问题。第一种解决的方式是略去模板蛋白存在的残基，留下一个必须补上的缺口。另一种情况是将主链截断，插入缺少的残基。

   • 侧链建模：构造各种可能的构象体，并利用基于能量的函数打分来实现侧链构象的选择的

4. 模型优化：模型产生初始坐标后，还要对分子结构进行进一步的优化，以此消除/缓解原子间的重叠或某些不合理构象（尤其是非保守区的构象）。优化一般采用分子力学或分子动力学。

   注意事项：

   • 对柔性区域的构象进行构象分析可以采用分子动力学方法

   • 在进行分子力学优化时，应尽量避免破坏模型的主链结构，一般要采取限制性优化方法

   • 在修正过程中，应考虑溶剂效应的影响

5. 模型评估： * 拉氏图：蛋白主链二面角的拉氏图（当落于核心区+允许区的氨基酸残基百分比>95%时，表明模型质量较高） * Profile-3D：Profile-3D 是UCLA 的David Eisenberg 教授开发的一种基于“穿线”（threading）法的模型评估程序。该方法采用3D-1D 的打分函数来检测所构建模型与自身氨基酸序列的匹配度关系。分数越高，说明同源模型的可信度越大。Verify Score 越接近或高于Verify Expected High Score，模型的质量越好。

实验：同源模建法构建蛋白质的三维结构

实验目的：

1. 掌握通过同源模建法进行未知跨膜蛋白 P2RY6 的三维结构的构建。

实验原理：

认为蛋白质三级结构的保守性远超过一级序列的保守性。所以根据待建模序列（目标）与一个或多个已知蛋白结构的序列（模板）间的同源性（序列一致性），直接由目标序列的一级结构预测其三级结构。模板蛋白和目标蛋白的序列一致性需要大于 30%，越大建模准确性越有保障。本实验使用 Discovery Studio 软件通过同源模建法进行未知跨膜蛋白P2RY6 的三维结构的构建。

本实验所用的软件环境：

• DS Version：19.1.0.18287

• PP Version：19.1.0.1963

• DS Client Version：19.1.0.18287

• OS Distribution：Windows

• OS Version：10.0.19044

实验步骤：

1. 查找目标蛋白序列：进入 Uniprot 网站（https://www.uniprot.org），搜索目标蛋白，本实验中搜索 P2RY6。点击 Q15077，点击 Download，再点击 FASTA(canonical)，下载 目标蛋白序列。

2. Blast 蛋白序列比对：

目标蛋白序列的填入：使用 Discovery Studio 软打开刚下载的名为Q15077.fasta的文件。结果如下图：

蛋白序列重命名：右击蛋白序列名称（如下图所示），点击 Rename Sequence，在弹出的输入框中输入如 P2RY6 等简单名称。（本例输入了 P2RY6，如果名称不简单不是英文，后面可能会出现错误，所以需要重命名）

Blast 蛋白序列比对：点击 Discovery Studio 软件上的 Macromolecules ——> Create Homology Models ——> BLAST Search (NCBI Server)来进行蛋白晶体数据库中相似蛋白的查找。设置参数如下

3. 选择合适蛋白模板：本实验选择了 4XNV_A（E-value 最小，E-value（E值）小于 1×10^-5，模板蛋白和目标蛋白的序列一致性大于30%），右键此选项，点击Load Selected Structures 进行蛋白模板的载入。

4. 创建序列-结构对齐矩阵：点击 Discovery Studio 软件上的 Macromolecules ——> Create Homology Models ——> Align Sequence to Templates 生成模板蛋白与目标蛋白的结构对齐矩阵。设置参数如下：

5. 蛋白建模：点击 Discovery Studio 软件上的 Macromolecules ——> Create HomologyModels ——> Build Homology Models 进行同源模型的建立。设置参数如下：

6. 模型评估：

拉氏图的绘制：点击 Discovery Studio 软件菜单栏上的 Chart ——> Ramachandran Plot 分别对每个同源模型进行拉氏图的绘制。

Profiles-3D：点击 Discovery Studio 软件上的 Macromolecules ——> CreateHomology Models ——> Verify Protein (Profiles-3D) 分别对每个同源模型进行 Profiles-3D 的绘制。设置参数如下：

以残基打分作图：选上全部新生成的模型，选择 AminoAcid 选项卡，点击 Discovery Studio 软件菜单栏上的 Chart ——> Line Plot，X 轴设为 name，Y 轴设为 verify score，Color by Data Series 设为Molecule 同时给每个同源模型进行残基打分作图，如下图所示。

实验结果：

Blast 蛋白序列比对的结果 , 创建序列-结构对齐矩阵的结果 , 蛋白建模的结果,

拉氏图的绘制的结果：

Profiles-3D的结果

以残基打分作图的结果：

基于人工智能的预测方法：Alphafold

2020年12月1日，国际蛋白质结构预测竞赛CASP公布的数据显示，谷歌旗下人工智能技术公司DeepMind开发的深度学习算法Alphafold的表现超过了大约100个其他团队，准确性达到了与实验室方法不分伯仲的水平，一举解决了困扰学界长达五十年之久的蛋白质折叠问题。

2021年7月，DeepMind发布了最新版本AlphaFold2的源代码，以及详细描述了它是如何开发的。https://alphafold.ebi.ac.uk/

基于结构的药物设计（SBDD）

基本概念与理论

靶点（受体）：能识别和结合生物活性物质（配体），并产生生物效应的结构。

生物活性物质（配体）：

• 内源性活性调解物：维持机体机能的基本生理机制，如多巴胺、激素等。

• 外源性药物：干预机体生理生化作用。

靶点学说（受体理论）

一种有效的药物必须符合以下三个要求： * 作用于靶点（与机体内的某一种或多种靶点发生相互作用）（药效学要求） * 暴露于靶点（药物到达靶点、达到适宜的浓度 C_max 并维持足够的时间 AUC ）（药代动力学要求） * 不与其他无关靶点作用（ 无 off target 效应） （安全性要求）

配体‐受体相互作用方式

“锁‐钥”原理

空间形状契合、相互作用力契合

诱导契合原理（“手‐手套” 原理）

配体构象变化、受体构象变化

分子识别：药物分子与靶蛋白发生相互作用的驱动力是分子识别。

定义： 药物在特定条件下通过分子间作用力（可逆 ）与靶点相结合的过程。其中，某些药物也可形成共价键 （不可逆 ）。

分子识别的两个决定性作用：互补性(形状和电性互补)和预组织作用（去溶剂化和构象扰动）。

范德华力（van der Waals Force, VDW）

• 相邻的电中性原子互相靠近，会有短暂的微弱吸引力，即范德华力

• 主要的分子间作用力

• 原子间距离为0.3～0.5 nm，范德华力的能量为1.9 kJ/mol

• 主要分为：色散力、诱导力、取向力

  色散力（Dispersion force）：所有分子或原子之间都存在，由于电子的运动产生的瞬时偶极键的作用力。

  • 分子的变形性越大（一般分子量越大，变形性越大），色散力越大

  • 分子的电离势越低（分子内所含的电子数越多），色散力越大。

  诱导力（Induction force）：在极性分子‐非极性分子之间，极性分子‐极性分子之间，都存在诱导力。

  • 极性分子与非极性分子接近时，极性分子的永久偶极产生的电场使非极性分子极化产生诱导偶极。永久偶极与诱导偶极间的吸引力称为诱导力。

  • 诱导力同样存在于极性分子之间，对极性分子而言，诱导力是附加的取向力。

  取向力（Orientation force）：取向力发生在极性分子‐极性分子之间。

  • 当极性分子互相接近时，分子的永久偶极之间同极相斥、异极相吸，使分子在空间按一定取向排列吸引，而处于较稳定的状态。这种永久偶极间的吸引力称为取向力。

  • 取向力与分子的偶极矩平方成正比，即分子的极性越大，取向力越大。

疏水键（Hydrophobic bond）

疏水键是两个不溶于水的分子间的相互作用。这种因能量效应和熵效应等热力学作用使疏水基团在水中的相互结合作用称为疏水键。

• 疏水作用与分子中疏水基团和片段的数目成正比，烷基越多，疏水性越强。

疏水键的本质

氢键 (Hydrogen bond)

氢键：氢原子与电负性大的杂原子共价结合，与另一具有未共用电子对的杂原子以氢为媒介，形成以X‐H…Y形式的一种弱的静电引力。

• 氢键具有饱和性和方向性，在药物和受体形成氢键时，双方的原子和基团在空间的取向和距离决定了结合强度

• 氢键在化学和生物学的多个领域中扮演着重要的角色，是驱动生物分子间可逆相互作用的主要力量

氢键的方向性

氢键供体能量强度

使用HCN作为受体探针计算34个氢键供体形成氢键的能量大小，排序结果如下：

Ref: J. Chem. Theory Comput. 2020, 16, 2846−2856

氢键受体能量强度

使用HCN作为受体探针计算67个氢键供体形成氢键的能量大小，排序结果如下：

静电作用 (Electrostatic interaction)

离子键（Ionic bond）

离子键：两个具有相反电荷的离子间的静电引力，也叫离子离子相互作用。

• 离子键无方向性，是随机转运的药物与受体的初始感受，对于两者间的趋近与识别有重要作用。

• 作用能较强，但在水中，去水合作用与静电引力对复合物形成的能量贡献相反，因而形成抵消效果。

离子‐偶极作用（Ion‐dipole interaction）：带电荷的原子或基团与含有偶极的基团之间的静电引力。对稳定药物受体复合物起到重要作用，但弱于离子键。

偶极‐偶极作用（Dipole‐dipole interaction）：两个偶极分子或偶极键之间的静电引力，弱于离子‐偶极作用。

π‐π相互作用（π‐π interaction）

π‐π相互作用：芳环含有环形π键，芳环之间的π键可发生π‐π相互作用，又称为π‐π堆积作用（π‐π stacking）。

• 作用力较弱，但在生物体内普遍存在，作用能量大约8~12 kJ/mol。

• 芳环堆积的三种类型

从静电吸引角度分析π体系的各种相互作用

阳离子‐π相互作用（Cation‐π interaction）

阳离子‐π相互作用：由有机或无机阳离子与π电子体系的负电势之间的静电引力，能力强度与氢键相近。

• 在水介质中，阳离子‐π相互作用强于离子键。

共价键（Covalent bond）：键能高，约140～800 KJ/mol，结合牢，驻留时间较长，治疗指数高。除了被特异性酶酶解而发生共价键断裂外，不易被破坏。因此产生的药物作用持久，为不可逆过程。

共价药物：共价键键合类型常见于作用于病原体或肿瘤细胞的靶点。药物的共价基团往往具有较高的化学活性而缺乏特异选择性。

配位键（Coordinate bond）

配位键：又称配位共价键，一种特殊的共价键，其共用的电子对是由其中一原子独自供应，另一原子提供空轨道。

• 过渡金属的核电荷高，半径小，有空的d轨道和自由的d电子，它们容易接受配位体的电子对，又容易将d电子反馈给配位体。因此，它们都能形成稳定的配合物。

水分子在药物设计中的作用

药物设计时需要考虑的溶剂化效应

• 药物和受体在结合前都处于水的溶剂化状态，为了双方的结合需去溶剂化，这是耗能过程，能量由结合能得到补偿。

• 如果去溶剂化能大于结合能，药物难以与受体稳定结合，所以分子设计时要确定形成氢键的因素是否是多余。

水分子既是氢键给体也是氢键受体，可参与药物与受体的结合

基于计算：WaterMap

• 由Schrödinger公司开发的一种较为流行的基于分子动力学模拟（MD）的水分子位点预测工具，可预测蛋白质结合位点中水分子的位置和能量。

• WaterMap算法步骤：1、进行约2ns的约束性MD模拟（溶质约束），根据密度分布对水分子位置和方向进行聚类；2、使用非均相溶剂理论计算每个水合位点的热力学性质（相互作用能、熵增和自由能），以此来决定靶向哪个特定的水分子。

基于计算：SZMAP (solvent‐zap‐map)

• SZMAP是OpenEye应用软件包中的一个模块，是一款水分子位置与稳定性预测软件。无需分子动力学模拟采样。

• SZMAP采用显式水分子探针在高介电连续溶剂(high‐dielectric continuumsolvent)中快速计算分子表面附近溶剂化自由能的大小与分布。

• 它可以预测特定水分子对配体结合的稳定化或不稳定化效应，还可以识别结合水的位置和哪些位置的水是无序的。在结合位点，更好地理解这些效应有助于先导化合物的优化和药物设计其他方面问题的研究。

基于计算：3D‐RISM (three‐dimensional reference interaction site model)

• 该程序集成在MOE、Flare和Amber中。

• 计算水H和O密度的时间平均分布，以及用于分析溶剂稳定性和溶剂化对结合自由能的贡献的自由能图。 3D‐RISM的独特之处在于它能够包含各种浓度的盐或疏水分子作为溶剂的一部分，从而生成离子，疏水基团以及水的密度图。

• 基于力场，无需分子动力学模拟采样。

热力学理论在药物设计中的应用

多数药物与受体作用是以非共价键形式结合，形成的复合物与游离的配体和受体之间呈动态平衡。这种结合作用是由热力学自由能(ΔG)变化所驱动，结合越强，配体‐蛋白复合物的离解常数(K_d越小，ΔG负值越大。ΔG与K_d的关系由吉布斯‐范特霍夫（Gibbs‐Van’t Hoff）方程表征:

ΔG = RTlnK_d R = 8.13 J∙mol^‐1∙K^‐1，T为绝对温度。当T = 300 K时，ΔG = 5.7logK_d (kJ/mol)

根据吉布斯‐亥姆霍兹（Gibbs‐Helmholtz）方程，ΔG是由焓（ΔH）和熵（‐TΔS ）构成，即：

ΔG = ΔH – TΔS

配体与受体结合产生自发过程的为负值，所以从上式看，ΔH为负值和‐TΔS为负值，有利于结合过程。ΔG中ΔH和TΔS有各自的贡献，参与并维持着游离的配体和受体与复合物之间的平衡作用。

配体与受体结合的焓‐熵变化，常常是有利的 ΔH 被不利的 ‐TΔS 所抵消，或者有利的 ‐TΔS 以不利的 ΔH 作负性补偿，这就是配体‐受体结合作用的焓‐熵补偿(enthalpy‐entropy compensation)。

如果│ΔH│＞ │TΔS│，则该结合为焓驱动结合，反之，则为熵驱动。

热力学量的实验测定——等温滴定量热法(isothermal titration calorimetry，ITC)

• 可直接测量生物分子结合过程中释放或吸收的热量。

• 准确地确定结合常数 (K_D)、反应化学量 (n)、焓 (∆H) 和熵 (‐TΔS)。

热力学量的计算方法——结合自由能(binding free energy)的计算

在分子模拟领域，准确计算结合自由能仍然是一个挑战。为此，人们发展了许多方法。如：

• 自由能微扰(FEP)与热力学积分(TI)：计算量很大, 不适用于生物大分子体系

• MM‐PBSA (Molecular Mechanics‐Poisson Bolzmann Surface Area, 分子力学泊松玻尔兹曼表面积)：准确度不如FEP，但计算量小, 在分子识别, 区分结合的强弱方面是一种有效的方法，可得到ΔG、ΔH 和‐TΔS的具体值。

分子对接技术 (Molecular docking)

分子对接的原理

定义：利用算法将配体小分子（如药物）放置到受体大分子（如靶标蛋白）的结合位点（Binding site），预测小分子与受体结合构象（Pose）及作用能的过程。

• 分子对接是一种广泛应用的研究分子间相互作用的分子模拟方法，其过程涉及小分子与受体之间的空间匹配和能量（性质，电性等）匹配。

• 利用分子对接方法研究蛋白‐配体作用模式在SBDD中的应用非常广泛。

分子对接的最初思想起源于1894年Fisher提出的受体学说。Fisher认为，药物与体内的蛋白质大分子即受体会发生类似药物与锁的识别关系，称为“锁钥模型”，这种识别关系主要依赖两者的空间匹配。

真正算法上的实现是在1982年，美国加州大学旧金山校区（UCSF）的Irwin Kuntz教授率先提出了“分子对接”的概念，并发表了第一个分子对接软件“DOCK”，至今仍被广泛应用。

分子对接中的两大关键问题：空间匹配和能量匹配

• 空间匹配是分子间发生相互作用的前提。

• 能量匹配是分子间保持稳定结合的基础。配体分子进入结合位点时，通过一定的程序计算它们之间的结合模式和结合能，并对结果进行评价，通过打分函数评判配体-受体的结合程度

分子对接的分类

根据对接过程中是否考虑研究体系的构象变化，将其分为

• 刚性对接：将配体和受体构象都当做刚性处理，以降低计算量

• 半柔性对接：受体保持刚性，配体则完全当做柔性分子进行处理，当前的大多数分子对接软件都采取这种模式

• 柔性对接：对接过程既考虑小分子柔性，又可考虑受体部分区域（如结合口袋中的残基）的柔性

常用的分子对接算法

刚性对接、半柔性对接：LibDock、CDOCKER、Glide、GOLD、AutoDock、DOCK、MOE_dock、FlexX、Surflex‐dock 柔性对接：Flexible docking、Induced fit docking（IFD） 片段对接：MCSS 共价对接：Covdock、GOLD、MOE 蛋白‐蛋白对接：ZDOCK、ICM‐Docking

分子对接的核心步骤

搜索算法：搜索配体与受体的可能结合模式。在分子对接中，配体小分子在受体结合口袋中会有不同的取向和构象，其中配体构象搜索方法起重要作用。对接问题的核心，就是在指定的构象空间内，采取合适的算法搜索合适的配体构象的过程。

• 系统性搜索（Systematic search）：系统搜索，就是在构象搜索时尽可能的去系统地遍历配体全部的自由度或结构参数的组合。但是，随着分子可旋转键数目的增多，配体构象空间会急剧膨胀，要做到完整采样并不现实。分子对接中常用的系统搜索算法有穷举搜索（Exhaustive Search）、基于片段的搜索（Fragment‐based Search）和构象系综搜索（Conformational Ensemble Search）等。

• 启发式搜索（Heuristic search）：启发式搜索（随机搜寻法），需对配体三维坐标进行随机改变产生新的构象，然后根据预先定义的评价函数对新构象进行取舍，通过不断迭代收敛于全局最优或近似最优构象。常用的启发式搜索算法包括蒙特卡洛（Monte Carlo, MC）算法、遗传算法（Genetic Algorithm, GA）、模拟退火（Simulated Annealing, SA）算法和群智能（Swarm intelligence, SI）算法等。

• 确定性搜索（Deterministic search）：确定性搜索，初始状态可以确定下一状态的改变趋势，即向能量减小的方向进行。如果体系初始状态完全相同，当使用相同的参数进行构象搜索，就会到达相同的终态。这种方法在进行构象搜索时一般很难越过势垒，因此容易陷入局部最小值。分子对接中常用的确定性搜索算法主要是分子动力学模拟，代表性程序为基于CHARMm力场的CDOCKER。

打分函数：打分函数是影响分子对接精度的另一个重要因素，其主要作用是评估配体和受体分子间的亲和力。用于评价同一分子的各种构象以及不同分子的最佳构象。

• 基于力场的打分函数（Forcefield）：基于力场的打分函数一般采用分子力场的非键作用能项来计算受体和配体之间的相互作用能。只考虑热焓对能量的贡献，不考虑熵的影响，一般情况下，采用标准力场的真空静电和范德华作用两项的加和来计算，有些打分函数还会考虑氢键项的贡献。

  • 基于力场的打分函数通式：\(ΔG_{binding}=ΔE_{vdw} + ΔE_{ele} + [ΔE_{Hbond}]\)

  • DOCK程序中采用AMBER的能量函数: \(E = \sum_{i=1}^{lig} \sum_{j=1}^{rec}(\frac{A_{ij}}{r_{ij}^{12}}-\frac{B_{ij}}{r^6_{ij}}+\frac{q_iq_j}{\epsilon(r_{ij})r_{ij}})\)

  • 其它常用的基于力场的打分函数：AutoDock的打分函数、Goldscore等

• 经验打分函数是受体和配体结合的自由能可以近似为不相关的不同能量项的贡献。通过计算并加和各个分能量项求得总的结合自由能。

  • 这些分能量项包括：氢键作用能、疏水作用能、离子化、熵的贡献、金属离子作用等

  • 常用的经验打分函数如PLP、ChemScore、GlideScore、LUDI等

  • ChemScore采用如下函数形式：\(ChemScore = S_{Hbond} + S_{metal}+ S_{lipophilic} + P_{rotor} + P_{strain} + P_{clash} + [P_{covalent} + P_{constraint}]\)

• 基于知识的打分函数，通过加和受体和配体之间的原子对统计势（Pairwise Statistical Potential）来计算配体和受体之间的结合能力。

  \[A = \sum_i^{lig}\sum_j^{pro}\omega_{ij}(r)\]

  • 距离依赖的原子对之间的作用势 \(ω_{ij}(r)\) 由下式计算得到：\(ω_{ij}(r) = -k_BTln[g_{ij}(r)] = -k_BTln[\frac{𝜌_{ij}(r)}{𝜌^∗_{ij}}]\); \(𝜌_{ij}(r)\) 是原子对i‐j在距离r时的密度；\(ρ^*_{ij}\) 是原子对在参考态时的密度

  • 常用的基于知识的打分函数：PMF、IT‐Score、DrugScore等

• 一致性打分（Consensus Score）：将几个不同打分函数的结果相结合，这种方法通常被称为“一致性评分”。 根据这种方法，使用多个打分函数来评分结合构象，相当于对对接结果进行二次打分（一致性打分）是非常有效的消除假阳性方法。

分子对接的一般流程

1. 受体结构的准备

• 靶点检验：X‐ray晶体结构要检查分辨率（分辨率越高越好）、种属、突变情况、结构的完整性等

• 受体结构处理（以晶体结构为例）：加氢、加电荷、带电残基的质子化状态、水分子的处理、助溶剂分子的删除、补充缺失的残基

2. 结合位点确定

• 复合物中的配体去除（最简单）

• 同源蛋白模型的结合位点信息

• 定点突变实验

• 程序算法：Cavity（DS）、SiteMap（Schrӧdinger）、SiteFinder（MOE）

3. 配体（小分子）结构准备

• 3D结构转化

• 加氢、加电荷、质子化状态（pH: 7.4）、手性和互变状态确定

• 如果是小分子化合物库，还需进行类药性分析和假阳性子结构过滤

4. 对接过程（参数设置）

• 对接算法的选择：使用重对接（native docking），将复合物晶体结构中的小分子重新对接到该复合物中，比较对接构象与晶体构象的均方根偏差值（RMSD），该值越小（通常要小于2 Å）说明对接构象与晶体构象越接近，对接算法对小分子结合构象的重现性越好，对接算法越可靠。

参考：J Mol Graph Model 2015 Vol. 60 Pages 142‐154

• 打分函数的评价（同一个受体结构）：

  • 同一配体的不同结合模式（pose）：对于重对接而言，比较各pose的打分值与其RMSD值的相关性，相关性越高，则说明打分函数越能反映对接构象的准确性。

  • 不同活性的配体：比较各配体打分最高的结合构象的活性值与打分值的相关性，相关性越高，说明该打分函数能较为准确的反映配体的活性（注：选择不同的配体活性值之间要至少相差一个数量级）。

• 其它参数设置：可先使用软件预设值，然后再根据初步对接结果（如重对接结果）进行调整

5. 对接结果分析

• RMSD值计算（重对接）

• 配体‐受体非键相互作用分析

• 配体‐受体相互作用2D图

• 受体结合口袋表面

• 残基相互作用分析

• 多种打分函数分析

• 结合能分析带溶剂化模型

分子对接的应用

确定配体和受体可能的结合构象，及其与靶标结合口袋内残基的相互作用细节，从而阐明活性化合物的作用机制，为分子的结构优化打下基础

对靶标进行基于结构的虚拟筛选，从而辅助发现对该靶标具有潜在活性的候选化合物

对同系列化合物进行分子对接，可以得到分子叠合构象，从而进行基于结构的3D‐QSAR研究

通过反向对接找到小分子可能作用的靶标

Discovery Studio中的分子对接算法简介

LibDock：LibDock根据小分子构象与受体相互作用的热区（Hotspot）匹配的原理将这些构象刚性对接到受体的结合口袋中，其最大的优势在于速度快，可以并行运算，适合于对大规模数据库进行快速精确的虚拟筛选

极性格点：受体中能与配体分子中可以形成氢键的原子相匹配的部分

非极性格点：受体中能与配体分子的非极性原子相匹配的部分

LibDock计算过程

CDOCKER：CDOCKER是基于CHARMm力场的半柔性对接程序，采用高温动力学来搜索配体分子的柔性构象空间，并采用模拟退火的方法将各个构象在受体活性位点进行优化，从而使对接结果更加精确。

CDOCKER计算过程

Flexible Docking：Flexible Docking可以实现受体‐配体的双柔性对接，模拟受体与配体的诱导契合效应，可以精细地研究受体‐配体相互作用信息，适用于作用机理的研究。

Flexible Docking计算过程

Discovery Studio中的非键相互作用分析

分子对接计算的注意点

小分子问题

• 小分子的起始构象、电荷和质子化状态对对接结果有一定影响

蛋白受体问题

• 如何选择靶标晶体结构？

• 定义活性位点的位置和大小？

• 活性位点氨基酸残基的结构和质子化状态是否正确？

对接问题

• 根据研究目的和研究体系，选择合适的对接算法

分子对接尚未解决的问题

1. 溶剂化效应

2. 熵效应

3. 蛋白的柔性

4. 打分函数的准确性

分子重对接实验

实验目的：

1. 学会使用 Discover studio 进行半柔性和柔性重对接。

2. 学会使用 CDOCKER, LIBDOCK, Flexible Docking 程序重对接并比较。

实验原理：

本实验使用 Discovery Studio 软件通过 CDOCKER, LIBDOCK 和 Flexible Docking 程序进行 3GEN 蛋白与自带的配体的重对接。

半柔性对接：受体保持刚性，配体则完全当作柔性分子进行处理，当前的大多数分子对接软件都采取这种模式。

柔性对接：对接过程既考虑小分子柔性，又可考虑受体部分区域（如结合口袋中的残基）的柔性。

基于 CHARMm 力场的 CDOCKER 半柔性对接程序使用确定性搜索算法。确定性搜索是通过初始状态确定下一状态的改变趋势，即向能量减小的方向进行。如果体系初始状态完全相同，当使用相同的参数进行构象搜索，就会到达相同的终态。这种方法在进行构象搜索时一般很难越过势垒，因此容易陷入局部最小值。CDOCKER 模拟退火的方法将各个构象在受体活性位点进行优化，从而使对接结果更加精确。

CDOCKER计算过程

LIBDOCK 程序根据小分子构象与受体相互作用的热区（Hotspot）匹配的原理将这些构象刚性对接到受体的结合口袋中，其最大的优势在于速度快，可以并行运算，适合于对大规模数据库进行快速精确的虚拟筛选。

LibDock计算过程

Flexible Docking 可以实现受体‐配体的双柔性对接，模拟受体与配体的诱导契合效应，可以精细地研究受体‐配体相互作用信息，适用于作用机理的研究。

Flexible Docking 计算过程：

本实验所用软件环境：

• DS Version：19.1.0.18287

• PP Version：19.1.0.1963

• DS Client Version：19.1.0.18287

• OS Distribution：Windows

• OS Version：10.0.22000

一般流程：

大分子结构的准备：通过实验方法所测得的生物大分子三维结构或多或少存在一些缺陷，比如分辨率低导致模型质量不佳、晶体结构无 H 原子和电荷信息、在模型构建步骤中相对粗糙的力场带来的结构误差等，这就需要在做 SBDD 前对生物大分子（如蛋白）的三维结构进行准备和能量最小化。

实验步骤：

◆ 靶点结构的获取及准备

1. 查找靶点结构：本实验中搜索 3GEN 蛋白作为靶点。进入 PDB 网站(https://www.rcsb.org/），搜索靶蛋白。本实验中搜索 3GEN，选择 3GEN，点击 Download Files，再点击 PDB Format， 下载靶蛋白结构文件 ，并用 Discovery Studio 打开下载好的文件。（如图所示）

2. 靶蛋白的准备：

   去除晶胞：点击Discovery Studio软件上菜单栏上的Structure –> Crystal Cell –> Remove Cell进行晶胞的去除。

   定义小分子配体：Discovery Studio软件上，选上小分子配体，也就是单击Ligand Groups， 再点击Receptor-Ligand Interactions –> View Interactions –> Define Ligand: <undefined>来进行小分子配体的定义。

   蛋白结构准备：Discovery Studio 软件上，点击 Macromolecules –> Prepare Protein，设置参数如下，来进行蛋白结构的准备得到新窗口3GEN_prep。接下来的操作都是在新的窗口当中进行。

   

   定义结合位点：Discovery Studio软件上，选上小分子配体，也就是单击Ligand Groups， 再点击Receptor-Ligand Interactions –> Define and Edit Binding Site –> From Current Selection来进行结合位点的定义。

◆ 小分子配体的获取及准备

1. 配体小分子的获取：本实验中使用 3GEN 自带的配体。Discovery Studio 软件上，点击菜单栏上的 File → New → Molecule Window 创建新窗口，剪切 3GEN 自带的配体粘贴。

2. 小分子的准备和能量优化：

   点击 Discovery Studio 软件上的 Small Molecules → Prepare or Filter Ligands → Prepare Ligands 进行小分子的准备，设置参数如下。

   

   点击 Discovery Studio 软件上的 Simulation → Change Forcefield → Apply Forcefield；

   

   再点击 Discovery Studio 软件上的 Simulation → Run Simulations → Minimization 进行能量优化。设置参数如下：

   

◆ 半柔性分子对接

1. LIBDOCK 分子对接：确保当前窗口为准备好的大分子，即 3GEN_prep，删除配体，点击 Discovery Studio 软件上的 Receptor-Ligand Interactions → Dock Ligands → DockLigands (LIBDOCK) 进行 LIBDOCK 半柔性分子对接。设置参数如下：

2. CDOCKER 分子对接：确保当前窗口为准备好的大分子，即 3GEN_prep，删除配体，点击 Discovery Studio 软件上的 Receptor-Ligand Interactions → Dock Ligands → Dock Ligands (CDOCKER) 进行 CDOCKER 半柔性分子对接。设置参数如下：

◆ 柔性分子对接

确保当前窗口为准备好的大分子，即 3GEN_prep，删除配体。点击选上 Discovery Studio 软件菜单栏上的 View → Explorers → Protocols，再然后双击 Discovery Studio软件上的 Protocols → Discovery Studio → Receptor-Ligand Interactions → Docking → FlexibleDocking 进行 Flexible Docking 柔性分子对接。设置参数如下：

◆ 结果分析

1. RMSD 值计算：新建窗口，即点击 File→New→Molecule Window，将对接结果中的对接 pose 和原晶体结构中的小分子配体复制粘贴到新窗口中，把原晶体结构中的小分子配体设为参考分子。选中参考分子，依次选择菜单栏上的 Structure → RMSD → SetReference 来设定参考分子。在需要计算的 pose 边上打勾（这里假设计算全部），并确保呈选中状态，依次选择菜单栏上的 Structure → RMSD → Heavy Atoms 进行 RMSD 值计算。

2. 对对接结果进行重打分：在 Discovery Studio 中打开对接结果窗口（应该是 3GEN窗口），点击 Receptor-Ligand Interactions → Dock Ligands → Score Ligand Poses 进行对配体poses 打分，设置如下图，会生成新窗口。

3. 在新窗口中，点击选上 Discovery Studio 软件菜单栏上的 View → Explorers → Protocols。再然后从 Protocols 窗口中双击 Discovery Studio→ Receptor-Ligand Interactions → Scoring and Analysis → Analyze Ligand Poses 进行配体poses 分析，设置参数如下图。

实验结果：

1. 靶蛋白的准备

2. 小分子配体的准备

3. 小分子配体的能量优化

4. LIBDOCK 分子对接

   • LIBDOCK 分子对接

   • RMSD 值计算

   Heavy Atom RMSD to 3GEN_ligand_origInal 32

   Name	Reference	RMSD (A)
   3GEN_ligand 1	3GEN_ligand_origInal 32	0.1735
   3GEN_ligand 7	3GEN_ligand_origInal 32	0.6876
   3GEN_ligand 38	3GEN_ligand_origInal 32	7.3890
   3GEN_ligand 39	3GEN_ligand_origInal 32	0.9175
   3GEN_ligand 40	3GEN_ligand_origInal 32	3.4741
   3GEN_ligand 41	3GEN_ligand_origInal 32	8.9469
   3GEN_ligand 42	3GEN_ligand_origInal 32	3.2732
   3GEN_ligand 43	3GEN_ligand_origInal 32	3.2962
   3GEN_ligand 50	3GEN_ligand_origInal 32	8.4331
   3GEN_ligand 52	3GEN_ligand_origInal 32	3.4042
   3GEN_ligand 53	3GEN_ligand_origInal 32	9.6773
   3GEN_ligand 59	3GEN_ligand_origInal 32	11.3117
   3GEN_ligand 60	3GEN_ligand_origInal 32	3.5681
   3GEN_ligand 63	3GEN_ligand_origInal 32	3.2734
   3GEN_ligand 64	3GEN_ligand_origInal 32	3.3751
   3GEN_ligand 65	3GEN_ligand_origInal 32	4.4569
   3GEN_ligand 68	3GEN_ligand_origInal 32	3.3365
   3GEN_ligand 69	3GEN_ligand_origInal 32	0.6875
   3GEN_ligand 74	3GEN_ligand_origInal 32	8.1215
   3GEN_ligand 76	3GEN_ligand_origInal 32	8.9682
   3GEN_ligand 77	3GEN_ligand_origInal 32	11.0938
   3GEN_ligand 79	3GEN_ligand_origInal 32	11.0721
   3GEN_ligand 80	3GEN_ligand_origInal 32	3.3874
   3GEN_ligand 85	3GEN_ligand_origInal 32	3.3914
   3GEN_ligand 91	3GEN_ligand_origInal 32	9.0667
   3GEN_ligand 93	3GEN_ligand_origInal 32	3.3918
   3GEN_ligand 95	3GEN_ligand_origInal 32	11.0424
   3GEN_ligand 97	3GEN_ligand_origInal 32	9.5115
   3GEN_ligand 98	3GEN_ligand_origInal 32	10.3211
   3GEN_ligand 99	3GEN_ligand_origInal 32	3.4550
   3GEN_ligand 100	3GEN_ligand_origInal 32	8.2923

   • 对 LIBDOCK 对接结果进行重打分 Score Ligand Poses

   • 对 LIBDOCK 对接结果进行重打分 Analyze Ligand Poses

5. CDOCKER 分子对接

   • CDOCKER 分子对接

   • RMSD 值计算

   Heavy Atom RMSD to 3GEN_ligand_origInal 11

   Name	Reference	RMSD (A)
   3GEN_ligand 1	3GEN_ligand_origInal 11	0.2238
   3GEN_ligand 2	3GEN_ligand_origInal 11	0.2856
   3GEN_ligand 3	3GEN_ligand_origInal 11	0.3219
   3GEN_ligand 4	3GEN_ligand_origInal 11	0.5222
   3GEN_ligand 5	3GEN_ligand_origInal 11	0.2150
   3GEN_ligand 6	3GEN_ligand_origInal 11	0.5501
   3GEN_ligand 7	3GEN_ligand_origInal 11	0.4920
   3GEN_ligand 8	3GEN_ligand_origInal 11	0.6795
   3GEN_ligand 9	3GEN_ligand_origInal 11	0.5764
   3GEN_ligand 10	3GEN_ligand_origInal 11	0.5110

   • 对 CDOCKER 对接结果进行重打分 Score Ligand Poses

   • 对 CDOCKER 对接结果进行重打分 Analyze Ligand Poses

6. FlexibleDocking 分子对接

   • FlexibleDocking 分子对接

   • RMSD 值计算

   Heavy Atom RMSD to 3GEN_ligand_origInal 7

   Name	Reference	RMSD (A)
   3GEN_ligand 1	3GEN_ligand_origInal 7	0.5484
   3GEN_ligand 2	3GEN_ligand_origInal 7	9.8311
   3GEN_ligand 3	3GEN_ligand_origInal 7	11.6764
   3GEN_ligand 4	3GEN_ligand_origInal 7	7.3178
   3GEN_ligand 5	3GEN_ligand_origInal 7	6.2607
   3GEN_ligand 6	3GEN_ligand_origInal 7	11.0410

   • 对 FlexibleDocking 对接结果进行重打分 Score Ligand Poses

   • 对 FlexibleDocking 对接结果进行重打分 Analyze Ligand Poses

讨论：

1. LibDock

   • 各pose的打分值与其RMSD值的相关性并不高，几乎没有相关性。

   

2. CDOCKER

   • 各pose的打分值与其RMSD值的相关性并不高，几乎没有相关性。

   

3. FlexibleDocking

   • 各pose的打分值与其RMSD值的相关性并不高，几乎没有相关性。

   

基于片段的药物设计(Fragment‐based drug design)(FBDD)

引言

结构优化的困惑

• 许多药物靶标的活性位点是由多个口袋组成。高通量筛选的化合物过于“成熟”，通常不能与靶蛋白很好地结合，而对于其中单个片段的优化往往会影响整个分子，甚至导致结合位置的改变高通量筛选或虚拟筛选得到的苗头化合物与靶标蛋白的结合模式。

• 苗头化合物或先导化合物的优化，大多以加入或变换基团，以增加与靶标结合的机会和强度。一般“不敢”轻易除去基团或片段，以免丢失参与结合的原子或基团（即药效团）。

换一种思路：从分子量较小的片段分子出发！

由高通量筛选获得的苗头化合物与先导化合物在分子大小、亲酯性和极性表面积等理化参数上几乎一样，但是活性却比先导物差了很多，这就说明了HTS苗头化合物的成药质量差（药效和理化性质的不平衡）。相反，片段分子的活性虽弱，但可覆盖更全 面的化学空间，可在优化过程中保持分子的成药质量。

基于片段的药物设计原理

基于片段的药物设计（Fragment‐based drug design, FBDD）是一种将高通量筛选和基于结构的药物设计相结合的药物发现新方法。由于筛选的起始物是较小的分子，与受体结合的自由度与可能性较大，这是有利的一面，但因为只与受体部分位点结合，显示的活性强度要弱于较复杂的分子量较大的化合物，这是有待提高的一面。低活性的片段分子，可借助与受体结合的结构特征，指导添加或连接有助于提高活性的基团或片段，同时控制分子的大小，以形成有成药前景的化合物。

设计理念

• 药物靶标蛋白的活性位点是由多个口袋（子口袋）组成（前提假设）

• 寻找能与靶标蛋白各个口袋特异性结合的片段（目标）

• 将上述各个片段用合适的连接子连接起来，组成高活性的化合物。片段连接后会引起结合自由能跳跃式下降，从而导致亲和力大幅度提高。

筛选得到的片段化合物通过结构优化成为先导物，以及最终成为候选药物，在提高活性的同时，尽量不要过多增加分子尺寸和亲脂性，为此需要在操作过程中用参数加以表征或“监视”活性与分子大小之间的关系。有如下参数：

• 配体效率（Ligand Efficiency, LE）配体效率是指配体中每个原子对结合能的贡献：

\[ΔG = 1.37pK_d; LE = ΔG/N_{非氢原子}\]

LE是将配体和受体的结合能力同分子大小相结合的参数，用以衡量苗头物或先导物，以及优化的化合物的质量如某化合物含41个非氢原子（MW：540），当 Kd = 10 nM，LE = 0.27 kcal/mol；而含有30个非氢原子的化合物，若 Kd = 10 nM，LE = 0.37 kcal/mol。Reynolds等分析了8000个配体同20个靶蛋白的结合常数，计算了配体效率，发现相对分子量差别较大的化合物活性与（非氢）原子数之间并不呈线性关系

• 配体‐亲脂性效率（Ligand‐Lipophilicity Efficiency, LLE）:LLE用来表征先导物和优化的质量：LLE = pIC50 (或pKi) – clogP (或clogD)。由于亲脂性对药物的溶解、溶出和代谢稳定性不利，所以高LLE值为结构优化的标志。

契合质量：Reynolds等提出了配体的契合质量的概念（fit quality, FQ），以消除因分子量增大，LE的变化被掩饰和拉平的效应。FQ的计算方法是将化合物的非氢原子数进行归一和标度化，得到相应的LE‐Scale：\(LE‐Scale = ‐0.064 + 0.873e^{(‐0.026×HA)}\) 其中HA表示Heavy Atoms（非氢原子）, 配体的契合质量FQ按照下式计算：FQ = LE/LE‐Scale

基于片段的药物设计研究方法

片段库的构建：用于高通量筛选的化合物多为类药性分子（drug‐like），大都遵循类药5规则。而用于FBDD筛选的化合物是分子小结构而简单的类先导物（lead‐like），其一般特征为：

• 分子多样性

• 有水溶性

• 合成可修饰性

• 药化中的优势结构或骨架

• 遵循类药3规则（分子量低于300，clogP低于3，氢键供体、受体和可旋转键分别不多于3）

片段筛选技术:

• NMR技术

• X‐射线单晶衍射

• 高浓度筛选

• 表面等离子共振技术（SPR）

• 热漂移检测

• 恒温滴定量热法（ITC）

发现活性的片段分子仅是研究的第一步，将活性片段转化为先导化合物甚至是候选药物才是基于片段的药物设计的最终目的。

从片段到先导物的设计方法：

• 片段生长（fragment growth）：以受体结合的第一个片段为核心，经过理性设计，在邻近处逐渐生长成活性强的较大分子。片段生长实例——FBDD首个上市药物 Vemurafenib 的发现历程。

参考：Nat Rev Drug Discov, 2013, 12: 5

• 片段连接与融合（fragment linking and fusion）:片段连接：与受体结合的相邻的两个片段经连接基团连接成活性强的较大的分子。片段融合：与受体结合的相互交盖或接近的两个片段合并成活性强的较大的分子。片段连接与融合实例——结核分枝杆菌泛酸合成酶(MtPS)命中片段的连接。

参考：https://doi.org/10.1021/ja100595u

通过片段融合设计——WDR5−MYC相互作用抑制剂。

参考：https://doi.org/10.1021/acs.jmedchem.0c00224

• 片段自组装（fragmentself‐assembly）

分别结合在活性位点中相邻口袋的两个活性片段含有可相互反应的基团，这两个片段可自发地连接成为高活性的化合物。

片段自组装实例：

优点：

1. 片段库仅需几百‐几千个分子

2. 化学结构新颖且多样，命中率高

3. 从片段优化至候选药物，分子量和活性同步增长，更加符合新药发现的一般规律，可行性强

缺点：晶体结构较难获得

适合于片段筛选的算法

MCSS：Multiple Copy Simultaneous Search, MCSS是一种基于片段的分子对接技术，由Miranker与Karplus开发（1991）。片段先被随机分布在结合位点中，然后程序采用CHARMm对这些随机片段进行能量优化，通过采用独立的MCSS_Score来打分和排序，以找到最适的片段位置。

• \(Score_{MCSS} = ‐(Energy_{Ligand} + Energy_{Ligand‐Protein})\)

• 计算效率高（多片段同时优化，片段之间互不影响）

MCSS的计算流程：

Ludi是应用最为广泛的基于片段的药物设计方法之一，它帮助我们找到新的分子骨架或者修改已有的分子骨架来提高小分子的活性。

Discovery Studio中的De Novo Protocol使用的就是Ludi算法，包括：

• De Novo Receptor：可以用来寻找与受体相互作用的全新片段，这种方法的优点是可以快速有效的对数据库中大量的片段进行筛选

• De Novo Link：可以将找到的活性片段连接或在已有骨架的基础上添加新的片段分子

• De Novo Evolution：也称为AutoLudi，主要用于进行Me better 药物设计工作，不仅可以搜索活性片段分子并且可以将片段和骨架自动连接直接产生新的分子。

基于Ludi进行FBDD的基本流程：

优点：

• 配体结构可能是全新的，不受现有知识的约束，也不受人的思维束缚

• 该方法既可用于先导化合物的发现，也可用于对先导化合物的结构优化

缺点：

• 可能只是理论分子，仅对接有效，且可能不易合成。

FBDD筛选和优化流程:

基于配体的药物设计(Ligand‐based drug design, LBDD)

基本概念

基于结构的药物设计（SBDD）：在已知生物大分子靶点结构的情况下，直接考虑药物与靶点的相互作用来进行药物设计。量体裁衣。

基于配体的药物设计（LBDD）：在生物靶点结构未知的情况下，通过研究与靶点具有特异性结合配体的结构信息，发现先导化合物的方法。量衣裁衣。

基于配体的药物设计（LBDD）：

• 定量构效关系（QSAR）：

  • 二维定量构效关系（2D‐QSAR） Hansch方法

  • 三维定量构效关系（3D‐QSAR） CoMFA方法

• 相似性搜索

• 药效团模型（Pharmacophore）

相似性搜索 Similarity search

相似性搜索的概念

分子相似性（molecularsimilarity）的概念最早于十九世纪八九十年代的早期，是现在许多虚拟筛选技术的理论基础。相似性原理指出，总体相似的分子应具有相似的生物活性。

相似性搜索主要有以下三类：

• 子结构搜寻

• 基于分子指纹的相似性搜索

• 基于分子形状的相似性搜索

分子相似性方法中的基本要素

描述符(descriptors)是分子相似性方法中的基本要素。目前已经有大量的描述符应用于分子相似性研究中，从描述符性质特征上看，描述符主要分为以下3个方面：

• 1D 描述符：由化合物本身属性衍生而来，如表示分子物化属性的 lgP，摩尔折射率等；

• 2D 描述符：由 2D 分子图形或者结构片断计算得来的，如拓扑指数，2D 分子指纹，连接表，图(或子图)，(子)结构描述符等；

• 3D 描述符：分子形状，3D 分子指纹、分子总表面积和电压等。

子结构搜寻（Substructure search）

子结构搜寻专注于分子拓扑结构上是否具有特定的官能团以及是否具有某些特定的对于分子识别至关重要的原子排布。

分子指纹（Molecular Fingerprint）方法

分子指纹是一种分子的抽象表征，它将分子转化（编码）为一系列比特串（即比特向量，bit vector），然后可以很容易地在分子之间进行比较。

• 分子指纹上的每个比特位对应于一种分子片段，假设相似的分子之间必然有许多公共的片段，那么具有相似指纹的分子具有很大的概率在 2D 结构上也是相似的。

相似性评估方法

• 评价两个分子指纹（向量）之间最常用的评估方法称为“ Tanimoto系数 ”。使用以下等式从两个分子A和B的位阵列指纹计算Tanimoto系数：

\[T=\frac{N_{A\&B}}{N_{A}+N_{B}-N_{A\&B}}=0.80\]

\(N_A、N_B\) 分别等于A、B中比特值为1的位数，\(N_{A\&B}\) 等于A和B中公共的比特值为1的位数

分子指纹的类型

根据将分子表征转换成位串方法的不同，分子指纹可以分为几种类型。 大多数分子指纹的编码方法仅使用了2D的分子图形信息，因此称为2D指纹;一些方法能够存储3D信息，比如药效团指纹。常见的2D分子指纹方法包括

参考：http://blog.molcalx.com.cn/2019/01/29/fingerprint.html

• 基于子结构的指纹(substructure key‐based fingerprint)

• 基于拓扑或路径的指纹(topological or path‐based fingerprint)

• 圆形指纹(circular fingerprint)。

基于子结构的指纹根据给定结构列表中某些子结构或特征的存在与否来设置位串。如MACCS和PubChem指纹。

一个10 BIT长度的指纹图谱示意图：每个位置编码一个子结构, 有3个位置所代表的子结构出现在上图的分子中（即被圆圈圈中的那些子结构），这些位置被标为1；其它的位置代表的子结构没有出现在分子中，则被标为0。

基于拓扑或路径的指纹通过分析从一个原子开始直至到达指定数量键的路径（通常为线性）上所有的分子片段，然后对每一个路径中的片段产生指纹。OpenEye的Tree指纹（Tree fingerprint）。

此类指纹适用于任意一个分子，并可以调整其长度，可以用于快速的子结构搜索与分子过滤。

圆形指纹是记录每个从原子出发直到指定半径内的环境。因此，此类指纹不适用于子结构查询（因为相同的片段可能具有不同的环境），但广泛用于完整结构的相似性搜索。如ECEP、FCFP、Molprint2D等。

以一个重原子为中心，搜寻在特定半径范围（有的实现用直径）内的结构特征，比如力场的原子类型，官能团，片段等等各种信息。

药效团指纹通常以类似于基于子结构的指纹方式编码分子的结构特征，但同时考虑了这些特征之间的距离，通常按距离范围对其进行分类生成比特位串。比如SYBYL软件包中的Tuplet。

Tuplet药效团指纹图谱示意图：一个分子不同的构象，将每个构象的药效团特征按聚类分类、编码为比特串。

计算分子指纹并用于虚拟筛选的常见软件：

开源、免费：Open Babel、RDKit、CDK、Indigo、jCompoundMapper

商业：OEChem、 TK/OpenEye、JChem/ChemAxon、Canvas/Schrodinger、MOE、Pipeline Pilot、SYBYL‐X

二维（2D）分子指纹方法的特点

采用2D指纹方法进行基于2D相似性的化合物搜索

• 指纹表征的是二进制格式的2D分子结构。

• 不需要生成用于筛选的分子3D构象，通过分子指纹之间的重叠来计算分子之间的结构相似性，2D指纹方法由于其计算效率和有效性而成为基于分子2D相似性筛选的高效方法。

• 指纹方法的有用性也被提出用于骨架跃迁研究。

分子形状相似性的方法

分子形状（Molecular Shape）是药物化学中的基本概念，是分子的3D形状的比对，其基于如果两个分子具有相似的形状，那么它们也可能具有相似的性质的原理。

• 形状相似性通常通过比对方法来评估，该方法寻求两种形状三维重叠的最大化，或者通过特征向量方法，将形状变换成能够有效地比较的特征的低维向量。

• 作为相似性计算的一部分，可以用静电或药效团特征进一步注释分子形状。常见的基于形状的相似性搜索算法有ROCS、Phase Shape等

形状相似性搜索涉及分子比对，分子间以最大化重叠体积或特征点进行比对，以找出两个分子的最佳叠合。

相似性搜索的应用

靶标预测：基于计算的靶标预测方法在药物发现中起重要作用。其中，分子相似性方法作为一种有希望的方法，其基于结构相似的分子往往具有与靶标相似的反应的假设，已成功应用于许多药物相关研究。

基于配体3D形状的结合构象预测：为了提高分子结合构象预测的准确性，研究人员开发了一种使用形状相似度进行构象预测的方法。

分子形状相似性应用于虚拟筛选：形状筛选只需要一个已知活性查询化合物及其构象（最好为药效构象），相对于基于分子对接和药效团的筛选更加节省时间和资源。

参考：Bioorg. Med. Chem. Lett. 25 (2015) 3442–3446

药效团模型 Pharmacophore model

药效团的基本概念

一个特定受体的药物小分子都具有某种特定的特征结构，换句话说，只要具备这些特征，一个化合物就有可能是这个受体潜在药物分子。这些特征点的集合，我们就称之为药效团。

药效团的由来

• 药物分子与受体靶点发生作用时，要与靶点产生几何匹配和能量匹配，此时药物分子的构象称为活性构象（药效构象）。

• 药物分子中的基团对于活性的影响不同，有些基团的变化对药物与靶点之间的相互作用影响很大，有些则影响不大。

• 对于这种差异的研究，人们发现具有相同活性的分子往往具有相同的某些特征。

• 只要具备这些特征，一个化合物就有可能是这个受体的潜在药物分子。

药效团（Pharmacophore），是指药物分子中可以与受体结合位点形成氢键、静电、疏水等非键相互作用的原子或原子团，及其相互之间的空间关系。具体官能团

IUPAC对药效团的定义：药效团是和特异的靶标结合并引发（或阻断）靶标生物响应所必需的一组立体和电子特征。抽象特征

药效团从结构层面上揭示了药物与受体结合并启动药理活性的微观特征，也就是说，一组活性分子共同呈现特定的药理活性的基础是由于有相同的药效团特征和分布

药效团模型方法

• 药效团模型方法一般包括两个层面的内容：即药效团模型的建立和基于药效团模型的数据库搜索。

• 药效团模型的建立仅仅是得到药效团模型。

• 通过药效团模型，希望能够找到新的先导化合物，这就需要采用基于药效团模型的数据库搜索。

• 通过数据库搜索，来寻找包含特定药效团特征的化合物，这些具有特定药效团特征的化合物可能具有相应的生物活性。

药效特征元素

进行药效团模建的前提是发现和确证配体含有的药效特征元素（药效团特征）

• 典型的药效团特征一般用球体来表示，球心为准确位置，半径为可偏离的范围。具体特征包括：氢键供体、氢键受体、正电中心、负电中心、疏水中心、芳环中心、排除体积、形状约束等。此外，还包括有自定义的其它特征。

氢键供体（hydrogen bond donor）

主要包括H以及与之相连的O和N，一般包括：

• 非酸性羟基（非羧酸或酚类中的羟基）。羧酸、酚类中的羟基之所以显酸性是由于羟基中氧的未成对p电子与苯环、羰基形成p‐π共轭体系，致使氧上电子云密度降低，氢很容易解离。

• 氨基。

• 次氨基，但不包括三氟甲基硫酰胺和四唑中的次氨基。

氢键受体（hydrogen bond accepter）

广义来讲，任何带有孤对电子的原子，如N、O、F、S。在一般的药效团模型方法中，仅考虑药物分子中最常见的氢键受体形式，包括：

• sp3或sp2杂化的O（与C等原子以单键或双键相连的O）

• 与C以双键形式相连的S

• 与C以双键或三键相连的N

氢键的方向性

由于配体与受体之间的氢键相互作用一般具有明确的方向性，因此一般都采用两个点来描述氢键特征，一个点表示氢键特征中重原子的空间位置，而另一个点表示氢键给体或受体的矢量方向。

对于氢键受体，矢量方向一般为重原子和其上孤对电子连线的方向。

对于氢键给体，矢量方向为重原子和与之相连的氢原子形成键的方向。

疏水中心（hydrophobic area）

配体与受体上的疏水基团总是倾向于形成紧密的疏水堆积作用，形成疏水性内核。

• 疏水相互作用本质上包含了熵效应和范德华相互作用。

• 疏水基团一般由非极性原子组成，有疏水相互作用的片段很多，如甲基、乙基、苯环等。

• 在药效团模型方法中，只需要用一个点表示。

芳环中心（aromatic ring）

芳环可以参与药物分子和蛋白受体之间的π电子离域系统的π‐π相互作用。

• 芳环中心主要包括五元和六元芳环，如噻吩、苯环等。

• 芳环需要由两个参量来定义：一个参量是芳环的空间位置，即芳环中所有原子的几何中心；另一个参量是芳环平面矢量方向，一般用垂直于芳环平面的矢量来描述。

电荷中心（charge interaction）

配体上的电荷中心是指配体上的带电基团，由于具有较多的部分电荷，这些基团往往可以和受体形成盐桥或较强的静电吸引作用。

• 电荷中心既可以是带有电荷的原子，也可以是在生理pH下会发生电离的中性基团，如氨基、羧基。

• π电子离域系统，如羧酸盐、胍基、脒基等也可能形成电荷中心。

正电荷中心

• 带正电荷的原子

• 伯、仲、叔脂肪胺中的N原子

• 氮‐氮双取代的脒基中的亚胺N原子或四氮取代的胍基中的亚胺N原子

• 至少含有一个未取代氢原子的脒基中的N原子中心或至少含有一个未取代氢原子的胍基中的N原子中心。

负电荷中心

• 带负电荷的原子

• 三氟甲基磺酰胺中的N原子

• 羧酸、亚磺酸或磷酸中羟基氧和氧代氧的原子中心

• 磷酸二酯和磷酸酯中羟基和氧代氧的原子中心

• 硫酸和磺酸中羟基氧和两个氧代氧的原子中心

• 磷酸单酯和磷酸中氧代氧和两个羟基氧的原子中心

• 四唑中的氨基氮原子

几何约束

药效特征元素可以抽象为点（如：疏水中心、电荷中心）、线（如：氢键）、面（如：芳环）的形式。一个完整的药效团模型中除了必须包含药效特征元素之外，还需要包括药效特征元素之间的空间约束，这些约束是指各特征元素的位置约束，各特征元素之间的距离、角度、取向等。

排除体积（excluded volumes）

• 在配体和受体相互作用时，在配体的某些取代位置上存在某些原子或原子团可能会和受体产生不利的原子碰撞，这些位置上的原子或原子团占有的位置就构成排除体积。

• 在排除体积中存在原子或原子团会大大降低化合物的活性。

药效团模型的应用

• 虚拟筛选

• 全新药物设计

• 先导物优化

• 多靶点药物设计

• 骨架跃迁

• 靶点预测（反向找靶）

常用的药效团模型构建工具

• Catalyst DS

• DISCO（distance comparison） Sybyl

• GASP（genetic algorithm similarity program） Sybyl

• Phase Schrodinger

• LigandScout

• Pharmacophore in MOE

Discovery studio中的药效团模型

药效团模型

药效团模型	药效团模型
基于配体	HipHop Hypogen
基于受体	Structure‐based Pharmacophore(SBP)
基于复合物	Complex‐based Pharmacophore(CBP)

基于配体的药效团模型构建与筛选的基本步骤

1. 活性化合物的选择及药效特征元素的定义

• 对于用来搭建药效团模型的化合物，其结构刚性要大一些，即分子的构象要少一些，这样叠合得到的药效团模型会较为明确。

• 所选分子应该具有尽量多的结构多样性。如果药效团识别方法只能识别分子公共的药效特征，则选择的分子要具有尽量高的活性。

• 如果药效团识别方法可以根据药效团模型和分子之间的叠合情况对分子的活性进行定量的预测，则选择的分子需要具有准确的生物活性。

2. 构象分析

• 对每个化合物都要进行构象分析，得到在某一能量范围内的构象。

• 这个能量范围一般是根据经验来指定，比如在CATALYST中，默认条件下，构象分析会保留和最低能量构象相差20kcal/mol之内的所有构象。

• 构象分析保留的构象要尽量少一些。

• 在一般的药效团识别软件中，都会对保留的构象数有明确的限制。

3. 分子叠合和药效团映射

• 得到药效团特征元素以及构象以后，就需要对这些构象进行叠合以找到共同的药效团模型。

• 在叠合时，药效特征元素作为分子间叠合的叠合点。

• 分子叠合往往会提供多个药效团，或者说分子存在多种公共的药效特征元素空间分布形式。

• 要想得到较为明确的药效团模型，则要尽量选择结构差别较大的分子。

4. 药效团模型的修正

• 分子叠合得到的药效团往往不是最优的，可能还需要根据自己的经验和其他实验或计算的结果对药效团模型做进一步的修正。比如加入排除体积和形状限制。

基于配体的药效团

• HipHop：基于分子共同特征的药效团模型。用于发现一系列配体小分子所共有的化学特征，并基于这些共同特性结构的比对叠合自动生成药效团模型，用户可以使用共有的特征药效团去搜索化合物数据库来寻找可能的先导分子

HipHop对于训练集的要求

• 输入的分子结构具有多样性

• 化合物数目值2～32个，6个左右比较理想

• 输入的结构需知活性强弱（强、中、弱）

• 需要包含principal和MaxOmitFeat性质参数

principal和MaxOmitFeat性质参数意义

• Principal 属性定义了分子的活性水平：

数值	活性水平	描述内容
2	有活性	参考分子，分子中所有化学特征在构建药效团模型时都要考虑。
1	中等活性	定位药效团特征元素时需要考虑该化合物的构象空间。
0	非活性	该分子在定位药效团特征元素时不考虑，用于选择性地定义排除体积。

• MaxOmitFeat 属性定义了每个分子中允许不与药效团模型匹配的特征元素的个数：

数值	描述内容
0	构建的药效团模型中所有特征元素都必需与化合物匹配上。
1	构建的药效团模型中允许有1 个特征元素不与化合物匹配。
2	构建的药效团模型中所有特征元素都无需与化合物匹配。
‐1	考虑所有药效团特征元素的子集

• HypoGen：具有活性预测能力的药效团。可以基于一系列针对特定生物靶标具有明确活性数值的化合物构建出具有活性预测能力的药效团模型。该算法首先构建得到活性分子能共享而非活性分子则不能共享的初始药效团模型，然后再经过模拟退火进行模型的进一步优化。最终构建得到的模型可以预测化合物的活性，以及指导化合物的优化以提高其活性。

HypoGen对训练集分子的要求

• 分子结构兼具多样性

• 分子的活性值至少跨越4个数量级

• 每个数量级的活性分子至少3个，化合物总数在16～31个

• 结构类似的化合物之间活性相差至少一个数量级

• 活性相似的化合物之间结构不同

• 需要包含Activ和Uncert性质参数

HypoGen算法流程

1. “构建”步骤

• 识别活性分子

• 识别活性最大的两个分子的所有药效团

• 用活性分子筛选药效团

2. “删除”步骤

• 由“构建”步骤产生的药效团如果出现在大部分非活性分子中，删除

3. 优化步骤

HipHop和HypoGen的区别

	HipHop	HypoGen
训练集数量	2～32	〉16
训练集是否有活性	全是活性分子	活性分子与非活性分子
活性值	不需要	需要
模型的类型	定性，只基于公共特征	定量，可预测活性值
模型的比较	训练集分子结构多样性	训练集分子结构和活性多样性
药效特征元素的个数	最多10个	最多5个

基于受体结构的药效团——Structure Based Pharmacophore (SBP)

SBP：基于经典的LUDI算法，从已知或假设的蛋白活性位点的特性直接得到相互作用位点图（feature interaction map），并且将这个信息转化成适用于快速三维数据库检索的药效团模型。

SBP 方法首先通过对活性部位的分析产生特征相互作用图，特征相互作用图用于反映可能的配体与蛋白之间的合理的相互作用，接着通过特征相互作用图产生药效团模型，进一步利用药效团模型发现潜在的活性化合物。

SBP操作流程

1. 导入蛋白结构及定义活性位点

2. 产生相互作用模型及药效团模型

3. 对feature聚类

4. 产生最终模型

Complex Based Pharmacophore (CBP)

CBP：基于LigandScout算法，通过识别受体‐配体相互作用的关键特性构建药效团模型进行先导物的发现与优化

该类药效团从受体配体间的相互作用出发，能更加准确直接地反应出它们之间的药效团特征

基于CBP的重要应用——反向找靶

一个分子与多个代表各靶标蛋白的药效团模型相互匹配，按照打分高低来判定可能与哪个靶标相互作用

PharmaDB：市场上最大的受体‐配体复合物药效团数据库

• 共25万多个药效团模型

• 基于scPDB中16000多个复合物晶体结构

药效团模型的验证

药效团模型的好坏直接影响经虚拟筛选后得到的候选化合物质量。为此，有必要进行药效团模型的验证，其目的就是为了评价所构建的药效团模型能否准确的将活性化合物与非活性化合物区分开。目前主要的验证方法有 3 种：

• Fischer’s randomization 验证方法

• 测试集验证方法

• Decoy 验证方法

Fischer’s randomization验证方法是用于检测我们所构建的模型是否具有较好的统计关联。在该方法中，原有训练集中所有化合物与各自的活性关系将会被自动打乱，并进行随机的重新分配，从而组成新的训练集。在参数设置和方法与之前相同的前提下，利用新的训练集构建药效团。统计显著性利用以下公式计算：Significance = 100(1‐(1+x)/y) 其中，x表示随机产生的药效团total cost值低于原始训练集的数目；y表示包括原始药效团在内的所有药效团数目。

如右图，19组随机得到的药效团模型其total cost值均明显高于初始模型（Hypo1），即说明该模型的统计显著性远远优于其它所有随机得到的药效团模型。

测试集验证：一个成功的药效团不仅要求能准确预测训练集中所有分子的活性，而且还必须具备准确预测训练集以外活性化合物活性的能力。

• 故需要选取一定数量（多于10个），具有一定活性梯度且不存在于训练集中的活性分子组成测试集

• 使用测试集分子预测模型的筛选能力，q^2越接近于1，预测能力越强

Decoy Set 诱饵分子数据集验证：是指通过已知化合物与非活性化合物建立decoy set，该数据集可验证该药效团模型对于活性化合物和非活性化合物的区分能力。模型筛选性能参数：命中率(hit rate)、富集因子(Enrichment Factor，EF)与ROC曲线(Receiver Operating Characteristic Curve， ROC Curve)

Decoy有一个现成的网站，可直接下载测试文件 http://dude.docking.org/

药效团模型方法vs定量构效关系方法

药效团模型可以基于不同类的化合物，它得到的是与生物活性有关的重要的药效团特征，这组药效团特征是对配体小分子活性特征的抽象与简化，因此药效团模型方法可以用于寻找结构全新的先导化合物。

定量构效关系方法一般是基于母体相似的同系列化合物，它得到的是化合物结构和活性之间的定量关系，这个定量关系一般只能用于先导化合物的改造。

总结

药效团模型的构建和验证实验

实验目的：

1. 掌握基于分子共同特征的药效团模型（HipHop）的构建和验证。

2. 掌握基于复合物的药效团（CBP）模型的构建和结果分析。

实验原理：

使用 Discovery Studio 软件进行基于分子共同特征的药效团模型（HipHop）和基于复合物的药效团（CBP）模型的构建，验证和结果分析。

HipHop：基于分子共同特征的药效团模型。用于发现一系列配体小分子所共有的化学特征，并基于这些共同特性结构的比对叠合自动生成药效团模型，用户可以使用共有的特征药效团去搜索化合物数据库来寻找可能的先导分子。

CBP：基于 LigandScout 算法，通过识别受体‐配体相互作用的关键特性构建药效团模型进行先导物的发现与优化。该类药效团从受体配体间的相互作用出发，能更加准确直接地反应出它们之间的药效团特征。

本实验所用软件环境：

• DS Version：19.1.0.18287

• PP Version：19.1.0.1963

• DS Client Version：19.1.0.18287

• OS Distribution：Windows

• OS Version：10.0.22000

实验步骤：

◆ 基于分子共同特征的药效团模型的构建和验证

1. 训练集分子的准备：本实验中使用老师提供的 5HT2c_ligands.sd 训练集。

2. 药效团特征元素的选取：

   药效特征元素的定义：点击 Discovery Studio 软件上的 Pharmacophores → Edit and Cluster Features → Feature Mapping 进行药效特征元素的定义。设置参数如下图。

   

   药效特征元素的查看：点击 Discovery Studio 软件上的 Pharmacophores → Edit and Cluster Features → Current Feature → All Features 进行药效特征元素的查看。

3. Common Feature Pharmacophore 的构建：点击 Discovery Studio 软件上的Pharmacophores → Create Pharmacophores Automatically → Common Feature Pharmacophore Generation 进行 Common Feature Pharmacophore 的构建。设置参数如下图，使用老师提供的 Decoy 数据集 active_ligands 和 Decoy 数据集 inactive_ligands 进行药效团验证。

4. 外部测试集分析：合并老师提供的 Decoy 数据集 active_ligands 与 Decoy 数据集 inactive_ligands（也就是复制 Decoy 数据集 active_ligands 中的全部分子，然后粘贴到 Decoy 数据集 inactive_ligands 中，或者直接 下载已经合并好的文件 ），点击 Discovery Studio 软件上的 Pharmacophores → Search, Screen and Profile → Ligand Profile 进行外部测试集的分析，设置参数如下。Input Ligands 选择前面合并的数据集，Input File Pharmacophores 选择之前运行 Common Feature Pharmacophore 构建流程所得到的 Output 文件夹中的所有模型。（按 SHIFT-选取）

◆ 基于复合物的药效团模型的构建和结果分析

1. 蛋白的准备：本实验中使用老师提供的 2irz.pdb 蛋白 。点击 Discovery Studio 软件上菜单栏上的 Structure→ Crystal Cell→ Remove Cell 进行晶胞的去除。点击Macromolecules → Prepare Protein → Prepare Protein，设置参数如下，来进行蛋白结构的准备得到新窗口 2IRZ_prep 。接下来的操作都是在新的窗口当中进行。

2. 配体的准备：从 2IRZ_prep 窗口中剪切 2IRZ 自带的配体（名称为Ligand1，在Ligand Groups中）粘贴到新窗口中，并重命名为 Ligand。（重命名的操作：点击选中 2irz 并单击鼠标右键，选择最后一项 Attribute of 2irz…，出现下图对话框，将 Name 改为 Ligand）

3. 药效团模型的构建：点击 Discovery Studio 软件上的 Pharmacophores → Create Pharmacophores Automatically → Receptor-Ligand Pharmacophore Generation 进行药效团模型的构建，设置参数如下。再把配体小分子再复制粘贴到结果窗口，观察药效团与配体小分子的匹配情况。

实验结果：

基于分子共同特征的药效团模型的构建和验证：

1. 药效特征元素的定义的结果

2. Common Feature Pharmacophore 的构建的结果

3. 外部测试集分析的结果

基于复合物的药效团模型的构建和结果分析：

1. 蛋白的准备

2. 药效团模型的构建

讨论：

参与基于分子共同特征的药效团模型的构建的六个分子：

这些分子结构具有多样性，包含principal和MaxOmitFeat性质参数，都有强活性。

从 HipHop 模型验证结果可以看出，HipHop 模型的敏感度很好，但是特异性很差，用于药物筛选时，可能会错误的把无活性的化合物认为潜在的Hit，但是应该不会错过潜在的Hit。其中名为 Pharmacophore_1 的药效团在验证中的表现最好，因为 AUC 值最大。

虚拟筛选 (Virtual Screening, VS)

虚拟筛选在药物发现中的意义 The significance of virtual screening in drug discovery

筛选途径的药物发现方法

实验筛选：是使用生物物理、生物化学等方法鉴别小分子和蛋白质的相互作用，从而筛选出对蛋白质功能有调节作用的先导化合物。基于大规模化合物库的实验筛选也称为高通量筛选（HTS）。

虚拟筛选：是指在进行生物活性筛选之前，根据预先设定的条件，在计算机上对化合物分子进行预筛选，以识别出最可能与靶标结合的小分子，从而大大降低实际筛选化合物数目，同时提高先导化合物的发现效率。

高通量筛选VS虚拟筛选

高通量筛选	虚拟筛选
体外活性测试（大量测试）	计算机模拟分析
命中率：0.01‐0.001% 假阴性较高	命中率：2‐24%
需要先有大量化合物	仅需测试少量化合物

虚拟筛选文章统计

虚拟筛选的基本流程 The basic process of virtual screening

虚拟筛选的通用流程

虚拟筛选化合物数据库

常用数据库：

ZINC、PubChem、DrugBank、ChEMBL、ChemDB、HMDB、BindingDB、SMPDB。

常用化合物库介绍：

1） 筛选化合物库分类：多样性化合物库、上市药物分子库、已知活性库、靶点化合物库，天然产物库，片段库等等；

2） 常用化合物库品牌：Chemdiv、Enamine、Lifechemicals、Specs、Chembridge、Maybridge、Microsource、Vitas‐M、Interbioscreen等

3） 各品牌化合物的优劣势介绍：

1. ChemDiv：140万库存分子，多样性被认为是商业库中最好的，每年新增20万以上的分子，每年新增的量就和Specs全部的分子数一样多了，而且有很多虚拟筛选的针对各靶点的子库，被大量的高通量筛选的客户选用。

2. Enamine：100万以上，乌克兰的公司，多样性也不错，但价格要高于chemdiv；

3. Lifechemicals：60多万种分子，库存量充足，如部分产品库存不足可很快完成合成及补充工作，进入中国市场较晚，价格中等；

4. SPECS：20多万种分子，价格最便宜，最早进入中国市场，已经被各大研究单位无数次的反复筛选过。化合物更新较慢，溶解性，多样性都相对比较差；

5. Chembridge：有70万种小分子，将化合物分成1，2，3 或4类，价格都不一样，是依赖于化合物合成的难易来分的。相对价格较贵

化合物数据库的预处理

1、类药性分析

目的：排除理化参数不佳的化合物, 降低筛选成本。保留具有一定类药性(drug‐likeness) 特点的化合物

参数	最小值	最大值
logP	‐2	5
分子量	200	500
氢键给体	0	5
氢键受体	0	10
摩尔折射率	40	130
旋转键数目	0	10
重原子数目	20	70
极性表面积（tPSA）	0 A²	120 A²
电荷	‐2	2

2、假阳性分子排除

目的：一些化合物易于与生物大分子发生化学反应, 在基于受体、酶或细胞检测实验中总是表现为阳性, 而实际上为假阳性。

虚拟筛选的不同方法

基于分子对接的虚拟筛选步骤

参考：http://dx.doi.org/10.1038/nchembio.2334

作者通过同源模建构建了MRGPRX2蛋白结构，在此基础上通过分子对接成功找到了一个选择性的激动剂，并得到实验验证，可作为药物开发的先导化合物。

基于片段设计的虚拟筛选步骤

基于片段的药物发现方法的主要内容是设计并建立由片段分子组成的化合物库，并对库中的片段进行生物活性的筛选从而找到苗头片段分子。再利用X射线晶体学、核磁共振等技术对这些片段与靶蛋白的结合模式与结合强度进行分析,根据这些结构信息对片段分子进行结构优化得到先导化合物。

研究案例

参考：doi:10.1021/jm501982k

凝血因子VIIa是一种新的抗血栓形成靶标，与S1口结合的大多数因子VIIa抑制剂包含疏水和阳离子基团，导致不好的膜渗透性和口服吸收。为了发现中性S1靶向抑制剂并改善药代动力学特征，Cheney等研究人员进行了基于片段的虚拟筛选研究。

代表性命中片段（K_i = 8.9mM）与凝血因子VIIa S1口袋结合并与氨基酸残基Gly218形成氢键。基于命中片段内酰胺的结合模式，发现了许多具有良好抑制效能和改善渗透性的类似物。如化合物5是细胞可渗透的，并且具有130 nM的K_i值。

基于药效团模型的虚拟筛选步骤

参考：doi:10.1021/ci500762z

研究案例：人组蛋白去乙酰化酶（HDACs）作为表观遗传调控的关键因子，已被确定为治疗多种癌症的药物靶标。总结现有的HDAC抑制剂发现，都存在一个锌结合位点（zinc‐binding group，ZBG）。

药效团的构建

模型验证：构建诱饵分子库：包含134个已报道的HDAC8抑制剂和10000个无活性的诱饵分子

筛选策略

基于相似性搜索的虚拟筛选

首先需要将已知的活性化合物和数据库中化合物进行相似性计算打分，得分越高的分子相似性越大；然后根据设定的相似性阈值，将满足相似性要求的分子从数据库中挑选出来。

研究案例1

参考：DOI: 10.1016/j.ejmech.2012.04.041

雌激素受体（ERs）是ER阳性乳腺癌和骨质疏松症的有效靶标。亚型选择性配体的数量的增加有助于ERa和ERb之间的差异的研究。

研究方法：

1.基于分子相似性使用参考化合物对SPECS数据库进行相似性搜索，相似系数设定为0.45‐0.99；所有命中化合物进行基于对接的虚拟筛选，最终获得具有合适结合构象的化合物；

2.购买纯度合适的命中化合物；

3.运用不同生物技术对化合物进行活性评估。

通过相似搜索和分子对接的组合策略的虚拟筛选，发现了21种化合物作为雌激素受体的有效配体，其中两个具有高的抑制率（> 70％）。

研究案例2

参考：https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4871132/

使用夫西地酸作为提问结构进行基于2D和3D相似性的虚拟筛选得到新型抗疟疾药。

不同虚拟筛选方法的联合应用

SBDD以锁匙模型为基础，寻找可以与靶标结合的配体分子；而SBDD以相似性原理为基础，寻找与已知抑制剂具有潜在相同作用方式的配体分子。SBDD结果准确，可直接预测结合能，缺点是需要受体结构，速度慢；LBDD速度快，无需蛋白结构，但是有时由于缺乏具体结合方式的指导，结果外推时常常误差很大。这两种方法、理论是独立发展，各有其适用条件和利弊的。所以一个虚拟筛选项目通常会整合多种CADD技术。

案例学习 Case studies

药效团模型为主的虚拟筛选方法

案例1. HIV‐1整合酶抑制剂的发现

参考：https://doi.org/10.1021/jm050837a

训练集：

药效团模型：

2个氢键给体、1个氢键受体、1个芳环中心

数据库搜索：

案例2. HIV‐1整合酶抑制剂的改造

参考：https://doi.org/10.1021/jm050549e

训练集：

药效团模型（具有活性预测能力的药效团）：

2个氢键受体，1个氢键给体，1个芳环中心

活性预测：

药效团活性预测结果与最终的实验活性测定结果进行比较，结果非常一致，线性相关常数r达0.85。

化合物改造

16个化合物中，有一个化合物与药效团模型的匹配方式如图A所示，该化合物只能与药效团中的三个特征相匹配，芳香疏水中心没有匹配上。依据这一信息，研究者在该化合物上添加了一个疏水的芳环侧链，然后进行了活性测试。发现改造之后的化合物（如图B）所示，与之前相比，活性有了很大的提高。

案例3. 基于CBP模型发现新型Lyp抑制剂

参考：https://doi.org/10.1021/jm500692u

药效团模型构建案例：

对接及其打分函数评价（诱饵分子库包括52个活性分子和2950个无活性的诱饵分子）

虚拟筛选方案：

量子力学 (Quantum Mechanics, QM)

量子力学基础 Quantum Mechanics Foundation

1900年Planck提出量子概念，是表现某物质或物理量特性的最小基本单位。

1913年Bohr提出：原子中的电子只能处于包含基态在内的定态上，电子在两个定态之间跃迁而改变它的能量，同时辐射出一定波长的光，光的波长取决于定态之间的能量差——量子论

量子力学：研究微观粒子运动规律的理论。用波函数描写粒子的运动状态，以薛定谔方程确定波函数的变化规律，并对各物理量进行计算。

量子化学：以量子力学的基本原理和方法来研究化学问题的学科。可计算分子的各种参数，如分子结构、成键特征、系统总能量、各个轨道的分子信息、分子间相互作用、各种光谱、振动谱、电磁谱性质和反应过渡态等

量子力学用波函数描写粒子的运动状态，以薛定谔方程(Schrödinger Equation) 确定波函数的变化规律，并对各物理量进行计算。

求解方程主要方法：

• 从头算（Ab initio）

• 密度泛函方法(DFT)

• 半经验计算法(Semi‐imperical)

从头计算法（ab initio methods）：

• 以基本物理常数以及元素的原子序数，不借助于任何经验参数，求解薛定谔方程(Schrödinger Equation)：Ĥψ(r) = E ψ(r)

• 计算结果精度高，可靠性大，但是计算量极大，消耗计算机时太多

• 从头计算方法包括Hartree‐Fork、微扰、耦合簇、组态相互作用等

密度泛函方法（Density Function Theory, DFT）：

• 用电子密度取代波函数做为研究的基本量，其将电子能量分为几个部分，动能、电子‐核相互作用、库仑排斥，以及其余部分的交换相关项，所有项只是电子密度的函数。

• DFT方法之所以计算效率更高，是因为它使用一个单变量积分来代替双电子的Coulomb积分以及用一个交换‐相关泛函(VXC) 来代替双电子的电子交换‐耦合积分。

1994年，最重要的一种交换‐相关泛函B3LYP被提出后，DFT就日趋流行，B3LYP几乎成为了计算各种问题的默认方法。在此之后比较重要的发展有以下这些：

• 以DFT‐D为代表的色散校正方法：解决B3LYP、PBE等传统泛函难以或完全不能描述色散作用的问题。

• 双杂化泛函：将MP2的二阶校正引入交换‐相关泛函。精读比普通泛函高一个级别，但计算量也高一个级别。

半经验计算法（semiempirical methods）

采用实验值拟合的经验参数，使薛定谔方程得到简化后再求解，提高计算速度

• 但计算精度较差，造成结构差异较大的体系进行性质比较时可靠性不高

• 可用于处理较大分子

• 半经验计算法的软件有：MOPAC和AMPAC等

量子化学计算的优缺点

优点：可计算出分子的理化参数、分子结构的电子结构，几何构型和易与亲电试剂或亲核试剂反应的部位

缺点：只适用于计算分子量较小分子，计算时间长

量子化学计算的常规流程：

常用的量子化学计算软件:

以计算原理分：

基于从头算或第一性原理方法：Gaussian、ORCA、PSI4、ADF、Dalton、Gamess、Crystal、VASP、Jaguar、Dmol3等

基于半经验或分子力学方法：MOPAC、EHMO、Hyperchem等

以研究对象分：

有限尺度体系(分子、簇合物等)：Gaussian、ORCA、PSI4、Jaguar、ADF、Dalton、Gamess、MOPAC、EHMO等

无限周期重复体系(晶体、固体表面、链状聚合物等)：Crystal、NNEW3、VASP、Wien、Dmol3等

“量子化学已经发展成为广大化学家使用的工具，将化学带入了一个新的时代。在这个新时代里，实验和理论能够共同协力探讨分子体系的性质。化学因此不再是单纯的实验科学了。” ——1998年瑞典皇家科学院诺贝尔化学奖发布公告

量子化学计算在药物化学研究中的应用 The Application of Quantum Chemistry in Medicinal Chemistry

1. 量子化学参数与定量构效关系

量子化学描述符可以成为定量构效关系的描述符，比如：总能量(Et)、电离能(Ee)、分子生成热(Ef)、前线轨道能量HOMO和LUMO等能量参数，酯水分配系数(LogP)、偶极距(Dipo)、分子极化率(a)、摩尔折射率(MR)、分子量(w)、分子范德华体积(y)和一些电性参数等

2. 分子力场参数的来源之一

3. 分子中各原子的部分电荷（partial charge）：主要用于分子模拟中，描述静电相互作用

4. 分子静电势（Electrostatic potential, ESP）

物理定义：空间某点的静电势是指从无穷远处移动单位正电荷至该点时所需做的功，换句话说，静电势实质上是静电相互作用力存在的根源之一。简单地讲，如果不考虑其他因素，当分子周围空间某一点静电势为正值时，意味着这一点对正电荷是排斥的；相反，如果某一点静电势为负值时，则意味着这一点对正电荷是吸引的。由于分子间的静电相互作用力是分子间主要远程相互作用，因此，静电势对于考察分子间静电相互作用、预测反应位点、预测分子性质等方面有重要意义，被广为使用。

表面静电势：不同表面区域静电势大小通过不同颜色展现，使分子表面上静电势的分布一目了然。分子表面一般都用Bader定义的范德华表面，即电子密度为0.001 e/Bohr^3的等值面。

静电势为负值——即带正电荷的微粒与其有较强的相互作用，容易与之靠近；相反，带负电荷的微粒则容易与H靠拢。据此，可对分子反应部位、分子间相互作用和分子识别做出初步判断。

5. 弱相互作用分析

参考：doi:10.1002/jcc.26068

为了揭示氢键本质、讨论如何估计氢键强度，此文首先构造了一批氢键二聚体

对上述42个氢键二聚体分别给出了结合能当中的静电作用(Elst)、色散作用(Dis)、诱导(Ind)、交换互斥(Rep)部分各自的贡献，其中前三项对于氢键的形成起到正贡献

6. 分子扭转能/内能（Strain Energy）

小分子构象:

• 低能构象/优势构象

• 活性构象/药效构象

配体分子与靶点结合的构象变化：从溶剂非结合构象系综变为到生物活性构象。

参考：doi:10.1021/acs.jmedchem.9b01720

量子化学计算在药物合成中的应用 The Application of Quantum Chemistry in Drug Synthesis

1. 前线轨道理论的应用

在亲电反应中，亲电试剂会进攻底物上HOMO值最大的地方，即电子能量最活泼的原子；相反，在亲核反应中，亲核试剂会进攻底物上LUMO值最小的，即可以接受电子的原子。

在前线轨道理论（Frontier Molecular Orbital）中，化学反应的发生，是因为两个反应物之间的最高占有轨道HOMO (Highest Occupied Molecular Orbital)和最低未占轨道LUMO (Lowest Unoccupied Molecular Orbital) 相互吸引，进而发生电子转移，导致反应发生。

预测卤代反应的位点

预测亲核取代反应的位点

立体选择性

2. 静电势的应用

质子酸性预测

静电势能图是将分子各个位点的静电势能映射到分子形状表面

3. 混合的量子力学/分子力学（QM/MM）方法

分子力学（MM）：计算精度不高、无法考虑化学反应

量子力学（QM）：计算量大、能处理体系小

量子力学/分子力学（QM/MM）联用应用：复合物结构优化，酶催化反应研究等。

如何给此项目贡献

第一步：安装Github并克隆此项目

1.从 https://desktop.github.com/ 下载安装好github desktop。

2.打开GitHub desktop。如果没有账号请点击Create your free account，如下图所示。

会弹出一个网站，填写网站所要求的信息,如下图所示。

输入邮箱收到的验证码，如下图所示。

后面的调查问卷可以不管。点击continue for free完成GitHub账号的注册。

4.注册完账户后，打开安装好的github desktop软件，点击sign in github。 会弹出一个网站，点击authorize desktop。然后会自动回到github desktop软件，点击finish。

5.接下来，点击Clone a repository from the Internet… 在弹出的框里选择URL。Repository URL中填写 https://github.com/abdusemiabduweli/CADD-tutorial Local path选择一个你想要保存此项目的文件夹，点击Clone，如下图所示。

就这样在电脑上安装配置好了Github，也下载好了一份此项目，项目所在文件夹地址就是你刚才在Local path输入的地址。

第二步：安装python

1.从 https://www.python.org/downloads/ 下载安装好python。

第三步：安装Visual Studio Code并搭建好环境

1.从https://code.visualstudio.com/ 下载安装好visual studio code。

2.打开GitHub Desktop，在如下图1所示的地方选择此项目 CADD-tutorial， 并点击2所示的按钮Open in visual studio code，如下图所示，打开Visual studio code。

如果出现如下图的提示，打勾1所示的框，点击2所示的按钮。

3.在visual studio code中点击1所示的extensions按钮，搜索ext:rst，点击2和3所示的插件的install按钮，安装这两个插件。

4.同样的，点击extensions按钮，搜索python，安装第一个插件。安装完后会出现如下图所示的界面，点击1所示的按钮，再点击选择2所示的文件。

5.点击Visual studio code菜单栏上的Terminal–>New Terminal，下图中的1所示的按钮，在下方窗口中输入pip install sphinx sphinx-autobuild，点击Enter键，等待安装完成。一般再次出现PS D:等字样时，表明安装完成。同样的操作，输入pip install furo，等待安装完成。

第五步：如何给此项目贡献

如果完成了上面的操作，那么你的电脑上就保存着一份此项目。项目就在上面所示的Local path中的地址。如果打开文件夹，你能看到如下图一样的一堆文件和文件夹。

其中红色方框中显示的就是本网站（ https://cadd-tutorial.readthedocs.io） 的相关文件。改变我的Github账号上的CADD-tutorial项目中的这些文件，就等于改变本网站。浅绿色方框中的文件夹你可能看不到，这是因为我在深绿色方框中的文件（.gitignore）里列举了这个文件夹（用记事本打开此文件就能看到其内容）。.gitignore文件的详细介绍可以参考这个文章：https://www.jianshu.com/p/699ed86028c2 。剩余的文件都是我提供给你们的参考文献。这些文献大部分来自中国药科大学计算机辅助药物设计课程。因为我没有得到此课程负责老师们的同意（我没问，我怕他们不同意我上传这些文件到Github），所以我们的第一个目标是把这些文件都转换成rst文件。rst文件都在此项目docs文件夹里，后缀为.rst，这些rst文件都是遵循rst文件的语法规则，rst语法规则可以参考这个文章 https://www.sphinx-doc.org/zh_CN/master/usage/restructuredtext/basics.html#source-encoding。这些rst文件编辑完后上传到我的github账号中此项目时，会被自动 https://readthedocs.org/ 网站转换成html等文件保存在 readthedocs团队提供的服务器中。我们的此项目成果 https://cadd-tutorial.readthedocs.io （CADD教程网站）就是由它们免费提供的技术支持实现的，这里特地感谢他们。readthedocs不仅把rst文件转换成html文件，还会转换成PDF，Epub文件，这些文件可以点击下图所示的按钮下载。

我提供的文献都是PDF格式的，除非在电脑上阅读，但是其他设备上阅读起来非常差，也不能编辑。如果我们把我们的文献全部转换成了rst文件，我就能把这些pdf文件都删掉，这样也不会有任何版权问题。我们的此项目遵循GNU协议，包括我们的rst文件，详情请阅读项目中的LICENSE文件。

那么怎么编辑我们的项目成果 https://cadd-tutorial.readthedocs.io （CADD教程网站）？ 前面提到了只要我们编辑此项目中docs文件夹中的文件和.readthedocs.yaml并上传到我的GitHub账号中的CADD-tutorial项目，CADD教程网站就会自动地被改变。.readthedocs.yaml文件是CADD教程网站的配置文件，可以先不管，若感兴趣，详情请阅读这个文章 https://docs.readthedocs.io/en/stable/config-file/index.html 。打开docs文件夹，你能看到以下文件

其中images保存此网站中用到的图片。Conf.py和requirements也是此网站的配置文件，也可以先不管，详情请阅读这些文章 https://www.sphinx-doc.org/en/master/usage/configuration.html https://pip.pypa.io/en/latest/user_guide/#requirements-files。其他的rst文件就是CADD-tutorial网站中各网页的内容，需要用rst语法规则编辑。rst语法规则可以参考这个文章 https://www.sphinx-doc.org/zh_CN/master/usage/restructuredtext/basics.html#source-encoding。打开我们的此项目成果-CADD教程网站 https://cadd-tutorial.readthedocs.io 就能看到首页，首页对应的rst文件是index.rst。在左侧能看到分子力学、分子动力学、定量构效关系、人工智能、基于结构的药物设计、基于配体的药物设计、如何给此项目贡献等网页。这些网页所对应的rst文件分别是 MolecularMechanics.rst MolecularDynamics.rst QSAR.rst AI.rst SBDD.rst LBDD.rst contribute.rst。编辑相应的rst文件就相当于编辑了相对应的网站。最后只需上传到我的github账号中CADD-tutorial项目中。

下面我给你们演示编辑“如何给此项目贡献”网页。 打开Github Desktop，选择CADD-tutorial，如下图所示

点击open in visual studio code打开visual studio code 点击docs,双击打开contribute.rst，如下图所示

编辑完此文件后要保存。如果想要看网页效果的话，点击右上角的如下图所示的按钮就能显示网页效果。

最终打开github desktop就能发现跟如下图相似的界面。

点击1所示的commit to master，再点击如下图所示的push origin

如果你账号没有此项目的fork的话，就会出现如下图所示的提示，需要点击fork this repository。这操作意味着你在你的账号中复制粘贴了我github中此项目的文件

如果你想要给此项目贡献，那你需要选择如下图中的第一个选项to continue to the parent project并点击continue。如果你想要独立建立此项目的克隆项目，那你选择第二个选项for my own purpuses并点击continue。 当然这里我希望你选择第一个。

然后你需要重新点击push origin。 然后，打开浏览器，打开github.com，登录你的账号，你账号里就能看到此项目的克隆项目，如下图所示

点击打开此项目，点击contribute，再点击open pull requesat，如下图所示。

最后点击create pull request，如下图所示，

按照如下图所示的填写表格 点击create pull request，你就能对此项目做出贡献了。当然你做出的贡献需要我的审核，等你擅长使用github等软件技术后，我就给你权限，这时你不需要等待我的审核就能给此项目做出贡献。

如果想要添加新的网页或者删除网页，只需添加或删掉对应的rst文件并在index.rst文件中.. toctree::的下面添加或删除相应rst文件名称。

非常感谢大家对此项目的支持，如果有任何问题可以联系我。我推荐使用github的isuues来提出问题，大家也可以帮你解决问题，我也会用我的第一时间去回答你的问题。

如果想要进一步学习相关技术软件，请阅读以下文章，教程视频。 关于git： 关于vs code: 关于sphinx：https://www.sphinx-doc.org/en/master/index.html。 关于rst：