AiPaper
论文状态:已完成

De Novo Design of Protein-Binding Peptides by Quantum Computing

发表:2025/03/07
原文链接PDF 下载
价格:0.10
价格:0.10
已有 4 人读过
本分析由 AI 生成,可能不完全准确,请以原文为准。

TL;DR 精炼摘要

本研究提出结合经典与量子计算的多尺度框架,利用量子退火器实现蛋白质结合肽段的从头设计,成功同时探索化学与构象空间。实验结果显示,设计的肽段化学多样且结合位点预测准确,体现量子计算在物理驱动药物设计中的潜力。

摘要

In silico de novo design can drastically cut the costs and time of drug development. In particular, a key advantage of bottom-up physics-based approaches is their independence from training datasets, unlike generative models. However, they require the simultaneous exploration of chemical and conformational space. In this study, we address this formidable challenge leveraging quantum annealers. Focusing on peptide de novo design, we introduce a multi-scale framework that integrates classical and quantum computing for atomically resolved predictions. We assess this scheme by designing binders for several protein targets. The D-Wave quantum annealer rapidly generates a chemically diverse set of binders with primary structures and binding poses that correlate well with experiments. These results demonstrate that, even in their current early stages, quantum technologies can already empower physics-based drug design.

思维导图

论文精读

中文精读

1. 论文基本信息

1.1. 标题

量子计算 De Novo 设计蛋白质结合肽段 (De Novo Design of Protein-Binding Peptides by Quantum Computing)

1.2. 作者

Lars Meuser, Alexandros Patsilinakos, Pietro Faccioli 所属机构:米兰比可卡大学 (University of Milano-Bicocca),Sa Biotech S.p.A.,意大利国家核物理研究院米兰比可卡分部 (INFN Sezione di Milano-Bicocca)。

1.3. 发表期刊/会议

该论文尚未在正式期刊或会议上发表,目前是发布在预印本服务器 arXiv 上的一篇预印本 (Preprint)。

1.4. 发表年份

2025年。

1.5. 摘要

In silico (计算机模拟) de novo (从头开始) 设计可以显著降低药物开发成本和时间。特别是,自下而上的 physics-based (基于物理) 方法的一个关键优势在于它们不依赖于训练数据集,这与 generative models (生成模型) 不同。然而,这些方法需要同时探索化学和构象空间。在本研究中,我们利用 quantum annealers (量子退火器) 解决了这一艰巨挑战。我们专注于肽段的 de novo 设计,引入了一个多尺度框架,该框架整合了经典计算和量子计算,用于原子级分辨率的预测。我们通过设计针对多种蛋白质靶点的结合剂来评估该方案。D-Wave 量子退火器快速生成了一组化学多样性丰富的结合剂,其一级结构和结合姿态与实验结果高度相关。这些结果表明,即使在当前早期阶段,量子技术也已经能够赋能基于物理的药物设计。

1.6. 原文链接

原文链接: https://arxiv.org/abs/2503.05458v1 PDF 链接: https://arxiv.org/pdf/2503.05458v1.pdf 发布状态:预印本 (Preprint)。

2. 整体概括

2.1. 研究背景与动机

论文试图解决的核心问题: 传统药物发现过程耗时且成本高昂,尤其是在寻找高亲和力、高选择性的药物分子方面。In silico de novo 设计作为一种有前景的替代方案,旨在从头开始生成新的药物分子。然而,现有的 de novo 设计方法面临两大挑战:

  1. 基于深度学习的生成模型 (DL-based generative models): 依赖于大量训练数据集,并且倾向于生成与训练集化学结构相似的分子,缺乏对化学空间的广泛探索,同时在生成高亲和力分子方面仍存在困难,且合成性不佳。
  2. 自下而上的基于物理的方法 (Bottom-up physics-based approaches): 虽然不依赖训练数据,但需要显式地考虑分子内和分子间的相互作用。最重要的是,这类方法需要同时探索巨大的化学空间(可能的分子组成)和构象空间(分子的三维结构),这是一个计算上极其复杂的优化问题。

为什么这个问题在当前领域是重要的: 药物发现的效率直接影响新药上市的速度和成本。解决 de novo 设计中的挑战,尤其是实现对化学和构象空间的有效探索,对于加速药物研发、发现全新作用机制的药物分子具有重大意义。肽段作为一种新兴的药物类别,具有毒性低、易合成等优点,但在计算设计上仍面临挑战。

现有研究存在哪些具体的挑战或空白:

  • 深度学习模型的数据依赖性和化学多样性不足。
  • 传统基于物理的方法在同时探索化学和构象空间时的巨大计算复杂性。
  • 如何有效整合新兴的量子计算技术与传统的物理模型,以解决药物设计中的优化难题。

这篇论文的切入点或创新思路: 论文提出利用 quantum annealers (量子退火器) 来解决基于物理的 de novo 肽段设计中同时探索化学和构象空间的难题。它构建了一个多尺度框架,将经典计算与量子计算相结合,以发挥两者的优势:量子计算用于粗粒度(简化)模型的优化,快速探索大规模搜索空间;经典计算用于精细化预测和原子级分辨率的结合姿态确定。这种混合方法旨在克服纯经典方法或纯量子方法在当前技术阶段的局限性。

2.2. 核心贡献/主要发现

论文最主要的贡献:

  1. 提出多尺度混合计算框架: 首次将经典计算与量子计算整合,用于基于物理的 de novo 肽段设计,实现原子级分辨率的预测。该框架利用量子退火器在粗粒度层面同时优化肽段的序列和构象,然后通过经典计算进行精细化。
  2. 量子编码 QUBO 问题: 将肽段设计问题严格编码为 Quadratic Unconstrained Binary Optimization (QUBO) 问题,使其可以在 D-Wave 量子退火器上求解。
  3. 实验验证框架的有效性: 通过对多个蛋白质靶点(包括3BRL、4DS1、3BFW和LC8蛋白)设计结合剂,并通过结构和序列两方面的统计分析,验证了该方法能够生成与实验结果高度相关的结合剂。
  4. 量子与经典优化器性能对比: 比较了 D-Wave 量子退火器与行业级经典优化器 Gurobi 以及传统模拟退火 (Simulated Annealing) 在解决 QUBO 问题上的性能。结果表明,D-Wave 在短时间内就能生成与 Gurobi 相似的最低能量值 (MEV),并且在保持低能量解的同时,能够产生更高化学多样性的结合剂。

论文得出了哪些关键的结论或发现:

  • 量子退火器能够有效应对 de novo 药物设计中同时探索化学和构象空间的挑战。
  • 即使在量子技术早期发展阶段,它也已能赋能基于物理的药物设计,提供与现有最先进经典方法竞争甚至互补的性能。
  • 所设计的肽段在结合姿态和一级结构上与实验数据表现出良好的相关性。
  • 量子优化器在生成化学多样性方面具有优势,这对于药物发现至关重要。

3. 预备知识与相关工作

