药物探索就如一场数字竞赛

 在药物研究中，可能存在的有机分子的数目估计超过10^60，但高通量筛选(HTS)方法的极限为10^6到10^9，其中通常小于10^3分子具有所需的生物活性。为了使得这种方法能发挥作用，我们必须确保我们的10^6个子集的分子能够代表比10^60更大的集合。然而以制药公司的传统方式收集到的分子通常不能代表更大的“化学空间”。是否有方法使收集的数据更具代表性是挑战之一。接下来的挑战是设计的新药既具有最大限度的疗效又能使昂贵的合成过程最少。通常统计学家可能会考虑设计实验(DOE)来实现这一点，但在多维世界的化学空间，这是一项具有挑战性的任务。JMP的独特的数据交互式可视化功能， DOE能力和数据整理工具，使我们能够结合化学原理的方法，系统地探索和评估大型复杂的数据集。通过这种方式，我们分析现有的数据，确定下一步的工作，大大加快了药物研究的进程。

药物研究是一个漫长的旅程，从数以百万计的潜在的开始点一直到一个单一的可能有一天会成为一种市场化药物的候选分子。纵观整个医药行业，从项目开始到一个药物到达市场约需要11到14年的时间。每一个成功药物的平均成本是50亿美元(福布斯，2014)。耗费如此长的时间和巨大成本的原因有很多。从统计角度的一个原因是我们在寻找一个极不可能的结果。需要给出可以把原子放在一起组成分子的各种各样的方法，并考虑自然规律与药物分子的限制，我们可以估计药物分子可能存在的总数。参与研究的原子总数多达17个，这就形成了1660亿个组合。药品经常包含40个或更多原子，估计就达到10^60个是潜在的药物分子了。

传统的方法从化合物的大集合中通过蛋白化验来寻找对靶蛋白具有某些活性的药品分子(亦称hits)，这种方法也被称为高通量筛选(HTS)，但是这种方式可以支持筛选的化合物的最大数量约是一百万。最近我们关注到基因编码库的创新方法，可以同时进行几个化合物的测试。这些方法能够相对容易上网的筛数十亿化合物，然而这仅仅是可能的化合物的小部分。如何能够成功地运用这些方法，是巨大的挑战。

在下一阶段的药物研究中我们面临类似的问题。有效复合物(或一系列类似的化合物组合)的识别使药物研究人员能够聚焦在单一分子的“结构支架”。紧接着的任务是人工合成和测试类似的化合物，寻找核心结构支架上不同取代基组，分析局部结构活性的关系(SAR)，见图1。这里的挑战是，即使限制在对少量的原子分析，依然有数百个潜在的子结构组选择，如果考虑原子越多选择组合就越多。作为药物研究人员迫切需要知道如何选择化合物才能快速、高效地探索结构活性关系。

为了简化问题，我们运用了类似化学结构具有相似性能的理论。这意味着一个多个维度的“化学空间”的存在，其中涵盖所有潜在有效的化学分子结构。因此一个理想的HTS筛选集合是具有均匀化学空间属性的一组化合物。从概念上讲，这是一个很好的方案，但我们没有有效的化学空间界定可以使用。我们有多种计算化学结构性能的方法，但发现经常出现化学结构的描述不完全，以及潜在结构的范围太大的问题。在现实中，虽然药物公司定期收集新化合物来补充筛选的需要，但是许多HTS的结论只有非常微弱的“hits” ，有的甚至没有“hits”，这是因为化合物的筛选集不能代表充分的化学空间。

                                                             图1 典型的药品(Iressa)分子结构显示                     图2 GPCR有效的化合物结构支架结构支架(黄色)和取代基(蓝色) X是分析对象组的位置

当化学反应锁定在分子的单个结构支架时，我们就可以确定可用的范围，开始应用统计方法使化合物的多样性最大化。图1显示了一个具有三个取代基组围绕在固定结构支架的化合物。在第一轮的探索中，我们倾向依次改变每一个结构组(在随后的几轮进行最好的组合)，得到仅有一个取代基位置不同的类似的化合物组成的合成库。

