基于HNP模型的强化学习状态空间表示方法

吴宏杰 韩佳妍 杨 茹 陆卫忠 傅启明*

1(苏州科技大学电子与信息工程学院 江苏 苏州 215009)2(苏州大学江苏省计算机信息处理技术重点实验室 江苏 苏州 215006)

0 引 言

蛋白质折叠构形预测问题是现如今生物信息学领域中具有高挑战难度的问题之一，清楚蛋白质折叠构形对生物学发展和人们生活都有重大意义[1]。生物信息学假说认为自然界存在的蛋白质形态是最稳定的形态(即自由能最低)[2]。蛋白质分子的折叠过程是指从蛋白质分子从一般状态变化到基态的复杂过程，它能使我们了解氨基酸序列是如何决定蛋白质分子结构，预测其结构及结构所表现出来的蛋白质分子的性能[3]。目前，国内外的相关研究提出了许多可计算的理论模型,如HP模型[4-6]、AB模型[7-8]和HNP模型[9-12]等。为了能够达到简化对序列空间搜索的目的，很多基于蛋白质设计方法应用于最简单的氨基酸两态分类的HP模型上并且取得了成功[5]。但是进一步研究表明，在很多情况下，简单HP模型存在一些问题，Shakhnovich[13]指出，当蛋白质理论模型设计中考虑更多的氨基酸类型时，这种情况就会得到改善。为此本文将基于强化学习搜索蛋白质序列空间的方法应用到HNP模型上。

如今，强化学习已经在多个领域得到广泛应用。例如车辆定位与识别、游戏自动化检测和机器人仿真等[14]；在生物科学领域也发挥着重要作用，例如生物序列对比、基因组测序等[15]。在强化学习中，Agent可以收到环境的状态反馈，并且与环境产生互动，通过得到的经验来进行自主学习，根据获得的奖赏反馈来找到整体最优解，而不会陷于局部最优。但是存在“维数灾难”(即状态空间的大小随着特征数量的增加而发生指数级的增长)和收敛速度迟缓这两个严重且普遍的问题[16]。本文选用强化学习中经典的Q-learning算法实现HNP晶格模型的最优结构预测。但Q-learning方法及其相关演变算法只在小的状态空间范围内比较有效，一旦训练的状态空间变大，则对计算机内存的要求较高，其存储和计算相应的状态动作对也会变得更加困难，这就是在利用强化学习算法计算蛋白质序列最优结构所遇到的“维数灾难”问题。

目前对于“维数灾难问题”，有研究者提出函数逼近的方法，它无须对状态空间进行预处理，可以用较小的存储空间来描述状态空间到动作的映射关系，但收敛性无法得到保证，理论基础还相当薄弱[16]。有研究者提出了分层强化学习，将复杂问题分解为多个小问题，通过解决多个小问题从而达到能够解决原问题的目的[17]。如果利用分层强化学习思想解决蛋白质序列结构预测问题则会出现序列特征摆放重复的问题，不符合真实的实际蛋白质序列结构。因此，基于HNP晶格模型，本文利用强化学习算法计算序列的最低能量的基础上研究强化学习的状态空间。通过比较全状态空间表示法与简单状态空间表示法，结合两者的优势，提出一种半状态空间表示法，有效解决蛋白质序列特征增多导致状态空间过大的问题，达到改善“维数灾难”问题的目的。在相同的实验环境下，比较这三种不同状态空间表示方法的实验结果，对比得出三种状态空间表示方法的优缺点：1) 全状态空间表示法在短序列的实验上可以准确计算最低能量，但其状态空间受计算机内存限制，无法进行长序列的能量计算；2) 简单状态空间表示法可以对长短序列都进行计算能量，但由于其状态空间设置过小，在计算最低能量的准确度上大打折扣；3) 半状态空间表示法有效解决了全状态空间表示法和简单状态空间表示法的缺点，能够对长序列的能量进行计算且比简单状态空间的准确度高。

1 HNP模型

1.1 概 念

很多学者不断提出针对于蛋白质折叠结构问题的简化模型，如HP模型[4-6]和AB模型[7-8]。而Miyazawa等[18]和Zhang等[19]采用统计方法计算了这20种氨基酸之间的相互作用能，根据相互作用能的大小，提出把20种氨基酸分为3类。根据相互作用能的大小，分别将氨基酸Cys、Met、Phe、Ile、Val、Trp、Tyr称为疏水(H)基团，氨基酸Ala、Gly、Thr、His称为中性基团(N)和Ser、Asn、Gln、Asp、Glu、Arg、Lys、Pro称为亲水(P)基团，即HNP模型[9]。

