1. 蛋白质折叠是什么

蛋白质折叠（Protein folding）是蛋白质获得其功能性结构和构象的物理过程。通过这一物理过程，蛋白质从无规则卷曲折叠成特定的功能性三维结构。在从mRNA序列翻译成线性的氨基酸链时，蛋白质都是以去折叠多肽或无规则卷曲的形式存在。

蛋白质的基本单位为氨基酸，而蛋白质的一级结构指的就是其氨基酸序列。蛋白质会由所含氨基酸残基的亲水性、疏水性、带正电、带负电等特性通过残基间的相互作用而折叠成一立体的三级结构。

2. 研究蛋白质折叠的目的

因为蛋白质的功能取决于其立体结构，而目前根据已知某基因序列可翻译获得对应蛋白质的氨基酸序列，即蛋白质的一级结构；如果从蛋白质的一级结构就能知道立体结构，那么即可直接从基因推测其编码蛋白质所对应的生物学功能。目前的问题在于，虽然蛋白质可在短时间从一级结构折叠至立体结构，研究者却无法在短时间中从氨基酸序列计算出蛋白质结构，甚至无法得到准确的三维结构。因此，研究蛋白质折叠的过程，可以说是破译折叠密码的过程。

3. “皇冠上的明珠”——蛋白质折叠问题的难度探讨

人体和其他生物体内的蛋白质，都由多种折叠而成。数千个氨基酸组成的长链能自发地折叠成一个稳定的三维结构。理论上讲，使用计算机我们可以推算出一个蛋白质的氨基酸序列折叠后形成的三维结构。因为氨基酸折叠成蛋白质的力学原理很明确，包括氢键、范德华力、疏水作用等相互作用，上千个氨基酸折叠后形成的三维结构，达到了力学最稳态。

不过实际上，蛋白质折叠问题的难度非常大。举个简单的例子，假设每个氨基酸都有2种状态——展开态和折叠态，如果一个蛋白质由100个氨基酸组成，那么它可能的三维结构数量就是2的100次方，这是个非常巨大的天文数字，而其中只有一个结构是稳定的三维结构。

100个氨基酸其实是非常小的蛋白，人体内大多数蛋白质都由数千个氨基酸组成，所以光靠超级计算机的“暴力计算”，是无法根据氨基酸序列预测出蛋白质结构的。这也是蛋白质折叠问题为什么被称为是现代分子生物学的“皇冠上的明珠”。

4. 传统物理和数学方向上的方法和成果

当前传统的研究方法中，诞生了许多分析蛋白质折叠问题的方法，其中效率最高，并最广为人知的就是AlphaFold2折叠。

AlphaFold2是基于AI技术发展而诞生的人工智能系统，由深度思维（Deep
Mind）公司研发，深度思维就是此前研发出“阿尔法狗”（AlphaGo）的那家公司。AlphaFold2最近在国际蛋白质结构预测大赛中夺冠，它的准确性均分达到了92.4/100，而过去的几十年中，其他传统方法只能在40分左右徘徊。首先让我们来了解一下AlphaFold2的工作原理：

AlphaFold2主要通过预测蛋白质中每对氨基酸之间的距离分布，以及连接它们的化学键之间的角度将所有存在的氨基酸对的测量结果汇总成2D形式的距离直方图，然后让卷积神经网络对这些图片进行学习，从而构建出蛋白质的3D结构。AlphaFold2主要架构如下图：

AlphaFold2主要架构

4.1 神经网络EvoFormer

具体来看，AlphaFold2主要利用多序列比对（MSA）把蛋白质的结构和生物信息整合到了深度学习算法中。主要包括两个部分：神经网络EvoFormer和结构模块（Structure
module）。

在EvoFormer中，主要是将图网络和多序列比对结合完成结构预测。图网络可以很好的表示出事物之间的相关性，它可以将蛋白质的相关信息构建出一个图表，以此表示不同氨基酸之间的距离。

神经网络和结构模块

AlphaFold2的研究人员利用卷积神经网络Attention机制构建出一个特殊的“三重自注意力机制”来处理计算氨基酸之间的关系图，他们会将得到的信息与多序列比对结合。多序列比对主要是使相同的残基位点位于同一列，暴露出不同序列之间的相似部分，从而推断出不同蛋白质在结构和功能上的相似关系。操作方式如下图：

“三重自注意力机制”

4.2 结构模块

结构模块是AlphaFold2架构的第二部分，它的主要工作是将EvoFormer得到的信息转换为蛋白质的3D结构。结构模块工作原理如下：

结构模块工作原理

在结构模块中，研究人员同样使用了Attention机制，它可以单独计算蛋白质的各个部分，称为“不变点注意力机制”。是以某个原子为原点，构建出一个3D参考场，根据预测信息进行旋转和平移，得到一个结构框架。具体步骤如下：

不变点注意力机制