传统的设计化学合成库的方法是选择一个易于处理的合成方法，并利用已有的试剂合成尽可能多的化合物。虽然在成本方面是比较经济，但是单一的合成路线将导致化合物缺乏多样性。

我们希望利用量化的因子来准确地描述代表特定的化合物集合，因此采用了实验设计的(DOE)原理，实现了使用更少的化合物也能确保更多的多样性。

通过筛选确认了一个A类G蛋白偶联受体有活性的化合物，其分子结构支架包含取代苯基环。我们希望探索取代基位置的结构活性关系。基于行业经验和分子结构不同的基本属性的原理，我们选择三个属性代表94个小的化学取代基。

这些属性是：

l Pi(一种亲脂性度量，一个区分水和有机溶剂指标);

l Molecular Refractivity(MR，衡量分子的大小);

l Sigmap(测量分子的电特性，如吸电子或供电子)

从主成分分析(PCA)看MR和Pi有一定的相关性，然而对于模型都是重要的，所以保留这两项在模型中，以示区分。

在使用连续变量作为DOE的输入时，有些组合在化合物集中是不存在的，因此我们通过将原有变量分成高中低三类转换成分类变量(“高”与“低”为上下四分位数，“中”是的中位数)。基于这些输入信息，定制设计给出了最小的9个化合物来表示这个集合，从而进行合成。(包含原先的一个)共有10种化合物可用于建模。

图3 10个初始化合物的JMP多元回归的                       图4 训练集的化合物3D结构图显示了大

结果显示Pi，MR，Sigmap与 pIC50 型取代基在分子内形成57度的夹角

无明确的联系

运用简单的多元线性回归(MLR)模型对10个化合物训练集的数据进行拟合，令人失望的是测量活性(pIC50)没有构建出理想的模型，参见图3。

这时我们需要考虑用来描述这些化合物其他属性，特别是可能描述的离群化合物和解释它们的显然反常行为的属性。我们意识到离群化合物有明显的不同于大多数集合的形状。我们使用分子力学计算每个分子的最小能量构造，选定了苯基环和相邻的羰基之间的面夹角(见图4)代表形状差异。在MLR模型中引入这个因素后极大地优化了模型的拟合度。Pi和Sigmap是不显著的(在95%置信度)，因此从模型中移除。最终的模型如图5所示，MR和苯基环和相邻的羰基之间的面夹角解释了～63%的pIC50a的方差。

图 5 10个初始化合物的JMP多元回归的                图 6 测试集的17个化合物的测量数据

结果显示MR和pIC50存在相关性 显示了多元回归模型的预测效果

我们也容易解读一个比较简单的模型，从而理解观察到的内在关系。图5参数估计表显示pIC50与MR有负相关，与面夹角正相关，即最有效的化合物是那些具有最大的扭曲度和最小尺寸的化合物。在此基础上又合成了17种化合物，其中大部分很好地预测了pIC50(同时加入了一些负面对照组)。图6显示了17种化合物训练集pIC50的实际测量值与预测值的关系。虽然不是每一个化合物都预测得很好，但其中10个化合物(图6中所示的空心方块)的预测是有效的，并被证实。

令人满意的是该模型的预测(告诉我们无效的化合物特性)和解释(告诉我们哪些属性对于效力是重要的)使我们确定了现有集合中可以实现最佳效力的取代基组，即使考虑更大的化合物集合也成为可能。我们也能确定不值得扩大的范围，比如使用具有更多原子的更大的取代基，因为我们已经观察到效力与分子大小的负相关性。

在这项工作中，我们已经表明，实验设计(DOE)的原理可用于药物研究，但仍需然仔细考虑构建SAR的问题，并尽可能地避免变异，从而使得通过少量因子代表化学多样性成为可能。

我们还发现，因子选择需要一定数量的试错法，很难说哪个因子是最重要的，是第一优先的。DOE在这个案例的成功运用，使得深入的探索结构活性的全貌，仅通过四分之一的可能化合物就聚焦在最有价值的属性空间。