1.2 能量函数表示

根据氨基酸相互作用能的大小，将氨基酸分为3类：疏水氨基酸(用H表示)、中性氨基酸(用N表示)、亲水氨基酸(用P表示)。通过氨基酸的划分把蛋白质链转化为以H、N、P构成的序列。

根据蛋白质分子的Hamiltonian 能量函数[20-21]：

(1)

式中：Si、Sj表示氨基酸分类H、N和P；ESi-Sj表示氨基酸与氨基酸之间的相互作用能量，EH-H=-5.59 RT，EH-N=-3.68 RT，EH-P=-3.07 RT，EN-N=-2.23 RT,EN-P=-1.78 RT,EP-P=-1.37 RT。rij指的是氨基酸与氨基酸之间的距离；l为假定的蛋白质分子的键长；δ(rij,l)是判断氨基酸是否形成一个紧密接触对的条件，如果rij=l，则δ(rij,l)=1，即形成一个紧密接触对，如果rij≠l，则δ(rij,l)=0，即不能形成一个紧密接触对。

2 基于HNP模型强化学习状态空间表示方法

2.1 基于不同状态空间的统一强化学习框架

强化学习(Reinforcement Learning, RL)是指从环境状态到行为映射的学习，目的是使Agent从环境中获得的积累奖赏值达到最大[22]。Q-learning算法是强化学习中一种独立于模型的TD控制算法，它可以用来为任何给定的马尔可夫决策过程(Markov Decision Process,MDP)找到最佳的行动选择策略。Q-learning算法的最简单形式被定义为：

(2)

Q-learning算法最大的优点就是无需环境模型，在学习过程中，该算法先建立一个Q值表，设其初值为0。当Agent采取动作at转移至下一状态st+1，并获得立即奖赏rt+1，根据式(2)，更新Q值表，经过不断迭代更新，Q值表会最终收敛到最优，根据Q值表选择最优动作值函数可以选取到最优动作。

本文通过强化学习中的Q学习算法来实现HNP模型结构的预测。其流程框架如图1所示。

图1 基于HNP模型的强化学习框架

(1) 序列输入与预处理。本文将氨基酸序列转化为HNP序列的形式，在运算时，为了方便考虑，通过-1、0和1来对应简化代替H(疏水)、N(中性)和P(亲水)三类氨基酸。

(2) 基本元素。在强化学习系统中，除了Agent和环境之外，还包含有动作、状态值函数和奖赏(回报函数)这3类基本元素。而马尔可夫决策过程(MDP)则包含5个模型要素，分别为状态(state)、动作(action)、策略(policy)、奖励(reward)和回报(return)[23]。

(3) 结果输出。通过理论可以得知，如果所有的状态动作对能够不断地被更新，那么Q值表能够被收敛到最佳结果，Agent可以根据收敛到的Q值更新表选择最优动作，从而得出HNP模型的最优结构预测。

2.2 状态空间表示方法

在强化学习中，状态是指学习器发现自身当前所处的情况：对于Agent而言，状态可以是其所有执行器的位置和角度，以及其地图位置、电池状态、其试图发现的目标等[24]。在马尔可夫系统中，所有这些变量都被赋予学习器。状态也可以被定义为学习器可用于发现其当前状况的所有信息的总结。这些信息包含在它所经历动作和取样的历史中，意味着完整的历史构成完美的状态，这也是衡量一个状态表示的标准[25]。

强化学习中的Q-learning算法是一种以马尔可夫决策过程为理论基础的算法。马尔可夫性是指Agent在当前状态s下采取当前动作a达到下一状态st+1以及得到的立即奖赏只与当前状态s和当前动作a有关，与之前的状态和动作都无关。也就是说，Agent和环境在一个离散时间序列(t=0,1,2,…)的每一步中都进行交互，在每个时间步t，Agent都能得到若干环境状态st∈S，其中S是所有可能状态的集合。

