新型苯并呋喃衍生物作为LSD1抑制剂的QSAR研究

QSAR Study of Novel Benzofuran Derivatives as Potent LSD1 Inhibitors

摘要: 组蛋白赖氨酸特异性去甲基酶1 (LSD1)是一种黄素腺嘌呤二核苷酸(FAD)依赖性胺氧化酶，可特异性识别H3K4和H3K9底物，并去除其单甲基或二甲基修饰。它介导许多细胞内信号通路，与肿瘤的发生和发展密切相关。因此，开发高效、特异的LSD1抑制剂不仅有利于研究LSD1的生物学功能，而且对抗肿瘤药物的开发具有重要的科学意义。建立定量构效关系(QSAR)模型可以预测分子的物理和化学性质。在本研究中，利用基因表达编程(GEP)建立了一个具有描述符的非线性定量构效关系(QSAR)模型，并预测了一系列新型苯并呋喃化疗药物的活性。这些描述符是在CODESSA软件中计算的，并基于启发式算法从描述符池中选择。选择4个描述符来建立多元线性回归模型。获得了训练集和测试集的最佳非线性QSAR模型，相关系数分别为0.92和0.80，平均误差分别为0.07和0.60。显然，基于GEP的QSAR模型具有更好的抑制剂疗效预测稳定性。这些发现对LSD1抑制剂作为高选择性的一流临床候选药物的设计提供了新的价值。

Abstract: Histone lysine specific demethylase 1 (LSD1) is a flavin adenine dinucleotide (FAD) dependent amine oxidase, which can specifically recognize H3K4 and H3K9 substrates and remove their monomethyl or dimethyl modifications. It mediates many intracellular signal pathways and is closely related to the occurrence and development of tumors. Therefore, the development of effi-cient and specific LSD1 inhibitors is not only conducive to the study of the biological function of LSD1, but also has important scientific significance for the development of anti-tumor drugs. Estab-lishing a quantitative structure-activity relationship (QSAR) model can predict the physical and chemical properties of molecules. In this study, gene expression programming (GEP) was used to build a nonlinear quantitative structure activity relationship (QSAR) model with descriptors and to predict the activity of a series of novel DNA-targeted chemotherapeutic agents. These descriptors were calculated in CODESSA software and selected from the descriptor pool based on heuristics. Four descriptors were selected to establish a multiple linear regression model. The best nonlinear QSAR model with a correlation coefficient of 0.92 and 0.80 and mean error of 0.07 and 0.60 for the training and test sets were obtained. It is apparent that the QSAR model based on GEP has better forecasting stability of inhibitor efficacy. These findings should be useful for the design of LSD1 in-hibitors as highly selective first-in-class clinical candidate.

1. 背景

肿瘤的发生和发展与遗传物质密切相关，即使核苷酸序列没有改变，仍然存在基因表达差异的可能性，从而导致疾病的发生。这种调控机制属于表观遗传学的范畴，主要包括DNA甲基化、组蛋白修饰和微小RNA的转录后调控等。DNA附着在组蛋白核小体上，其转录受组蛋白修饰的调节。组蛋白特异性去甲基酶1 (LSD1)可以调节组蛋白修饰，其异常表达与许多肿瘤疾病的发生和预后密切相关 [1] 。

组蛋白去甲基化酶1 (LSD1)，又被称作KDM1A/AOF2，是一种高度保守的腺嘌呤二核苷酸(FAD)依赖性氧化酶 [2] ，LSD1可以特异性地去除组蛋白H3第4位赖氨酸(H3K4)和第9位赖氨酸(H3K9)的单、二甲基化，从而调节基因转录(包括激活或抑制)，而这种表达调控在癌症的发生中起着关键作用，因此LSD1的过度表达可能导致抑癌基因的异常沉默 [3] [4] [5] [6] 。此外，在细胞和动物模型中对LSD1的敲除或抑制的研究表明，LSD1的减少可以导致H3K4me2的水平增加，并重新激活肿瘤抑制基因的表达 [7] [8] [9] 。这些发现表明抑制LSD1是治疗癌症的有效策略。目前相关研究已经发现，LSD1在肺癌、胃癌、结直肠癌、膀胱癌、食管癌、前列腺癌、胶质瘤、急性髓系白血病等癌症疾病中高表达，其表达程度与癌症分级、恶性程度和预后有关 [2] [10] [11] [12] [13] 。

