CN108875298B - 基于分子形状匹配的药物筛选方法 - Google Patents

Abstract

本发明提供了一种基于分子形状匹配的药物筛选方法，其包括：根据目标形状，利用三维Zernike函数计算目标形状的展开系数C_nlm标和形状描述因子F_nl标；利用三维Zernike函数，计算化合物库中每个化合物的展开系数C_nlm分子和形状描述因子F_nl分子；通过比对每个化合物的F_nl分子与F_nl标之间的差异性，对化合物库进行初步筛选；以及再将初次筛选后得到的备选化合物的C_nlm分子与C_nlm标进行对比，对备选化合物进行二次筛选。该筛选方法，通过分步衡量化合物库中每个化合物与目标形状的相似性，极大的降低了筛选过程中的计算量，大幅缩减了备选分子的数目，筛选效率高、且增加药物筛选的成功率。

Description

基于分子形状匹配的药物筛选方法

技术领域

本发明涉及药物领域，具体而言，涉及一种基于分子形状匹配的药物筛选方法。

背景技术

药物筛选是现代药物开发过程中获取具有特定生理活性化合物的一个重要步骤，其是指从大量化合物库中选择对某一特定作用靶点具有较高活性的化合物的过程。随着药物开发技术的发展，药物研发已经从早期的对系列化合物的生理活性进行验证性实验，逐渐转变为对化合物库进行大规模筛选，极大地缩短了药物研发的进程。

目前,普遍使用的筛选方法是基于计算机辅助的活性位点分析法和数据库搜寻。该类方法通常需要先分析靶标生物大分子的活性位点、建立一系列活性分子的药效构象，抽提出共有的药效基团，进而在现有的化合物数据库中寻找符合药效基团模型的化合物。然而，由于化合物数据库中化合物数量繁多，通过这种基于药效团的筛选，得到符合标准的化合物依旧数量庞大，难以较为准确的筛选到先导化合物，筛选效率低。

发明内容

本发明的目的在于提供一种基于分子形状匹配的药物筛选方法，这种药物筛选方法的计算量小、筛选效率高，筛选后的备选分子数目大幅减小，且备选分子有效性高。

为了实现本发明的上述目的，特采用以下技术方案：

一种基于分子形状匹配的药物筛选方法，其包括：

根据与靶标生物大分子的结合部位相符的目标形状，利用三维Zernike函数计算所述目标形状的展开系数Cnlm标和形状描述因子Fnl标；

利用三维Zernike函数，计算化合物库中每个化合物的展开系数Cnlm分子和形状描述因子Fnl分子；

通过比对每个化合物的形状描述因子Fnl分子与所述目标形状的形状描述因子Fnl标之间的差异性，对所述化合物库进行初步筛选得到备选化合物；以及

再将所述备选化合物的展开系数Cnlm分子与所述目标形状的展开系数Cnlm标进行对比，对所述备选化合物进行二次筛选。

优选的，所述化合物的形状描述因子Fnl分子与所述目标形状的形状描述因子Fnl标之间的差异性是通过所述化合物的形状描述因子Fnl分子与所述目标形状的形状描述因子Fnl标之间差值的二阶矩来进行评价。

优选的，所述化合物的形状描述因子Fnl分子与所述目标形状的形状描述因子Fnl标之间的差异性采用由式(3)计算的卡方评分X2来表示：

式(3)中，n代表Zernike展开的主阶数，从0到设定值nmax,l代表依赖于主阶数n的副阶数，从0到n,且需要满足(n-l)是偶数。

优选的，所述初步筛选包括：将所述化合物库中每个化合物按照卡方评分X2从小到大的顺序进行排序，取排名靠前的10～20％的化合物作为备选化合物；或者根据卡方分布，选取分数低于平均值一个标准差的化合物作为备选化合物。

优选的，所述备选化合物的展开系数Cnlm分子与所述目标形状的展开系数Cnlm标之间的相似度通过皮尔森相关系数cc进行评价，取皮尔森相关系数cc超过预设阈值的化合物作为备选分子。

优选的，所述皮尔森相关系数cc采用式(4)计算：

