基于分子对接技术对维生素D₃羟化酶底物谱的筛选及验证

Screening and Validation of Vitamin D₃ Hydroxylase Substrate Profile Based on Molecular Docking Technique

摘要: 维生素D₃羟化酶(Vdh)是生物法合成维生素D₃ (VD₃)活性形式25(OH)VD₃和1α,25(OH)₂VD₃的关键酶之一。然而，至今对Vdh的天然底物了解甚少，限制了其工业上的应用。为了提高Vdh的商业价值，找到Vdh的最佳底物，探究Vdh与底物的作用机制，深入挖掘Vdh的应用潜力。本研究以VD₃与Vdh的结合模型为基础，利用Autodock分子对接技术对ZINC化合物数据库进行筛选，结合打分，性质筛选以及广义伯恩表面积法(MM/GBSA)计算结合自由能，最终从中挑选出4个小分子进行了实验验证。实验结果与筛选结果一致，4个小分子均能与Vdh发生反应，范德华力是Vdh与小分子结合的主要驱动力，推测Vdh是生产类固醇类药物的潜在生物催化剂。我们的研究为生物合成类固醇药物提供了新思路，也为深入了解Vdh提供了理论基础。

Abstract: Vitamin D₃ hydroxylase (Vdh) is one of the vital enzymes in the biosynthesis of 25(OH)VD₃ and 1α,25(OH)₂VD₃. However, little is known about the natural substrates of Vdh so far, which limits its industrial application. Therefore, in order to improve the commercial value of Vdh and to deeply explore the application potential of Vdh, it is necessary to find the best substrate for Vdh and to ex-plore the mechanism of Vdh-substrate interaction. This study is based on the VD₃-Vdh binding model, using molecular docking technology to screen the ZINC compound database, combining scoring, property screening, and generalized Berne surface area method (MM/GBSA) to calculate the binding free energy. Finally, 4 small molecules are selected for the analysis of experimental ver-ification. The experimental results were consistent with the screening results, four compounds from ZINC have the potential to be hydroxylated by Vdh, and van der Waals forces were the main driving force for the binding of Vdh to small molecules. The screening of Vdh substrate spectrum was com-pleted in this study. It was speculated that Vdh was a potential biocatalyst for the production of steroids. Our study provides new ideas for the biosynthesis of steroid drugs and a theoretical basis for a deeper understanding of Vdh.

1. 引言

从自养假诺卡氏菌中分离得到的维生素D₃羟化酶(Vdh)是一种细胞色素P450酶(CYP107BR1)，能将VD₃转化为其生理活性形式，即25(OH)VD₃或1α,25(OH)₂VD₃ [1]。研究表明，Vdh可以对VD₃进行持续羟化，第一步把初底物VD₃羟化为25(OH)VD₃，再以25(OH)VD₃为底物，进一步反应生成1α,25(OH)₂VD₃。Vdh现已在大肠杆菌，酵母和红球菌等工程菌中成功表达 [2]，它操作简单，生长速度快，转化率也进一步提高，是目前微生物转化VD₃中活性最高的酶，已被逐渐运用在工业上，具有很大的商业价值。但VD₃并不是Vdh的天然底物，Vdh对它的催化活力并不高，因此探寻Vdh的最佳催化底物是提高Vdh催化效率的一种策略。VD₃属于类固醇类，类固醇又称甾体类激素，主要由1个五元环和3个六元环组成，对生物的内分泌具有重要的调节作用。目前，以类固醇为基础的药物制剂被广泛使用，是仅次于抗生素市场化程度最高的一类药品，年产量超过100万吨 [3]。大多数类固醇物质都具有消炎和抗过敏的功能，可作为抗炎、抗肿瘤、抗过敏的药物，也可作为性激素治疗生殖系统疾病等 [4] [5] [6]。另外，在预防和治疗许多严重疾病，如癌症、肥胖症、糖尿病、类风湿性关节炎、哮喘、激素代谢综合征、神经退行性疾病等方面类固醇也发挥着重要的作用。羟化的类固醇是类固醇发挥作用的重要形式，与极性较小的非羟基化类似物相比，羟基化的类固醇通常表现出更高的生物活性 [7]。虽然目前大量运用化学法得到羟化类固醇，但研究人员依旧在寻找对类固醇工业羟化的新型生物催化剂，生产出具有更强治疗性的新型羟基类固醇。

