基于TIES_PM分子动力学模拟的结核病耐药性从头预测：快速、精准与可重复性新策略

　　本综述重点介绍了TIES_PM（Thermodynamic Integration with Enhanced Sampling for Protein Mutations）这一基于分子动力学（MD）模拟与自由能计算（ΔΔG）的创新方法，其在预测结核分枝杆菌（Mycobacterium tuberculosis）对利福平（RIF）耐药性方面展现出卓越性能。文章系统评估了61种RNA聚合酶（RNAP）突变，结果与WHO临床分类高度一致，并揭示部分突变存在非典型耐药机制。该方法具备快速（单突变仅需约5小时）、低成本、可扩展等优势，为当前结核病诊断流程提供了有力补充，并对耐药性筛查、药物研发及全球抗微生物耐药性（AMR）监测具有重要推动意义。

　　

ABSTRACT

结核病作为全球最致命的传染病之一，每年导致数百万新发病例和死亡。耐药结核病，特别是对利福平等一线药物的耐药性，对全球健康构成严峻挑战。本研究应用基于集成学习的分子动力学计算机模拟方法TIES_PM，通过自由能计算估算结合亲和力，预测RNA聚合酶中的利福平耐药性。通过分析61种突变（包括利福平耐药决定区内的突变），TIES_PM得出可靠结果，与临床参考数据高度一致，并识别出表明替代耐药机制的异常数据点。未来，TIES_PM能够识别和选择耐药进化风险较低的先导化合物，对于较小蛋白质，它可能通过分析所有可能的密码子排列来系统预测抗生素耐药性。此外，其灵活性允许将预测扩展到其他一线药物和耐药性疾病。TIES_PM为当前诊断流程提供了快速、准确、低成本和可扩展的补充，特别是在研究和临床领域的耐药性筛查方面。

IMPORTANCE

抗微生物耐药性（AMR）作为全球性威胁，挑战结核病（TB）的早期诊断和治疗。本研究采用自由能计算方法TIES_PM，通过量化细菌RNA聚合酶（RNAP）突变如何影响利福平（RIF）结合，有效预测AMR。在模拟61种临床观察到的突变后，结果与WHO分类一致，并揭示出模糊案例，表明存在替代耐药机制。每个突变仅需约5小时，提供快速、经济高效的预测。集成方法确保了统计稳健性。TIES_PM可扩展到较小蛋白质进行系统密码子排列分析，实现全面抗生素耐药性预测，或适用于识别低耐药风险的药物先导化合物。它还可应用于其他TB药物和耐药病原体，支持个性化治疗和全球AMR监测。这项工作为完善耐药突变数据库和表型分类标准提供了新工具，增强早期诊断，同时推进转化研究和传染病控制。

INTRODUCTION

抗微生物耐药性（AMR）被列为全球公共卫生和发展的十大威胁之一。2019年，细菌AMR直接导致全球127万人死亡，并在495万死亡案例中起到推动作用。除了对死亡和残疾的影响，AMR还带来巨大的经济负担。根据世界银行估计，到2050年，AMR可能导致额外的1万亿美元医疗费用，到2030年，全球年度国内生产总值（GDP）损失在1万亿至3.4万亿美元之间。

AMR发生在细菌、病毒、真菌和寄生虫对 antimicrobial 药物无反应时。因此，它威胁到现代医学的许多进展。肺炎、结核病、败血症、淋病和食源性疾病等感染变得越来越具有挑战性，在某些情况下几乎无法有效治疗。其中，结核病（TB）是全球关注的问题。2022年，TB导致约130万人死亡，包括16.7万HIV感染者。在全球范围内，TB是第二大传染性杀手，仅次于COVID-19，领先于HIV/AIDS。

由需氧细胞内专性病原体结核分枝杆菌（Mycobacterium tuberculosis）引起，TB是一种 endemic 细菌感染，通过空气传播，在全球范围内构成持续的公共卫生挑战。除了肺部，TB还可以传播到其他器官，如胸膜、淋巴结和中枢神经系统。

