下（应用篇）| 量子计算加速蛋白质折叠

本文将延续上篇文章，通过应用VQE算法模拟解决蛋白质折叠问题的实验，解决使用传统方法耗时长、准确率低的问题，从而极大提升现代分子生物学的研究效率，为破解蛋白质折叠谜题带来新希望，进一步推动科学界前进。

1. 实施原理

将蛋白质折叠中氨基酸的数量与量子计算机中的量子比特数对应上，然后用这些量子比特的哈密顿量的本征态演化过程来模拟蛋白质折叠过程，这是一个类比过程，当量子从基态演变为一个稳定的本征态时，我们认为蛋白质折叠也达到了一个稳定状态，此时形成的蛋白质三维构型与自然形成的蛋白质三维构型就应该非常接近。为了验证构型，我们会按照以下实施流程图运行：

图一：实施流程图

为了解决经典计算机难以完成长链指数级增长的构象问题，我们考虑在量子框架中，用量子比特与氨基酸的数量进行线性扩展并相互对应，目的是通过量子的随机配置开始经过不断优化以降低能量直到最终确定蛋白质的最小能量构象，具体实施可以通过将蛋白质折叠问题编码为量子位算子并确保满足所有物理约束来实现。

对于问题编码，我们使用：

配置量子位：用于描述配置和不同磁珠相对位置的量子位

相互作用量子比特：编码不同氨基酸之间相互作用的量子比特

对于我们的实验案例，我们使用如下图所示的立方晶体结构，通过配置量子位对蛋白质折叠进行编码。

图二：立方晶体结构

一组量子位的系统哈密顿量为

H(Q)=Hgc(Qcf)+Hch(Qcf)+Hin(Qcf,Qin)

Q代表量子比特，Qcf是所需量子比特的总数，Qin是量子比特的寄存器，Hin(Qcf,Qin)代表纠缠态的下的能量项，Hgc是几何约束下的哈密顿量（控制一级序列的生长而无分岔的氨基酸初级序列的生长），Hch是手性约束下的哈密顿量（为系统强制执行正确的立体化学），Hin是系统的相互作用能量项。在这次实验案例中，我们只考虑最近的邻居相互作用。

2. 准备步骤

2.1 设定蛋白质主链

蛋白质由一条主链组成，该主链是氨基酸的线性链。对于不同残基的命名，我们使用氨基酸定义的单字母代码。在此实施案例中，证明了神经肽中量子比特算子的产生，我们分别采用了6个主链（WHLGEL），7个主链（AWHLGEL）和8个主链（AWHLGELV）代表的字母作为氨基酸的组成。其中每个字母代表一个氨基酸的英文缩写，如下图：

图三：氨基酸的英文缩写

2.2 侧链的生成

在蛋白质的主链之外，可能存在附着在主链残基上的氨基酸，我们称之为侧链。本实验方案中，侧链需要在物理约束下进行作用。

2.3 氨基酸之间的相互作用

对于蛋白质残基间接触相互作用的大小，我们使用了由Miyazawa.S.和Jernigan.R.L通过准化学近似推导的电位统计模型。除了这个模型之外，我们还允许运行随机接触交互，并引入了自定义的相互作用图谱，以增强蛋白质的某些构型（例如α螺旋，β片等），最终通过某个函数来实现测量接触相互作用大小的目的。

‍

2.4 物理约束

为了确保遵守所有物理约束，我们引入了三种惩罚函数：

（1） Chiral：用于施加正确手性的惩罚参数。

（2） Back：用于惩罚同一轴转弯的惩罚参数，此处用于消除连续两次选择同一轴的序列。通过这种方式，我们不允许链条折回自身。

（3） Punishment：用于惩罚最近邻接触内，磁珠之间的局部重叠的惩罚参数。

‍

2.5 肽定义

多个氨基酸脱水缩合形成的化学键叫肽。基于主链和可能存在的侧链，我们定义了肽对象。其中包括建模系统的所有结构信息，对象包括主链，侧链，氨基酸之间的相互作用以及三个物理约束。

‍

2.6 蛋白质折叠问题

基于定义的肽、相互作用和我们为模型定义的惩罚项，我们定义了返回量子位算子的蛋白质折叠问题，也就是返回最终的蛋白质折叠构型。为了使蛋白质折叠问题的解对应于哈密顿H（q）基态，我们首先准备了一个带参数的变分电路，包括配置寄存器。它由一个初始化块组成，该块具有Hadamard门和参数化的单量子比特RZ门，然后是一个纠缠块和另一组单量子比特旋转。见下图：

图四：参数化量子电路

其中θ=(θcf,θin)表示量子比特的角度集合，θcf代表所需量子比特的θ，θin是量子比特寄存器。与量子力学情况不一样的是，对于蛋白质折叠问题的解决方案，我们不需要估计哈密顿期望值，我们只需要对能量分布的低能量尾部进行采样。也就是说，最小化目标函数的条件是：Egrd≈O(θ*)，其中Egrd是蛋白质折叠损耗的最小能量，O(θ*)是迭代后所有量子处于本征态的能量，θ*是量子处于本征态的θ，如下图：

图五：公式

因此，角度θ的优化是使用变分量子本征求解器（VQE）算法执行。同时，我们用CVaR定义了一个目标函数，它基于由α值界定的分布尾部的平均值，见下图：