虚拟筛选(Virtual Screening, VS)是借助计算机技术对特定的疾病或靶标，以药物设计理论为基础，从化合物小分子数据中搜寻出与指定受体结合较好的配体的技术。该方法可以在少量的实验基础上，探究酶可能催化的未知底物。有助于提高Vdh的利用范围，降低实验操作的工作量，节省大量资金成本。目前已经成为了新药发现的实用工具，是各大药物公司必不可少的一项技术 [8]。虽然虚拟筛选技术已经日趋成熟，但目前采用此方法拓宽Vdh底物谱的研究还比较缺乏。本研究以拓宽Vdh底物谱为目的，利用分子对接，分子动力学等虚拟筛选化合物的方法，开展了Vdh潜在未知底物的研究。

2. 材料与方法

2.1. 材料

含有pET28a-Vdh表达质粒的大肠杆菌BL21 (DE3)菌株由实验室构建并保藏。卡那霉素(KAN)、六水氯化铁(FeCl₃-6H₂O)、异丙基-β-D-硫代半乳糖苷(IPTG)和6-磷酸葡萄糖(G6P)购自Solarbio公司。5-氨基乙酰丙酸(5-ALA)、菠菜铁氧还蛋白(Fdx)和铁氧还蛋白还原酶(Fdr)、葡萄糖脱氢酶(GDH)、还原型烟酰胺腺嘌呤二核苷酸磷酸(NADPH)均购自Sigma公司。LB培养基购自青岛高科技工业园海博生物技术有限公司。SDS-PAGE变性丙烯酰胺凝胶快速制备试剂盒、蛋白预染Marker、二氯甲烷、甲醇等均购自生工生物工程(上海)股份有限公司，β-环糊精(PMβCD)购自国药(上海)国际医药卫生有限公司。

2.2. 蛋白诱导表达

将实验室构建的包含目的基因pET28a-Vdh的大肠杆菌BL21 (DE3)涂布在含有50 μg/mL卡那霉素的LB琼脂培养基上，37℃培养过夜。从LB琼脂培养基上随机挑取单克隆菌落加入到含有50 μg/mL卡那霉素，0.475 mM FeCl₃-6H₂O，0.5 mM 5-ALA的LB液体培养基中，37℃培养至菌液OD600 ≈ 0.6。添加 0.5 mM诱导剂IPTG，37℃条件下继续培养6 h，未添加诱导剂的为阴性对照。4000 ×g离心10 min，弃去上清，收集菌体沉淀。使用PBS (pH7.4)使菌体沉淀重悬，在冰上使用超声破碎仪充分破碎菌体。离心后分别收集上清和沉淀，沉淀使用8 M尿素，50 mM Tris HCl，300 mM NaCl (pH 8.0)进行溶解，以备下一步实验使用。

2.3. 配体的准备

从PDB数据库中下载Vdh酶与VD₃的复合物结构(PDB ID: 3A50) [9]，以其复合物的结构为模板进行虚拟筛选。该晶体使用PyMol保留同源二聚体中的A链，并分别保存Vdh蛋白和VD₃的结构。使用Chimera对蛋白进行处理 [10]，只保留A链并补全缺失的结构，去掉水分子，添加氢原子，对酸性和碱性氨基酸进行质子化处理。将蛋白导入AutoDock 4.2软件中，合并非极性氢，最后保存为.pdbqt格式的文件备用。

2.4. 小分子库的准备

以VD₃为结构基础在ZINC库中进行筛选，随后用OpenBabel软件将.sdf格式转换为.mol2格式 [11]。接着对小分子化合物进行预处理：ChemBio3D Ultra13.0中的calculations-MM2 minimize进行第一次能量最小化，最后再使用calculations-GAMESS interfere-minimize (Energy/Geometry)进行能量最小化，保存能量最小化结构的小分子，最后利用脚本将小分子生成.pdbqt格式文件进行下一步操作。