本文将蛋白质序列中的氨基酸可能摆放的位置看作一个状态，想要计算出蛋白质序列的最低能量也就是需要罗列出Agent采取每一步动作所有可能的状态的集合。而当状态空间不断增大时，Agent与环境交互的时间越来越长，其存储和计算相应的状态动作对也会变得更加困难，对计算机内存要求更高。需要简化这个历史，而不丢失信息。因此，基于HNP蛋白质模型，对以下三种状态空间设置方法进行研究。

(1) 全状态空间设置方法。对于一个长度为n的蛋白质序列，在二维折叠空间中，有四个可能的折叠方向，分别为：L(左)、U(上)、R(右)和D(下)，即A=[L,U,R,D]。如图1(a)所示，首先随机固定第一个氨基酸的状态，后续的每个氨基酸可能的状态是上一个氨基酸4个可能方向L(左)、U(上)、R(右)和D(下)的集合，也就是说后一个氨基酸的所有可能状态个数是前一个氨基酸状态个数的4倍，状态转移表如表1所示。

表1 全状态空间转移表

续表1

随机固定第一个氨基酸的状态s1，第一个氨基酸状态个数为1。s1采取四个可能的折叠方向，转移到的第二个氨基酸的可能状态{s2,s3,s4,s5}，则第二个氨基酸所有可能状态个数为4。如果s2变为当前状态时再转移到第三个氨基酸的可能状态{s6,s7,s8,s9}；同理，如果s3、s4、s5变为当前状态时也转移到第三个氨基酸的可能状态{s10,s11,s12,s13}、{s14,s15,s16,s17}、{s18,s19,s20,s21}，则第三个氨基酸所有可能状态个数为4×4=16。以此类推，将整个状态集S的状态总数视为一个初值为1、公比为4的等比数列求和：

(2) 半状态空间设置方法。对于全状态空间过大，简单状态空间略小的问题，本文提出一种半状态空间设置方法，如图1(b)所示，如果当前氨基酸的上一状态采取的动作相同，则采取所有可能的动作转移到相同的下一状态，相当于固定了16个状态空间，状态转移表如表2所示。

表2 半状态空间转移表

续表2

随机固定第一个氨基酸的状态s1，则第一个氨基酸状态个数为1。s1采取四个可能的折叠方向，转移到的第二个氨基酸的可能状态{s2,s3,s4,s5}，则第二个氨基酸所有可能状态个数为4。如果s2变为当前状态时再转移到第三个氨基酸的可能状态{s6,s7,s8,s9}；同理，如果s3、s4、s5变为当前状态时也转移到第三个氨基酸的可能状态{s10,s11,s12,s13}、{s14,s15,s16,s17}、{s18,s19,s20,s21}，则第三个氨基酸所有可能状态个数为16。从第三个氨基酸开始，考虑当前动作和上一动作转移到下一氨基酸的可能状态。即如果s6作为当前状态，上一动作为L，采取当前动作为L，则转移到s22；采取当前动作为R，则转移到s23；采取当前动作为U，则转移到s24；采取当前动作为D，则转移到s25。如果s7作为当前状态，上一动作为R，采取当前动作为L，则转移到s26；采取当前动作为R，则转移到s27；采取当前动作为U，则转移到s28；采取当前动作为D，则转移到s29。如果s8作为当前状态，上一动作为U，采取当前动作为L，则转移到s30；采取当前动作为R，则转移到s31；采取当前动作为U，则转移到s32；采取当前动作为D，则转移到s33。如果s9作为当前状态，上一动作为D，采取当前动作为L，则转移到s34；采取当前动作为R，则转移到s35；采取当前动作为U，则转移到s36；采取当前动作为D，则转移到s37。如果s10作为当前状态，上一动作为L，采取当前动作为L，则转移到s22；采取当前动作为R，则转移到s23；采取当前动作为U，则转移到s24；采取当前动作为D，则转移到s25。以此类推，从第三个氨基酸开始，每个氨基酸的所有可能状态都为16，则整个状态集S的状态总数为：

1+41+42+42…+42=1+4+16×(n-2)=16n-27

所以一个蛋白质分子的半状态集S={s1,s2,…,s16n-27}。