式中，<ρ(r)>为分子或目标形状的密度函数；<ρ₁(r)ρ₂(r)>为分子与目标形状的重叠程度；σ(ρ₁(r))σ(ρ2(r))为分子密度的涨落程度；r为空间向量。

优选的，式(4)中所述密度函数ρ(r)为近似密度函数所述近似密度函数采用式(5)计算：

式中，Cnlm为展开系数Cnlm分子或Cnlm标，Znlm是化合物或目标形状的Zernike函数；最高阶数nmax为10～30。

优选的，还包括：将经过二次筛选后得到的备选分子与所述靶标生物大分子进行结构对接，并采用分子动力学的方法进行优化以对所述备选分子进行三次筛选，得到用于实验检测的目标分子。

优选的，在对所述化合物库进行初步筛选之前，还包括：根据靶标生物大分子的药效团对所述化合物库进行预筛选，取符合所述药效团的化合物集合进入所述初步筛选。

优选的，所述目标形状是通过模拟所述靶标生物大分子的活性部位来定义的三维空间形状；或者，所述目标形状是通过抽提能够与所述靶标生物大分子结合的活性分子的晶体结构来确定的三维空间形状。

与现有技术相比，本发明的有益效果如下：

本发明提供的这种药物筛选方法，是基于目标形状来实现的，其通过利用正交的三维Zernike函数，得出每个化合物或目标形状在三维Zernike函数空间的展开系数{C_nlm}作为该形状的数值描述，并进一步得出旋转不变的形状描述因子{F_nl}，分步衡量化合物库中每个化合物与目标形状的相似性，对化合物库进行筛选，极大的降低了筛选过程中的计算量，大幅缩减了备选分子的数目，筛选效率高、且增加药物筛选的成功率。

同时，该药物筛选方法中的目标形状与靶标生物大分子的结合部位的形状相符。基于这种目标形状筛选到的备选分子与靶标生物大分子的匹配程度高，提示其有很大可能性进入靶标生物大分子的结合部位并与之稳定结合产生相互作用，因此筛选到的备选分子的有效性高。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，以下将对实施例或现有技术描述中所需要使用的附图作简单地介绍。

图1为用三维Zernike函数近似得到的ZINC化合物数据库中标号为ZINC10383266的小分子三维形状；(a)为小分子的化学分子式，(b)为小分子的原子结构(用化学键的形式表示)和用最高阶数为20描述的该分子的外形(白色透明的外表面)。

图2为本发明实施例1提供的筛选方法的流程图。

图3为本发明实施例1提供的大麻素蛋白受体的药效团分布(a)以及环绕抑制剂分子的二维描述(b)。

图4为本发明实施例1中基于大麻素蛋白受体构建的目标形状。

具体实施方式

下面将结合实施例对本发明的实施方案进行详细描述，但是本领域技术人员将会理解，下列实施例仅用于说明本发明，而不应视为限制本发明的范围。实施例中未注明具体条件者，按照常规条件或制造商建议的条件进行。所用试剂或仪器未注明生产厂商者，均为可以通过市售购买获得的常规产品。

药物分子的形状和表面电荷分布是其与靶标生物大分子(简称受体分子)相互作用的基础。药物分子的形状与受体分子的结合位点的匹配程度直接影响结合的亲和力和稳定性：药物分子的尺寸和形状必须满足特定的要求才能进入到受体分子的结合部位并与之产生相互作用、实现其功能。药物分子尺寸太大，会导致其与受体分子之间产生不必要的碰撞，甚至不能进入到受体分子真正的结合位置；如果药物分子尺寸过小，则会导致其不能与受体分子的全部功能位点进行有效结合。药物分子的形状对其与受体分子的结合也会有类似的影响。

鉴于此，本发明提供一种基于分子形状匹配的药物筛选方法，其在药物筛选的前期阶段，根据受体分子的结构信息确定药物分子的外形特征(即得到目标形状)，再与数据库里的小分子化合物进行形状匹配，筛选潜在药物，从而增加药物筛选成功率。

该药物筛选方法包括以下步骤：