2.5. 分子对接

如图1所示，对接中心设置以VD₃为中心，向四周扩展5.0 Å的距离，盒子应尽量包括周围的活性残基，但又不能太大，才能保证对接结果的准确性。设置盒子x，y，z大小分别为60 Å × 48 Å × 66 Å，中心设置为3.837，−18.67，−34.066，Spacing设置为默认值0.375。使用脚本依次将小分子对接到蛋白的活性中心。输出结果以打分高低自动排序，从中选出能量和结构较好的小分子进行下一步分析。结合能越低，表明小分子与蛋白的亲和力越高。

Figure 1. Amino acids at 5.0 Å near VD₃

图1. VD₃附近5.0 Å的氨基酸

2.6. 分子动力学

经过以上步骤，挑选出结合能最佳的16个小分子进行分子动力学模拟。小分子在GaussView 5.0中优化 [12]，选用B3LYP/6-31G (d, p)来进行计算RESP电荷替代原有的电荷信息。

模拟过程是在AMBER14中进行的 [13]。使用tleap模块给蛋白质加载ff14SB力场 [14] [15]，antechamber模块给小分子加载gaff立场 [16]。采用显示溶剂化模型TIP3P [17]，溶剂盒子边缘离蛋白质最近的距离为10 Å。在体系添加24个钠离子使整个体系呈中性环境。最后生成蛋白质和小分子的拓扑文件和坐标文件。

分子动力学分为能量最小化，加热，平衡，和动力学模拟四个过程，氢原子键长采用SHAKE算法约束 [18]，用PME (partical mesh ewald)方法处理非键相互作用 [19]。能量最小化分两个部分来完成，各进行2500步的计算。第一次能量最小化采用1000步的最陡下降法，1500步的共轭梯度法对体系的水分子和离子进行能量最小化，约束蛋白质。第二次对整个体系进行了2500步的能量最小化，前1000步采用最陡下降法，后1500步使用共轭梯度法。

加热过程采用拉格朗日(Langevin)恒温器控制，体系温度将从0 k加热到300 k，用时60 ps。随后在恒温恒压(NPT)系综下进行50,000步，100 ps的缓慢升压和密度处理使压强达到一个大气压(1 atm)，再次对体系进行250,000步的模拟 [20]。最后为了消除边界效应采用周期性边界条件的方法，将平衡后的体系在300 k温度下，非键相互作用的截断半径为10 Å [19]，进行了100 ns无约束的动力学模拟。

使用Amber自带的程序计算蛋白与小分子复合物的RMSD。

2.7. 结合自由能计算

Hou等人研究表明在计算蛋白质–小分子复合物体系的结合自由能时，广义伯恩表面积法(MM/GBSA)方法比分子力学泊松–玻尔兹曼表面积法(Molecular Mechanics/Poisson-Boltzmann Surface Area, MM/PBSA)方法更加准确 [21]，因此本文采用MM/GBSA方法来探究小分子与Vdh的结合情况。对于小分子与Vdh形成的复合物结构，我们采用MM/GBSA计算结合自由能。结合自由能可以简单的由如下公式(1)表示：

公式(1)

$\begin{array}{c}\Delta {\text{G}}_{\text{bind}}=\Delta {\text{G}}_{\text{complex}}-\Delta {\text{G}}_{\text{receptor}}-\Delta {\text{G}}_{\text{ligand}}\\ =\Delta {\text{G}}_{\text{vdw}}+\Delta {\text{G}}_{\text{ele}}+\Delta {\text{G}}_{\text{GB}}+\Delta {\text{G}}_{\text{SA}}\end{array}$ (1)

其中，ΔG_complex代表复合物，ΔG_receptor代表受体，ΔG_ligand代表配体。每一项的总能(G)由分子力学(G_MM)和溶剂化效应(Gsol)两部分组成，即G = G_MM + G_sol，其中分子力学(G_MM)等于范德华作用能(E_vdw)和静电作用能(G_ele)的和(G_GMM = E_vdw + G_ele)，溶剂化效应(G_sol)等于极性溶剂化自由能(G_GB)和非极性溶剂化自由能(G_SA)的和(G_sol = G_GB + G_SA)。