迄今为止，在文献中已经发现了许多LSD1抑制剂，并根据其分类作用机制。可以考虑两种不同类型的抑制剂：1) Tranylcypromine (TCP)及其类似物作为不可逆抑制剂的代表，其通过与辅因子FAD共价结合来抑制LSD1的活性；2) 可逆性抑制剂可以与底物竞争结合活性位点 [14] - [26] 。据相关研究所知，苯并呋喃衍生物首次被描述为LSD1抑制剂 [27] 。因此，可以分析化合物的结构，为新药的开发提供帮助。

定量构效关系(Quantitative Structure-Activity Relationship，简称QSAR)即通过数学方法建立化合物的分子描述符与其生物活性/毒性之间的线性或非线性关系模型，在分子水平阐明结构与生物学及物理化学特性之间的关系。近十几年来，随着分子生物学和计算机科学的快速发展，使QSAR研究提高到了一个新的水平，计算机辅助药物设计(CADD)在新药的研发中起着越来越重要的作用。QSAR研究对于设计和筛选生物活性显著的药物，以及阐明药物的作用机理等均具有指导作用 [28] [29] [30] [31] 。

建立QSAR模型的数学算法有很多，如逐步回归(SMR)、多元线性回归(MLR)、径向基函数神经网络(RBFNN)、人工神经网络(ANN)和支持向量机(SVM)等。近年来的研究表明，基因表达编程(GEP)算法具有简单、灵活和高效的优点，因此GEP方法被广泛应用于QSAR模型的建立。

本研究旨在建立了一个可靠的QSAR模型，使用GEP方法揭示靶向LSD1的新型苯并呋喃衍生物的分子结构与其肿瘤抑制活性之间的关系。与此同时，该模型可以为进一步探索和预测癌症新靶向药物提供重要的理论基础。

2. 计算与研究方法

2.1. 数据收集及其划分

QSAR研究的第一步是选择合适的数据和建立数据集，化学物的活性应该在活性范围内均匀分布，为训练集选择的化学物质应具有足够的结构多样性，以跨越与所研究的生物活性相关的化学空间范围。从文献 [27] 获得38种化合物和IC50值，并列于表1中。为了获得良好且有意义的结果，我们使用系统时间作为随机种子；我们将数据集分为30个复合训练集和8个复合测试集。训练集用于构建、训练和优化模型，测试集用于评估模型的预测能力。

Table 1. IC50 of 38 compounds (HM)

表1. 38种化合物的IC50值(HM)

2.2. 描述符的选择与计算

在QSAR研究中，良好的模型的构建必须使用准确合理的分子描述符，分子描述符测量分子性质的某些方面，无论是其物理和化学性质，还是通过基于分子结构的各种算法导出的数值指数。已经报道了各种各样的描述符用于QSAR分析，例如静电、拓扑、几何和量子化学描述符。

在这项研究中，使用Chemdraw软件绘制了化合物的化学结构，然后将分子结构导入Hyperchem [32] 软件中进行优化，实验的最低能量稳态是化合物的最佳状态。目前最流行的半经验量子化学计算方法中使用范围最广泛的方法，也是精确度、准确度最高的方法是AM1 (Austin Model 1)和PM3 (Parametric Method 3)，为了实现化合物的最佳实验状态，使用半经验AM1或PM3方法 [33] 进行精确优化。用Polak-Ribiere算法对分子结构进行优化，直至均方根梯度达到0.01，然后得到三组描述子。

2.3. 启发式线性模型