耐药结核病（DR-TB）进一步加剧了问题。DR-TB源于TB药物的不当使用，如错误处方、劣质药物或患者过早停止治疗。它占全球AMR死亡的13%，由新发展的耐药性和人与人传播驱动，突出了其严重的公共卫生影响。多药耐药结核病（MDR-TB）是DR-TB的一种形式，由抵抗至少异烟肼（INH）和利福平（RIF）的细菌引起，这两种是最有效的一线TB药物。MDR-TB需要超过18个月的治疗方案，并需要二线药物，这些药物通常效果较差、毒性更大且更昂贵。此外，由抵抗最有效二线TB药物的细菌引起的TB使患者治疗选择非常有限。

TB耐药性基于Mtb基因组（结核分枝杆菌基因组）中的遗传突变，这些突变影响药物靶点或激活酶功能。此类机制包括外排泵、药物修饰和细胞包膜的不渗透性。耐药TB的诊断和指导适当治疗依赖于使用分子诊断和基因组测序检测此类突变。当前对TB耐药性的研究基本侧重于阐明这一现象的分子机制、发现新药物以及重新利用现有药物以对抗耐药菌株。它还意味着开发更快、更可靠的诊断技术，能够尽早识别耐药TB菌株，以有效控制疾病传播。

通过 universal 药物敏感性测试和系统筛查接触者及高风险人群进行早期诊断，对TB的有效治疗和管理至关重要。2022年，有750万新诊断的TB病例，这是WHO自1995年开始监测以来的最高记录，而估计病例为1060万。估计TB发病率与诊断之间的差距表明，在改进诊断工具方面还需要进一步工作。在估计的41万例多药耐药TB病例中，仅175,650例接受了治疗，更加严重地突出了诊断不足。

TB诊断的主要方法包括显微镜检查、细菌培养和分子测试。Ziehl–Neelsen（ZN）显微镜检查是最常见的技术，特别是在资源有限地区，但其敏感性在22%至80%之间，变化很大，取决于细菌负荷，使其对细菌计数低的患者效果较差。Lowenstein–Jensen等细菌培养方法高度特异和敏感，敏感性在80%至85%之间，但需要4-6周才能提供结果，这与早期诊断的要求相冲突。在分子测试中，GeneXpert MTB/RIF具有89%的敏感性和99%的特异性，是一种快速的TB/DR-TB测试（少于2小时）。然而，它昂贵且无法区分活菌和非活菌。这些限制要求进一步研究抗TB药物耐药性的分子机制和诊断技术的发展。

传统诊断方法的局限性，特别是在AMR方面，需要更敏感的方法。更复杂的耐药形式，特别是在TB等疾病中，需要先进的诊断功能和更具体的耐药机制分子方面信息。正是在这一领域，计算方法可能在AMR研究中取得重大进展。

在AMR研究中，由于分子生物学、生物信息学、化学信息学以及高性能计算（HPC）能力的进步，计算方法取得了显著进展。超越蛋白质与配体相互作用的初始阶段，使用分子对接，该方法已扩展到分子动力学（MD）模拟，能够模拟如TB相关蛋白质等复杂分子系统中发生的耐药机制。基于物理（PB）的计算方法的发展里程碑首先包括1980年代早期引入分子对接用于预测药物结合，随后1990年代MD模拟研究分子时间依赖性行为。2000年代，引入了自由能计算方法以更精确地预测结合亲和力。最近，HPC与这些方法的集成使得能够进行大规模模拟以生成更准确的预测，特别是使用基于集成的方法。AutoDock、Glide、GROMACS和AMBER等工具已广泛用于药物发现和分子模拟，而机器学习（ML）的最新进展也有助于使用大型数据集预测耐药突变。最近的 efforts 已经前瞻性地评估了使用盲数据的耐药突变，使计算预测更接近临床应用。

结合亲和力提供了计算途径以支持临床决策。突变通常通过降低抗生素与其靶点之间的结合亲和力导致耐药性。为了量化这一点，我们可以计算野生型和突变型蛋白质之间的结合自由能差异（ΔΔG），使用相对结合自由能（RBFE）方法。最常用的RBFE方法是一种化学方法，其中野生型氨基酸在模拟中逐渐且非物理地变为突变型，并在每一步测量能量差异。这些计算进而告知突变如何影响药物结合，从而预测可能的耐药性。

