谈谈分子动力学模拟蛋白-配体复合物过程中配体发生脱离的原因
谈谈分子动力学模拟蛋白-配体复合物过程中配体发生脱离的原因
On the reasons why ligand detaches during molecular dynamics simulation of protein-ligand complex
文/Sobereva@北京科音2022-Jan-13
有不少人问怎么他做基于经典力场的蛋白质-配体复合物的分子动力学模拟的过程中配体跑出去了。鉴于这个问题被问得很频繁，笔者遂写个小文把常见的可能因素说明一下，便于读者排查原因。
首先要凭基本化学常识和直觉判断配体是否真有可能稳定与蛋白结合，可以在VMD程序中把配体周围一定距离的蛋白质原子都显示出来观察，这需要利用within语句来定义显示范围，详见《VMD里原子选择语句的语法和例子》（//www.umsyar.com/504）。配体和蛋白结合常见作用机制有氢键（通常主要本质是静电吸引作用，见//www.umsyar.com/513）、盐桥（配体显离子性部分与氨基酸侧链带相反电荷基团的静电吸引）、范德华作用（具体来说是指色散吸引效应，pi-pi堆积本质上与此也相同）、疏水作用（本质上来自于溶剂的熵效应，不懂的话Google）。如果你对弱相互作用类型不太了解，可以看看此文里相关科普：《谈谈“计算时是否需要加DFT-D3色散校正？”》（//www.umsyar.com/413）。基本上，对弱相互作用有一定了解、接触过较多蛋白-配体复合物体系研究的人都能凭感觉初步估计在一个给定的复合物结构中，一个小分子和蛋白质是否有可能稳定结合上。如果在初始结构下，上述这些导致配体和蛋白相互结合的作用非常不显著，显然小分子可能由于热运动跑到溶剂里去，或者跑到蛋白质上其它更容易结合当前配体的区域去。值得一提的是，有一些程序可以帮助你大致了解一下配体和蛋白质之间的相互作用，比如LigPlot+程序（https://www.ebi.ac.uk/thornton-srv/software/LigPlus/）和https://proteins.plus在线服务器里的PoseView。此外，proteins.plus里的DoGSiteScorer还可以找出蛋白质中的口袋并图形化展现出来，并给出Drug Score体现口袋的类药性（结合药物分子的可能性），有时能帮助你做一些粗略的判断、辅助确定对接需要考虑的空间范围。
初始的复合物结构是怎么得到的，显然直接关系到跑动力学期间配体能否在初始位置附近呆住。有这么几个情况，应当分情况考虑
(1)用的是X光衍射测的复合物结构：这样的初始结构通常质量较高，小分子理应能在动力学过程中结合在初始的位置。但是也要注意，晶体状态的复合物结构只是实验测定条件下的最稳定结构，而在现实条件的动力学模拟过程中，由于热运动和溶剂效应，小分子未必就一定能始终稳定呆在晶体中原本的位置。另外，有些晶体结构的解析度较低，比如三点几埃，此时配体、氨基酸残基位置的不确定性很高，可能此结构下侧链与配体的相互作用情况和实际有明显偏差，进而可能导致模拟开始后没多久配体就跑掉。
(2)用的是分子对接得到的复合物结构：如果没有复合物晶体结构但是有蛋白质晶体结构，通常就是用对接得到的打分最高到的复合物结构当初始结构。受制于采样充分程度、打分函数的准确度、对受体的柔性考虑等因素，对接产生的初始结构的可靠性明显不如较高解析度的晶体衍射测出来的结构。很多人做过分子对接后都是再专门跑一下分子动力学检验配体能否基本维持在对接产生的位置，从而检验对接产生的结构的合理性。使用这种结构当初始结构时，如果配体在动力学模拟过程中跑掉了，很可能就是因为对接产生的打分最高的结构并不合理。此时可以考虑用你觉得看上去也靠谱但打分不是最高的复合物结构试试，或者换其它对接程序得到的结构试试。
(3)用的是对接得到的结构，而且蛋白质结构是预测出来的：这无疑是可靠性最低的结构，因为蛋白质本身结构的可靠性都成问题。这种情况配体如果在模拟中跑出来，除了考虑对接的问题外，当然也得考虑蛋白质三维模型是否真的合理。蛋白质预测程序/方法很多，其中一种不理想时可以再考虑别的，选择余地较大，比如基于同源模建的Modeller、Robetta等，基于深度学习预测蛋白质结构的AlphaFold 2和RoseTTAFold。
用上面(2)、(3)方式得到的初始结构，在跑动力学之前记得仔细观察配体和受体结合区域的特征，凭常识估计是否真有可能靠非共价的方式稳定结合，明显没什么希望的结构直接pass掉。
用的力场不合理也是可能导致小分子在动力学过程中跑掉的原因，但记住这绝对不是唯一原因，别一发现小分子跑掉就总是赖力场。目前模拟蛋白质+小分子体系的常见组合是AMBER（蛋白质）+GAFF（小分子）、CHARMM（蛋白质）+CGenFF（小分子）、GROMOS（直接描述蛋白质，并结合ATB创建的配体拓扑文件）、OPLS-AA（直接描述蛋白质，并结合LigParGen创建的配体拓扑文件）。其实以上这些蛋白质+小分子描述的组合，对于前述的范德华作用的描述没太大差异（而且通常范德华作用对蛋白质-配体结合的贡献相对来说较次要），而且也不直接影响疏水作用（疏水作用是在动力学模拟过程中自然而然体现的，不是在力场层面上描述的），也没法说哪个一定更好。对于GROMACS用户，我建议优先考虑AMBER14SB结合sobtop（//www.umsyar.com/soft/Sobtop）创建的基于GAFF力场的小分子拓扑文件，主要是因为这俩力场都比较鲁棒，AMBER14SB有GROMACS现成的力场包而且描述蛋白质构象较好（虽然目前也有更新的AMBER19SB，但截止到撰文时网上还没有现成的GROMACS力场包，对于描述蛋白-配体作用方面它比AMBER14SB没有额外优势），而且sobtop程序用起来很方便又特别可靠。对小分子拓扑文件创建不了解的话看《几种生成有机分子GROMACS拓扑文件的工具》（//www.umsyar.com/266）。静电作用对于蛋白-配体结合的稳定程度起至关重要的影响，而且它直接决定了对氢键和盐桥的描述合理性。以上力场中描述静电作用都是靠原子电荷基于库仑公式来计算的，而这些力场对蛋白质的标准氨基酸残基都已经明确定义了原子电荷，因此在特定的蛋白&配体力场组合下，影响结合稳定性的关键在于配体的原子电荷的确定，这是非常有讲究的。对于AMBER+GAFF的组合，给配体用RESP2(0.5)原子电荷即便不是最理想至少也是最理想之一，可以用Multiwfn程序（//www.umsyar.com/multiwfn）基于Gaussian或ORCA等 程序优化的结构和产生的波函数直接算出来，见《RESP2原子电荷的思想以及在Multiwfn中的计算》（//www.umsyar.com/531）。RESP2(0.5)对于其它蛋白-配体力场组合不一定是最理想的，因为牵扯到与范德华参数、二面角参数的耦合等问题，这里不打算展开细说，但至少用RESP2(0.5)仍是合理选择之一。使用上述力场组合（假设蛋白质力场部分都是较新版本），对于蛋白质结构的描述都没有什么问题，但对于成键方式不够典型的配体来说，有可能出现在某个力场下它的结构能正确维持，而在另一个力场下则明显不合理的情况。如果配体自身结构都没法被正确描述（如某个局部本应是平面的而在动力学过程中却严重弯曲了），自然也就不能指望配体和蛋白口袋的相互作用肯定能合理描述。因此，做动力学的时候要观察轨迹，留意小分子在模拟过程中的结构是否和期望一致，必须确保没有明显异常。如果你没有化学直觉，不知道一个比较特殊的配体分子（假设是有机的）理应是什么样，可以用比如B3LYP-D3(BJ)/6-311G*这种常见也靠谱的几何优化级别用 程序做优化，还可以用molclus（http://www.keinsci.com/research/molclus.html）做构象搜索。但也应注意，对于柔性较大的配体，蛋白质与它的相互作用往往会明显影响其构象，可能在结合时它的构象明显和孤立状态下自由能最低构象不同，这是很正常的。
如果在AMBER14SB+GAFF结合RESP2(0.5)原子电荷下模拟蛋白-配体复合物时出现配体中途跑掉的情况（而且看上去小分子自身结构的描述没啥问题），换力场也大概率解决不了。虽然也可以再试一个别的力场组合，但如果发现还不行，一般也没必要再试了，而应该从本文提及的其它角度考虑跑掉的原因。
有些特殊的蛋白-配体相互作用不是靠以上常见的力场能直接描述的。比如碘代苯和蛋白质上负电性原子间可以形成显著的卤键，这光靠原子电荷是没法展现的，而必须表现出卤原子周围电荷分布的各向异性，这需要修改原先的模型，考虑额外点电荷，或考虑原子多极矩。例如在JCTC, 14, 5383 (2018)有人提出给GROMOS力场的卤原子的σ-hole区域增加额外的正电性的点电荷使得卤键能够定性正确表现出来。PS：在GROMACS中如果想用此文的方式，需要依靠虚拟点来实现。另外，如果蛋白-配体之间有涉及过渡金属的配位键出现，光靠原子电荷描述的静电吸引作用一般也不足矣正确描述和维持之，此时往往需要加上限制势人为强行维持住。而对于涉及到高价主族金属离子的蛋白-配体结合，光靠原子电荷表现静电作用可能也描述不充分，因为明显的极化作用对结合会产生稳定化，这一般需要 方法才能准确描述，因此这种情况下若配体在动力学过程中跑飞了的话也得考虑加上限制势来强行维持住金属离子与周围原子的作用。
配体、残基侧链在结合时所处的质子化状态也是要恰当考虑的，否则无法正确表现蛋白-配体间的实际相互作用。有的配体、氨基酸残基在不同pH下有不同的质子化状态，而且蛋白-配体的相互作用还可能影响质子化状态。对于组氨酸尤其要注意考虑，它在模型体系中的pKa接近7，在生理环境下模拟时其侧链的质子化态有较大不确定性，受环境影响明显，既可以处于中性（侧链上的一个氮被质子化），也可以带一个正电荷（侧链上两个氮都被质子化）。当配体与组氨酸可能存在氢键作用时，组氨酸的质子化状态应当是有利于氢键形成的。你缺乏相关知识和判断能力的话也不要紧，有一些程序能自动给你判断蛋白和配体合理的质子化状态，例如https://proteins.plus在线服务器里的Protoss工具（介绍见https://www.zbh.uni-hamburg.de/forschung/amd/software/protoss.html）。如果你是用GROMACS的话，pdb2gmx产生的蛋白拓扑文件时可以加上-his选项来人工选择各个组氨酸的质子化态成为Protoss给出的情况，而使用小分子拓扑文件产生工具时应当提供合理的质子化态的结构文件来产生。
动力学模拟设置不合理也是可能导致小分子跑掉的原因。做蛋白-配体的正式动力学模拟之前一般都是在给蛋白质和配体加上位置限制势的情况下跑一段短暂的动力学，以使得溶剂分子在这个过程中能充分弛豫、迎合复合物结构。如果没做这一步而直接跑常规的动力学的话，有可能会由于初始加的水的位置不够理想，不仅没有描述充分实际的水环境，甚至造成出现一些不合理的互斥，导致配体跑出去。另外，做生物体系的动力学我一般建议让参考温度从0 K开始缓慢线性上升，这样最为稳妥。如果让程序随机产生初速度，而且初速度对应的温度又不太低且没加限制势，也没准会导致在模拟初期配体和蛋白的相互作用尚未充分稳定的时候配体由于拥有过高的动能而跑走。
如果复合物结构是通过对接方式产生的，发现动力学模拟还没有跑多久配体就脱离了，可以考虑给小分子和与之相互作用的残基之间恰当地加上一个或多个较弱的距离限制势（怎么加根据具体结构随机应变，比如如果配体和小分子之间通过氢键维持，可以给氢键作用的原子间加上限制势维持住氢键距离），避免小分子跑掉。等动力学跑一阵子后，蛋白质的口袋区域的构象就弛豫了，蛋白-配体相互作用会比一开始的状态有所增强，之后再去掉限制势进行正式的模拟，说不定配体就能一直呆住了。注意，除非是由于前述的力场的局限性而不得不加限制势，否则正式模拟时限制势必须去掉，否则相当于给模拟引入了明显的人为因素，发文章的话容易被审稿人抨击。
不要以为配体在现实中就应该始终和蛋白是维持结合状态的。配体L和蛋白P结合成复合物LP的速率为k_on[L][P]，其中k_on为结合速率常数；而LP解离为L和P的速率为k_off[LP]，其中k_off为解离速率常数。当达到稳态时满足k_on[L][P]=k_off[LP]，并且有Kb=k_on/k_off=[LP]/([L][P])=1/Kd的关系，其中Kb称为结合常数，Kd称为解离常数。显然，复合物的形成和复合物的解离是同时进行的，不是说配体结合到蛋白后就一直不出来了，也因此Kb不是个无穷大的值。k_off越大体现配体越容易解离、能结合的平均时间越短，这和解离过程的自由能垒有密切联系。k_off较大的配体说明它结合得就是不很稳定，就是容易在动力学过程中跑出去，所以不能一味地认为模拟中出现解离的现象就一定是错误的。有的人跑复合物的模拟并非才跑了没多久配体就脱离了，而是跑了好几百ns才出现配体脱离现象，很可能正是体现出真实的配体解离现象。
还应当考虑自己做的模拟是否有真的意义，期望出现的现象是否真的是真实的。比如之前我看到有人研究的是酶，问怎么小分子在里面没法一直呆住。实际上要求小分子呆住本来也没意义，本身底物小分子偶然进去之后被催化完就走掉了，何故期望小分子能稳定呆在里面？