步骤S1:根据与靶标生物大分子的结合部位相符的目标形状，利用三维Zernike函数计算目标形状的展开系数和形状描述因子

步骤S2:利用三维Zernike函数，计算化合物库中每个化合物的展开系数和形状描述因子

步骤S3：通过比对每个化合物的形状描述因子与形状描述因子 之间的差异性，对化合物库进行初步筛选；以及

步骤S4：再将初次筛选后得到的备选化合物的展开系数与展开系数进行对比，计算出备选化合物与目标性状的密度相关系数，并以此对备选化合物进行二次筛选。

其中，在步骤S1中的目标形状须与靶标生物大分子的结合部位相匹配，即与该目标形状相似度高的药物分子能够进入该靶标生物大分子的结合部位并与之结合。

进一步地，该目标形状是通过模拟靶标生物大分子的活性部位来定义的三维空间形状，即利用计算机辅助模拟并探测靶标生物大分子的活性位点，由活性位点分析得到有关靶标生物大分子结合的信息，并由此定义出结合部位(即目标形状)。目前，能够实现活性位点分析的软件有DRID、GREEN、HSITE等。

或者，该目标形状是通过抽提能够与靶标生物大分子结合的活性分子的晶体结构来确定的三维空间形状。即，通过分析靶标生物大分子的一系列活性分子的空间结构的共有特征，确定出能够与靶标生物大分子结合的药物分子的外形特征(即目标形状)。若有多个活性分子的结构，则可以分别作为目标形状，也可以通过平均操作形成一个更加平滑的目标形状。单一目标性状可以筛选出更有针对性的化合物，但是有可能被目标性状局限。平均多个活性分子结构之后形成的平滑目标形状会降低针对性，但是会增加备选分子的多样性。若计算资源允许，建议进行组合筛选，得出备选分子的交集(或者合集)进行实验检测。

在确定目标形状后，发明人利用三维Zernike函数对化合物库中的小分子与目标形状进行形状识别：

三维Zernike函数，由两部分组成，在球坐标的三维空间中，球表面函数用正交归一的球谐函数(spherical harmonics)标识(为了清晰，阶数l,m未标注)；在径向用特殊设计的在[0,1]区间正交归一函数R(r)标示(阶数为n,l)。这两个函数的乘积就组成了三维Zernike函数最终得出的阶数由(n,l,m)三个参数确定。简言之，三维Zernike系列函数在单位球体内是正交归一的，也就是说，可以把任意在单位球体内定义的密度函数ρ(r)在Zernike函数上做投影展开，这些展开系数{C_nlm}是在Zernike空间内描述密度函数ρ(r)的完备且唯一的参数集合。

即，在三维空间内的任一三维物体皆有与之对应的一组Zernike函数的展开系数(也被称为Zernike矩)，将该三维物体的密度函数ρ(r)缩放到单位球体内，并在Zernike函数上做投影展开，即可得到能够在Zernike空间内完备且唯一地描述该物体的展开系数{C_nlm}。因此，在该步骤中，发明人利用三维Zernike函数计算得到目标形状的展开系数以及化合物库中待筛选的化合物的展开系数用以后续的对比、筛选。

其中，展开系数和按下式计算：

式中，n、l、m为阶数，|r|为径向距离，黑体r标示三维空间任意位置，Z^*是Zernike函数的复数共轭，ρ(r)为三维空间密度函数。

公式(1)通过积分可以得到密度函数ρ(r)在阶数为(n,l,m)的Zernike函数的展开系数C_nlm(也被称为Zernike矩)，由其对目标形状以及待选化合物形状进行描述。一系列的展开系数构成在Zernike函数空间内描述原来密度函数的另一种表达方式。

目标形状与待选化合物形状的相似性需要进行量化比对才能做出合理的评估。在评估过程中，发明人采用两种度量来分布描述分子的相似性：

(1)形状描述因子{F_nl}：

该形状描述因子{F_nl}是通过将展开系数{C_nlm}中具有旋转依赖性的阶数m求和所得的一维向量。也就是说，将Zernike矩{C_nlm}中具有相同阶数(n,l)的参数排列成一个向量，用该向量的长度来获得旋转不变的参数F_nl。形状描述因子 按式(2)计算：