图六：函数图表

表示为CVaRα(θ)=<ψ(θ)∣H(q)∣ψ(θ)>α。其中ψ(θ)是与θ关联的波函数。通过将由U(θ)表示的参数化电路应用于任意起始状态ψ，通过改变参数θ的经典控制器迭代优化估计，最小化<ψ(θ)∣H(q)∣ψ(θ)的期望值。

该算法我们名为条件风险值（CVaR）VQE或简称为CVaR-VQE。与传统的VQE相比，CVaR-VQE极大地加速了哈密顿量的优化，如图一所示。同时，逻辑门参数的经典优化是使用差分进化优化器执行的，它模拟角度θ在希尔伯特空间中的自然选择，优化过程也总结在图一实施流程图中。

2.7 使用具有CVaR期望值的VQE来解决问题

实验方案中要解决的问题现在已经实现了所有的物理约束，并且有一个哈密顿量，方案中我们针对的是单比特字符串，它为我们提供了最小的能量。因此，我们可以使用具有风险条件值（CVaR）期望值的变分量子特征求解器来解决问题并找到最小配置能量。这里我们使用的是VQE作为优化部分。

3. 实验步骤

首先找到一个由13个氨基酸（DCAWHLGELVWCT）经过自然折叠而形成的蛋白质折叠构型rep8_β片，见下图：

其对应的ribbon图是可以看清折叠位置的，如下图所示：

考虑到实际能够使用的量子比特数及模拟器的运行耗时，我们从要模拟的13个氨基酸中挑出带有折叠部分的6个、7个和8个氨基酸来作为我们的算法验证。

3.1 步骤

首先我们把量子位的系统哈密顿量H(Q)=Hgc(Qcf)+Hch(Qcf)+Hin(Qcf,Qin)转换为实际可编程的哈密顿量H（q）。

Hy是相关系数，qi=（1-σz）/2 ,σz是泡利矩阵,Yi ∈ {0, 1},Nt是项的总数。

系数Hy是肽珠间的一种相互作用，这些作用是通过强制执行蛋白质肽对象的惩罚参数来对蛋白质折叠进行约束，相互相互作用的大小是通过调用蛋白质折叠的化学函数来执行得到的，此处我们只要设定好惩罚参数就可以得到不同的系数。系数hy基本决定了哈密顿量H(q)是如何找到最低能量的。计算结果如下图（列出了部分数值）：

图七：系数Hy与蛋白质折叠构型

其中I代表单位矩阵，Z代表泡利矩阵。每6个I和Z组合而成的系列就是一种蛋白质折叠构型。

3.2 对哈密顿量求解

我们用CVaR-VQE对每一个肽珠间的哈密顿量进行求解，得到哈密顿量的值（肽珠序列一开始时被设定为WHLGEL）。见下图（只列出了部分值）

图八：折叠构型与哈密顿量

当两个肽珠之间占据相邻位点距离L<1，即两个量子比特同为0或者同为1时，物理相互作用表现为吸引；当两个肽珠之间占据相邻位点距离L>1时，即两个量子比特一个同为0，另一个为1，物理相互作用表现为排斥，排斥意味着折叠。具体折叠方式取决于最小哈密顿量对应的蛋白质折叠构型。

图八中，每个哈密顿量乘上图六中对应的蛋白质折叠构型系数Hy就是该蛋白质折叠的最小能量哈密顿量，也就是蛋白质折叠的最小能量Egrd。我们分别把6个氨基酸（WHLGEL），7个氨基酸（WHLGELG）和8个氨基酸（WHLGELGE）形成蛋白质折叠的输出结果对比如下图：

从上图的Ribbon图可以看出，其中折叠的位置基本保持发生在LG之间，与13个氨基酸自然形成的折叠构型非常接近，说明算法拟合的构型与自然折叠的构型相符，实验结果是成功的。最后让我们来看一下CVaR-VQE梯度下降的效果：

在20次迭代后目标函数逐渐收敛，说明我们可以利用这一方法来解决蛋白质折叠这一难题，并且速度更快。

4. 另一种方法QAOA算法

我们还可以使用QAOA算法来实现上面的VQE算法功能，现在将两种方式结果进行一个对比（代码省略）：

7个氨基酸WHLGELG算法和QAOA算法计算比较

8个氨基酸WHLGELGE的VQE算法和QAOA算法比较

从效果来看VQE算法和QAOA算法比较接近，这至少证明QAOA算法是VQE算法的一种替代方法。

本次的实验说明了量子计算减少了蛋白质折叠的检索空间，使NP问题减低为O(n)的问题。同时量子计算加速了蛋白质折叠问题的研究，提升了蛋白质折叠问题的效率，之前使用其他的方法解决蛋白质折叠问题时，即便一个典型的蛋白质也有10^300种可能的构型，使用VQE算法后，蛋白质折叠问题得到了指数级的简化。

比如在本次案例中，如果采用传统经典计算，7个蛋白质有20^7=1280000000种折叠方式，剔除强制惩罚项后的折叠方式也还需要 1.8*10^5次迭代,，而本案采用CVaR-VQE算法后我们只用了50次迭代就接近目标函数了，效率提升肉眼可见。