2.8. 筛选结果验证

将16个化合物按照相似性分为4组，最后从中选出4个可以化合物进行催化验证。按照表1配置酶催化反应体系。在30℃下加入NADPH启动反应。48 h后，停止反应，用乙酸乙酯萃取反应混合液后旋转蒸发浓缩，用薄层色谱(TLC)分析法进行验证，展开剂为甲醇:二氯甲烷 = 1:10。

Table 1. Catalytic reaction system

表1. 催化反应体系

3. 结果与讨论

3.1. 初步筛选

用VD₃骨架在小分子化合库ZINC15中搜寻潜在的化合物，最终获得176个候选化合物。利用RDKits对176个小分子进行性质分析。目前已经上市的VD₃的分子量为384.63，logP为7.6，供受体均为1，环的数量为3个。分析可知，从Zinc数据库所选的化合物与VD₃有着较高的相似性，大部分含有三个六元环并一个五元环，且都含有一条侧链骨架，为类固醇类物质，此类小分子均存在发生羟化的活性位点的可能性。另外，结合图2及图3可以发现所选取的化合物的理化性质均与化合物VD₃相似，进一步说明Vdh可能具有羟化类固醇的潜力。

Figure 2. Statistical distribution of relative molecular masses of compounds

图2. 化合物的相对分子质量的统计分布

Figure 3. Statistical distribution of physicochemical properties of compounds

图3. 化合物的物理化学性质的统计分布

3.2. 对接结果

为进一步验证对接的准确性，我们将蛋白质晶体中的小分子VD₃作为配体，与处理好的蛋白受体对接，对接后两者构象重叠很好，RMSD = 0.538 < 2.0 Å，说明本实验中设置的对接参数合理，对接模型可靠(图4)。将小分子库与Vdh的活性位点进行对接，对接结合能越低，说明小分子与酶结合得越好。按照对接得分，对潜在的底物小分子进行排序，挑选对接能最低的16个小分子进行验证(表2)。这些小分子的对接结合能均小于Vdh与VD₃的对接结合能(−7.80 kcal/mol)，说明这16个小分子可能会与Vdh稳定结合。

Figure 4. The result of docking of Vdh with VD₃. The green small molecules in the figure indicate the conformation of VD₃ after docking; the blue color indicates the conformation of VD₃ in the original crystal

图4. Vdh与VD₃对接结果。图中绿色小分子表示对接后VD₃的构型；蓝色表示原晶体中VD₃的构型

Table 2. Docking energy and binding free energy of Vdh and small molecules

表2. Vdh与小分子的对接能和结合自由能

3.3. 分子动力学结果

经过两次能量最小化，加热，两次平衡和动力学模拟之后，计算小分子及VD₃与Vdh结合的结合自由能，结合自由能越小，表明受体与配体的亲和力越高，反之亦然 [22]。最终筛选出来16个小分子的结合自由能都与原底物VD₃相似，部分甚至比VD₃有更好的结合自由能打分(表2)。将16个化合物按照相似性分为4组，最后从中选出4个小分子进行进一步的分析。Vdh与VD₃的结合自由能为−42.5313 kcal/mol，挑选的小分子ZINC000248240127为−56.0776 kcal/mol，ZINC000252462409为−55.387 kcal/mol，ZINC000253497611为−51.2147 kcal/mol，ZINC000257349893为−48.3397 kcal/mol，结合自由能均低于Vdh，说明挑选出来小分子均有可能与Vdh发生反应，这说明最终筛选出的潜在底物分子有一定的可靠性。我们分析了四种小分子与Vdh的结合作用力。范德华力和表面溶剂作用有利于化合物小分子的结合，极性溶剂化作用不利于小分子与Vdh的结合(表3)。实验结果提示，在Vdh与小分子的体系中，范德华力起主导作用，可能是与活性中心周围的非极性氨基酸如M184，V283，P287，L387等形成的疏水环境有关。