其中，n、l为阶数，而第三个阶数m的范围从-l到l。公式(2)每一个元素都是Zernike矩，{C_nl,-l,...,C_nl,l}为密度函数ρ(r)在阶数为(n，l)的Zernike函数的展开系数的集合，形状因子F_nl则是对应的集合的几何长度。

该形状描述因子{F_nl}有两个特征，第一是把旋转变量依赖的阶数m通过求和给去掉了，降低了数据量；第二是这个因子可以直接用于计算分子形状的差异，不必再进行昂贵的模型旋转操作。即计算{F_nl}的差异时，无需将分子在三维空间内进行旋转以搜索最大的重合度，因此极大的降低了计算的数据量，缩短了计算时间。

发明人利用三维Zernike函数，通过式(2)分别得到化合物库中每个化合物的形状描述因子和目标形状的形状描述因子通过比对形状描述因子与形状描述因子之间的差异性，对化合物库进行初步筛选。

进一步地，形状描述因子与形状描述因子之间的差异性是通过与差值的二阶矩来进行评价。优选地，该差异性采用由式(3)计算的卡方评分X²来表示：

式(3)中，n代表Zernike展开的主阶数，从0到设定值n_max,l代表依赖于主阶数n的副阶数，从0到n_max,且需要满足(n-l)是偶数。这一步的计算是简单的向量相减，速度非常快。

进一步地，初步筛选包括：将化合物库中每个化合物按照卡方评分X²从小到大的顺序进行排序，取排名靠前的10～20％的化合物作为备选化合物，进入后续的二次筛选中；或者取卡方分布中分数低于平均值一个标准差的化合物作为备选分子，进入后续的二次筛选中。

利用旋转不变量的形状描述因子{F_nl}对化合物库进行初筛，能够去除形状与目标形状差异很大的化合物，将化合物库中待选分子的样本量降低至10～20％，便于后续利用Zernike展开系数{C_nlm}对备选分子形状与目标形状进行更细致的比较，极大的降低了计算量。

(2)皮尔森相关系数cc：

展开系数与所述展开系数之间的相似度通过皮尔森相关系数(Pearson Correlation Coefficient)进行评价。

利用Zernike展开系数可以对分子进行旋转得出最优的空间方位使得参与比较的分子与目标形状的重合最大，进而得出公式(4)

所描述的皮尔森相关系数cc：

式中，<ρ(r)>为分子或目标形状的密度函数；<ρ₁(r)ρ₂(r)>为分子与目标形状的重叠程度；σ(ρ₁(r))σ(ρ₂(r))为分子密度的涨落程度，r为空间向量，描述三维空间内任意一点的位置信息。

皮尔森相关系数cc的计算需要满足：(1)目标形状与分子的中心是重合的；(2)目标形状与分子的相对空间取向给出二者最大程度的重叠。注意，在具体匹配操作中，目标形状的空间取向固定不变，化合物分子应该遍历SO(3)旋转空间的全部取向，选取匹配最佳的方位得出皮尔森相关系数cc作为该化合物与目标形状的匹配度衡量。

进一步地，二次筛选包括：将初次筛选后得到的备选化合物按照皮尔森相关系数cc由大到小的顺序进行排序，取排名靠前的10～20％的化合物作为备选分子；或者，取皮尔森相关系数cc超过预设阈值的化合物作为备选分子，该预设阈值依赖于最高阶数n_max，对于n_max＝20、cc>0.7都在备选分子的可接受范围，对于cc>0.8的化合物优先考虑作为备选分子。