之后Attention机制会对所有原子进行预测，最终汇总得到一个高度准确的蛋白质结构。AlphaFold2的研究人员还强调AlphaFold2是一个"端到端"的神经网络，他们会反复把最终损失应用于输出结果。然后再对输出结果进行递归，不断逼近正确结果，这样做既能减少额外的训练，还能大幅提高预测结构的准确性。

4.3 预测结果

下面我们可以来看看AlphaFold2对于蛋白质折叠预测的效果：

图中显示的是327aa的蛋白与同源结构最高identities = 30%的结构预测结果（青色为预测的，绿色的为解析出来的结构）。两个结构RMSD=0.86埃，预测结果非常好。说明AlphaFold2的预测准确性是很高的。

4.4 AlphaFold2的部分代码

为了对AlphaFold2的运行原理理解更透彻，我们也可以看看它进行蛋白质折叠预测的几个算法运行的例子：

4.4.1 折叠单体

假设有一个具有序列的单体。输入fasta应该是：

    >sequence_name

    <SEQUENCE>

然后运行以下命令：

    python3 docker/run_docker.py \

    --fasta_paths=monomer.fasta \

    --max_template_date=2022-05-10 \

    --model_preset=monomer \

    --data_dir= $DOWNLOAD_DIR

4.4.2. 折叠homomer

假设有一个来自原核生物的同源异构体，它有3个相同序列的拷贝。输入fasta就是：

    >sequence_1

    <SEQUENCE>

    >sequence_2

    <SEQUENCE>

    >sequence_3

    <SEQUENCE>

之后运行以下命令：

    python3 docker/run_docker.py \

    --fasta_paths=homomer.fasta \

    --is_prokaryote_list=true \

    --max_template_date=2022-05-10 \

    --model_preset=multimer \

    --data_dir= \$DOWNLOAD_DIR

4.4.3. 折叠异聚体

假设有一个未知来源的异聚体A2B3，即具有的2个副本和3个副本，输入fasta：

    >sequence_1

    <SEQUENCE A>

    >sequence_2

    <SEQUENCE A>

    >sequence_3

    <SEQUENCE B>

    >sequence_4

    <SEQUENCE B>

    >sequence_5

    <SEQUENCE B>

运行以下命令：

    python3 docker/run_docker.py \

    --fasta_paths=heteromer.fasta \

    --is_prokaryote_list=false \

    --max_template_date=2022-05-10 \

    --model_preset=multimer \

    --data_dir= $DOWNLOAD_DIR

4.4.4. 一个接一个的折叠多个单体

假设有两个单体，monomer1.fasta和monomer2.fasta。我们可以使用以下命令按顺序折叠两者：

    python3 docker/run_docker.py \

    --fasta_paths=monomer1.fasta,monomer2.fasta \

    --max_template_date=2022-05-10 \

    --model_preset=monomer \

    --data_dir= $DOWNLOAD_DIR

4.4.5. 一个接一个的折叠多个多聚体

假设有两个多聚体，multimer1.fasta和multimer2.fasta。两者都来自原核生物，可以使用以下命令按顺序折叠两者：

    python3 docker/run_docker.py \

    --fasta_paths=multimer1.fasta,multimer2.fasta \

    --is_prokaryote_list=true,true \

    --max_template_date=2022-05-10 \

    --model_preset=multimer \

    --data_dir= $DOWNLOAD_DIR

4.4.6. AlphaFold2输出

根绝需要预测的氨基酸形式的不同，我们做好每一段的前置代码之后，运行最终的输出程序。输出包括计算的MSA、非松弛结构、松弛结构、排序结构、原始模型输出、预测元数据和部分时间：

    <target_name>/

    features.pkl

    ranked_{0,1,2,3,4}.pdb

    ranking_debug.json

    relaxed_model_{1,2,3,4,5}.pdb

    result_model_{1,2,3,4,5}.pkl

    timings.json

    unrelaxed_model_{1,2,3,4,5}.pdb

    msas/

    bfd_uniclust_hits.a3m

    mgnify_hits.sto

    uniref90_hits.sto

5. 传统方法的局限性及量子计算的优势

虽然目前以AlphaFold2为代表的人工智能技术，凭借之前通过生物学方法积累的大量蛋白质序列和结构信息，在蛋白质折叠问题上能够对蛋白质结构进行比较精准的预测，但是它们的预测速度还是收到很大的限制，以及无法解决问题内在的NP-hard复杂性和蛋白质序列与结构间高维到高维的映射关系。

而凭借量子态叠加和量子纠缠带来的强大并行计算的能力，量子计算能够加速分析过程，解决上述难题。启科量子也在相关实验过程中，总结了自己的量子算法，针对蛋白质折叠问题卓有成效。

6. VQE量子计算：

6.1. 实施原理；

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

图一：实施流程图

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

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

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

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

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

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

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

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