(3) 简单状态空间设置方法。全状态空间的设置方法全部采取四个可能的折叠方向，全部罗列导致状态空间过大，其存储和计算序列的状态动作对会比较困难。因此本文提出第一步随机固定第一个氨基酸的状态，从第二步开始每一个氨基酸可能的状态都固定为四个状态空间，图形表示如图1(c)所示，状态转移表如表3。随机固定第一个氨基酸的状态s1，则第一个氨基酸状态个数为1。s1选择四个可能的折叠方向转移到下一氨基酸的可能状态{s2,s3,s4,s5}，则第二个氨基酸所有可能状态个数为4。如果s2变为当前状态时再转移到第三个氨基酸的所有可能状态{s6,s7,s8,s9}；同理，如果s3、s4、s5变为当前状态时也转移到第三个氨基酸的所有可能状态{s6,s7,s8,s9}，则第三个氨基酸所有可能状态个数也为4。如果s6变为当前状态时再转移到下一可能状态{s10,s11,s12,s13}；同理，如果s7、s8、s9变为当前状态时也转移到下一可能状态{s10,s11,s12,s13}，则第四个氨基酸所有可能状态个数也为4。由此类推，把整个状态集S视为第一个氨基酸的状态加上后面n-1个氨基酸状态所有可能状态的集合；状态总数为：

1+41+41+…+41=1+4×(n-1)=4n-3

所以一个蛋白质分子的简单状态集S={s1,s2,s3,…,s4n-3}。

表3 简单状态空间转移表

续表3

综上所述，可以得出三种不同的状态空间设置表示的总个数。为了便于观察以及方便进行对比，将三个方法的状态空间集记录列于表4中。

表4 三种状态集表示

2.3 奖赏函数设置与策略选择

在强化学习中，奖赏函数把环境中感知到的状态动作对映射为一个单独的数字，即奖赏，表示该状态的内在需求。Agent的唯一目标是在长期运行过程中实现奖赏最大化。本文将强化学习思想与HNP模型相结合，为了能够找到能量最低的结构模型，也为了能够让Agent更好地训练出所需要的理想化结构，本文选择根据上述所说的Hamiltonian能量函数计算公式，将氨基酸与氨基酸之间的相互能量作为奖赏函数。因此，奖赏函数的设置如下：

(1) 当序列未完全摆放好时，其奖赏设置为0。

(2) 当序列摆放到最后一个氨基酸时，判断氨基酸与氨基酸之间是否满足形成紧密接触对的条件。如果满足，则相应的奖赏值为氨基酸之间的相互作用能的绝对值|E|；如果不满足，则奖赏值为0。

MDP策略是按状态给出的，动作的条件概率分布在强化学习下属于随机性策略。在算法的训练阶段，本文采用“探索-贪心(ε-greedy)策略”进行动作选择。ε-greedy策略既考虑当前回报，又兼顾将来回报，算法中随机采取概率ε(0<ε<1)选取下一个动作，利用剩下1-ε的概率去采取最优动作，从而不断计算更新Q值表。

2.4 训练算法

流程如算法1所示，其中：S表示状态设置，Q表示状态动作Q值矩阵，at表示动作，R为奖赏函数，st为当前状态，st+1为下一状态，terminal为终点状态,n为当前状态数，N表示状态总数。

算法1HNP模型训练算法

1. //设置全状态空间

2. state_transfer=0

3. for n in range(np.int((4 ** N-1)/3)):

4. state_transfer+=1

5. Or //设置半状态空间

6. state_transfer=list(range(2，22))

7. for n in range(4，N+1):

8. state_transfer+=list(range((n*16-42),n*16-26))*4

9. Or //设置简单状态空间

10. state_transfer=list(range(2，6))

11. for n in range(2，N):

12. state_transfer+=list(range((n*4-2),n*4+2))*4

13. Initialize:st,Q(s,a)

14. Repeat (for each episode){

15. Initializest

16. While(st!=ternimal){

17.at← Q(ε-greedy)

18. Take actionat, returnR,st+1

20.s←st+1

21. }

22. }

3 实验与结果分析

3.1 实验参数选择

为了保证在实验过程中的客观性以及值函数的准确性和收敛性，需要设置一定参数来进行实验。实验中参数有学习参数γ、步长参数α和探索概率ε。学习参数γ又称折扣率，决定着未来奖赏的当前价值。在强化学习问题中，一般取学习参数γ=0.9，使得Agent在学习探索过程中能够找到全局最优解而不陷于局部最优解。步长参数和探索概率共同影响着值函数的准确性和收敛性。步长参数过大，会引起值函数震荡，难以收敛；过小则使得值函数收敛速度缓慢。探索概率过大，Agent为寻求更优解会不断更新值函数，不易收敛；过小则缺少探索信息，值函数收敛缓慢。因此需要选择最优参数以达到收敛效果最佳的目的。