RBFE在药物先导优化中广泛应用，以预测药物微小变化对其与靶点结合的影响，达到约1 kcal/mol的准确性。研究表明，RBFE可有效预测各种疾病中的耐药突变，包括与金黄色葡萄球菌和嗜血杆菌中甲氧苄啶耐药相关的小蛋白质中的那些。尽管RBFE在小型蛋白质系统中预测耐药性非常有效，但由于准确性和计算限制，将其应用于更大的大分子复合物仍然具有挑战性。然而，随着详细的大分子结构3D模型的不断增加以及HPC的进步，使得药物耐药性预测更加有希望，正如本文报道的工作所证明的那样。

通过充分利用MD模拟结合相对结合自由能（RBFE）计算，我们能够以高准确性预测RNA聚合酶（RNAP）蛋白中的RIF耐药突变，该蛋白由rpoB基因编码；RIF通过结合到RNAP的β亚基，阻止RNA延伸（图3显示了它们的位置）。

在本研究中，我们应用了一种称为TIES_PM（Thermodynamic Integration with Enhanced Sampling for Protein Mutations）的基于集成的方法来预测RNAP-RIF系统中的耐药性。这项工作的创新之处不在于计算框架本身，而在于其大规模应用（RNAP-RIF系统包括超过44万个原子）、与 curated 表型数据的集成以及跨61种临床观察到的突变的系统评估。论文结构如下：在第2节中，我们详细介绍了方法。在第3节中，我们展示了耐药性预测结果，包括许多落在利福平耐药决定区（RRDR）内或外的突变。在第4节中，我们讨论了第3节的结果，强调了该方法对潜在生物医学和临床应用的更广泛影响。

MATERIALS AND METHODS

本节介绍了自由能计算的基础理论、核心计算方法（TIES_PM）、蛋白质和配体准备、模拟过程以及模拟生成数据的后处理。

Free energy theory and TIES_PM

药物与其靶蛋白之间的结合亲和力是药物设计和耐药性预测中的重要属性，表明药物与靶蛋白相互作用的能力。主要思想是配体结合由吉布斯自由能?G的变化驱动。?G越负，蛋白质与药物的复合物越稳定，表明结合越强。因此，吉布斯自由能的变化计算为：

ΔG_binding = (G_A - G_B)

其中ΔG_binding是吉布斯自由能的变化，G_A是复合物的吉布斯自由能，G_B是蛋白质和药物的吉布斯自由能之和。

对于耐药性预测，如果突变发生后?G变得不那么负，表明药物与突变靶蛋白之间的亲和力降低，则表明发生耐药性。相反，自由能降低（即?G更负）对应于结合亲和力增加，表明突变仍然敏感。

TIES_PM（Thermodynamic Integration with Enhanced Sampling for Protein Mutations）是一种开源计算方法，旨在计算由于蛋白质突变引起的结合自由能变化。这个标准TIES方法的扩展版本在模拟中纳入了蛋白质突变与单一抑制剂，从而能够详细分析特定氨基酸替换如何改变蛋白质-配体相互作用的特性。

TIES_PM的基本思想是通过改变问题来节省计算资源。最初，我们想比较突变前后的结合亲和力。因此，直接解决方案是计算原始蛋白质-配体对ΔG_binding^wt和突变蛋白质-配体对ΔG_binding^mut之间的差异。然而，当我们分解每个能量项并进行等效变换时，我们可以将公式改为，其中两个蛋白质-配体对（B-A和B′-A′）之间的差异变为过程A-A′和B-B′之间的差异。后者被称为化学过程，因为这些中间步骤是非物理的，尽管最终状态是真实的——最终状态之间的差异给出ΔΔG，因为自由能与路径无关。

ΔΔG = ΔG_binding^mut - ΔG_binding^wt = (G_A′ - G_B′) - (G_A - G_B) = (G_A′ - G_A) - (G_B′ - G_B)

ΔΔG = ΔG_com - ΔG_prot

当ΔΔG < 0时，突变蛋白表现出更好的亲和力，意味着敏感；相反，当ΔΔG > 0时，原始蛋白表现出更好的亲和力，意味着耐药。由于化学过程的系统变化远小于蛋白质-配体对的变化，相对结合自由能（RBFE）计算将减少采样空间，节省计算资源，并提高计算的准确性。