在进行二次筛选的过程中，需要计算每个备选化合物与目标形状之间的皮尔森相关系数cc。计算皮尔森相关系数cc需要将备选化合物在三维空间进行旋转以搜索最大的重合度，所以计算量非常大。通常计算一个备选化合物与目标形状的相关系数需要至少1秒钟的CPU时间(以目前主流计算机配置为例，CPU主频2.0GHz，不考虑内存的限制)，所以如果对1000万个分子进行同样的计算，所需的时间大概2700CPU小时，用100个CPU同时计算也要至少一天的时间。而通过本实施方式提供的这种药物筛选方法，经过首轮筛选(即利用化合物与目标形状的形状描述因子F_nl之间的差异进行的筛选)，已经剔除了80～90％的化合物，仅需对剩余的10～20％的化合物进行皮尔森相关系数cc的计算，因此，极大的降低了计算量，提高了筛选效率。

进一步的，为了更加快速准确的评估化合物与目标形状之间的匹配程度，在进行皮尔森相关系数cc的计算过程中，采用近似密度函数取代公式(4)中的密度函数，近似密度函数采用式(5)计算：

式中，C_nlm为展开系数或Z_nlm是化合物或目标形状的Zernike函数。

式(5)表示原始密度函数ρ(r)可以通过以展开系数作为权重的Zernike函数线性叠加得到近似密度函数其近似度由最高阶数n_max决定，如果是n_max无穷大，则原始的密度函数可以完美描述。然而由于计算量的限制，在实际应用中，最高阶数n_max的数值要根据实际情况来选择，以平衡计算量和近似误差。因为Zernike函数展开是在另外一组正交系统里描述三维空间的形状，低阶展开描述的是最显著的特征，阶数越高描述的内容越精细。因此，在基于形状匹配进行的药物筛选的应用中，发明人选择最高阶数n_max为10～30，优选为18～25，预设值为20。阶数越高，对原有函数的描述越精确，同时，计算量也相应上升，推荐最高阶数为n_max＝20。图1展示了n_max设为20的情况下，用三维Zernike函数近似得到的小分子的三维形状。

为了提高计算速度，结合Zernike展开系数{C_nlm}的性质，相关系数的计算利用快速傅里叶变换进行全局扫描，再针对重点空间取向利用单纯形方法(Simplex Method)进行局部优化。这种方法既提高了计算速度，也对计算结果的准确性有很好的保证。

具体方法描述如下：

三维Zernike函数的旋转操作函数与球谐函数一样，对三个欧拉角(α，β，γ)确定的旋转操作可以分解为式(6)：

其中，系数和旋转函数表达式如下：

注意，式(7)中的和为三维Zernike矩，分别对应阶数(n,l,m₁)和(n,l,m₂)。星号*为复数共轭。