Table 3. Calculation results of each energy term of MM/GBSA of the composite system of Vdh and small molecules

表3. 复合体系的MM/GBSA各能量项的计算结果

RMSD的计算可用于研究蛋白质–小分子复合物的变化和稳定性。从图5可以看出，四个小分子的RMSD基本保持稳定，前30 ns出现了明显波动，之后基本维持稳定。对于Vdh-VD₃体系，初始位置至30 ns的波动值为1.8 Å，ZINC000248240127初始位置至30ns的波动值为1.75 Å，ZINC000253497611初始位置至30 ns的波动值为2.5 Å，ZINC000257349893初始位置至30 ns的波动值为2.2 Å。其中，Vdh的收敛值为2.0 Å，ZINC000248240127为1.5 Å，ZINC000252462409为2.0 Å，ZINC000253497611为2.0 Å，ZINC000257349893为2.0 Å，说明ZINC000252462409与Vdh的结合更稳定，与动力学结果基本一致。四个小分子与Vdh的结合与Vdh-VD3相互作用的特征保持一致，说明四个小分子对蛋白整体构象的影响与VD3基本一致。

Figure 5. RMSD results in molecular dynamics. A: ZINC000248240127; B: ZINC000252462409; C: ZINC000253497611; D: ZINC000257349893

图5. 分子动力学中RMSD结果。A：ZINC000248240127；B：ZINC000252462409；C：ZINC000253497611；D：ZINC000257349893

3.4. Vdh蛋白表达及筛选结果验证

对含有重组表达载体的大肠杆菌BL21使用0.5 mM IPTG在37℃下诱导6个小时，在冰上破碎后用Mag-Beads His-Tag蛋白纯化磁珠进行纯化。使用SDS-PAGE鉴定目的蛋白。SDS-PAGE结果显示在目的蛋白大小处均有明显的条带(45 kDa) (图6)。

通过薄层色谱TLC对催化产物检测，由图7可知，四个化物反应标准品与反应物均不在同一个位置，说明发生了催化反应。其中结合最强的ZINC000248240127小分子的结果出现了两个产物点，说明可能发生了两个不同位点的羟化，结合力最差的ZINC000257349893产物中出现了底物，说明转化率不高，也进一步验证了结合自由能低。说明Vdh可以催化这些小分子，进一步证明利用虚拟筛选扩宽底物谱具有较高的准确性。

Figure 6. Vdh protein SDS-PAGE electrophoresis results

图6. Vdh蛋白SDS-PAGE电泳结果

(a)(b)(c)(d)

Figure 7. Thin layer chromatography of the results of Vdh-catalyzed small molecule reactions. (a): ZINC000257349893; (b): ZINC000252462409; (c): ZINC000253497611; (d): ZINC000248240127

图7. Vdh催化小分子反应结果的薄层层析。(a)：ZINC000257349893；(b)：ZINC000252462409；(c)：ZINC000253497611；(d)：ZINC000248240127

4. 讨论与总结

本文通过分子模拟对接技术与实验结合的方法从ZINIC数据库中筛选Vdh的底物谱。运用AutoDock 4.2分子对接程序包从ZINC库中筛选到176个潜在底物小分子。筛选的结果表明，小分子和VD₃的结构具有相似性，说明对接倾向于把具有相似结构的小分子确定为最佳配对。之后又结合了AMBER中的MM/GBSA程序计算了结合自由能进一步筛选。研究结果表明，范德华力是Vdh与小分子结合的主要作用力。实验结果显示从四个分组中挑选出4种化合物小分子都能与Vdh发生反应催化，证明了筛选出的底物库与Vdh有较好的结合效果，也说明了分子对接技术对于新底物的预测是可取的。结果表明，范德华力是Vdh与小分子结合的主要作用力。Vdh具有羟化类固醇物质的潜力，是分子对接技术操作简单、筛选潜在底物、发现新的可催化物质的一个好方法，现在已经有越来越多的研究者开始运用在蛋白质工程中了。我们的结果拓展了Vdh的底物谱，为Vdh更加广泛地应用于工业生产提供了理论基础。

NOTES

^*通讯作者。