3.1. 基础概念

  • In silico 设计 (In silico design): 指的是在计算机上进行模拟、建模和分析的设计过程。在药物发现中,它通常涉及虚拟筛选、分子对接、分子动力学模拟以及 de novo 设计等计算方法,旨在加速药物分子的发现和优化。
  • De novo 设计 (De novo design): 意为“从头开始设计”。与传统的虚拟筛选(从现有化合物库中寻找候选分子)不同,de novo 设计旨在根据目标靶点的特性,从基本的化学构件开始,逐步构建出全新的、具有预期药理活性的分子,而非仅仅优化现有分子。
  • Physics-based 方法 (Physics-based approaches): 基于物理原理和化学定律来描述分子间相互作用的方法。这些方法通常通过分子力场(如范德华力、静电力、氢键等)来计算分子构象和结合能量,不依赖于预先训练的数据集,因此在理论上具有更强的泛化能力和对新化学空间的探索能力。
  • Generative models (生成模型): 在机器学习中,生成模型旨在学习数据的内在分布,并能够生成与训练数据相似的新数据样本。在药物设计中,generative models 通常指基于深度学习的模型(如 GANVAETransformer 等),它们通过学习已知药物分子的结构和活性数据,来生成具有潜在药用价值的新分子结构。
  • Quantum annealers (量子退火器): 一种专门用于解决组合优化问题的量子计算机。它利用量子隧穿和叠加等效应,通过模拟物理退火过程,将系统逐渐演化到其基态(最低能量态),从而找到优化问题的近似解或精确解。D-Wave 是当前最著名的量子退火器供应商。
  • QUBO (Quadratic Unconstrained Binary Optimization) (二次无约束二值优化): 一种数学优化问题,其目标函数是二值变量的二次多项式,且没有额外的约束条件。许多 NP-hard 问题(如旅行商问题、图着色问题等)都可以被转化成 QUBO 形式,这使得它成为量子退火器和经典优化算法(如模拟退火)的理想输入格式。QUBO 问题的形式通常为 miniQiixi+i<jQijxixj\min \sum_i Q_{ii} x_i + \sum_{i<j} Q_{ij} x_i x_j,其中 xi{0,1}x_i \in \{0, 1\}
  • Ising model (伊辛模型): 统计物理学中的一个经典模型,用于描述磁性材料中自旋之间的相互作用。它与 QUBO 问题密切相关,因为通过简单的变量替换(si=2xi1s_i = 2x_i - 1,其中 si{1,1}s_i \in \{-1, 1\} 是伊辛自旋变量),QUBO 问题可以转化为伊辛模型问题。量子退火器通常通过寻找伊辛模型的基态来解决优化问题。
  • 肽段 (Peptides): 由短链氨基酸通过肽键连接而成的小分子。它们在体内发挥多种生物学功能,并作为一类重要的药物分子(peptide drugs)被开发,具有高特异性、低毒性、易于合成等优点,但也面临膜渗透性差、易被代谢等挑战。
  • 结合亲和力 (Binding affinity): 衡量两个分子(例如药物分子和靶蛋白)之间结合强度或倾向性的物理量。通常用结合常数 (KaK_a) 或解离常数 (KdK_d) 来表示。高结合亲和力意味着分子能够更紧密、更稳定地与靶点结合,这对于药物的有效性至关重要。
  • 粗粒化模型 (Coarse-grained model): 一种简化分子表示的方法,旨在降低计算复杂度。在粗粒化模型中,多个原子或分子基团被表示为一个“珠子”或一个“粗粒化位点”。例如,本文中将氨基酸表示为单个珠子,并将其分组为化学家族,从而减少了需要处理的变量数量,使得在量子退火器上进行模拟成为可能。
  • Lennard-Jones (LJ) potential (伦纳德-琼斯势): 描述原子或非极性分子之间非键相互作用的经验势函数,包括短程排斥和长程吸引。其形式通常为 V(r)=4ϵ[(σr)12(σr)6]V(r) = 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right],其中 rr 是距离,ϵ\epsilon 是势阱深度,σ\sigma 是相互作用为零时的距离。
  • Miyazawa-Jernigan knowledge-based potential (Miyazawa-Jernigan 知识势): 一种基于统计学的蛋白质相互作用势。它通过分析大量已知蛋白质结构中氨基酸对的接触频率,推导出氨基酸之间相互作用的倾向性。这种势函数反映了氨基酸在蛋白质折叠和结合中的偏好,常用于蛋白质结构预测和设计。
  • Autodock CrankPep (ADCP) (Autodock CrankPep): Autodock 软件的一个专门版本,设计用于肽段的分子对接。分子对接 (molecular docking) 是一种计算方法,用于预测配体(如肽段)与受体(如蛋白质)结合时,配体的最佳结合姿态 (binding pose) 和结合亲和力。
  • 精确率-召回率分析 (Precision-recall analysis): 评估二分类模型性能的常用方法,特别是在类别不平衡的数据集中。Precision (精确率) 表示被模型预测为正例的样本中,真正是正例的比例;Recall (召回率) 表示所有真正的正例中,被模型正确识别为正例的比例。Precision-recall curve 通过绘制不同阈值下的精确率和召回率来展示模型的性能。
  • AUC (Area Under the Curve) (曲线下面积): 通常指 ROC (Receiver Operating Characteristic) 曲线下的面积,或 Precision-recall curve 下的面积 (AUCPr)。AUC 值通常用于量化分类器性能的综合指标,值越接近1表示模型性能越好。在本文中,特指 AUCPr

3.2. 前人工作

  • 传统虚拟筛选 (Virtual Screening): 药物发现中常用的方法,通过计算机筛选大规模化合物库(多达 101010^{10} 个化合物)以识别潜在的“命中”候选分子。然而,筛选出的分子仍需进一步优化以提高亲和力、溶解度、可递送性和降低毒性,这个过程依然耗时且昂贵。
  • De Novo 设计算法:
    • 片段组装 (Fragment assembly): 将预先选定的分子片段组装成新的分子(例如 [2, 3])。
    • 随机突变和交叉 (Random mutations and crossovers): 通过对现有分子结构进行随机改变和组合来生成新分子(例如 [4, 5])。
  • 深度学习 (DL) 方法: 近年来,深度学习在 de novo 药物设计中取得了显著进展,使用了各种神经网络架构(如 GANVAETransformer 等)。然而,它们面临的主要挑战是:
    • 高亲和力分子生成困难 (Difficulty in generating high-affinity molecules): 尽管可以生成大量分子,但其中具有高结合亲和力的比例不高 [21]。
    • 合成性问题 (Synthesizability issues): 生成的分子可能难以或无法通过化学方法合成 [22]。
    • 数据依赖性和偏差 (Data dependency and bias): 大多数 DL 方案需要大量的目标特异性数据库进行训练,并倾向于生成与训练集化学邻域相似的分子。虽然主动学习 (active learning) 方案可以缓解这一问题 [6-8],但计算成本更高。
  • 量子计算在药物发现领域的应用: 随着量子硬件的快速发展,量子计算也被探索用于解决药物发现中的各种任务,包括:
    • 分子对接 (Molecular docking): [44-47]
    • 溶剂配置预测 (Solvent configuration prediction): [48]
    • 采样稀有构象蛋白质转变 (Sampling rare conformational protein transitions): [49]
    • 蛋白质折叠 (Protein folding): [50-54]
    • 蛋白质设计 (Protein design): [55-57]
    • 结合经典和量子神经网络识别小分子抑制剂 [58]。

3.3. 技术演进

药物设计领域的技术演进经历了从早期的经验法则和试错法,到基于结构的药物设计(SBDD)中的虚拟筛选,再到计算能力增强后的 de novo 设计。DL-based de novo 设计通过学习大数据模式,极大地扩展了化学空间探索的效率,但其固有的数据依赖性和潜在偏差限制了其在发现真正新颖分子方面的能力。与此并行,physics-based 方法始终追求从第一性原理出发,但受限于复杂的优化问题。本文的工作正是在这一背景下,将新兴的量子计算技术引入 physics-based de novo 设计,旨在克服纯经典方法在探索复杂化学和构象空间时的局限性,从而在不依赖大量训练数据的情况下,实现更高效、更具多样性的分子发现。

3.4. 差异化分析

本文的方法与相关工作中的主要方法相比,核心区别和创新点在于:

  • DL-based 方法的区别:
    • 独立于训练数据集: 本文的 physics-based 方法不需要大量的目标特异性训练数据集,避免了 DL 模型的数据依赖性、训练偏差和生成分子缺乏新颖性的问题。
    • 同时探索化学和构象空间: DL 模型通常需要分开处理序列和构象,或者在生成阶段就固定一些构象信息。而本文的方法利用量子退火器在粗粒度层面同时优化肽段的序列和其在结合口袋中的构象,这是其独特优势。
  • 与传统 physics-based 方法的区别:
    • 利用量子计算解决优化难题: 传统的 physics-based 方法虽然原理上可靠,但在同时探索化学和构象空间时会遇到巨大的计算复杂度,导致优化问题难以求解。本文通过将问题编码为 QUBO 并在量子退火器上求解,有效利用了量子隧穿和叠加等量子效应,以更高效的方式探索大规模组合搜索空间。
  • 多尺度混合框架: 本文不只依赖量子计算,而是提出了一个结合经典和量子计算的多尺度框架。量子计算用于粗粒度搜索以获得全局最优解,而经典计算(如 Autodock CrankPep)则用于将粗粒度解精细化为原子级分辨率的结合姿态,从而平衡了计算效率和预测精度。
  • 化学多样性: 实验结果表明,量子优化器倾向于生成具有连续能量谱的解,这暗示了其在生成具有相似低结合能但不同序列的分子方面具有优势,从而提高了设计结果的化学多样性,这对于后续的药物优化和开发至关重要。

