关于HP模型蛋白质折叠问题的序参量

《Journal of Computational Science》：Towards order parameters for the HP-model protein folding problem

　　

• H-ratio序参量用于识别最大计算难度中的峰值。
• 将GEV分布拟合到每个H-ratio的最大值，以估计其临界值。
• 最难实例需要大量回溯；中位数难度的实例可能逐步折叠。
• 已求解140524个长度为10至30的独特HP模型蛋白，并以开源形式提供。

**引言**

蛋白是地球上生命的根本构成单元。其分子组成的结构造就了一个极其困难的优化问题：蛋白质结构预测，也称为蛋白质折叠。本文将对该问题的离散化版本进行实验研究。

蛋白是由氨基酸组成的长链。对于天然存在的蛋白，这些链会折叠成称为天然态的三维结构。这源于疏水效应和分子间相互作用，包括氢键、离子键和范德华力[1]、[2]。在稳定的天然态构象中，蛋白执行各种重要功能，如催化酶促反应、调控基因表达，或为细胞结构提供支架[3]。另一方面，不稳定或错误折叠的蛋白是导致多种严重疾病的主要原因，如阿尔茨海默病、帕金森病、亨廷顿病或囊性纤维化，这些统称为"蛋白质错误折叠疾病"[4]、[5]。目前对抗错误折叠疾病仍然困难，但可以通过多种方式进行。例如通过所谓的"构象校正剂"或阻断错误折叠诱导的聚集[6]。了解是什么导致蛋白以稳定或不稳定的方式折叠，将有助于更好地理解错误折叠疾病，从而提供替代干预策略。

寻找最稳定构象是一个具有挑战性的连续优化问题，其空间维度与蛋白长度成正比[3]。Folding@home[7]和Rosetta@home[8]项目旨在通过利用分布式志愿者计算来克服真实蛋白的巨大计算需求挑战。然而，这些模拟的高度自由度仍然是其可扩展性的制约因素。虽然天然蛋白的折叠时间随长度呈指数增长，但折叠模拟的计算时间大致呈二次方增长。

进展虽然持续但缓慢，需要大量计算资源[9]、[10]。值得注意的是，从头蛋白质结构预测问题仍然是计算困难的，即使在通过HP模型[11]大幅降低自由度后也是如此。尽管HP模型蛋白折叠存在一些生物学局限性，如处理冷变性[12]或无法模拟热量法两态协同性[13]，但它仍作为计算困难的优化问题被广泛研究。该问题一个有趣的特点是需要在优化解决方案时同时考虑可行性，这与旅行商问题（TSP）[14]和划分问题（PP）[15]等其他问题不同。此外，能量景观粗糙，局部极小值之间由高能壁垒分隔[16]。这些特征加上简化的约束集合，使得HP模型蛋白折叠成为此类实证优化研究的理想问题。

HP模型是动态真实世界蛋白质折叠的静态组合优化抽象，因此大幅降低了复杂性[11]。它将蛋白表示简化为疏水（）或极性（）的氨基酸，并通过仅允许在二维或三维笛卡尔网格的相邻单元格中放置来离散化空间。因此，由??个氨基酸组成的蛋白的二维构象可以通过一个包含???1个移动步骤（↑、↓、←、→）的序列来描述，表示链中下一个氨基酸相对于上一个应放置的方向。在本研究中，折叠是相对于二维网格进行编码的。前一步移动的反向移动总是导致无效构象，因此从选项中移除。这使每个氨基酸放置只剩下3个可选移动。每个未在链上相连但相邻放置的一对疏水氨基酸之间存在"接触"[1]，使构象的自由能降低1。因此，更低的能量值对应更稳定的蛋白。由氨基酸链??=(???,???,…,????)组成的折叠链的自由能由公式(1)给出，其中??是网格上的曼哈顿距离[参见19]。一个在糖尿病治疗中有应用的短蛋白是Exendin-4，发现于毒蜥的唾液中[20, Section 2.3.3]。其衍生迷你蛋白Trp-Cage TC5b（图1左），氨基酸链为NLYIQWLKDGGPSSGRPPPS[21]、[22]，在HP模型中可表示为字符串[23]。TC5b可实现的自由能难以先验估计，但存在一些约束：任何疏水氨基酸在二维笛卡尔网格上最多只能有2+1个接触[2]，且一对疏水氨基酸之间要有接触可能，它们在链上的距离必须为3+2???，其中??∈?。

一种（精确的）深度优先分支定界（DFS-BnB）算法表明，TC5b在二维网格上的最优构象为←←↑↑←←←↓→→↓←↓→↓→↑→→（图1右），其稳定性为?8。这是在总共2.63?8???2个可能构象中的一个，其中前两个氨基酸固定以简单消除沿一个轴线的对称性。2.638的增长因子源于自避免行走的枚举，其展现相对于??的渐近指数扩展率[24]。因此，短于10个氨基酸的蛋白可能具有略小的有效状态空间。TC5b大约有3.8×10?个二维构象。虽然这种构象空间已经极其庞大，但与天然存在的蛋白相比仍远为逊色，后者不仅体积大得多，而且存在于连续状态空间中。