启发式算法(heuristic method)是二维定量结构–活性关系研究常用的方法之一，在计算机药物辅助设计中也得到广泛应用。通过CODESSA软件中的启发式算法可对药物分子结构进行分析，对计算得到的大量分子描述符进行筛选，能够快速建立预测药物活性的最佳多元线性方程。对方程中描述符的深入研究，可获得影响药物活性的理化因素，为随后的药物设计提供思路。HM建立线性关系的步骤模型是：首先，选择单参数描述符，使用R²、F检验和T检验作为标准选择，并删除相关性低的描述符。其次，选择两个参数描述符，并使用R²和F检验作为标准选择。最后，选择n个参数描述符。在获得具有最佳统计特征的双参数相关系数之后，添加在先前选择过程中未使用的描述符。重复此步骤，直到建立的相关方程包含最多的参数。模型评估的标准为R²、F检验值、t检验值和R²cv。因此，HM [34] 开发了一个具有4个描述符的线性模型。

2.4. GEP的非线性模型基因表达式编程

遗传表达式编程(Genetic Expression Programming, GEP)是一种用于解决机器学习、数据分析、优化和模式识别等领域问题的进化算法 [35] 。开发了GEP算法来自遗传算法(Genetic algorithm, GAs)和遗传规划算法(Genetic programming, GP) [36] ，但GEP与GA和GP使用的个体有本质区别。GEP有两个基本元素，染色体遗传操作子会作用于染色体而不是表达树，这样保证每一次基因操作都是正确的，不需要检查其是否符合句法 [37] 。近年来，GEP算法遗传算法编码具有简单方便的特点操作和灵活性，保留了GP表达式搜索能力强大的优势，通过使用不同的基因组和表现型巧妙地避免了GA和GP的缺点，提高了解决

Figure 1. GEP flow chart

图1. GEP算法的计算步骤

问题的效率 [38] [39] 。

在GEP中，基因由头部和尾部组成。标题由终端和功能符号组成，尾部仅由终端符号组成。假设头部长度为h，尾部长度为t，关系如下：

$t=hn-1+1$

在该方程中，n表示需要函数符号集中最多变量的函数。显然，基因的长度可以通过以下公式计算：

$l=h+t=n\ast h+1$

GEP算法的主要步骤如图1所示。首先，一定数量的染色体个体是随机的作为初始种群生成并表达。接下来，从一组适合度样本计算每个个体的适合度。对群体循环执行遗传操作，直至满足终止准则为止 [40] 。否则，将根据适合度值选择这些个体进行基因操作。最后，产生具有新特征的后代。重复该过程，直到出现良好结果。

3. 结果及讨论

3.1. HM的结果

CODESSA软件为每个计算了487个描述符合物。为了获得与LSD1抑制剂活性最相关的一组描述符，首先建立了描述符数量范围为1到9的线性回归模型。不同数量的描述符对R2和S的影响如图2所示。结果表明，随着描述符数量的增加，R2增加，S2减少。在描述符数量增加到5之后，R2的增加速率和S2的降低速率略有下降，并且模型使用过多的描述符的概括能力很差。综合考虑后，模型具有4个描述符被选为预测抑制剂活性的最佳线性模型。

Figure 2. The effects of different numbers of descriptors on R2 and S2

图2. 不同数量的描述符对R2和S2的影响

表2显示了这些描述符的详细信息。如前所述，模型的实用性是最终目的，但要应用模型就必须保证模型有可靠高效的预测能力。因此本论文从建立QSAR模型的各个步骤考虑，试图解决目前QSAR研究中还有待完善的问题，特别是在化合物活性构象的选择、引入新的建模方法以及新的建模策略等方面，旨在尽可能考虑QSAR的不同角度提高模型的可靠性以及外部预测能力。同时，把建立的具有很好预测能力的定量构效关系模型用于活性化合物的设计和筛选。

HM模型的图形如图3所示。

Table 2. The selected molecular descriptors and their physical-chemical meaning, coefficient, and T-test

表2. 选择的分子描述符及其物理–化学意义，系数和T检验

Figure 3. Plot of measured and calculated IC50 by HM