从式(8)可见，旋转角((α和γ)可以通过傅里叶变换来加速，而β可以通过对简化的wigner操作算子计算出来。

通过采用近似密度函数能够更加准确的描述化合物的形状，进而提高评估备选化合物与目标形状之间相似度的准确性，并最终提高药物筛选的成功率。

进一步地，在步骤S4之后，该药物筛选方法还包括：

将经过二次筛选后得到的备选分子与靶标生物大分子进行结构对接，并采用分子动力学的方法(比如Gromacs，Amber，Charmm，NAMD等软件)进行优化以对所述备选分子进行三次筛选，得到用于实验检测的目标分子。

利用三维Zernike函数展开得到的形状匹配结果已经提供了最优的空间取向，因此可以直接用于与受体蛋白结合位点的对接。具体操作是，把二次筛选后得到的备选分子固定到最优的空间取向，然后平移到靶标生物大分子的结合位点，就可以得到一个初步的对接模型。接下来可以用分子模型来进行进一步优化，缓解化学分子与靶标生物大分子之间的挤压，增强有利的相互作用。优化之后的化学分子可以用分子动力学方法进行平衡状态下的仿真模拟，来进一步改善化学分子与靶标生物大分子在结合位点的相互作用。另一方面，数值模拟的结果也可以用来评估该化学分子是否能够保持与靶标生物大分子的稳定结合。

进一步地，在步骤S2和/或步骤S3之前，该药物筛选方法还包括：

根据靶标生物大分子的药效团对化合物库进行预筛选，取符合药效团的化合物集合进入初步筛选。

药效团(Pharmacophore)是针对任一给定的受体分子的结构而定义的功能集团的分布。药效团能够在三维空间内有效描述功能集团(比如芳香环、疏水集团、正负电荷富集区域)在受体分子结合位置的分布情况，具有很高的敏感度和准确性。药效团的选择可以参考Molecular Operating Environment(MOE)商用软件。

尽管基于药效团的分子筛选在技术上比较成熟，但直接采用这种方法，却不能准确区别分子的大小和形状，导致大量化学复合物被错误地选为备选药物分子，且筛选到的备选分子数量依旧庞大。而在本实施方式中，发明人在进行初步筛选前，先基于药效团对化合物库中的化合物进行预筛选，将符合药效团的化合物作为备选分子，进入基于形状描述因子{F_nl}的初步筛选，以降低初步筛选的样本量，提高筛选效率，最后再进行精确的形状匹配，利用皮尔森相关系数进行筛选。

下面结合实施例对本发明的特征和性能作进一步的详细描述：

实施例1

本实施例以大麻素受体蛋白(CBI)为靶标生物大分子，示例本发明提供的这种基于分子形状匹配的药物筛选方法，流程图如图2所示：

a.基于药效团的预筛选：

利用商用软件Molecular Operating Environment(MOE)分析大麻素受体蛋白的结构，并选择和设计药效团，如图3所示：图3中a图为大麻素受体蛋白的药效团分布；b图为大麻素受体蛋白环绕抑制剂分子的二维描述，其中，HYD为疏水基团；ARO-ARO为连接的芳香环；ARO+为含有电荷的芳香环；H-BOND为氢键；ARO为芳香环。

再根据药效团对ZINC数据库的2100多万个化合物进行预筛选，筛选后得到符合药效团排布的化合物约为24万个，进入下一轮的筛选。

b.目标形状的构建及的计算：

采用大麻素受体蛋白的活性底物(PubChem库中的ZDG分子，PubChem CID:122198481)的晶体结构，对该大麻素受体蛋白的活性部位的腔室进行模拟，构建得到目标形状，如图4所示。

再根据图4中的目标形状，利用三维Zernike函数，按照式(1)计算展开系数

计算得到的为下列zernike矩(此处仅示意性的列举前10个Zernike矩)，其对应的阶数{(n,l,m)}为{(0,0,0),(1,1,-1),(1,1,0),(1,1,1),(2,0,0),(2,2,-2),(2,2,-1),(2,2,0),(2,2,1),(2,2,2)}。

{C_nlm}＝{(0.00417339628552081+0j),

(0.00021995729734553622+0.00021250111777450047j),

(0.00031106659304898837+0j),

(-0.00021995729734553622+0.00021250111777450047j),

(-0.005572994736218416+0j),

(-0.001092886706908194-0.0011740970007127437j),

(0.000143022177348094+0.0001265828466184282j),

(-0.0010662958985742223+0j),

(-0.000143022177348094+0.0001265828466184282j),

(-0.001092886706908194+0.0011740970007127437j),...}

再根据目标形状和展开系数按照式(2)计算形状描述因子

计算得到的{F_nl标}。

与上述{C_nlm标}对应的Zernike形状因子{F_nl标}为：{F00,F11,F20,F22}＝{8.708618277999448e-06，1.4191915036651722e-07,2.3258315531360964e-05,1.8707012318107897e-05}

c.按照步骤b的方法，分别计算步骤a中筛选得到的化合物的展开系数{C_nlm分子}和形状描述因子{F_nl分子}，计算结果略，此处仅实例性的列出文献中研究过的两个化合物所对应的{C_nlm分子}和{F_nl分子}，给出一个直观的比较。

(1)化合物1(Otenabant)：

{C_nlm分子1}＝{(0.006519774637355752+0j),

(0.00038327102106809803+0.00042291974738548673j),

0.0005980988425479588+0j),

(-0.00038327102106809803+0.00042291974738548673j),

(-0.008142121004171972+0j),

(0.0021276007610844315+0.0001089353107761309j),

(-0.00025149263105106755-0.00021787062155226164j),

(-0.0022291231709351926+0j),

(0.00025149263105106755-0.00021787062155226164j),

(0.0021276007610844315-0.0001089353107761309j)}