为了在原子水平上化学地将蛋白质从野生型变为突变状态，在TIES_PM中定义了一些中间状态或窗口，介于两个端点（野生型和突变型）之间，如图2A所示。每个窗口都必须在自由能差异计算过程中进行模拟。在我们的工作中，设置了13个窗口，每个窗口是两种物理状态在不同比例下的混合物，受突变参数λ的程度控制。野生和突变蛋白分别由λ=0和λ=1状态表示，中间混合状态由λ ∈ (0,1)表示。

吉布斯自由能的变化使用热力学积分（TI）确定，基于以下方程：

ΔG = ∫₀¹ (?G(λ)/?λ) dλ = ∫₀¹ ??V(λ,x)/?λ?_λ dλ

其中λ是突变程度，V(λ,x)是势能，??V(λ,x)/?λ?_λ是状态λ中势能导数的系综平均。

在每个这些窗口进行五个副本的模拟，以便可以采样广泛的构象。这给出了一个集成（图2B），其中预测通过捕获通常被一次性模拟遗漏的构象多样性而更加精确。这种集成方法至关重要，因为分子动力学本质上是混沌的——单次模拟无法可靠地捕获这种多样性或确保可重复性。通过使用多个副本，我们控制错误并有效提高预测准确性。TIES_PM将通过计算配体结合和脱辅基蛋白状态之间的自由能差异，精确量化这些突变如何影响药物结合和激活。均值和标准误（SE）使用 bootstrap 重采样（B = 10000）估计，而95%置信区间根据Student’s t分布计算，为X? ± t?SE，其中t_{0.975, 4} = 2.776，适用于n = 5个副本（即自由度=4）。这种方法与经典TIES框架一致，该框架已在多个靶点（包括TYK2）上经过全面基准测试，表现出强大的准确性。TIES-PM已成功应用于多种系统中的蛋白质突变，如FGFR3激酶和雌激素受体，其中收敛分析确认4 ns生产运行对大多数情况足够。

System preparation and simulation

Protein and ligand preparation

准备蛋白质和突变的计算模型需要多个步骤。

首先，我们使用结核分枝杆菌RNA聚合酶全酶与利福平复合物的晶体结构（PDB ID：5UH6）作为起始模型。与晶体结构常见的情况一样，从PDB文件检索的模型可能包含缺失的环或原子。为了补全缺失部分，使用AlphaFold2预测的结构用于补全缺失片段，基于与天然X射线结构的共同选定区域进行比对。

其次，从具有多个构象体的蛋白质中去除冗余构象体，以保持结构的一致性。

接下来，运行AmberTools中的“reduce”以添加氨基酸的正确质子化状态，从而正确表示蛋白质的生物学正确状态。对于像RIF这样的药物分子，进行键校正和添加氢原子，以便为MD模拟做好准备。然后，基于标准力场进行参数化，如用于蛋白质的ff14SB和用于药物分子的GAFF2。使用双拓扑RBFE引入突变。野生型和突变型侧链共存，它们与环境的相互作用相应缩放，以再现野生型和突变型状态之间的逐渐过渡。然后通过溶剂化结构、添加离子和应用周期性边界条件来设置系统，控制模拟盒大小、缓冲距离、离子键和原子数等参数。所有这些都使用Amber的tleap工具完成，以准备必要的输入文件。然后在提交到超级计算机运行模拟之前，彻底检查系统的正确性。

Simulation process

在生产运行之前进行能量最小化和1 ns平衡，生产运行涉及每个副本4 ns的MD模拟，集成中使用五个副本，使用NAMD。轨迹和能量导数每10 ps保存一次以供进一步分析。模拟在Argonne国家实验室的Polaris超级计算机上执行，该计算机配备AMD EPYC处理器和NVIDIA A100 GPU，峰值容量为44 petaflops。我们系统的模拟速率约为24 ns/天，平衡阶段的挂钟时间为1.3小时，模拟阶段为一个副本4.1小时。