图3. 通过HM测量和计算IC50的绘图

线性模型方程如下：

IC50 = 2.3176 + 58.489*氢键供体原子的表面电荷面积加权表分数 + 1.351*环的相对数量 + 0.0098*补充信息内容 + 5.2309*H原子的最小部分电荷

根据方程中系数的绝对值，描述符对LDS1抑制剂活性的影响如下：补充信息内容 > 氢键供体原子的表面电荷面积加权表分数 > H原子的最小部分电荷 > 环的相对数量。所选的四个描述符是在下文中进行了解释，以更深入地了解可能影响LSD1抑制剂的活性。

HDCA-2/TMSA表示氢键供体原子HDCA-2的面积加权表面电荷；HDCA-2的计算如下：

$HDCA2={\displaystyle \underset{D}{\sum }\frac{q_D\sqrt{s_D}}{\sqrt{S_{tot}}}}D\in H_{H-donor}$

它对分子的构象变化具有重要意义 [41] 。最大HDCA-2/TMSA系数为正值，表明HDCA-2/TSA值越高，LDS1抑制剂药物的活性越强。同时，可以设计具有更高HDCA-2/TMSA值的分子来提高LDS1抑制剂药物的活性。

Complementary Information content表示补充信息内容，计算如下：

${}^{k}CIC={\log }_{2}n-{}^{k}IC$

${}^{k}IC=-\underset{i=1}{\overset{k}{{\displaystyle \sum }}}\frac{n_i}{n}{\log }_2\frac{n_i}{n}$

Min partial charge for a H atom和Relative number of rings分别表示H原子的最小部分电荷，以及环的相对数量，他们反映了化合物的分子组成 [42] 。这三种描述符的正系数表明他们三者值越高，LDS1抑制剂化疗药物的活性越高。

3.2. GEP的结果

数据集被随机分成30个化合物的训练集和8个化合物的测试集，然后是非线性模型是使用软件自动问题求解器(APS)建立的以集成GEP算法的实现。为了获得更准确和普遍的模型，我们在APS中介绍了六个描述符，并通过GEP功能集建立了非线性模型。表3列出进化过程中使用的参数。

Table 3. Parameters for the simple symbolic regression problem

表3. 简单符号回归问题的参数

最后，训练集和测试集的相关系数分别为0.92和0.80，并且在830rd代中获得了0.07和0.60的标准误差。图4显示了GEP模型图。

Figure 4. Plot of measured and calculated log (IC50) by GEP

图4. 通过GEP测量和计算的图IC50

4. 结论

任何QSAR建模的主要目标是开发的模型应该足够强大，可以对新化合物的生物活性进行准确和可靠的预测。在本研究中，由HM，GEP建立的非线性模型在预测38种LDS1抑制剂的活性方面显示出更好的稳定性和预测能力。为了深入了解描述符的重要性，通过基于系数的绝对值对4个描述符进行排序，得到描述符对IC50的影响力大小。所建立的模型揭示了影响LDS1抑制剂活性的因素，并为进一步设计用于抗肿瘤的高效药物提供了方向和指导。

致谢

三年时间转瞬即逝，我的肿瘤学硕士研究生阶段的学习即将结束。本论文是在导师邱文生教授和郭婧老师的悉心指导和督促下完成的。导师渊博的专业临床医学知识，雄厚的科研能力，严谨的治学态度，精益求精的工作作风，诲人不倦的高尚师德，严以律己、宽以待人的崇高风范，朴实无华、平易近人的人格魅力对我影响深远。不仅使我树立了远大的学术目标还教会我怎样成为一名优秀的临床医生。

衷心感谢邱文生主任、丁爱萍老师、吕静老师、齐卫卫老师、赵淑芬老师、郭婧老师和我的同门师兄妹三年来对我工作学习上指导与鼓励。

衷心感谢青岛大学附属医院提供的良好科研平台。

基金项目

北京市希思科临床肿瘤学研究基金会(Y-HR2018-185, Y-2019AZZD-0471)。