4. 方法论

4.1. 方法原理

本文的核心思想是开发一个基于物理的多尺度框架,用于 de novo 肽段设计,该框架能够整合经典计算和量子计算的优势。其方法原理是:

  1. 统计物理学问题转化: 将寻找最优配体的问题形式化为统计物理学的最小化问题,即最小化配体与目标蛋白的结合能量,同时考虑配体与普通蛋白质表面的平均相互作用,以确保设计的选择性。
  2. 粗粒化简化: 为了在量子硬件上实现,首先对肽段的化学和构象空间进行粗粒化。氨基酸被分组为更少的化学家族,肽段的构象则被离散化为在晶格上的路径。相互作用也采用简化的知识势进行描述。
  3. 量子编码 QUBO 将粗粒化后的优化问题编码成 Quadratic Unconstrained Binary Optimization (QUBO) 问题,这是一种适合量子退火器求解的数学形式。QUBO 的低能量态对应于具有高结合亲和力的肽段序列和构象。
  4. 量子退火求解: 使用 D-Wave 量子退火器来快速探索巨大的化学和构象空间,找到 QUBO 问题的基态(或接近基态的解),从而获得优化的肽段序列和粗粒度结合姿态。
  5. 经典计算精细化: 量子计算得到的粗粒度结果(优化的序列和大致构象)随后通过经典的分子对接软件(如 Autodock CrankPep)进行精细化,以获得原子级分辨率的 off-lattice (脱晶格) 结合姿态。 通过这种混合方法,论文旨在利用量子计算在探索复杂高维空间中的独特优势,同时结合经典计算的精确建模能力,共同解决药物设计中长期存在的挑战。

4.2. 核心方法详解

4.2.1. 配体设计问题的统计物理学公式

从统计物理学的角度来看,为给定靶蛋白 PP 识别最佳配体的普遍问题可以被形式化为最大化以下表达式:

maxΓ,Σexp(U(Γ,Σ;P)kBT)PΓexp(U(Γ,Σ;P)kBT)(1) \operatorname* { m a x } _ { \Gamma , \Sigma } \frac { \exp \left( - \frac { U ( \Gamma , \Sigma ; P ) } { k _ { \mathrm { B } } T } \right) } { \sum _ { P ^ { \prime } } \sum _ { \Gamma ^ { \prime } } \exp \left( - \frac { U ( \Gamma ^ { \prime } , \Sigma ; P ^ { \prime } ) } { k _ { \mathrm { B } } T } \right) } \quad (1)

其中:

  • U(Γ,Σ;P)U ( \Gamma , \Sigma ; P ) 表示由化学组成为 Σ\Sigma、构象状态为 Γ\Gamma 的配体与处于天然构象的蛋白质 PP 组成的系统的相互作用能。

  • kBk _ { \mathrm { B } } 是玻尔兹曼常数。

  • TT 是温度。

  • 分母中的求和项 PΓ\sum _ { P ^ { \prime } } \sum _ { \Gamma ^ { \prime } } 遍历所有可能的蛋白质靶点 PP' 和配体 Σ\Sigma 的所有构象状态 Γ\Gamma'。这个分母确保了设计的配体能够选择性地最大化与给定靶点 PP 的亲和力。

    这个问题 (1) 可以等价地重新表述为一个最小化问题:

minΓ,Σ(U(Γ,Σ;P)G(Σ))(2) \operatorname* { m i n } _ { \boldsymbol { \Gamma } , \boldsymbol { \Sigma } } \left( U ( \boldsymbol { \Gamma } , \boldsymbol { \Sigma } ; P ) - G ( \boldsymbol { \Sigma } ) \right) \quad (2)

其中 G(Σ)G ( \Sigma ) 被解释为与给定化学结构 Σ\Sigma 相关的自由能:

G(Σ)kBTlnPΓexp(U(Γ,Σ;P)kBT)(3) G ( \Sigma ) \equiv - k _ { \mathrm { B } } T \ln \sum _ { P } \sum _ { \Gamma } \exp \left( - { \frac { U ( \Gamma , \Sigma ; P ) } { k _ { \mathrm { B } } T } } \right) \quad (3)

计算 G(Σ)G ( \Sigma ) 是一个极其艰巨的任务,因为它需要考虑所有可能的蛋白质靶点,并且对于每个靶点,还需要对所有配体构象状态的玻尔兹曼因子进行求和。为了降低计算成本,作者对 G(Σ)G ( \Sigma ) 进行了截断到最低阶的累积量展开 (cumulant expansion),得到简化后的优化条件:

minΓ,Σ(U(Γ,Σ;P)U(Σ)0)(4) \operatorname* { m i n } _ { \Gamma , \Sigma } \left( U ( \Gamma , \Sigma ; P ) - \langle U ( \Sigma ) \rangle _ { 0 } \right) \quad (4)

其中 U(Σ)0\langle U ( \Sigma ) \rangle _ { \mathrm { 0 } } 是配体 Σ\Sigma 与蛋白质的平均相互作用:

U(Σ)01NSPΓU(Σ,Γ;P)NsPΓ1 \begin{array} { r l r } { \langle U ( \Sigma ) \rangle _ { \mathrm { 0 } } } & { \equiv } & { \frac { 1 } { \mathcal { N } _ { S } } \sum _ { P } \sum _ { \Gamma } U ( \Sigma , \Gamma ; P ) } \end{array} \quad \text{和} \quad \begin{array} { r l } { \mathcal { N } _ { s } } & { { } \equiv } \end{array} \textstyle \sum _ { P } \sum _ { \Gamma } 1

根据条件 (4),最优配体是使其与给定靶点 PP 的结合能相对于其与蛋白质的平均相互作用能最小化的配体。这意味着,配体不仅要与目标蛋白结合紧密,还要尽量避免与非目标蛋白发生强结合,以提高选择性。

4.2.2. 粗粒化模型

为了将上述问题映射到量子退火器,作者将配体限制为肽段,并开发了一个粗粒化的数学表示。

化学和构象空间粗粒化:

  • 氨基酸表示: 氨基酸被表示为单个珠子 (single beads)。
  • 化学家族: 将氨基酸分组为 DD 个不同的化学家族,其中 D<20D < 20(自然界有20种常见氨基酸)。例如,D=2D=2 可以将氨基酸分为疏水和亲水两组。
  • 构象离散化: 肽段的构象空间通过一个填充靶蛋白口袋 PP 的方格晶格 (square lattice) 进行离散化。
  • 晶格间距: 晶格间距设置为与肽键长度 (0.38 nm) 匹配。
  • 肽段构象: 肽段在口袋中的每个构象状态都对应于晶格上的一个自避路径 (self-avoiding path)。
  • 蛋白质结构: 蛋白质的三维结构采用 off-lattice (脱晶格) 连续模型表示。

相互作用粗粒化: 为了推导相互作用能 U(Σ,Γ;P)U ( \Sigma , \Gamma ; P ) 的表达式,作者采用了 Miyazawa-Jernigan knowledge-based potential (Miyazawa-Jernigan 知识势),特别是 Kim and Hummer 提出的包含了不同氨基酸之间 Lennard-Jones (LJ) 对势相互作用的公式 [27]。 具体定义了两个 20×2020 \times 20 的参数矩阵,其条目由 i,j{120}i , j \in \{ 1 \dots 20 \} 索引:

  1. 能量矩阵 ϵ^\hat { \epsilon } εij=λ(eije0)(5) \varepsilon _ { i j } = \lambda ( e _ { i j } - e _ { 0 } ) \quad (5) 其中:

    • λ=0.159\lambda = 0 . 1 5 9 提供了一个整体尺度因子。
    • e _ { i j } 是原始 Miyazawa-Jernigan 矩阵 [26] 中的条目。
    • e0=2.27kBTe _ { 0 } = - 2 . 2 7 k _ { \mathrm { B } } T 是一个整体能量偏移。
  2. 相互作用范围矩阵 σ^\hat { \sigma } σij = σi+σj2(6) \sigma _ { i j } ~ = ~ \frac { \sigma _ { i } + \sigma _ { j } } { 2 } \quad (6) 其中 σi\sigma _ { i } 表示类型为 ii 的氨基酸的 van der Waals (vdW) 直径。