此外，将TIES_PM模拟与临床数据进行比较以验证预测。例如，来自WHO耐药目录的数据，可作为CSV文件获取，包括突变分类和一些关联的临床证据，在EVIDENCE列中格式化为JSON。这些数据在突变的临床意义方面提供信息，并有助于将突变分类为耐药（R）或敏感（S），以最好地代表临床现实，特别是考虑到缺乏RNAP-RIF突变的实验ΔΔG值。包括这些证据对于确保TIES_PM模拟最好地代表临床现实非常重要，即使承认定量值，例如最小抑制浓度或滴定热法，最适合用于药物耐药性和敏感性确定。

RRDR and the choice of mutations

结核分枝杆菌中的利福平耐药决定区（RRDR）通常跨越rpoB基因的428-452位氨基酸，大多数临床观察到的利福平（RIF）耐药突变发生在这里。这些残基形成药物结合口袋的一部分，位于RNA聚合酶（RNAP）的β亚基上；因此，这里的突变通常破坏RIF的抑制作用。

然而，并非RRDR中的所有非同义突变都 confer 耐药性（例如，rpoB L443F）；一些 confer 耐药性的突变（例如，rpoB I491F, V170F）位于其外。因此，尽管RRDR是RIF耐药性检测的关键靶点，更全面的突变筛查仍然重要。在我们的研究中，我们选择了61种突变，大约一半在RRDR内，一半在外（图3）。

Resistance classification and data validation

Clinical threshold

结合自由能的 positive 变化（ΔΔG > 0）意味着突变后抗生素与其靶点的结合效果较差， thus 表明耐药性。临床上，如果样品的最小抑制浓度（MIC）超过临界值，通常为流行病学截断值（ECOFF/ECV）——野生型样品的MIC的99百分位数，则被视为“耐药”。遵循从酶抑制动力学 adapted 的热力学框架，耐药突变的MIC分布可以转换为ΔΔG阈值。对于RNAP–RIF系统，根据CRyPTIC Consortium的ECOFF/ECV值和耐药变体的几何平均MIC，这 yield 预期的ΔΔG约为~1.2 kcal/mol。

Classification of SOLOs and phenotype ratio

一个关键概念是将突变定义为SOLO（单次出现，单独观察）事件——突变在样品中 exclusively 发生而没有共存变体的实例。通过比较样品中SOLO突变的R/S频率，分配它们的耐药效应分类。例如，C138R突变总是产生耐药性，因此可以可靠地分类为R，而像I133T这样的突变产生可变的R/S结果，因此它们的耐药分类未知。

为了更好地描述参考临床数据的原始结构，我们定义“表型比率”为R或S频率除以总SOLOs。这提供了一个基础，通过比较样品中的R/S结果来发现突变的临床重要性。一些突变，如D435G，其SOLOs具有50% R和50% S分布， yet 是耐药的，可能因为它们位于RRDR内且没有明确的表型数据。因此，临床数据需要基于其数据结构进行批判性解释。

Statistical assessment of consistency between data sets

为了将我们的计算与以前的研究进行比较，我们使用了配对t检验和F检验。t检验发现均值差异，F检验发现方差。我们还测试了模型性能的敏感性（真阳性率，正确预测的耐药突变比例）和特异性（真阴性率，正确预测的敏感突变比例）。此外，使用F1分数（精确率和召回率的调和平均值）来评估模型在正确检测耐药突变的同时避免错误分类敏感突变方面的表现。

RESULTS

在本节中，我们首先对七种重要突变进行了计算。其中三种是临床耐药的，而四种是敏感的。我们还与Fowler等人以前的研究结果进行了比较。基于这些预测的信心，总共计算了61种突变，增加了研究结果的说服力。

Comparison with the previous computational study

第一组七种突变是根据它们的位置和对利福平耐药性的影响选择的。S450L位于rpoB的RRDR内。V170F和I491F，虽然位于RRDR外，但靠近S450L突变和抗生素结合位点，用于 confer 耐药性；I491F以具有可变MIC为特征。阴性对照是L443F，它在RRDR内且靠近结合位点但不产生耐药性，S388L和T585A，它们离结合位点较远，出现在临床样本中但不 confer 耐药性。S428C，一种RRDR突变，没有耐药关联且对侧链影响最小，被选为人工阴性对照，因为它未在临床样本中观察到，并且其侧链朝向远离药物。预测结果与以前的工作进行了比较，如图4和表1所示。