本文在设置学习参数γ=0.9的前提下，选择不同的步长参数α和探索概率ε进行实验。以序列HHPNHHNNHH(长度为10)作为实验对象，进行5轮实验，训练迭代次数为30万次，改变步长参数α和探索概率ε，多次实验计算蛋白质序列在不同参数下达到最低能量所需的平均次数，从而比较得出最优参数。得出数据如表5所示。

表5 三种状态空间设置方法在不同参数下序列

由表5所示，在探索概率ε数值相同的情况下，随着步长参数α的增大，序列达到最低能量所需平均次数逐渐下降，当α=0.9的时候，实验数据趋于收敛；在步长参数α数值相同的情况下，探索概率ε=0.5测试出的所需次数较ε=0.2和ε=0.8更优。综上所述，本文选取α=0.9、ε=0.5作为最优参数。

3.2 对比实验

根据三种状态空间表示方法设置的不同对蛋白质序列进行模型训练。本文在PIR数据库中随机选取了11条长度不同的序列作为实验对象，将其转化为HNP序列，将序列集信息记录于表6。为避免偶然性，按照三种状态空间设置的方法对每条序列进行5轮测试实验，每轮训练迭代次数设置为30万次，保证其他参数条件以及实验环境不变，记录在迭代次数内每条序列可以计算出的最低能量并统计于表6，其中未能测出序列最低能量的结果以“-”表示。

表6 三种状态空间方法测得HNP序列的最低能量 单位：RT

由表6可以看出，当序列长度较短时，以前面5条序列为例，基于全状态空间方法和半状态空间方法的实验都可以测出序列的最低能量；而基于简单状态空间方法的实验则不能测出序列3和序列5的最低能量，相比于另外两种方法在计算序列最低能量的准确度上降低了40百分点。此外，基于简单状态空间方法的实验在训练序列2和序列4的实验中不能每次都计算出最低能量，相比于基于全状态空间方法的实验在五轮训练实验的准确度上平均降低了12百分点；基于半状态空间方法的实验在训练序列3的实验中不能每次都计算出最低能量，相比于基于全状态空间方法的实验在五轮训练实验的准确度上平均降低了4百分点，相比于基于简单状态空间方法的实验在五轮训练实验的准确度上平均提高了8百分点。当序列长度逐渐增大时，以后面6条序列为例，基于全状态空间方法的实验则会出现状态空间过大的情况，导致在相同实验环境下无法计算后面6条长序列的最低能量；基于半状态空间方法的实验相比于基于简单状态空间方法的实验能计算出更低的能量，尤其对于序列6、7、9、10和11，基于半状态空间方法的实验比基于简单状态空间方法的实验在每轮训练实验中都计算出更低的能量。由此可见，对于长序列而言，半状态空间方法比全状态方法和简单状态空间方法能找到更符合实际蛋白质序列的最低能量。

综上所述，半状态空间方法和简单状态空间方法有效解决全状态空间无法计算长序列的缺点。简单状态空间较全状态空间有3条序列预测出更低能量，半状态空间较全状态空间方法全部6条长序列都预测出更低能量，且半状态空间预测的能量平均值较简单状态空间降低了9.83百分点。

4 结 语

本文以强化学习算法为核心，对蛋白质HNP模型结构进行训练模型，计算最接近真实的蛋白质序列结构。针对强化学习中的状态空间问题，对比研究了三种不同的状态空间设置方法。全状态空间表示方法全面地列出了序列所有的可能状态空间，能够较为准确地计算出最接近真实的蛋白质序列结构，但每当序列长度增加时，状态空间也会以指数级增长，导致状态空间过大，其存储和计算状态动作对较为困难，对计算机内存的要求也变高。半状态空间表示方法和简单状态空间表示方法在全状态空间表示方法的基础上将其状态空间进行简化，缩小了状态空间，可以计算长度更长的序列，但也降低了一定的准确度。此外，半状态空间表示方法相较于简单状态空间表示方法在计算蛋白质序列最低能量的准确率上更具优势。

本文是强化学习在生物信息领域新的尝试。目前还存在需要解决的问题，如无法在确保准确度的前提下合理缩小氨基酸序列的状态空间。所以，在研究HNP模型的蛋白质分子构象和折叠问题中，需要提出更进一步的算法改进。