类型为 ii 的氨基酸和类型为 jj 的氨基酸在相对距离 rr 处的 Lennard-Jones (LJ) 相互作用由下式给出:

uij(r)={4εij[(σijr)12(σijr)6],εij<04εij[(σijr)12(σijr)6]+2εij,ϵij>0,r<rij04εij[(σijr)12(σijr)6],ϵij>0,rrij0.(7) u _ { i j } ( r ) = \left\{ \begin{array} { l l } { 4 | \varepsilon _ { i j } | \left[ \left( \frac { \sigma _ { i j } } { r } \right) ^ { 12 } - \left( \frac { \sigma _ { i j } } { r } \right) ^ { 6 } \right] , } & { \varepsilon _ { i j } < 0 } \\ { 4 \varepsilon _ { i j } \left[ \left( \frac { \sigma _ { i j } } { r } \right) ^ { 12 } - \left( \frac { \sigma _ { i j } } { r } \right) ^ { 6 } \right] + 2 \varepsilon _ { i j } , } & { \epsilon _ { i j } > 0 , r < r _ { i j } ^ { 0 } } \\ { - 4 \varepsilon _ { i j } \left[ \left( \frac { \sigma _ { i j } } { r } \right) ^ { 12 } - \left( \frac { \sigma _ { i j } } { r } \right) ^ { 6 } \right] , } & { \epsilon _ { i j } > 0 , r \geq r _ { i j } ^ { 0 } . } \end{array} \right. \quad (7)

其中 rij021/6σijr _ { i j } ^ { 0 } \equiv 2 ^ { 1 / 6 } \sigma _ { i j }。所有 LJ 相互作用的截止距离设置为 8.5A˚8.5 \mathrm { \AA }

精简化学字母表: 为了进一步粗粒化,氨基酸被分组为 D<20D < 20 个家族。这个分组过程通过寻找映射 a(i){1,,D}a ( i ) \in \{ 1 , \dots , D \} 来最小化以下损失函数:

L=i,j=120(eijea(i)a(j))2(8) L = \sum _ { i , j = 1 } ^ { 20 } ( e _ { i j } - e _ { a ( i ) a ( j ) } ^ { \prime } ) ^ { 2 } \quad (8)

其中,相应的有效(聚类后)相互作用矩阵的条目定义为聚类元素上的平均值:

eij=1Nk,l=120eklδa(k),iδa(l),j(9) e _ { i j } ^ { \prime } = \frac { 1 } { \mathcal { N } } \sum _ { k , l = 1 } ^ { 20 } e _ { k l } \delta _ { a ( k ) , i } \delta _ { a ( l ) , j } \quad (9)

并且 N\mathcal { N } 是归一化因子,对包含在两个集合中的元素数量进行归一化。同样,有效相互作用范围矩阵通过对聚类元素进行平均获得:

σI=1Nk=120δa(k),Iσk(10) \sigma _ { _ { I } } ^ { \prime } = \frac { 1 } { \mathcal { N } } \sum _ { k = 1 } ^ { 20 } \delta _ { a ( k ) , I } \sigma _ { k } \quad (10)

粗粒化模型中的设计优化条件: 为了解决设计问题 (4),需要评估肽段 Σ\Sigma 与蛋白质的平均相互作用 U(Σ)0\langle U ( \Sigma ) \rangle _ { 0 }。为了在有限的量子比特资源下估计这一项,作者引入了平均场近似 (mean-field approximation):

U(Σ)0Ncn=1lΣj=120fjεi(n)j(11) \langle U ( \Sigma ) \rangle _ { \scriptscriptstyle 0 } \approx { \mathcal N } _ { c } \sum _ { n = 1 } ^ { l _ { \Sigma } } \sum _ { j = 1 } ^ { 20 } f _ { j } \varepsilon _ { i ( n ) j } \quad (11)

其中:

  • lΣl _ { \Sigma } 是肽段长度。
  • i ( n ) 是链上位置 nn 处的氨基酸类型。
  • f _ { j } 是氨基酸类型 jj 在典型蛋白质表面上的相对频率。
  • Nc\mathcal { N } _ { c } 是肽段残基与典型蛋白质表面上的氨基酸形成的平均接触数。

4.2.3. 设计优化问题的量子编码 (QUBO)

为了使用 D-Wave 量子退火器解决肽段结合剂设计问题,条件 (4) 被编码为一个 Quadratic Unconstrained Binary Optimization (QUBO) 问题。这需要将有利的肽段序列和结合姿态映射到适当定义的量子哈密顿量 H^\hat { H } 的低能量态。

二值变量:

  • qi(k){0,1}q _ { i } ^ { ( k ) } \in \{ 0 , 1 \}:在每个晶格点 ii 处,如果该位置被类型为 k{1,,D}k \in \{ 1 , \ldots , D \} 的残基占据,则设置为1。
  • qij{0,1}q _ { i j } \in \{ 0 , 1 \}:表示相邻晶格点 iijj 之间的化学键的形成。
  • qij(k)q _ { i j } ^ { ( k ) }:辅助变量,最多是二值变量的二次项。

整体 QUBO 哈密顿量: 整体哈密顿量 HH 由几个项组成:

H=Hint+Hext+Hanc+Hocc+Hpath+Hchain(12) H = H _ { \mathrm { i n t } } + H _ { \mathrm { e x t } } + H _ { \mathrm { a n c } } + H _ { \mathrm { o c c } } + H _ { \mathrm { p a t h } } + H _ { \mathrm { c h a i n } } \quad (12)

各项解释如下:

  1. 外部相互作用项 HextH _ { \mathrm { e x t } } 代表肽段与口袋中蛋白质残基的相互作用。

    Hext=ik=1D(Ei(k)E0(k))qi(k)(13) H _ { \mathrm { e x t } } = \sum _ { i } ^ { \prime } \sum _ { k = 1 } ^ { D } ( E _ { i } ^ { ( k ) } - E _ { 0 } ^ { ( k ) } ) q _ { i } ^ { ( k ) } \quad (13) 其中:

    • i\sum _ { i } ^ { \prime } 表示对晶格点的求和。
    • Ei(k)E _ { i } ^ { ( k ) } 是类型为 kk 的孤立氨基酸在晶格点 ii 处与靶蛋白口袋中氨基酸相互作用的(预先计算的)能量。
    • E0(k)E _ { 0 } ^ { ( k ) } 根据方程 (11) 计算,表示该孤立氨基酸与一般蛋白质表面的平均相互作用,用于体现条件 (4) 中的选择性。
  2. 链内相互作用项 HintH _ { \mathrm { i n t } } 代表肽段内部的非键相互作用。

    Hint=i,jk,l=1Dukl(rij)(qi(k)qij(k))qj(l)(14) H _ { \mathrm { i n t } } = \sum _ { i , j } ^ { \prime } \sum _ { k , l = 1 } ^ { D } u _ { k l } \big ( r _ { i j } \big ) \Big ( q _ { i } ^ { ( k ) } - q _ { i j } ^ { ( k ) } \Big ) q _ { j } ^ { ( l ) } \quad (14) 其中:

    • r _ { i j } 表示晶格点 iijj 之间的欧几里得距离。
    • u _ { k l } ( r _ { i j } ) 是类型 kk 和类型 ll 氨基酸之间的 LJ 相互作用(见方程 (7))。
    • 因子 (qi(k)qij(k))( q _ { i } ^ { ( k ) } - q _ { i j } ^ { ( k ) } ) 确保排除共价键连接的氨基酸之间的相互作用。
  3. 辅助变量约束 HancH _ { \mathrm { a n c } } 确保辅助变量 qij(k)q _ { i j } ^ { ( k ) } 的定义一致性。它们仅在相邻的晶格点 iijj 上定义,并且当 qij=qi(k)=1q _ { i j } = q _ { i } ^ { ( k ) } = 1 时(即类型为 kk 的残基在 ii 点,且 i,j 成键),qij(k)q _ { i j } ^ { ( k ) } 被设置为1。

    Hanc=Ai,jk=1D(3qij(k)+qi(k)qij2qi(k)qij(k)2qijqij(k)),(15) \begin{array} { c } { { \displaystyle H _ { \mathrm { a n c } } = A \sum _ { \langle i , j \rangle } ^ { \prime } \sum _ { k = 1 } ^ { D } \left( 3 q _ { i j } ^ { ( k ) } + q _ { i } ^ { ( k ) } q _ { i j } \right. } } \\ { { \left. - 2 q _ { i } ^ { ( k ) } q _ { i j } ^ { ( k ) } - 2 q _ { i j } q _ { i j } ^ { ( k ) } \right) , } } \end{array} \quad (15) 其中:

    • i,j\sum _ { \langle i , j \rangle } ^ { \prime } 表示对相邻晶格点对的求和。
    • AA 是一个正的常数,设定了违反约束的能量惩罚。
  4. 占据约束 HoccH _ { \mathrm { o c c } } 确保每个晶格点最多被一个氨基酸占据。

    Hocc=AiklDqi(k)qi(l)(16) H _ { \mathrm { o c c } } = A \sum _ { i } ^ { ' } \sum _ { k \ne l } ^ { D } q _ { i } ^ { ( k ) } q _ { i } ^ { ( l ) } \quad (16) 此项通过惩罚同一晶格点被两种不同类型氨基酸占据的情况来强制实现这一约束。

  5. 路径约束 HpathH _ { \mathrm { p a t h } } 强制肽段的线性拓扑结构,即链内部残基与两个邻居键合,末端残基与一个邻居键合。

    Hpath = A(ht+hs+hr) hs = (1k=1Dqs(k))2+(k=1Dqs(k)j(s,j)qsj)2 ht = (1k=1Dqt(k))2+(k=1Dqt(k)j(t,j)qtj)2 hr = js,t(2k=1Dqr(k)j(s,j)qrj)2(1719) \begin{array} { r l } { { \displaystyle H _ { \mathrm { p a t h } } ~ = ~ { \cal A } \left( h _ { t } + h _ { s } + h _ { r } \right) } } & { } \\ { { ~ h _ { s } ~ = ~ - \left( 1 - \displaystyle \sum _ { k = 1 } ^ { D } q _ { s } ^ { ( k ) } \right) ^ { 2 } + \left( \displaystyle \sum _ { k = 1 } ^ { D } q _ { s } ^ { ( k ) } - \displaystyle \sum _ { j \in ( s , j ) } ^ { \prime } q _ { s j } \right) ^ { 2 } } } & { } \\ { { ~ h _ { t } ~ = ~ - \displaystyle \left( 1 - \displaystyle \sum _ { k = 1 } ^ { D } q _ { t } ^ { ( k ) } \right) ^ { 2 } + \left( \displaystyle \sum _ { k = 1 } ^ { D } q _ { t } ^ { ( k ) } - \displaystyle \sum _ { j \in ( t , j ) } ^ { \prime } q _ { t j } \right) ^ { 2 } } } & { } \\ { { ~ h _ { r } ~ = ~ \displaystyle \sum _ { j \ne s , t } ^ { \prime } \left( 2 \displaystyle \sum _ { k = 1 } ^ { D } q _ { r } ^ { ( k ) } - \displaystyle \sum _ { j \in ( s , j ) } ^ { \prime } q _ { r j } \right) ^ { 2 } } } \end{array} \quad (17-19) 其中:

    • sstt 表示肽段的起始和终止晶格点。
    • k=1Dqs(t)(k)\sum _ { k = 1 } ^ { D } q _ { s ( t ) } ^ { ( k ) } 表示 ss (或 tt) 点是否被任何氨基酸占据。
    • j(s,j)qsj\sum _ { j \in ( s , j ) } ^ { \prime } q _ { s j } 表示 ss 点与其邻居有多少个键。
    • 对于起始点 ss 和终止点 tt,表达式中的项确保这些点被一个氨基酸占据且只形成一个键(作为肽段的末端)。
    • 对于中间的残基 rr (js,t\sum _ { j \ne s , t } ^ { \prime }),表达式确保它们被一个氨基酸占据并且形成两个键(作为链的内部)。
    • A\mathcal{A} 是一个大的正惩罚常数。 为了使这些哈密顿量项作为硬约束 (hard constraints) 起作用,能量尺度 AA (以及 A\mathcal{A}) 必须相对于所有软相互作用足够大。
  6. 长度约束 HlengthH _ { \mathrm { l e n g t h } } 允许设定肽段的期望长度 L0L_0

    Hlength=w(L0i,jqij)2(20) H _ { \mathrm { l e n g t h } } = w \Big ( L _ { 0 } - \sum _ { \langle i , j \rangle } ^ { \prime } q _ { i j } \Big ) ^ { 2 } \quad (20) 其中:

    • ww 是权重因子。
    • i,jqij\sum _ { \langle i , j \rangle } ^ { \prime } q _ { i j } 是肽段实际形成的键的总数,即其长度。
    • 通过调整 ww,可以严格固定长度,或允许一定范围的长度波动。

QUBO到Ising模型的转换: 为了在 D-Wave 量子退火器上实现,需要将 QUBO 问题转换成 Ising model (伊辛模型) 的形式。这通过变量替换 σιz=2qι1\sigma _ { \iota } ^ { z } = 2 q _ { \iota } - 1 来实现,其中 ι\iota 遍历 QUBO 哈密顿量中的所有二值变量 qιq_\iota。这将二值变量 qι{0,1}q_\iota \in \{0,1\} 映射到自旋变量 σιz{1,1}\sigma_\iota^z \in \{-1,1\}。 生成的广义 Ising 哈密顿量包含二次项和线性项:

HIsing=lhlσlz+l>mJlmσlzσmz H _ { \mathrm { I s i n g } } = \sum _ { l } h _ { l } \sigma _ { l } ^ { z } + \sum _ { l > m } J _ { l m } \sigma _ { l } ^ { z } \sigma _ { m } ^ { z } 其中 hlh_l 是局部磁场系数,JlmJ_{lm} 是自旋之间的耦合系数。 经典哈密顿量通过将每个自旋变量替换为泡利 zz 算符 σ^ιz\hat { \sigma } _ { \iota } ^ { z },提升到量子 Ising 哈密顿量 H^Ising\hat { H } _ { \mathrm { I s i n g } }。量子计算机的量子比特被识别为 σ^lz\hat { \sigma } _ { l } ^ { z } 算符的本征态。

量子退火过程: 肽段设计问题被映射到寻找广义量子 Ising 哈密顿量 H^Ising\hat { H } _ { \mathrm { I s i n g } } 的基态。这通过绝热切换 (adiabatic switching) 过程实现:

  1. 初始化: 量子计算机的波函数在易于求解的哈密顿量 H^in\hat { H } _ { \mathrm { i n } } 的基态中初始化,例如 H^in=hlσ^lx\hat { H } _ { \mathrm { i n } } = - h \sum _ { l } \hat { \sigma } _ { l } ^ { x },其中 hh 是任意实常数。
  2. 绝热演化: 量子退火器的波函数根据时间相关的哈密顿量 H^(t)=a(t)H^in+b(t)H^Ising\hat { H } ( t ) = a ( t ) \hat { H } _ { \mathrm { i n } } + b ( t ) \hat { H } _ { \mathrm { I s i n g } } 演化一段时间 t _ { f }。调度函数 a ( t )b ( t ) 被定义为在时间间隔 t _ { f } 内,H^(t)\hat { H } ( t )HinH _ { \mathrm { i n } } 切换到 H^Ising\hat { H } _ { \mathrm { I s i n g } },即 a(0)b(0)a ( 0 ) \gg b ( 0 )a(tf)b(tf)a ( t _ { f } ) \ll b ( t _ { f } )
  3. 结果: 绝热定理保证,如果扫描过程足够慢 (tfhΔEt _ { f } \gg \frac { h } { _ { \Delta E } },其中 ΔE\Delta E 是最小能隙),系统将保持在其瞬时基态。最终测量量子比特的状态即可得到 QUBO 问题的解。 本文中使用的是 D-Wave 量子退火机的混合求解器,它结合了量子退火与经典的预处理和后处理步骤。

4.2.4. 肽段设计算法 (多步流程)

该算法遵循以下多步程序,如图1所示:

Figure 1

FIG. 1: 流程示意图。I:通过问题哈密顿量(公式 (12))的(量子)最小化,生成连接区域 A 和 B 的粗粒度肽段。配置被冻结,序列优化在原子级分辨率下进行。II:经典分子动力学模拟预测先前生成序列的脱晶格全原子表示。

  1. 粗粒度序列和构象优化 (Coarse-grained sequence and conformational optimization):

    • 使用第二节 C 节描述的二值编码和第二节 B 节定义的粗粒度模型。
    • 同时优化肽段的序列和其结合姿态。
    • 为满足现有量子计算硬件的限制,化学字母表被限制为 DD 个家族 (5D105 \lesssim D \lesssim 10)。
    • 此步骤利用量子退火器进行。
  2. 精细序列优化 (Refined sequence optimization):

    • 将上一步获得的氨基酸在晶格上的位置固定。
    • 这样可以释放更多的量子比特资源,用于对序列进行更精细的优化,此时可以使用所有 20 种天然氨基酸。
    • 固定构象后,可以移除辅助量子比特以及 HancH _ { \mathrm { a n c } }HpathH _ { \mathrm { p a t h } }HchainH _ { \mathrm { c h a i n } } 等项。
  3. 原子级结合姿态预测 (Atomistically detailed binding pose prediction):

    • 将选定的肽段序列 Σ\Sigma 输入到最先进的对接软件中。
    • 该软件(本文中使用 Autodock CrankPep (ADCP) [38])返回 off-lattice (脱晶格)、原子级详细的结合姿态。

5. 实验设置

5.1. 数据集

  • 结构验证数据集:
    • 使用了来自 LeadsPep 数据集 [41] 的三个不同的蛋白质-肽段复合物的 PDB 代码:3BRL, 4DS1, 3BFW。
    • 这些肽段分别包含 10 和 11 个氨基酸。
    • 实验设置:从 PDB 文件中移除肽段,并在相应的蛋白质口袋中设置一个晶格,将其放置在实验结合姿态中 Cα 原子 7.6 Å 范围内的区域。晶格的起始点 ss 和终止点 tt 被选择为最接近实验结合肽段末端的晶格点。
  • 序列验证数据集:
    • 使用了一个包含 111 个已知与 LC8 蛋白特定口袋结合的肽段的数据集 [42, 43]。LC8 是一种分子枢纽蛋白,参与细胞内稳态。
    • LC8 结合肽段的特征是存在一个 8 个氨基酸的识别基序。

5.2. 评估指标

对论文中出现的评估指标进行说明:

5.2.1. 原生接触分数 fnatf _ { \mathrm { n a t } } (Fraction of Native Contacts)

  • 概念定义: f_nat 是一种结构相似性指标,用于评估预测的结合姿态与实验参考姿态之间在关键相互作用上的吻合程度。它量化了预测姿态中与实验参考姿态共享的“原生接触”的比例。在这里,“原生接触”被假定为由 PDB 结构中的肽段与蛋白所形成的关键相互作用,这些相互作用被认为是所有良好结合剂的通用特征。高 f_nat 值表示预测的结合姿态在形成关键相互作用方面与实验参考姿态高度一致,从而暗示了良好的结合质量。
  • 数学公式: 论文中未直接给出 f_nat 的具体计算公式。但根据其描述“将 fnat>0.5f_nat > 0.5 的姿态标记为阳性”,可以推断 f_nat 衡量的是预测姿态中与实验参考姿态共享的接触数量占实验参考姿态总接触数量的比例。 假设 NnativeN_{native} 是实验结合姿态中形成的 native contacts 总数,而 Npredicted_nativeN_{predicted\_native} 是预测结合姿态中与实验 native contacts 重叠的接触数量。那么 f_nat 可以定义为: fnat=Npredicted_nativeNnative f_{\mathrm{nat}} = \frac{N_{\mathrm{predicted\_native}}}{N_{\mathrm{native}}} 这种定义方式在蛋白质-蛋白质相互作用预测和分子对接评估中是常见的。
  • 符号解释:
    • fnatf _ { \mathrm { n a t } }: 原生接触分数。
    • Npredicted_nativeN_{\mathrm{predicted\_native}}: 预测结合姿态中与实验参考姿态共享的原生接触数量。
    • NnativeN_{\mathrm{native}}: 实验参考结合姿态中形成的原生接触总数。

5.2.2. 精确率-召回率曲线下面积 AUC (Area Under the Precision-Recall Curve)

  • 概念定义: AUC 在本文中指的是精确率-召回率曲线下面积 (AUCPr)。这是一个衡量二分类模型性能的综合指标,特别适用于评估在类别不平衡数据集中(例如,只有少数“好的”结合姿态)的性能。精确率衡量预测为阳性的样本中有多少是真正的阳性,召回率衡量所有真正的阳性样本中有多少被正确识别。AUCPr 越高,表示模型在识别高 f_nat 的结合姿态时,能够同时保持高精确率和高召回率,即性能越好。
  • 数学公式: 精确率-召回率曲线下面积 (AUCPr) 通常通过计算曲线上所有点的精确率与召回率变化量的累加来估算。如果给定一系列按置信度排序的预测结果,可以得到一系列 (recall, precision) 点 (Rk,Pk)(R_k, P_k)AUCPr=k=1N(RkRk1)Pkinterp \mathrm{AUCPr} = \sum_{k=1}^{N} (R_k - R_{k-1}) \cdot P_k^{\text{interp}} 其中,RkR_kPkP_k 是第 kk 个阈值下的召回率和精确率,PkinterpP_k^{\text{interp}} 是在 RkR_k 处的插值精确率(通常取高于 RkR_k 的最大精确率以形成非递减曲线)。
    • 精确率 (Precision) 定义为:P=True PositivesTrue Positives+False PositivesP = \frac{\text{True Positives}}{\text{True Positives} + \text{False Positives}}
    • 召回率 (Recall) 定义为:R=True PositivesTrue Positives+False NegativesR = \frac{\text{True Positives}}{\text{True Positives} + \text{False Negatives}}
  • 符号解释:
    • AUCPr: 精确率-召回率曲线下面积。
    • True Positives (TP): 被正确识别为高 f_nat(阳性)的结合姿态。
    • False Positives (FP): 被错误识别为高 f_nat(阳性)的结合姿态。
    • False Negatives (FN): 真正的阳性(高 f_nat)但被错误识别为阴性(低 f_nat)的结合姿态。

5.2.3. 最低能量值 MEV (Minimum Energy Value)

  • 概念定义: MEV 是指在 QUBOIsing 模型优化问题中,所能达到的最低能量值。在本文中,它代表了设计算法所预测的肽段序列和结合姿态的结合能量。MEV 越低,通常意味着预测的肽段与蛋白质靶点的结合亲和力越强,也代表了优化器找到了一个“更好”的解。
  • 数学公式: MEV 对应于 QUBO 哈密顿量 HHIsing 哈密顿量 HIsingH_{Ising} 的最小值。 MEV=minqi{0,1}H(q1,,qN) \mathrm{MEV} = \min_{q_i \in \{0,1\}} H(q_1, \ldots, q_N) 或者在 Ising 模型中: MEV=minσi{1,1}HIsing(σ1,,σN) \mathrm{MEV} = \min_{\sigma_i \in \{-1,1\}} H_{\mathrm{Ising}}(\sigma_1, \ldots, \sigma_N)
  • 符号解释:
    • MEV: 最低能量值。
    • HH: QUBO 哈密顿量,包含各种相互作用和约束项。
    • qiq_i: QUBO 问题中的二值变量。
    • HIsingH_{\mathrm{Ising}}: Ising 哈密顿量。
    • σi\sigma_i: Ising 模型中的自旋变量。

5.2.4. 氨基酸家族频率分布 (Amino Acid Family Frequency Distribution)

  • 概念定义: 在序列验证中,为了比较设计肽段和实验已知结合肽段的序列相似性(而非精确匹配),作者将氨基酸聚类为 DD 个化学家族。然后,计算在结合基序的每个位置上,每个氨基酸家族出现的相对频率。通过比较设计肽段和实验肽段的频率分布,可以评估设计算法在捕捉氨基酸偏好性方面的准确性。
  • 数学公式: 假设一个结合基序有 LL 个位置,氨基酸被聚类成 DD 个家族。对于基序的每个位置 p{1,,L}p \in \{1, \ldots, L\} 和每个氨基酸家族 k{1,,D}k \in \{1, \ldots, D\},其频率 Fp,kF_{p,k} 可以计算为: Fp,k=在位置 p 处属于家族 k 的肽段数量所有肽段的总数量 F_{p,k} = \frac{\text{在位置 } p \text{ 处属于家族 } k \text{ 的肽段数量}}{\text{所有肽段的总数量}}
  • 符号解释:
    • Fp,kF_{p,k}: 在结合基序位置 pp 处,氨基酸家族 kk 的相对频率。
    • LL: 结合基序的长度。
    • DD: 氨基酸家族的数量。

5.3. 对比基线

论文将自己的方法与以下基线进行了比较:

  • 随机生成肽段: 用于评估算法性能的下限。通过随机生成 30 个肽段并进行对接,其 AUCPr 分数分布作为衡量设计算法有效性的基线。如果设计算法的性能优于随机肽段,则表明其具有一定的设计能力。
  • 重新对接实验原始肽段: 用于评估算法性能的上限。将 PDB 文件中原始的、已知结合的肽段重新进行对接模拟。其 AUCPr 分数代表了在当前对接软件(ADCP)和评估方法下所能达到的理想性能,作为设计算法的“理论上限”。
  • Gurobi (经典优化器): 一款行业级、最先进的商业优化求解器,广泛用于学术和工业研究。论文使用 Gurobi 求解相同的 QUBO 问题,以评估量子退火器在质量和效率上与经典优化器的竞争力。
  • 模拟退火 (Simulated Annealing, SA): 一种传统的启发式全局优化算法,常用于解决组合优化问题。论文使用 SA 进行长时间运行,以了解其在生成 MEV 谱和多样性方面的表现,并与 D-Wave 的结果进行对比。

6. 实验结果与分析

6.1. 核心结果分析

6.1.1. 结构验证

Figure 2

FIG. 2: PDB 条目 (a) 3BFW、(b) 4DS1 和 (c) 3BRL 对接模拟的精确率-召回率曲线下面积。每次模拟中,30 个随机肽段的得分(红色)与之前移除的肽段(橙色)、未考虑序列自由能的生成肽段(黑色)以及考虑了序列自由能的生成肽段(蓝色)进行比较。更高的得分表示在对接验证框架中(参见正文和 SM 节)更好地结合到肽段。3BRL 的重新对接姿态和生成分子的对接姿态显示在 (d) 中。

图2展示了针对三个蛋白质-肽段复合物(3BFW、4DS1、3BRL)的结构验证结果,通过精确率-召回率曲线下面积 (AUCPr) 来衡量。

  • 理想情况 (橙色线): 重新对接原始肽段(橙色垂直线)在所有三个案例中都获得了高 AUCPr 分数,这代表了理论上的最佳性能,因为它使用的是实际的结合分子。
  • 随机基线 (红色分布): 随机生成的肽段(红色分布)的 AUCPr 分数广泛分布,平均值接近0.4,表明随机肽段的结合质量很差。
  • 本文算法结果 (蓝色线): 使用本文算法设计的肽段(蓝色垂直线)的 AUCPr 分数显著高于随机生成肽段的平均值,并且更接近实验姿态。特别是对于蛋白质 3BRL (图2c),设计算法取得了非常好的结果,接近重新对接原始肽段的理想性能。
  • U(Σ)0U(Σ)0 项的重要性 (黑色线): 图中还展示了在设计过程中忽略平均相互作用项 U(Σ)0\langle U ( \Sigma ) \rangle _ { 0 }(黑色垂直线)的结果。这些结果显著差于随机生成的肽段,甚至对 3BFW (图2a) 和 4DS1 (图2b) 来说,性能几乎为零。这表明 U(Σ)0\langle U ( \Sigma ) \rangle _ { 0 } 项在确保设计的肽段具有选择性并匹配结合口袋的化学环境方面至关重要。如果缺少这一项,算法可能倾向于生成疏水性序列,而非特异性结合目标口袋。

6.1.2. 序列验证

Figure 3

FIG. 3: 顶部 50 个生成结合剂与 111 个实验已知 LC8 枢纽蛋白结合剂的比较。这些直方图对应于 LC8 锚定基序的八个位置。对于每个位置,显示了生成结合剂(蓝色)和 LC8 枢纽结合剂(橙色)的查询氨基酸类型。氨基酸已被聚类为第二节 B 部分所述的五个组。

图3展示了针对 LC8 枢纽蛋白的序列验证结果。作者将氨基酸聚类成 5 个家族,并比较了设计肽段(蓝色条)和实验已知结合肽段(橙色条)在 8 个结合基序位置上每个氨基酸家族的相对频率。

  • 整体相关性: 总体而言,设计肽段和实验肽段的氨基酸家族频率分布表现出良好的相关性。
  • 特定位置的强相关性: 在 8 个位置中的 6 个(即 1、2、4、6、7 和 8),实验数据集中出现频率最高的氨基酸家族,在设计数据集中也属于出现频率最高的两个家族之一。
  • 算法识别能力: 位置 2、6 和 8 尤其显示出数据集之间的强相关性,这表明设计代码能够识别在这些位置上促进结合所需的氨基酸类型。
  • 局限性讨论: 作者也指出,实验数据集可能不是无偏样本,不一定代表最优结合分子,但即便如此,这种比较仍提供了定性的评估。

6.1.3. 经典与量子优化对比

Figure 4

FIG. 4: 使用经典和量子优化为 PDB 条目 3BRL 生成结合剂所获得的 MEV。橙色点:1000 次 Gurobi 运行中获得的最低(圆圈)和平均(三角形)MEV,作为单个模拟运行时间的函数。水平蓝色线:300 次 D-Wave 运行中,5 秒混合退火时间获得的最低(虚线)和平均(点线)MEV。为了公平比较,当使用 Gurobi 进行经典最小化时,我们将 HancH _ { \mathrm { a n c } }HoccH _ { \mathrm { o c c } }HpathH _ { \mathrm { p a t h } } 所施加的条件编码为硬约束。

图4比较了 D-Wave 量子退火器(混合求解器)和经典优化器 Gurobi 在解决针对 3BRL 蛋白的 QUBO 问题时的 MEV (最低能量值)。

  • 最低 MEV 相似: 两种方法在达到最低 MEV 方面表现非常相似。

  • Gurobi 的收敛性: Gurobi 在运行约 10 秒后,其最低 MEV 不再显著改善(橙色圆圈趋于平稳)。

  • D-Wave 的效率: D-Wave 在 5 秒的混合退火时间内的 300 次运行,其平均 MEV(蓝色点线)低于 Gurobi 运行约 200 秒才能达到的平均 MEV(橙色三角形)。这表明 D-Wave 在较短时间内能够更有效地找到低能量的解。

    Figure 5

FIG. 5: 经典和量子优化中获得的 MEV 谱。左图:Gurobi 在不同最小化时间(行)和 D-Wave 混合求解器在 300 次 5 秒运行中获得的谱。右图:Gurobi 运行、经典模拟退火 (SA) 和 D-Wave 混合求解器获得的谱。

图5进一步比较了不同优化器生成的 MEV 谱(频率分布)。

  • 量子优化器的连续 MEV 谱和化学多样性 (左图和右图中的蓝色曲线): D-Wave 量子退火器生成了一个连续的 MEV 能量谱。作者通过直接检查发现,D-Wave 获得的 300 个 MEV 对应于 300 种不同的肽段序列。这表明量子优化器在找到低能量解的同时,能够提供丰富的化学多样性,这对于药物发现中寻找多种“命中”候选分子至关重要。

  • Gurobi 的离散 MEV 谱 (左图和右图中的橙色曲线): Gurobi 生成的 MEV 谱是离散的,表现为几个孤立的峰值。每个峰值通常对应于一个单一的设计序列。即使延长 Gurobi 的运行时间,也只是增加了这些已发现低 MEV 序列的相对出现频率,而没有发现更多的新序列。这表明 Gurobi 倾向于收敛到少数几个最优解。

  • 模拟退火 (SA) 的表现 (右图中的绿色曲线): 模拟退火也能产生连续的 MEV 谱,但其平均 MEV 显著高于 D-Wave 在仅 5 秒运行中获得的 MEV,尽管 SA 运行时间长达 1200 秒。这表明,虽然 SA 能够探索多样性,但在寻找低能量解方面的效率不如 D-Wave

    总结: 这些结果表明,即使在当前技术发展的早期阶段,量子退火器在肽段设计等组合优化问题上,已经能够与最先进的经典优化器竞争,并在提供高质量(低 MEV)解的同时,展现出在化学多样性方面的独特优势。

6.2. 数据呈现 (表格)

论文中没有提供具体的实验结果表格,所有实验数据和对比都是通过图表形式呈现和分析。

6.3. 消融实验/参数分析

  • 平均相互作用项 U(Σ)0\langle U ( \Sigma ) \rangle _ { 0 } 的重要性: 如图2所示,如果设计过程中忽略了 U(Σ)0\langle U ( \Sigma ) \rangle _ { 0 } 项,即仅最小化 U(Γ,Σ;P)U ( \Gamma , \Sigma ; P ),结果会显著变差(黑色垂直线),甚至低于随机生成的肽段。这表明该项对于确保设计的肽段具有选择性、避免与非目标蛋白结合以及匹配结合口袋的特定化学环境至关重要。缺乏此项会导致算法倾向于生成通用疏水序列,而非特异性结合剂。

  • 氨基酸聚类大小 DD 的影响: 在序列验证中,为了验证分析的稳健性,作者不仅使用了 D=5D=5 个氨基酸家族进行聚类分析(图3),还在补充材料中提供了 D=10D=10 个家族的类似分析(图S1)。

    Figure S1

    FIG.S1: 拟议的生成结合剂与 LC8 枢纽蛋白已知结合剂的比较,与图 3 相同。氨基酸已被聚类为 10 种氨基酸类型。

    即使在更精细的 D=10D=10 家族分辨率下,仍观察到设计的肽段序列与实验已知肽段序列之间存在整体的积极定性相关性。在 8 个位置中的 4 个(位置 1、2、7 和 8),实验数据集中出现频率最高的氨基酸类型仍然是设计数据集中出现频率最高的两个类型之一。这表明作者的 `coarse-grained` 模型在不同粒度下都能捕捉到蛋白质-肽段相互作用的关键特征。

7. 总结与思考

7.1. 结论总结

本文提出了一个创新的多尺度 de novo 肽段设计框架,该框架有效地整合了经典计算和量子计算的优势。通过将肽段设计问题严格地编码为 QUBO 问题,并利用 D-Wave 量子退火器进行粗粒度搜索,该方法能够同时且高效地探索肽段的化学序列和构象空间。随后,经典的分子对接软件对粗粒度结果进行精细化,生成原子级分辨率的结合姿态。

实验验证表明,所设计的肽段在结构上(与实验 PDB 结构的原生接触分数和 AUCPr 评估)和序列上(与已知结合肽段的氨基酸家族频率分布)均与实验结果高度相关。与经典优化器 Gurobi 和传统模拟退火的对比研究显示,D-Wave 量子退火器不仅能够快速找到与经典方法相媲美的低能量解,而且在生成化学多样性丰富的结合剂方面表现出独特的优势。这些发现强有力地证明了,即使在量子技术尚处于发展早期阶段,它也已能够成为 physics-based 药物设计的重要赋能工具。

7.2. 局限性与未来工作

作者指出的局限性:

  • 量子技术的早期阶段: 现有量子技术仍处于初期,硬件的规模和稳定性限制了当前的应用。随着系统规模的增大和分子表示的详细化,QUBO 哈密顿量形成的能量景观将类似于自旋玻璃,导致离散优化问题的复杂性呈指数级增长。
  • 经典计算资源: 虽然经典模拟是在桌面计算机上进行的,可以通过更强大的计算资源加速,但经典算法和硬件已高度优化,像 Gurobi 这样的最先进求解器从 GPU 加速中获益不大。
  • 定量评估的复杂性: 量子和经典优化器性能的定量评估并不直接,超出了本文的范围。优化器的效率会受到其内部参数(如退火时间、调度)和所用计算资源的影响。

未来可能的研究方向:

  • 直接实验验证结合亲和力: 未来的工作应包括对预测序列进行直接的实验验证,以评估其与靶蛋白的真实结合亲和力。
  • 推广到小分子设计: 将当前的肽段设计方案推广到更复杂的小分子设计。
  • 考虑受体柔性: 在优化过程中纳入受体蛋白质的柔性,因为蛋白质结合位点通常是动态的。
  • 整合 ADMET 性质: 在优化过程中考虑 ADMET (吸收、分布、代谢、排泄和毒性) 性质,以设计更具药物潜力的分子。
  • 扩展应用领域: 类似量子赋能的基于物理的设计方法可以扩展到药物发现之外的其他应用,例如有机半导体或分子纳米传感器的设计。

7.3. 个人启发与批判

个人启发: 这篇论文提供了一个非常清晰且有前景的路线图,展示了如何将新兴的量子计算技术融入到传统的计算化学和药物发现领域。最受启发的是其多尺度混合计算框架的设计哲学:

  1. 优势互补: 巧妙地将量子计算的全局搜索能力(在粗粒度空间探索复杂能量景观)与经典计算的精确计算能力(在原子级分辨率上精细化构象)结合起来,这在当前量子硬件仍不成熟的背景下是极其务实且高效的策略。
  2. 数据独立性: Physics-based 方法不依赖于大量的训练数据,这使得它在探索全新化学空间和处理数据稀缺的靶点时具有 DL-based 方法无法比拟的优势。结合量子计算,有望克服传统 physics-based 方法在计算复杂性上的瓶颈。
  3. 化学多样性: 量子优化器生成连续 MEV 谱的能力,意味着它能提供一组具有相似低能量但化学结构多样的候选分子。这对于药物发现至关重要,因为多样性增加了找到具有良好 ADMET 性质和可合成性的“命中”分子的概率。

批判:

  1. 缺乏直接的实验验证: 尽管论文通过 AUCPr 和序列分布进行了详尽的计算验证,但最终药物分子的价值在于其实验结合亲和力。论文明确指出未来工作需要进行直接实验验证,这在当前阶段仍是其主要局限。

  2. f_natAUCPr 公式未直接给出: 尽管这些是领域内常用指标,但为了论文的自洽性和对初学者的友好性,在方法部分或补充材料中提供其精确的数学定义会更好。我在本解析中补充了这些公式的常见定义。

  3. 粗粒化模型的精确性权衡: 粗粒化简化了计算,但也可能牺牲了某些重要的原子级细节,这可能影响最终预测的准确性。尽管通过经典对接进行了精细化,但粗粒度阶段的偏差是否会向下传递,仍需更深入的分析。

  4. 量子硬件的可扩展性: 论文提及了当系统规模增大时,复杂性呈指数级增长。当前 D-Wave 等量子退火器虽然拥有大量量子比特,但其连接性和噪声水平仍是挑战。未来的硬件改进将是该方法取得突破性进展的关键。

  5. QUBO 编码的复杂性: 将所有物理约束和相互作用编码为 QUBO 形式本身就是一个复杂且可能引入近似的过程。哈密顿量中的惩罚项 (A,A,wA, \mathcal{A}, w) 的选取也可能影响结果的质量。

    总的来说,这篇论文为 de novo 药物设计领域开辟了新的思路,将量子计算的潜力与成熟的物理模型结合起来,为未来更高效、更智能的药物发现提供了有力的工具。

相似论文推荐

基于向量语义检索推荐的相关论文。

暂时没有找到相似论文。