因此，即使使用HP模型的简化，蛋白质结构预测仍然是计算上不可行的。找到超过某些目标稳定性的构象是一个????-完全问题[25]、[26]，而找到最稳定构象是????-困难问题[27]。此优化变体可以被表述为一个不可逆的、确定性的、完全信息[28]的序贯决策过程，具有单一决策者。它等价于一个有限时域马尔可夫决策过程（MDP）[29]。不可逆性源于部分构象形成自避免行走[26]、[30]，因此早期移动会限制后续可能的移动。在这个意义上，该问题在结构上类似于国际象棋或围棋的博弈树，但没有对手。HP模型蛋白质折叠问题也可以建模为用于深度强化学习的MDP[31]。这种方法为长达50个氨基酸的链实现了已知最优解[32]。然而，它不能保证最优性，且深度学习的泛化能力受训练数据多样性的限制[33]。因此，需要精确方法来研究该问题的难度。

难以处理的????-完全和????-困难问题已被研究数十年，之前的研究表明，单个????-完全和????-困难问题实例的难度实际上具有一定可预测性[34]、[35]、[36]、[37]、[38]。研究人员研究了各种各样的问题，如??-着色问题（K-COL）[35]、??-可满足性问题（K-SAT）[35]、[39]、[40]、哈密顿圈问题（HC）[35]、[41]、完美矩形装箱问题（PRPP）[42]、二次分配问题（QAP-SAT）[43]、旅行商问题（TSP）[35]、[44]、非对称旅行商问题（ATSP）[36]和项目调度问题（PSP）[37]。一个共同点是，最难实例在一定程度上可以通过称为序参量的数据分析来识别。序参量描述了最大计算难度中的一个峰值，对于????-完全问题，这与可满足性的是/否相变相吻合。

本研究的目标是通过提出寻找????-困难HP模型蛋白折叠问题序参量的深入分析来扩展这一知识。这反过来可以用于开发新算法或算法比较的基准。在实际应用中，具体实例的难度比问题的理论最坏情况复杂度更为重要。因此，可以利用序参量，通过逐实例算法选择[45]来，例如，精确求解简单实例并近似求解困难实例。

**相关工作**

长期以来已经确立，????中许多"困难"问题的实例实际上很容易求解[34]。复杂度类别基于最坏情况实例，而实际难度在不同问题实例之间差异巨大。1991年，Cheeseman等人研究了可满足性问题的[3]

**方法论**

我们提议区分一级和二级序参量。一级序参量识别????-困难问题的最大计算难度峰值。二级序参量进一步聚焦，通过在一级序参量的一个区间内区分实例难度。在本研究中，三种不同特征被测试作为一级和二级序参量。它们的适用性有助于阐明HP模型[模型]计算难度的来源

**结果**

首先讨论关于-ratio的结果。Van Eck和Van den Berg于2023年的先前研究发现，probit函数是描述使用DFS-BnB算法折叠数据集时计算难度范围的良好描述符[19]。对??_random和[两个]数据集的初步调查显示，其实例求解所需的递归次数呈现出类似趋势（图3）。然而，将probit函数拟合到[数据上]……

**讨论**

据我们所知，本研究首次量化了HP模型蛋白质折叠问题计算难度的来源。它表明，控制序列中疏水氨基酸的比例（-ratio）会影响从数据集中采样实例时的难度分布。完全随机数据集的难度分布遵循probit曲线，而控制-ratio则使其大幅偏离。

更有趣的是……

**未来工作**

寻找HP模型蛋白质折叠的序参量及其临界值估计是探索实例难度来源的重要发现。然而，它们并不直接解释来源。-ratio仅描述了生成这些实例所需的条件，因此，有必要进一步研究其他能更详细说明难度来源的特征。我们尝试了……

**CRediT作者贡献声明**

Okke van Eck：原稿撰写、可视化、验证、软件、资源、项目管理、方法论、调查研究、形式分析、数据管理、概念构思。Anton Feenstra：审阅与编辑、原稿撰写、监督。Reitze Jansen：审阅与编辑。Kristian Verduin：审阅与编辑。Mario Acosta：审阅与编辑、获取资金。Daan van den Berg：审阅与编辑、监督。Ruben Horn：

**利益冲突声明**

作者声明他们没有已知的竞争性经济利益或个人关系，这些利益或关系可能被认为会影响本论文中报告的工作。

**致谢**

本研究使用了荷兰国家电子基础设施，由SURF合作社支持，资助编号为EINF-12999。Okke van Eck衷心感谢巴塞罗那超级计算中心在本研究最终阶段提供的支持与自由。最后，我们感谢ChatGPT在构思过程中以及在缺少人类同事的情况下对合理性的非正式评估方面的有用性。

Okke van Eck | Anton Feenstra | Reitze Jansen | Kristian Verduin | Mario Acosta | Daan van den Berg | Ruben Horn