建立案例所用模型的代码如下：

    from
    qiskit_nature.problems.sampling.protein_folding.interactions.random_interaction
    import (RandomInteraction,)

    from
    qiskit_nature.problems.sampling.protein_folding.interactions.miyazawa_jernigan_interaction
    import (MiyazawaJerniganInteraction,)

    from qiskit_nature.problems.sampling.protein_folding.peptide.peptide
    import Peptide

    from
    qiskit_nature.problems.sampling.protein_folding.protein_folding_problem
    import (ProteinFoldingProblem,)

    from qiskit_nature.problems.sampling.protein_folding.penalty_parameters
    import PenaltyParameters

    from qiskit.utils
    import algorithm_globals, QuantumInstance

    algorithm_globals.random_seed = 22

6.2. 准备步骤：

6.2.1. 设定蛋白质主链

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

图二：氨基酸的英文缩写

使用到的代码如下：

main_chain = "APRQ"

# main_chain = "DCAWHKGELVWCT"

6.2.2. 侧链的生成

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

对于侧链的编译，我们使用如下代码：

side_chains = [""] 4

print('side_chains', side_chains)

6.2.3. 氨基酸之间的相互作用

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

random_interaction = RandomInteraction()

mj_interaction = MiyazawaJerniganInteraction()

6.2.4. 物理约束

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

6.2.4.1. Chiral：用于施加正确手性的惩罚参数。

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

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

penalty_back = 10

penalty_chiral = 10

penalty_1 = 10

penalty_terms = PenaltyParameters(penalty_chiral, penalty_back, penalty_1)

6.2.5. 肽定义

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

peptide = Peptide(main_chain, side_chains)

print('peptide', peptid)

6.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极大的加速了哈密顿量的优化，如图一所示。同时，逻辑门参数的经典优化是使用差分进化优化器执行的，它模拟角度θ在希尔伯特空间中的自然选择，优化过程也总结在图一实施流程图中。

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

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

6.3. 实验步骤：

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

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

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

6.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组合而成的系列就是一种蛋白质折叠构型。

这部分的实验代码如下：

protein_folding_problem = ProteinFoldingProblem(peptide, mj_interaction,
penalty_terms)

qubit_op = protein_folding_problem.qubit_op()

print(qubit_op)

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

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

CVaR算法代码如下：

from qiskit.circuit.library import RealAmplitudes

from qiskit.algorithms.optimizers import COBYLA

from qiskit.algorithms import NumPyMinimumEigensolver, VQE,OAOA

from qiskit.opflow import PauliExpectation, CVaRExpectation

from qiskit import execute, Aer

# set classical optimizer

optimizer = COBYLA(maxiter = 50)

# set variational ansatz*

ansatz = RealAmplitudes(reps = 1)

# set the backend*

backend_name = "aer_simulator"

backend = QuantumInstance(

Aer.get_backend(backend_name),

shots = 8192,

seed_transpiler = algorithm_globals.random_seed,

seed_simulator = algorithm_globals.random_seed,

)

counts = []

values = []

def store_intermediate_result (eval_count, parameters, mean, std):

counts.append(eval_count)

values.append(mean)

# initialize CVaR_alpha objective with alpha = 0.1*

cvar_exp = CVaRExpectation(0.1, PauliExpectation())

# initialize VQE using CVaR*

vqe = VQE(

expectation = cvar_exp,

optimizer = optimizer,

ansatz = ansatz,

quantum_instance = backend,

callback = store_intermediate_result,

)

result = vqe.compute_minimum_eigenvalue(qubit_op)

print(result)

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

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

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

完成这部分内容，代码如下：

import matplotlib.pyplot as plt

fig = plt.figure()

plt.plot(counts, values)

plt.ylabel("Conformation Energy")

plt.xlabel("VQE Iterations")

fig.add_axes([0.44, 0.51, 0.44, 0.32])

plt.plot(counts[8:], values[8:])

plt.ylabel("Conformation Energy")

plt.xlabel("VQE Iterations")

plt.show()

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

6.4. 另一种方法

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

#另外一种方法

result = cvar_exp(qubit_op)

cobyla = COBYLA()

cobyla.set_options(maxiter=25)

cobyla = COBYLA()

cobyla.set_options(maxiter=50)

quantum_instance = QuantumInstance(backend=backend_name, seed_simulator=seed,
seed_transpiler=seed)

qaoa_mes = QAOA(optimizer=optimizer,

quantum_instance=backend,

callback=store_intermediate_result)

qaoa = MinimumEigenOptimizer(qaoa_mes)

result = qaoa.solve(qubit_op)

print(result)

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次迭代就接近目标函数了，效率提升肉眼可见。

7. 量子计算未来展望

启科量子本次应用VQE算法解决蛋白质折叠问题的实验，能够很好地解决使用传统方法耗时长、准确率低的问题，从而极大提升现代分子生物学的研究效率，为破解蛋白质折叠谜题带来新希望，进一步推动科学界前进。同时，随着量子计算进一步成熟及应用普及，相信在不久的将来各行各业都能够更多地享受到量子计算带来的成果，为人们的生活带来更大的福祉！