{F_nl分子1}＝{2.1253730660953664e-05,5.046189010477904e-07,5.308479401728197e-05,4.028083033059879e-05}

(2)化合物2(AM6545)：

{C_nlm分子2}＝{(0.004244131815783876+0j),

(0.00023114156670208785+0.0002460539258441589j),

(0.0003690620595496471+0j),

(-0.00023114156670208785+0.0002460539258441589j),

(-0.005352091885838367+0j),

(0.001675496587967557+0.00013809037812919435j),

(-6.444217646029072e-05+2.7618075625838835e-05j),

(-0.0017009191373026883+0j),

(6.444217646029072e-05+2.7618075625838835e-05j),

(0.001675496587967557-0.00013809037812919435j)}

{F_nl分子2}＝{9.006327434874472e-06,1.820723601803322e-07,2.271498345368534e-05,1.8600280054080706e-05}

d.基于{F_nl}的初步筛选：

按照式(3)来计算每个待筛选化合物的卡方评分X²，得出每个分子的形状描述因子F_nl分子与F_nl标之间的差异性。

得到每个化合物的卡方评分X²(以化合物1和化合物2为实例)：

目标形状与化合物1的卡方为1.347x10^-8；

目标形状与化合物2的卡方为1.125x10^-9；

由计算结果可知，化合物2与目标形状更为近似。这与实验事实完全一致，化合物2被认定为激活剂，而化合物1是抑制剂分子。

将待筛选的化合物库中(约24万个)化合物按照卡方评分X²从小到大的顺序进行排序，取排名靠前的15％的化合物作为备选化合物(约3万个备选分子)，进入下一轮的筛选。

经过药效团和形状因子的筛选之后，备选分子2100多万个化合物已经降低为3万多个。利用下一步的形状匹配方法，用单台计算机10个小时之内就可以完成评估并提出最可能的结合方位。

e.基于{C_nlm}的二次筛选：

将备选化合物的展开系数与目标形状的展开系数进行对比，利用如式(4)所示的皮尔森相关系数cc评价二者的相似度，对备选化合物进行二次筛选。

其中，式(4)中的化合物的密度函数ρ(r)为近似密度函数如式(5)所示：

经过计算后，按照皮尔森相关系数cc排序化合物作为备选分子，共筛选得到24个备选分子，其相关卡方和相关系数及结构如表1所示：

表1.筛选得到的备选分子的相关信息

由此可见，通过步骤a～e的分析和评估，发明人从ZINC数据库的2100多万个化学分子中，逐步把研究重点集中到24个备选分子的集合上，并对24个备选分子进行后续步骤的研究。

f.模型对接及分子动力学模拟：

再通过模型对接及分子动力学模拟对上述24个备选化合物进行优化，判别备选分子的化合位置和稳定性，对备选化合物进行优化筛选。应用Gromacs软件(开源软件)进行了500纳秒的全原子模型动力学模拟，发现这24个化合物都有很好的结合性，没有破坏受体蛋白的结构，也没有从结合位点脱离的迹象。

g.实验检测：

对f步骤筛选得到的备选小分子进行实验检测，得到大麻素受体蛋白的拮抗剂。

尽管已用具体实施例来说明和描述了本发明，然而应意识到，在不背离本发明的精神和范围的情况下可以作出许多其它的更改和修改，比如，最直接的一个改进是用同样的方法来描述电荷在三维空间的分布，从而实现化合物的电荷分布于结合位点区域的电荷呈互补的形式存在。因此，这意味着在所附权利要求中包括属于本发明范围内的所有这些变化和修改。

Claims (10)

1.一种基于分子形状匹配的药物筛选方法，其特征在于，其包括：

根据与靶标生物大分子的结合部位相符的目标形状，利用三维Zernike函数计算所述目标形状的展开系数C_nlm标和形状描述因子F_nl标；

利用三维Zernike函数，计算化合物库中每个化合物的展开系数C_nlm分子和形状描述因子F_nl分子；

