（凝聚态物理专业论文）蛋白质的并行分子动力学模拟及全电子结构计算.pdf

中文摘要 摘要 利用并行分子动力学的计算机模拟和蛋白质分子全电子结构计算，来研究蛋 白质分子的结构、动力学性质和生物活性。这里，主要包括了以下两个方面的工 作： 1 对肌红蛋白和铜锌超氧化歧化酶进行分子动力学模拟。通过对方均位移以及 一系列自由能非线性指标的分析，验证了准弹性中子散射实验的结果：温度 在约2 0 0k 以下时，两种蛋白质的行为都表现为类似于谐振固体的运动；当 高于约2 0 0k 时，出现了明显的动力学相变，进入非谐振运动。并且发现水 环境对此动力学相变起着重要作用，而这个相变温度很可能与酶的活性有很 大关联。 2 结合蛋白质全电子结构计算方法、并行计算分子动力学模拟和蛋白质表面的 拓扑结构计算，研究了b a m h i 限制性内切酶的热敏感性突变体p 1 7 3 l 的活性。 证明了实验中确认的4 个活性位点：7 7 、9 4 、1 l l 和1 1 4 号氨基酸残基；并 通过比较不同温度下这四个氨基酸残基所占有的前线轨道能级与最高占有轨 道能级差的变化，找出活性随温度变化的规律。 关键词：蛋白质，并行分子动力学，动力学相变，全电子结构计算，前线轨道， 活性 英文摘要 a b s t r a c t p a r a l l e l m o l e c u l a r d y n a m i c s s i m u l a t i o t i sa n dw h o l ee l e c t r o n i c s t r u c t u r e c a l c u l a t i o na r eu s e dt o i n v e s t i g a t es t r u c n l r e d y n a m i c s a n db i o l o g i c a l a c t i v i t y o f p r o t e i n s i nt h i st h e s i s ，t w ot o p i c s w i l lb ed i s c u s s e d 1 m y o g l o b i na n dc u z ns u p e r o x i d ed i s m u t a s e a r es i m u l a t e db yp a r a l l e lm o l e c u l a r d y n a m i c s b ya n a l y s i s o fm e a ns q u a r e d i s p l a c e m e n t s a n ds e v e r a ln o n l i n e a r i n d i c e s ，t h e s i m u l a t i o n s r e p r o d u c e t h e q u a s i e l a s t i c n e u t r o n s c a t t e r i n g e x p e r i m e n t a lr e s u l t s ：b e l o w a r o l l i l d2 0 0kt h e s et w op r o t e i n sb e h a v ea sh a r m o n i c s o l i d sw i t he s s e n t i a l l yo n l yv i b r a t i o n a lm o t i o n ，w h i l ea b o v et h i st e m p e r a t u r et h e r e i sa s t r i k i n gd y n a m i c t r a n s i t i o ni n t oa n h a r m o n i cm o t i o n m o r e o v e r ，t h e s i m u l a t i o n sf u r t h e rs h o wt h a tw a t e rm o l e c u l e sp l a ya ni m p o r t a n tr o l ef o rt h i s d y n a m i c a lt r a n s i t i o n ，a n dt h ed y n a m i c a lt r a n s i t i o nm i g h th a v eg r e a tr e l e v a n c e w i t l le n z y m e a c t i v i t y 2 t h eh e a t s e n s i t i v em u t a n to f b a m h i ，n a m e da sp 1 7 3 l ，w a ss i m u l a t e db yp a r a l l e l m o l e c u l a rd y n a m i c s c o m b i n e dw i t hw h o l ee l e c t r o n i cs t r u c t u r ec a l c u l a t i o na n d s u r f a c et o p o l o g yc a l c u l a t i o n ，i t sb i o l o g i c a la c t i v i t yi si n v e s t i g a t e d w ev a l i d a t e t h er e s u l ti nt h ee x p e r i m e n t s ：r e s i d u e7 7 ，9 4 ，111a n d113a r ea c t i v es i t e s i n a d d i t i o n ，t h ec h a n g eo f t h ee n e r g y g a p s b e t w e e nt h ef r o n t i e ro r b i t a le n e r g yl e v e l s o c c u p i e db yt h e s ef o u rr e s i d u e sa n dt h eh i g h e s to c c u p i e dm o l e c u l a ro r b i tr e v e a l s t h er u l eo fe n z y m e a e t i v i t yv a r i a t i o na l o n gw i t l lt h et c m p e r a t u r e k e yw o r d s ：p r o t e i n ，p a r a l l e l m o l e c u l a r d y n a m i c s ，d y n a m i c a lt r a n s i t i o n ，w h o l e e l e c t r o n i cs t r u c t u r ec a l c u l a t i o n ，f r o n t i e ro r b i t ，a c t i v i t y 2 - 引言 ”t h en e w p a r a d i g m , n o we m e r g i n g , i st h a ta l lt h e g e n e s tw i l lb e k n o w n ( i nt h e s e n s eo fb e i n gr e s i d e n t 加d a t a b a s e sa v a i l a b l e e l e c t r o n i c a l 幼，a n d t h a tt h e s t a r t i n gp o i n to f a b i o l o g i c a l i n v e s t i g a t i o nw i l lb et h e o r e t i c a l a ni n d i v i d u a ls c i e n f i s tw i l lb e g b 2 w i t hat h e o r e t i c a lc o n j e c t u r e , o n t yt h e nt u r n i n gt oe x p e r i m e n tt o f o l l o wo rt e s tt h a th y p o t h e s i s l l j 胁t e rg f l b e r t , n o b e ll a u r e a 船1 9 9 1 2 1 世纪是生命科学的世纪，生命物质和生命现象成为科学家研究的重要对 象。尤其是在历时1 3 年、耗资数十亿的著名的人类基因组计划( h u m a ng e n o m e p r o j e c t ，简称h g p ) 胜利完成后，后基因时代拉开了序幕。面对网络数据库中呈 几何级数增长的核酸和蛋白质序列、结构的信息，我们需要多学科工作者的共同 努力，于是一门新兴的交叉学科生物信息学( b i o i n f o r m a t i c s ) 应运而生。 生物学不再是仅仅基于观察和实验的科学。正如本文开头所引用的诺贝尔 奖得主w a l t e rg i l b e r t 在关于生物学研究模式的转变一文中所说的，生物学 研究的出发点将转变为理论的，科学家从理论推测出假定，然后再回到实验室去 验证这些理论和假定。现今，生物信息学研究者们以核酸、蛋白质等生物大分子 数据库中的资料为主要对象，以计算机的硬件、软件和网络为主要工具，以数学、 物理、计算机科学为主要手段，对庞大的原始数据和资料进行收集、分析、整理 和解释，使它们成为具有明确生物意义的生物信息。已构建的数据和信息系统， 例如基于关键词的数据库查询系统p d b ( p r o t e i n d a t ab a n k ) 1 2 ，基于序列相似性比 对的数据库搜索b l a s t ( b a s i cl o e a la l i g n m e n ts e a r c ht 0 0 1 ) 1 3 1 等，为进一步的研究 工作，如分子模拟和设计、药物设计和基因组序列外显子、内含子、启动予识别 等创造了良好的信息平台和交流环境。 生物信息学的研究重点主要体现在基因组学( g e n o m i c s ) 和蛋白组学 ( p r o t e o m i c s ) 两方面。在人类基因组图谱公布以后，研究重点从基因转向了蛋白质。 为什么人类的基因仅比果蝇这样的低等生物多两倍，却有如此复杂的生命活动 呢? 蛋白质就是答案。只有搞清楚蛋白质的结构和功能，揭示出其在生命活动中 的作用，才能进一步开发出诊断疾病的新方法和治疗疾病的新药物，这项研究对 于全人类的健康和制药产业的发展都有非凡的意义。 诚然，蛋白质的研究方法有很多，为什么如今分子动力学( m o l e c u l a r 3 引言 d y n a m i c s ，简称m d ) 模拟会成为最基础和重要的方法之一呢? 最根本的原因是， 蛋白质的功能主要决定于其结构，而分子动力学是研究物质结构的学科，为蛋白 质研究提供了有效的手段。其优势体现在如下三个方面。第一，随着更多更大更 复杂的生物大分子结构被测定和存入网络数据库，计算能力更强的超级计算机正 遵循着摩尔定理( m o o r e sl a w ) 不断发展，并且各种并行分子动力学软件也在不断 更新。这为研究蛋白质这样的复杂生物大分子提供了良好的硬件和软件环境，能 够计算的分子总数越来越大，方法越来越好( 势函数的改进) ，模拟的时间 f t i m e s c a l e ) 越来越长。第二，分子动力学模拟能够提供粒子运动随时间变化的最 详尽的信息，甚至能回答实验上所不能回答的问题。例如，使蛋白质正确折叠的 伴随体( c h a p e r o n ) g r o e l 在发挥功能时，有一个结合a t p ( 三磷酸腺甙) 的过程。 实验上通过x 射线晶体衍射只能测到“开”和“关”两个状态，但是整个进入的 过程是未知的【4 j 。然而，在分子动力学模拟的帮助下，我们可以观察整个过程【5 j 。 第三，模拟是一个人为控制的过程，不会受到实验上一些条件的限制。我们曾经 通过替换一种酶t c 2 3 0f t h e r m o s t a b l ec a t h e c h o l2 ，3 d i x o y g e n a s e ) o p 的金属离子和 若干残基( r e s i d u e ) 来证明它们对其热稳定性( t h e r m o s t a b i l i t y ) 的贡献，而其中有些 替换在实验上是很难实现的。这好比是一种“电脑炼金术( c o m p u t e ra l c h e m y ) ”， 可以根据需要改变结构，切割序列，实现低投资高效率的新分子搜寻，这也正是 分子药物设计的初衷。当然，并不是说分子动力学模拟可以完全脱离了实验，毕 竟分子动力学所用的势场还需要不断的改进。只有在和实验结果的对比下证明模 拟是可信的，其提供的详细的信息才是有意义的。 分子动力学帮助我们从原子分子水平去探索生命现象的奥秘，这列于生物 学的研究是跨出了一大步。然而，在大多数生物化学过程中，如光合作用、酶促 反应、生物体内信号传递，都与电子传递有关，这就需要从电子的层次来认识生 命活动。而且，生物大分子是高度组织的分子体系，在生命活动中是以整体来作 用的，每部分之间都互相影响，所以对蛋白质做整体量子化学计算就显得十分的 重要了。计算生物大分子的整体电子结构，这跨出了更大的一步。传统的量子化 学计算对象，往往仅由几十个原子构成，但是计算量已经非常大了。即使现在的 计算机发展速度很快，遇到蛋白质这样的高分子，常规的量子化学计算方法仍然 难以胜任。比如，猪胰岛素是最小的蛋白质分子之一，只有7 8 2 个原子和3 0 7 6 个电子，采用极小基集合做从头计算，需要2 4 1 8 个基函数，由此可见计算量是 惊人的。因此，对于蛋白质这样的生物大分子有必要利用新的计算方法，使这种 整体量子化学计算在当前的计算机上可以实现。 本文谈及的两个课题就用到了并行计算分子动力学模拟和全电子结构计算 ( w h o l ee l e c t r o n i cs t r u c t u r ec a l c u l a t i o n ) 。首先，第一章中简单介绍了蛋白质的基本 4 引言 概念、并行计算分子动力学的方法和全电子结构计算的方法。第二章中选取了蛋 白质的动力学相变( d y n a m i c a lt r a n s i t i o n ) 作为研究的首要对象，因为这是蛋白质的 一个普遍性质，并且对于它和蛋白质功能活性之间的关系还存在着争议。通过对 方均位移( m e a ns q u a r ed i s p l a c e m e n t ) 以及一系列自由能非线性指标( n o n l i n e a r i n d i c e s ) 的分析，模拟重现了肌红蛋f j ( m y o g l o b i n ，简称m b ) 和铜锌超氧化歧化酶 ( c u z ns u p e r o x i d ed i s m u t a s e ，简称s o d ) 准弹性中子散射实验的结果。而且，研 究还涉及了水环境对动力学相变的影响和相变与酶活性的关联。第三章选取了限 制性内切酶( r e s t r i c t i o ne n d o n u c l e a s e ) b a m h i 的一种热敏感性突变体( h e a t s e n s i t i v e m u t a n t ) p 1 7 3 l 作为研究对象，结合了并行分子动力学模拟、全电子结构计算方 法( 其中使用了溶液模拟和重叠二体近似) 和蛋白质表面拓扑结构计算，验证了 实验中指出的活性位点( a c t i v es i t e ) 。并且，通过4 个活性位点所占有前线轨道 ( f r o n t i e r o r b i t ) 的能级变化，不但验证了实验中给出的活性对温度变化的规律，还 给出了其它的预测。 5 第一章预备知识 生物信息学是一个高度交叉学科，需要运用生物、数学、物理、化学和计 算机等各方面的知识。因此，在进入我们的两个研究课题之前，有必要首先对课 题的研究对象和研究手段进行简单的认识，作为以后章节中将出现的一些重要概 念的预备知识。内容包括氨基酸和蛋白质的基本结构和种类、并行计算分子动力 学、蛋白质全电子结构计算方法。 1 1 氨基酸和蛋白质 蛋白质是细胞的重要组成部分。它在细胞中执行两种主要功能：一种是作 为在活体组织内催化生物反应的酶( e n z y m e ) ，另一种是作为细胞中特殊的结构， 例如细胞壁。若按照其外形和在生物组织里的作用，可粗略地分为纤维蛋白 ( f i b r o u s p r o t e i n ) ，膜蛋( m e m b r a n e p r o t e i n ) 和球蛋白( 9 1 0 b i n ) 三大类。其中球蛋白 的种类最多，包括了很多催化生物化学反应的酶，本文中涉及的两个主要蛋白质 都属于球蛋白。 氨基酸 氨基酸的结构如图1 1 ( a ) 所示，它的中心有一个碳原子，称为a 碳原子( c 。) 。 c 。的四个化学键，一个连着氢原子h ，一个连着氨基n h 2 ，另一个连着羧基 c o o h ，最后一个连着侧链r 。根据侧链r 的不同，形成了2 0 种不同的氨基酸 ( 见表1 1 ) 。虽然自然界和实验室里合成的氨基酸不止这个数目，但所有的蛋白 质都是由这2 0 种氨基酸组成的。 r 侧链 r ，r z l i 一 n h 。- - c r c o o h+ h 。卜阱c o n h + _ 鼯c 0 0 一 氨基 l 羧基i。 i hh 肚键 h ( a )( b ) 图1 i 结构示意图。 ( a ) 氨基酸：( b ) 肽键。 - 6 - 第章预各知识 表1 12 0 种氨基酸的符号和名称 符号( 简略符号)中文名称符号( 简略符号)中文名称 a l a ( a )丙氨酸l e u ( l )亮氨酸 a r g ( r ) 精氨酸 l y s ( k ) 赖氨酸 a s n ( n ) 天冬酰胺m e t ( m )甲硫氨酸 a s p ( d ) 天冬氨酸p h e ( n苯丙氨酸 c y s ( c ) 半胱氨酸 p r o ( p )脯氨酸 g i n ( q ) 谷氮酰胺 s e r ( s )丝氨酸 g l u ( e )谷氨酸t h r ( t )苏氨酸 g l y ( g ) 甘氨酸 t r p ( w ) 色氨酸 h i s ( h )组氮酸 t y r ( y ) 酪氨酸 n e ( i )异亮氨酸 v a l ( v )缬氨酸 蛋白质的结构 相邻氨基酸的氨基和羧基脱去一个水分子就形成了肽键( p e p t i d eb o n d ) ，见 图1 1 ( b ) 。由肽键连接而成的氨基酸链( 多肽链) 称为蛋白质的一级结构( p r i m a r y s t r u c t u r e ) ，所有的c 。原子构成了蛋白质的骨架( p r o t e i nb a c k b o n e ) 。 多肽链的片断在氢键( h y d r o g e nb o n d ) 的作用下，局部折叠成稳定的结构，称 为蛋白质的二级结构( s e c o n d a r ys t r u c t u r e ) 。最主要的二级结构包括o 【螺旋m h e l i x ) 和1 3 折叠片( d p l e a t e ds h e e t ) ，其它的都属于无规则卷l 益t ( r a n d o mc o i l ) 。 蛋白质的三级结构( t e r t i a r ys t r u c t u r e ) 目p 蛋白质分子多肽链的三维空间构象。 一般三级结构的形成，使蛋白质在热力学处于最低自由能的结构稳定状态。 四级结构( q u a t e r n a r ys t r u c t u r e ) 仅存在于由多条多肽链构成的蛋白质中，实际 上是三级结构的亚单位彼此在三维空间中的相对位置关系。 1 2 并行计算分子动力学 分子动力学方法的应用已有相当长的历史。其基本思想虽起源于v a nd e r w a a l s 和b o l t z m a n n ，但真正有意义的工作开始于上世纪5 0 年代a l d e r 和 w a i n w r i 曲t 以刚球模型相互作用为基础的工作【刚。1 9 6 4 年，r a h m a n 的文章将软 球模型( l e n n a r d - j o n e s ) 应用到分子动力学o p 7 1 ，是一个重要的突破。紧接着，是 关于复杂的液体系统的分子动力学模拟【8 1 和陆续展开的许多对同质系统 ( h o m o g e n e o u ss y s t e m ) 的平衡或非平衡行为的研究 g j o l 。这些前期工作为分子动力 学在生物分子上的应用做了很好的铺垫。1 9 7 7 年，m c c a m m o n 等人首次成功使 用分子动力学模拟了生物分子系统，模拟对象是牛胰岛素抑制n ( b o v i n e 7 第一章预各知识 p a n c r e a t i ct r y p s i ni n h i b i t o r ，简称b p t i ) 川。虽然用的势函数比较粗糙，计算时闯 也很短，但对蛋白质这样的复杂系统是第一次有意义的尝试。从首次对蛋白质的 分子动力学模拟到今天，已经经历了大约2 7 年的时间。分子动力学在生物分子 研究上的应用范围也愈来愈广，不但可以研究蛋白质构象变化、稳定性变化等动 力学过程，还可以辅助x 射线晶体衍射和核磁共振( n m r ) 来确定蛋白质结榭1 2 1 。 目前网络上可免费下载的分子动力学软件有很多，在我们的工作中使用的 是依里诺依大学( u f u c ) 的生物物理组开发的并行分子动力学软件n a m d l 】。所 以，以下对于分子动力学的原理、算法等的介绍，都是n a m d 中使用的。 特别要指出，用于n a m d 的输入文件和输出文件的制作、处理软件都是本 研究组开发或搜索的，n a m d 本身并没有带这些程序。 v e r l e t 算法 分子动力学的基础就是牛顿力学，简单地归纳为以下的公式 粤：f f ：v i - v , u 叫，2 其中，m ，表示第f 个粒子的质量，r 表示其位置，f i 表示其所受的力，u 是整个 系统的势场函数。由 两式相加，可得到 r ( f + 占r ) * r 。) + v ( r ) j f + 三a ( 0 8 f 2 r ( t - 万，) zr o ) 一v ( f ) j r + 三a 。) 万f 2 r ( t + 8 0 2 r ( t ) 一r ( t 一8 t ) + a ( t ) s t 2 ， ( 1 1 ) 这就是v e r l e t 算法的公式，通过粒子在t 时刻的位置、加速度和r 占t 时刻的位置， 就可以推出f + j f 时刻的位置。所以，一旦知道了所有原子的初始位置和势场函 数，就可以算出以后任意时刻的原子位置。 c h a r m m 势场 接下来简单介绍一下n a m d 软件所使用的势场c h a r m m 势【1 4 】。 分予动力学的模型将分子描述为由弹簧( 键) 连接的带电质点，用势场来描述 键长、键角、转矩、非键范德华作用和静电作用等对原子能量的影响，表示为： 一8 ，【 第一章预备知识 其中 u = u b 皤+ u 渊目+ ud i i l e + u 嘶口r + u 吣q b a 叫， u 唰= ( 只- 0 0 ) 2 u m 。= 群“( 1 + c o s ( n ，谚+ 一) ) 2 ( ，= 爿”( 。，一。) 2 = 莩利一删+ 莩善弩 前四项分别对应于图1 2 中的四种运动模式，最后一项表示非键的范德华作用、 氢键和静电作用。 并行计算 图1 2c h a r m m 势中包括的四种键作用 n a m d 采用以c + + 为源码的c h a r m + + 语言( 含m p i 库) 。并行计算的基 本单位是称为c h a r e 的类，在程序运行过程中，各种c h a f e 类之间互相通l 乳以实 现并行处理完成任务。 整个c h a l w + + 语言由s e q u e n t i a lo b j e c t s ，c h a r e s ( c o n c u r r e n to b j e c t s ) b r a n c h e dc h a r e s ( af o r mo fr e p l i c a t e do b j e c t s ) ，s h a r e do b j e c t s 和c o m m u n i c a t i o n 9 - 第一章预备知识 o b j e c t s 五个库目录组成。 具体的分子动力学模拟方法将在第二章中作介绍。 1 3 蛋白质全电子结构计算 为了进一步地了解生命活动，除了分子动力学提供的原子水平的认识，电 子结构的计算也十分重要，因为许多生物过程都与电子传递有关。 此处所采用的软件早期由l a d i k 教授领导的小组开发【l ”】并运用到生物大 分子体系，后来由美国康乃尔大学叶元杰教授改进并用于研究生物分子的活性和 电导率 1 8 , 19 】。此软件需经过整理和改进，并且配合我们编写的大量处理程序以后， 才成为一个完整的计算蛋白质整体电子结构的软件。若稍加改进，亦可开发成第 一原理的计算。它使用的计算方法是f o c k 矩阵方块对角化和e n f c 本征值计算 方法，计算速度较快，且用加电荷的方法解决了溶液环境的模拟问题。 数值计算方法 假设一个蛋白质分子由个氨基酸残基构成，可以把每一个残基作为一个 基本单元。每个残基使用m 个基函数参与整体计算，即第，z 个单元使用m 。个原 子轨道p ，( ”) ，( 1 n n ，1 l ，兰m 。) 。由l c a o 近似( 把分子轨道写成原予轨道 的线性组合，每个参与共轭的原子提供一个原子轨道) ，整个分子体系的分子轨 道可以写成所有单元使用的基函数的线性组合： n y = q ( 一) 竹( 阼) ， n = lj = l 其中，一是单元序号，表示原子轨道序号，纺( 功表示第疗个单元中的第，个原 子轨道。由线性变分法，可以导出h a t r e e f o c k - r o o t h a a n 方程： h c = s c e ， ( 1 2 ) 其中，c 是由线性组合系数q ( ) 构成的分子轨道组合系数矩阵，h 是f o c k 矩 阵，s 是重叠矩阵， j 乱，= j 虻( r ) 吼p ) d r = i 以，= 协( ，) 岛，( ，) 咖= e 是由所有分子轨道能级构成的对角矩阵。 1 0 第一章预备知识 由于蛋白质是非常大的分子体系，所以公式( 1 2 ) 的维数非常大，最小的蛋 白质也要达到数千阶，计算量非常的大。考虑到蛋白质由于分子体系较大，很多 原子之间的距离也非常的大，使得很多矩阵元的值接近于零。这样，就可以采取 近似的方法，只考虑空间距离比较近的单元之间的相互作用和轨道重叠：忽略空 间距离较远的原子之间的作用，矩阵元罱零。所以，最后保留的矩阵元包括单元 内的、序列上相邻的单元之间的、序列上不相邻但三维结构中距离小的单元之间 的相互作用和重叠积分块( 比如：二硫键相联的两个半胱氨酸之间的相互作用不 可忽略) 。如此得到的h 和s 近似为如下的形式： h = 日1 1日1 2 0 2 1日2 2日2 3 - s = 0 o 以h - 0 _ 墨_ 。 日v ； h m ： o - o 以i m h 一1 0h n n 4 h n q 0 s 2 ， 最，。i 砖。吒 只砘 0 s 一 o 。 0 s n 。n s n n ( 1 3 ) ( 1 4 ) 其中，序列上不相邻、但在三维结构中距离近而不能忽略的三个相互作用和重叠 积分矩阵块，用h 作为行标，疗：( f _ 1 ，) 作为列标来标示。 接着，对本征值方程进行数值求解。将方程( 1 2 ) 改写成对于单个本征值的 本征方程，变形后得： 第一章预备知识 ( 日一a s ) c = 0 ， ( 1 5 ) 兄为本征值，c 为本征向量。通过数值方法，希望找到一个近似值五，使 【( 片一旯苫) c r 【( 日一旯s ) c 】 ”2 而) ，使 n e ( x ：) 一n e ( x 。) = 1 。也就是找到小区间( x 1 ，z 2 】，只包含一个本征值。接着，再 压缩这个小区间使l x ：一z 。 占，最后取( 而+ x ：) ，2 作为本征值的近似值。 同样，利用e n f c 定理，很容易能得到近似态密度函数 p ( e ) = 壁訾， ( 1 6 ) q ( e ) 是方程( 1 5 ) 在区间( 一m ，e 】中的本征值的个数。 溶液环境模拟 多数生物分子需要在溶液环境下才能发挥功能，所以计算蛋白质全电子结 构的时候，必须考虑溶液效应。 在水溶液中，一部分氨基酸残基是带电荷的，例如赖氨酸( l y s ) 和精氨酸( a 唱) 带一个正电荷，谷氨酸( g t u ) 和门冬氨酸( a s p ) 带一个负电荷，所有肽链n 端和c 端都会带电荷。这些带电荷的残基会吸引异种电荷到它们周围，形成双电层。所 谓的溶液环境模拟，就是在这些带电荷的残基周围放置符号相反的若干点电荷， 这些点电荷的总电荷量和氨基酸残基所带总电荷量相同。 重叠二体近似 以氨基酸为单元，将整个蛋白质分子体系按每两个单元分割成一些二体 ( d i m m e r ) 。这些二体互相重叠，尾端的单元被计算了1 次，中间的单元至少被计 算2 次有些单元可能被计算3 4 次。图1 3 中，尾端单元1 被计算了1 次；单 元3 、4 被包括在2 个二体中，所以计算了2 次；单元2 由于被一个二硫键连接， 计算了3 次；单元5 由于与另外的单元距离非常近，也被计算了3 次。 对于二体的边界条件的处理，采用的是加赝原子的方法，这些原子主要是 与二体相连的单元的原子。例如，与单元相连的肽键需要保留，非单元内的c 。 原子用氢原子代替；与单元相连的二硫键也应该保留，且另外保留二硫键外的一 个甲基。 1 3 第一章预备知识 耍耍耍耍 s 耍 s k 画 查丛 图1 3 重叠二体近似示意图 接着，对这些二体作自洽场计算。如前面所说的，由于蛋白质的分子体系 特别大，计算时涉及大量的双电子积分，它们的数量理论上是4 的数量级( n 是基函数个数) ，因此需要寻找简化计算的方法。g a z d y 等人于1 9 8 4 年指出吲， 对于自洽场方法，采用局部自洽后再构建整体f o c k 矩阵所得到的结果和直接做 整体白洽场计算所得的结果非常接近，这样可以大幅减少计算量，是较好的近似。 那么，原先需要对整个生物大分子体系所作的自洽场计算就可分解为对这些二体 所作的自洽场量子化学计算。 完成对所有二体的自洽场计算以后，就得到了它们的f o c k 矩阵和重叠矩阵。 接着用它们建立整体的特征方程( 即公式f 1 3 - 1 5 1 ) ， 巧川( ”，h ) 2 去k = l 巧川( h ，”，n 聊= 1 ，2 3 ( 1 - 7 ) 巧川( n ，聆) = 巧即( 珂) ， n ( 1 8 ) 上标( m ) 表示整体方程的f o c k 矩阵，上标( d ) 表示二体计算所得的f o c k 矩阵。 公式( 1 7 ) 表示的是f o c k 矩阵中的对角方块矩阵元，m 表示某个单元被m 个二体 所包含，计算了m 次，需要对它们作加和平均后得到整体f o c k 矩阵的相应矩阵 元。公式( 1 8 ) 对应的是非对角方块矩阵元。 这样，再利用e n f c 计算方法对构建出的水溶液中的蛋白质整体特征方程 进行数值求解，可以得到整体电子结构，包括电子能级和分子轨道波函数。 1 4 第二章蛋白质的动力学相变研究 动力学相变是蛋白质的一个基本性质。本章通过对肌红蛋白和铜锌超氧化 歧化酶的分子动力学模拟，以方均位移随温度的变化验证了实验上中子散射实验 的结果，并且进一步验证了水环境对这个动力学相变的关键作用。通过不同残基 的方均位移的不同变化程度，讨论了相变和活性的关系。还采用了一系列非线性 的指标，对自由能数据处理，给出动力学相交的另一种定量描述。最后，给出双 阱势的单粒子模拟来验证动力学相变的原理。 2 1 背景和理论介绍 我们选取的模拟对象是肌红蛋白和铜锌超氧化歧化酶。肌红蛋白是较小的 蛋白质分子，由1 5 3 个残基和一个血红素( h e m e ) 组成，结构上由8 个a 螺旋组成。 由于它的分子数少，结构又足够复杂，具有一定的代表性，所以适合作为模式蛋 白质。铜锌超氧化歧化酶也是一种蛋白质，铜离子有稳定蛋白质结构的作用，而 锌离子能起到催化酶反应的作用。铜锌超氧化歧化酶能促使超氧化自由基o 。分 解为氧气和过氧化氢，而过氧化氢会进一步分解为氧气和水。所以它对于细胞的 存活和人类的健康、寿命都有很重要的生物学意义。 中子散射实验 d o s t e r 等人于1 9 8 9 年和1 9 9 0 年在自然杂志发表了关于肌红蛋白中子 散射实验的结果 2 3 2 4 1 ，发现在2 0 0k 附近，蛋白质的运动模式会从谐振变为非谐 振，并建议可以使用分子动力学的模拟来验证和研究这个现象。1 9 9 5 年，生物 物理杂志上发表了a n d r e a n i 等人对铜锌超氧化歧化酶所做的中子散射实验【2 5 1 ， 同样也发现了2 0 0k 的相变点。在2 0 0k 以下，纯弹性散射，2 0 0k 以上纯弹性 散射大幅下降。之后，陆续有很多关于其它的蛋白质的中子散射实验都证明了蛋 白质在1 8 0 2 2 0k 附近，其内部运动会经历一个所谓的“玻璃相变( g l a s s t r a n s i t i o n ) ”【2 3 2 “，表现为原子方均位移 随温度的突变。 中子散射实验使我们能研究蛋白质在o 1 1 0 0 皮秒时间范围内的运动，这 正是i v i d 模拟普遍能够达到的时间范围。由于氢原子核的中子散射截面异常的 大，所以这个方法主要观测到的是氢原子的运动。好在氢原子在蛋白质中非常的 多，而且分布很广且均匀，因此中子散射实验提供了关于蛋白质整体运动的信息。 第二章蛋白质的动力学相变研究 实验上的观测量是不连续动力学结构因子( i n c o h e r e n td y n a m i cs t r u c t u r ef a c t o r ) s ( q ，t o ) ，其中壳q 和a 分别对应于系统和入射中子之间交换的动量和能量。s ( q ，t o ) 是系统密度涨落( d e n s i t yf l u c t u a t i o n s ) 时间关联函数( t i m e c o r r e l a t i o nf u n c t i o n ) 的傅 里叶变换，可以直接从m d 模拟的结果中求出。当温度较低时( 低于2 0 0k ) ， 蛋白质的弹性强度( e l a s t i ci n t e n s i t y ) 呈g a u s s i o n 分布，它的谐振运动可以用 d e b y e - w 乱l e r 因子来描述，s ( q ，* o ) = e x p 卜 q 2 ) ，其中方均位移 与温 度呈线性关系。当温度升到2 0 0k 左右的时候，弹性强度急剧减小，表示有新的 自由度被激发 2 7 , 2 8 。 蛋白质动力学相变原理 不只是中子散射实验，还有很多光谱学的实验1 2 9 , 3 0 1 也发现了这个现象。那 究竟是什么原因造成这个动力学相变的呢? 1 9 9 8 年，d o s t e r 在关于“生物学中 水的作用”的专题讨论会上提出：蛋白质动力学相变来源于氢键的作用1 3 “。其 分别来自于分子内部的氢键，和蛋白质分子与水溶剂之间的氢键，后者的数目相 当庞大，因此占有主导地位。随着温度的升高，原来被氢键锁住的氨基酸侧链的 振动增强，使得蛋白质在亚能级构象上布居数变化，造成了相变；另一方面，蛋 白质表面的水分子层和蛋白质之间形成的氢键网络重组更是重要原因。在相变温 度以下，水分子的平动和转动都受到了限制，因此也限制了蛋白质分子内部运动： 但温度升高到一定的程度之后，它们的运动开始加剧，使得蛋白质的构象变化成 为可能。 实验上，改变溶剂浓度后，发现动力学相变后方均位移的变化幅度的确不 同口”。然而，以此证明“溶剂推动”造成了动力学相变还是不够的。2 0 0 0 年， v i t k u p 等人做了一个很有趣的m d 模拟 a 2 1 。将蛋白质和水溶液的温度分别设置 在3 0 0k 或1 8 0k ( 高于或低于动力学相变的温度2 2 0k ) ，这样就有了4 种组合 的结果( 表2 1 ) 。结果很令人惊讶，蛋白质本身的温度对方均位移的贡献很小， 反而溶剂的温度是主要决定因素。这样一个超乎自然规律的实验，给了“溶剂推 动”理论一个强有力的证明。 表2 1 平均方均位移 ( p 表示蛋白质，s 表示溶剂) 1 6 第二章蛋白质的动力学相变研究 至于蛋白质动力学相变和其活性的关系，在过去的近1 0 年内也备受关注， 但至今没有找到一个确切的答案。一些文章指出这个相变温度就是蛋白质活性开 始的标志 3 3 - 3 5 ，但是一些实验结果却和这个论点背道而驰，因为实验发现在相变 温度下，蛋白质的活性仍然存在【3 6 】。那么，这个动力学相变还值得进一步的研 究，所以我们选取了它作为模拟的首要对象。 2 2 动力学模拟过程 我们使用的分子动力学软件是n a m t ) ，但是n a m d 本身并没有带输入文件 制作和输出文件处理的程序。所有这些输入文件制作、数据处理软件都是本研究 组搜索或开发得到的。下面先介绍一下运行蛋白质分子动力学的具体过程。 蛋白质分子动力学运行过程 如图2 1 所示，整个运行蛋白质分子动力学的过程，简单地说分为三个步骤： 输入文件的制作、主程序n a m d 的控制运行和输出文件的数据处理。 首先，是输入文件的制作。输入文件包括p d b 文件、p s f 文件、c o n f 控 制文件、拓扑文件( + t o p ) 和参数文件( + p a r a m ) 。其中，拓扑文件和参数文件都可 以直接从网上下载，由c n s l 3 7 】提供的。需要制作的是p d b 文件和p s f 文件。先 从p d b 网络数据库【2 j 中取出原始的蛋白质p d b 文件，其中记录了所有原子的x 射线晶体衍射三维坐标。然后，用p s f g e n 软件补全晶体衍射实验测量中的遗 漏，产生完整的蛋白质三维坐标p d b 文件和记录所有键长键角等参数的p s f 文 件。接着，用d o w s e r 软件口9 1 和s o l v a t e 软件1 3 9 】在蛋白质分子周围加上水分 子包围层和若干n a + ，c 1 - 离子，调节p h 值。这样就完成了n a m d 所需要的输 入文件。 控制n a m d 的运行是通过编写c o n f 控制文件实现的，它由t c l 语言编 写。文件的内容包括一些参数的设定( 如范德华力的截断距离c u t o f f 、记录数 据的密度等) 和对整个分子动力学过程的控制。n a m d 的模拟过程也如图2 1 中所示。通常，先要经过能量最小化( e n e r g ym i n i m i z a t i o n ) 的过程使读入的蛋白质 结构达到能量尽可能低的状态，越接近0k 越好，但实际情况一般是降到几十k 。 随后，根据需要控制温度逐步上升，且每一次上升r 后，系统要经过足够长的 时间才能达到平衡，但是太长的时间，上升的温度会超过r 。所以，上升每个 r 运行多少时间，是需要多次试验观察才能确定的。另外，公式( 1 1 1 中占t 也就 是时间步长( t i m e s t e p ) 的值直接影响到计算的准确性和速度，是如何确定的昵? 分 子动力学受到分子内最高振动频率的限制，最理想的时间步长应该取最高振动频 率的十分之一。多数情况下，蛋白质分子里的c h 键振动频率最高，约为 1 7 第二章蛋白质的动力学相变研究 输出文件的处理 图2 1蛋白质分子动力学运行过程流程图。 1 8 件 第二章蛋白质的动力学相变研究 1 0 _ 1 4 秒，所以时间步长通常都取为l 飞秒( 1 0 _ ”s ) 。 最后是输出文件的数据处理和分析。n a m d 分子动力学软件产生的输出文 件包括d c d 文件、l o g 文件、v e l 文件( 记录速度) 和c o o r 文件等。其中， 最主要的是d c d 文件和l o g 文件。d c d 文件是二进制格式，记录了各个时刻 的所有原子的三维坐标；l o g 文件是文本格式，记录了各个时刻总体系的自由 能和温度。利用可视化软件( 如v m d ) 可以直观地观察整个动力学模拟过程中蛋 白质的结构变化；通过计算方均位移和自由能的数据的一些指标计算，可以了解 蛋白质内部分子振动状态的变化；当然也可以结合电子结构计算的软件，进一步 研究蛋白质的电子能级的变化。 本课题的模拟过程 首先从p d b 数据库取出蛋白质的初始三维坐标。我们选取了1 h j t ( 抹香 鲸的肌红蛋白) 1 4 0 作为模拟对象，其结构是利用x 射线散射在2 9 8k 温度下测 定的，分辨率为1 7 0a 。1 h j t 包含由1 5 3 个残基组成的蛋白质和血红素( 共1 5 4 0 个原子或原子集团) ，在加入了3 个n a + 、1 5 个c 1 - 和2 9 5 3 个水分子后，原子总 数达到了1 0 4 1 7 。铜锌超氧化歧化酶的初始三维坐标是从