所有四种阴性对照（S388L, S428C, L443F, T585A）被预测对RIF的有效性没有影响。两种 confer 耐药性的突变（V170F和S450L）不仅显示 positive ΔΔG值，而且超过ECOFF/ECV衍生的临床阈值，确认它们 confer 耐药性的事实。有争议的突变I491F被预测为敏感，与临床实验结果相反（如表1所示）。然而，重要的是要注意，I491F的临床评估表型比率仅为55.5%。临床数据的可靠性显著影响预测准确性的评估，如下一节所示。

此外，配对t检验 yield t = ?0.672，P值为0.527，表明我们的结果与以前研究之间的平均ΔΔG值没有显著差异。类似地，F检验产生F = 0.204，P值为0.660，表明两个数据集的方差在统计上也相似。这些结果 confirm 我们的计算值与以前的研究一致，支持TIES_PM方法的稳健性。

值得注意的是，我们的预测表现出更高水平的分类信心，因为与以前的研究相比，跨越耐药阈值的偏差更少且幅度更小。这表明TIES_PM提供了更精细和可靠的耐药性和敏感性评估，减少了边界案例的模糊性。

Predictions among the most common mutations

为了全面评估我们的TIES_PM方法，我们预测了总共61种突变的耐药性，整合了那些具有既定临床分类以及最初标记为未知类型的突变。选择突变以涵盖关键的利福平耐药相关残基，覆盖 diverse 表型类别（耐药、敏感、未分类）、热力学稳定性范围（ΔΔG: -1.6 至 18.2 kcal mol?1）和临床表型比率（0%–100%），从而严格评估TIES_PM方法在机械模糊和统计不确定场景中的稳健性。

如图5所示，我们将突变分为两组：组A包括WHO已提供明确分类为耐药（R）或敏感（S）的突变，而组B由最初由于 insufficient 或冲突的临床证据而分类为未知的突变组成。这种区分使我们能够评估我们的方法在 well-characterized 和模糊案例中的表现。详细的数值在表S1中提供。

从图5A中，十九种突变被分类为R，九种为S。 among 这些，我们的模型正确预测了15种耐药突变和所有九种敏感突变， even 考虑误差条。敏感性和特异性分别为78.95%和100%，而F1分数达到88.24%，反映了性能的强平衡。这些结果证明了我们的模型在准确识别耐药突变的同时最小化假阳性的可靠性，使其成为耐药性预测的强大工具。

在四种被错误分类或 inaccurate 预测的耐药突变（H445C, I491F, L452P, S441Q）中，两种表现出低临床表型比率：I491F，先前在 preceding 部分讨论过，在173次SOLO测试中发现96次耐药， yield 表型比率仅为55.5%；类似地，L452P在259次SOLO测试中临床鉴定为168次耐药，表型比率仅为64.9%。相比之下，H445C和S441Q表现出高于90%的正常表型比率。S441Q的平均ΔΔG超过决策阈值，但其误差范围相对较大。增加副本数量有望提高准确性，从而能够正确预测耐药性。至于H445C，差异的可能解释将在讨论部分进一步讨论。

比较图5A和B，在图5B中观察到误差条的显著增加，这对应于最初由WHO分类为未知的突变。这种分类源于 insufficient 或冲突的临床数据，而我们计算中较大的误差条进一步反映了这些突变的复杂性。因此，对于预测中误差条跨越耐药阈值两侧的情况，我们将其视为“未知类型”，标签为“U”。这并不罕见，特别是当估计值接近零时。这 indeed 是一些突变被分类为“Unknown”的原因之一。它们对药物敏感性的影响最小，导致临床观察可能走向任一方向。因此，在计算错误指示的不确定性方面，我们的模型预测与这些突变固有的分类挑战之间存在一致性。为了进一步调查这一点，我们评估了λ点重叠、副本间一致性和ΔΔG收敛性。大多数系统显示稳定的自由能估计和足够的采样，符合Zhang等人提出的标准。这些发现 confirm 观察到的不确定性不是由于采样缺陷，而是反映了突变固有的结构或物理化学复杂性，例如电荷变化、庞大侧链或局部重排。

此外，基于其位置分类为耐药的突变——在图5B中以浅红色突出显示——尽管误差较大，但表现出更高的平均ΔΔG，与它们预测的耐药趋势一致。

总体