通过比对每个化合物的形状描述因子F_nl分子与所述目标形状的形状描述因子F_nl标之间的差异性，对所述化合物库进行初步筛选得到备选化合物；以及

再将所述备选化合物的展开系数C_nlm分子与所述目标形状的展开系数C_nlm标进行对比，对所述备选化合物进行二次筛选；

所述备选化合物的展开系数C_nlm分子与所述目标形状的展开系数C_nlm标之间的相似度通过皮尔森相关系数进行评价；

利用三维Zernike函数的展开系数C_nlm标对备选化合物分子进行旋转得出最优的空间方位使得参与比较的分子与目标形状的重合最大，在具体匹配操作中，目标形状的空间取向固定不变，化合物分子遍历SO(3)旋转空间的全部取向，选取匹配最佳的方位得出皮尔森相关系数cc作为该化合物与目标形状的匹配度衡量；相关系数的计算利用快速傅里叶变换进行全局扫描，再针对重点空间取向利用单纯形方法进行局部优化。

2.根据权利要求1所述的基于分子形状匹配的药物筛选方法，其特征在于，所述化合物的形状描述因子F_nl分子与所述目标形状的形状描述因子F_nl标之间的差异性是通过所述化合物的形状描述因子F_nl分子与所述目标形状的形状描述因子F_nl标之间差值的二阶矩来进行评价。

3.根据权利要求2所述的基于分子形状匹配的药物筛选方法，其特征在于，所述化合物的形状描述因子F_nl分子与所述目标形状的形状描述因子F_nl标之间的差异性采用由式(3)计算的卡方评分χ²来表示：

式(3)中，n代表Zemike展开的主阶数，从0到设定值n_max，1代表依赖于主阶数n的副阶数，从0到n，且需要满足n-1是偶数。

4.根据权利要求3所述的基于分子形状匹配的药物筛选方法，其特征在于，所述初步筛选包括：将所述化合物库中每个化合物按照卡方评分χ²从小到大的顺序进行排序，取排名靠前的10～20％的化合物作为备选化合物；或者根据卡方分布，选取分数低于平均值一个标准差的化合物作为备选化合物。

5.根据权利要求1所述的基于分子形状匹配的药物筛选方法，其特征在于，所述备选化合物的展开系数C_nlm分子与所述目标形状的展开系数C_nlm标之间的相似度通过皮尔森相关系数cc进行评价，取皮尔森相关系数cc超过预设阈值的化合物作为备选分子。

6.根据权利要求5所述的基于分子形状匹配的药物筛选方法，其特征在于，所述皮尔森相关系数cc采用式(4)计算：

式中，＜ρ(r)＞为分子或目标形状的密度函数；＜ρ₁(r)ρ₂(r)＞为分子与目标形状的重叠程度；σ(ρ₁(r))σ(ρ₂(r))为分子密度的涨落程度；r为空间向量。

7.根据权利要求6所述的基于分子形状匹配的药物筛选方法，其特征在于，式(4)中所述密度函数ρ(r)为近似密度函数所述近似密度函数采用式(5)计算：

式中，C_nlm为展开系数C_nlm分子或C_nlm标，Z_nlm是化合物或目标形状的Zernike函数；最高阶数n_max为10～30。

8.根据权利要求1所述的基于分子形状匹配的药物筛选方法，其特征在于，还包括：将经过二次筛选后得到的备选分子与所述靶标生物大分子进行结构对接，并采用分子动力学的方法进行优化以对所述备选分子进行三次筛选，得到用于实验检测的目标分子。

9.根据权利要求1所述的基于分子形状匹配的药物筛选方法，其特征在于，在对所述化合物库进行初步筛选之前，还包括：根据靶标生物大分子的药效团对所述化合物库进行预筛选，取符合所述药效团的化合物集合进入所述初步筛选。

10.根据权利要求1所述的基于分子形状匹配的药物筛选方法，其特征在于，所述目标形状是通过模拟所述靶标生物大分子的活性部位来定义的三维空间形状；或者，所述目标形状是通过抽提能够与所述靶标生物大分子结合的活性分子的晶体结构来确定的三维空间形状。