文摘

分子动力学模拟必须足够长的得出可靠的结论。然而,没有方法的存在是为了证明一个模拟聚合。我们建议的方法“滞后RMSD-analysis”作为一种工具来判断如果一个MD模拟尚未运行足够长的时间。分析是基于RMSD值之间对配置变量Δ时间间隔隔开t。除非RMSD(Δt)已经到了一个固定形状,模拟尚未聚集。

1。介绍

数学建模和计算生物学已经被证明是非常成功的在努力了解免疫反应的机制,例如,hepatitisC后对t细胞增殖动力学建模,艾滋病毒感染(1),或一个免疫原性肿瘤的生长2]。统计力学的方法可以成功地应用为了定量理解免疫系统(3,4),和不同的建模方法是获得足够的免疫力,在调查Perelson和Weisbush开创性的论文5]。建模与仿真可以执行在不同的水平,从基于代理模型的顶层(6,7)合作的众多大型生物分子在免疫突触的形成8- - - - - -11)到更详细的抗原结合的模型(12)、识别和信号(13- - - - - -20.]。通过interleukin2 t细胞增殖调制模拟(21]。应用程序旨在预测抗原表位的反应22)提高疫苗设计(23),接种肿瘤(24),改善病人护理已经设计了(25- - - - - -27),也在个性化医疗(28]。特别是机制如何t细胞受体(tcr)检测抗原肽(p)提出的主要组织相容性复合体(MHC)分子可以通过分子动力学(MD)研究在分子甚至原子水平(29日- - - - - -32]。然而,TCRpMHC复合物是巨大的蛋白质复合物,必须研究在水中,添加了一个更大数量的原子模型和模拟。由于这些原因,足够长的模拟是强制性的为了获得现实的结果。这是更是如此,因为分子识别现象可能操作相当长的时间尺度比通常的MD轨迹的长度。因此,努力在这样的模拟评估抽样质量是强制性的。

考虑到实际的MD模拟轨迹,很明显从第一原则,没有任何正式的检查可以证明完成抽样取得。如果部分相空间尚未访问,没有这个可能性变得明显从只看数据从那些已经访问的其他部分。Zhang et al。33)提出一个正式的方法来分解相空间到泰森多边形法多面体(34,35]。

无论如何,正式的抽样检查质量只能利用仿真数据已经产生在最检测可能的缺点的轨迹。情况类似于测试的随机性伪随机数发生器(36]。一个永远不能证明随机性。只是可能检测具体偏离随机性,如串行相关性高、任何此类检测只能预选的错误率。

同样,MD轨迹可能追究nonrandomness的特定标记,最重要类型的能量和分子变形趋势。这项工作侧重于后者,偏差在形状量化通过均方根偏差(表示)37在时间 2 对于一个给定的参考结构 1 : R 米 年代 D 1 , 2 = 1 = 1 为 为 2 − 1 为 为 2 1 / 2 , ( 1 ) 在哪里 ( )是原子的位置 在时间 和 是分子中原子的总数。

通常,第一帧的轨迹 ( 1 ) 作为参考,值吗 R 米 年代 D ( 1 , 2 ) 计算连续所有 ( 2 > 1 ) 帧;看到黑色的曲线在图1。RMSD监控这种方式显示大型快速波动的长期变化和跳跃。一般来说,难以识别这种长期的变化,这可能与功能模式。最近的一项研究[38反映了视觉RMSD单独检验的不足。


图1
表示为时间的函数trajectory1 (10 ns)。黑色曲线:表示之间首先配置和连续的。蓝色曲线:对于一个给定的配置 ( ) ,表示每个其他配置相同的轨迹计算,平均收益率 R 米 年代 D ( ) 。这个值的 的 R 米 年代 D ( ) 然后采用最大作为参考吗 ( 1 ) 整个情节RMSD值轨迹。红色曲线:蓝色曲线,但最少 R 米 年代 D ( ) 。注意RMSD之间的参考配置和本身被定义为零。

广泛的工作已经发表在使用RMSD结构性变化的表征,飘,和趋势;参见[33)和参考引用。Grossfield和扎克曼39]给出了一个具有开创性的概念讨论的遍历性,绝对和相对收敛,提出了整体抽样检查质量。提出了块平均(40减少短期波动和获得更可靠的长期趋势指标。然而,平均RMSD在块内的所有配置涉及到许多不同的时间间隔 Δ = − 的特异性,从而减少由此产生的一个特定的时间间隔的平均水平。

因此我们认为RMSD双配置之间相隔一个常数时滞 Δ 节中描述2.2。1。以这种方式得到更稳定的估计(平均),然而完全具体的对于一个给定的时间间隔 Δ 。

2。材料和方法

2.1。分子动力学模拟
2.1.1。就业结构

我们提出的方法应用于共6 MD模拟。为此我们使用两种不同TCRpMHC复合物:免疫野生型肽细胞之间的绑定/ MHC和相同的细胞/免疫原性突变肽MHC的少。我们使用的晶体结构LC13 TCR绑定HLA-B * 08: 01和eb病毒肽FLRGRAYGL。这种结构可以从蛋白质数据库(PDB) [411军情五处[]通过PDB加入代码42),被称为野生型。突变的复杂我们酪氨酸的侧链取代位置7肽的丙氨酸(Y7A)。这是使用SCWRL[执行43)因为之前我们可以表明,这个工具是最适合突变pMHC复合物(44,45]。

复杂的细胞受体)LC13 FLRGRAYGL肽和HLA-B * 08: 01以来分子动力学模拟是一个理想的测试集描述的这个复杂的结晶,其部分和作为一个整体。最初Kjer-Nielsenetal。结晶细胞(46)和MHC (47)分别发表了一整个TCRpMHC系统结构时(42之后)。这些数据给大量的洞察LC13 / EBV / HLA-B * 08: 01系统,导致我们的野生型和突变系统的选择研究。整个系统见图2。


图2
在复杂的插图的LC13 TCR HLA-B * 08: 01。蓝色:MHC、红:β2-microglobulin,灰色:肽,橙色:TCR alpha-chain,黄色:TCR beta-chain。
2.1.2。分子动力学模拟协议

野生型和突变体复杂模拟独立运行10、50、200 ns的共有6模拟(见表1)。模拟采用以下协议。首先,我们最小的能量系统使用最速下降法。然后我们沉浸复合物在显式程控48)允许的最小距离的人工水浴锅2 nm盒子边界和蛋白质。接下来,我们温暖的复杂使用位置限制310 K。最后,我们进行了模拟使用GROMACS4 [49)和GROMOS96力场(50]。依照(所有进一步的参数设置51]。

表1
分子动力学模拟运行。
2.1.3。表示计算

模拟完成后,我们计算了每个配置RMSD值在给定轨迹对其他使用标准配置相同的轨迹g_rmsGROMACS的函数。这产生一个 × 矩阵表示的值 是配置的数量的轨迹。

2.2。RMSD配置之间的轨迹
2.2.1。平均滞后RMSD和建模

给定一个配置x( )的MD模拟作为参考,RMSD其他配置x( )可能被视为“距离测量时间间隔 Δ = − 。如果 Δ t是足够短,这可能是沿着整个转移轨迹和RMSD值采样。平均 R 米 年代 D ( Δ ) 特点是配置之间的区别吗 Δ t在特定的模拟运行。

小的值 Δ 描述配置接近对方。增加 Δ 意味着彼此比较配置更遥远,直觉上说——“不相关的”。这一事实应该反映在 R 米 年代 D ( Δ ) 增加的值 Δ 。

作为 Δ 增加,依赖性应该减少和方法的水平“不相关的”或“独立”配置。为了量化这种饱和趋势,我们应用了希尔方程(52]: R 米 年代 D ( Δ ) = ⋅ Δ + Δ , ( 2 ) 的参数 反映了最大值的函数渐近(“高原价值” R 米 年代 D ( Δ → ∞ ) ), 是时间滞后 Δ 的 R 米 年代 D ( Δ ) = / 2 (即。,the value of Δ 实现half-saturation),希尔系数 是一个参数,它决定了形状,也就是说,sigmoidicity水平,模型的功能。被拟合估计的参数希尔方程(2)的测量值 R 米 年代 D ( Δ ) 使用nlinfit函数实现的统计工具箱MATLAB (Mathworks纳蒂克,妈,美国)。的最长时间滞后 Δ 选择相应的轨迹仿真总时间的一半。

2.2.2。评估初始条件的影响

最初的阶段 ≤ o ff 年代 e t 每个MD模拟的强烈依赖于配置开始,通常是晶体结构,决不可以代表一个配置从一个轨迹。即使如此能量最小化和热身之后。因此,如果MD的初始阶段轨迹包含在 R 米 年代 D ( Δ ) ,则会导致偏见。这不仅是对个人的价值观 R 米 年代 D ( Δ ) 但也从合适的参数估计。特别是,限制高原价值 = R 米 年代 D ( Δ = ∞ ) 也会有偏见和依靠 o ff 年代 e t :事实上, = ( o ff 年代 e t ) 。我们模仿这种依赖 o ff 年代 e t = 0 + ⋅ e x p − ⋅ o ff 年代 e t , ( 3 ) 在哪里 0 是极限值, 0 = ( o ff 年代 e t = ∞ ) ,β和 缩放参数。请注意, 0 是一个外推估计 R 米 年代 D ( Δ = ∞ ) 和 o ff 年代 e t = ∞ 的估计,即表示两个完全不相关的配置之间的轨迹,独立于初始阶段的影响。

值 o ff 年代 e t 选择在间隔吗 0 ≤ o ff 年代 e t ≤ 米 一个 x / 2 ,在那里 米 一个 x 表示总各自的轨迹的仿真时间。

3所示。结果与讨论

3.1。收敛RMSD随着时间滞后

每个面板的图3显示一个代表性的情节(出于简洁性我们只显示结果的数据轨迹1 - 3)意思表示的值(红圈, 设在)配置之间获得一个MD轨迹(参见图标题和表1),由各自的时间滞后 Δ ( 设在)。之间的时间间隔配置的增加,意味着RMSD方法高原。

(一)
(一)
(b)
(b)
(c)
(c)
图3
依赖的意思表示的时间滞后 Δ 为 o ff 年代 e t = 0 。(一)10 ns trajectory1, (b) 50 ns trajectory2, (c) 200 ns trajectory3;见表1。注意,由于不同的模拟的总持续时间,最长时间滞后认为(等于一半的仿真长度)不同。竖线指示的时间间隔 ,half-saturation实现;参见(2)。水平线表示估计高原的意思表示时,相应的参数 在(2)。

适合的模型(2)的值 R 米 年代 D ( Δ ) 显示为实线。从适合获得参数随时可以解释如下。参数的估计 在(2)代表的极限值 R 米 年代 D ( Δ ) 并由水平线表示。参数的估计 对应于一个“特色”(“half-saturation”)的时间间隔 Δ 和显示为一条垂直线的阴谋。

请注意,每个轨道的初始阶段确定的形状 R 米 年代 D ( Δ ) 应正确安装所代表的模型。相比之下,其余的轨迹是长期趋势的特征 R 米 年代 D ( Δ ) 。虽然显示了更小的变化RMSD,其余部分包含更多的数据点相对于自然的初始阶段。达到一个适当的平衡,我们增加了滞后长度在小的步骤在最初的阶段(参见红圈图的最初浓密的继承3),在较大的步骤(更宽松的圆圈)。这个过程使体重增加数据点的初始阶段。

3.2。偏移量的影响

应用不同的时间偏移量 o ff 年代 e t 在开始分析各自的轨迹变化的依赖 R 米 年代 D ( Δ ) 作为时间的函数滞后 Δ ;看到典型的结果显示在图4。

(一)
(一)
(b)
(b)
图4
依赖的意思表示 o ff 年代 e t trajectory1 (10 ns)。(一): o ff 年代 e t = 0 ps, (b): o ff 年代 e t = 2 0 0 ps。竖线表示时间滞后 ,half-saturation实现;参见(2)。水平线表示估计高原的意思表示时,相应的参数一个在(2)。

更大的偏移量 o ff 年代 e t 从一开始的配置,时间间隔越小 Δ 必要的系统其平整 R 米 年代 D 高原。的模范显示拟合参数的依赖性 o ff 年代 e t 正如先前所显示的,系统地分析如下;参见图5。随着 o ff 年代 e t 高原值 和half-saturation参数 减少,而形状参数 是相当恒定的,几乎独立的 o ff 年代 e t 。

(一)
(一)
(b)
(b)
(c)
(c)
图5
增加 o ff 年代 e t 的变化拟合参数 R 米 年代 D ( Δ ) ,(2)。(一)10 ns轨迹,(b) 50 ns轨迹,(c) 200 ns轨迹;见表1。

在下一步的系统依赖推断高原 o ff 年代 e t 通过安装(3);参见图6。(monoexponential衰变模型3)代表最吝啬的选择,给出仿真结果,显示在图的形状6。

(一)
(一)
(b)
(b)
(c)
(c)
图6
适合和外推的 R 米 年代 D 高原的值增加 o ff 年代 e t 。(一)10 ns trajectory1, (b) 50 ns trajectory2, (c) 200 ns trajectory3;见表1。红圈表示高原值 ( o ff 年代 e t ) 误差线表示其渐近标准错误。固体蓝色曲线显示了一个非线性最小二乘的3)到高原的价值观 ( o ff 年代 e t ) 。从后者符合限制高原价值 0 = ( o ff 年代 e t = ∞ ) 是提取并显示为一个坚实的水平线及其渐近获得合适的(标准误差3)。

数据6(一)- - - - - -6 (c)说明了补偿的影响 o ff 年代 e t 从一开始就在各自的配置 R 米 年代 D 高原值作为参数的估计 的模型(2)。 R 米 年代 D 高原值增加而逐渐减少 o ff 年代 e t 。

4所示。结论

蛋白质的分子动力学模拟,“收敛”和抽样的问题是重要的问题如果明智的结论可以从这样的模拟。自收敛的一个特定的仿真(在某种意义上,统计抽样的相空间足够完成对特定现象研究)无法提前判断(39建议),各种技术证明模拟尚未聚合(37,53- - - - - -56]。

在这里,我们提出一个方法来评估如果模拟不然而运行足够长的时间,基于连续RMSD分析配置给定MD轨迹,相隔时间滞后 Δ 可变长度的模拟(总数的一半时间)。只要一个模拟尚未融合,函数的形状 R 米 年代 D ( Δ ) 仍相当依赖于 o ff 年代 e t ,沿着轨迹的时间点,分析时间间隔 Δ 就开始了。作为 o ff 年代 e t 足够大的形状吗 R 米 年代 D ( Δ ) 成为固定;也就是说,在图4意思表示曲线的形状是不同的在左边和右边面板中,但进一步增加 o ff 年代 e t 不会进一步改变其整体外观。这也反映在各自的模型参数,收敛于常量值大的值 o ff 年代 e t ;参见图5。

来描述这 Δ 依赖的意思表示值, R 米 年代 D ( Δ ) 对于给定的 o ff 年代 e t 使用了,希尔函数。这个函数类型是经常应用于酶动力学模型饱和现象结合在酶活性站点的数量。在目前的工作,希尔函数证明足够灵活模型的函数形式 R 米 年代 D ( Δ ) 并确定相应的参数(高原价值、形状参数和half-saturation时间)。一个简单的指数饱和函数相反,希尔函数能够不同程度的sigmoidicity模型。

为了量化初始条件的影响,平衡阶段 R 米 年代 D ( Δ ) ,我们系统地增加 o ff 年代 e t 在开始分析每个轨迹。除了形状参数 所有其他参数的函数描述 R 米 年代 D ( Δ ) 表现出一种独特的依赖初始阶段的轨迹,从而表明有偏见的估计结果如果初始阶段将包括在分析中。

在引言中指出,这里的方法结合了平均的平滑作用,同时保留之间的精确时间间隔配置的特异性进行比较。例如,外推的高度 R 米 年代 D 高原 0 (见图6构型)可能被解释为“距离”两个任意选择,完全不相关的特定部分的配置相空间目前访问。如果一个选择任意配置的参考,所有其他访问的轨迹将会平均距离 0 。相比之下,图1显示了两个极限情况的参考配置:(i)的配置和最大平均距离(即。,R米年代D)to一个ll others (blue curve), which represents the “maximum outlier” ever seen in this trajectory and (ii) the configuration with minimum average distance (i.e., RMSD) to all others (red curve), which is “most central” within the trajectory.

注意,推断高原值显著大于 R 米 年代 D ( Δ ) 实际达到的轨迹(见图3)。只有当我们认为推断高原的函数 o ff 年代 e t ,我们获得现实的(即。,lower) estimates (see Figure6在实际的模拟)瞄准。在这个意义上,我们可能属性提出拟合过程一些预测能力的程度RMSD最后预计如果进行实测轨迹。这推断意味着RMSD对应于双配置分离等时间间隔足够大来考虑配置大约不相关的,独立的配置空间的代表各自的分子。

同样,“half-saturation时间” 从适合获得随增加而减小 o ff 年代 e t 。因此,在配置足够大 o ff 年代 e t 作为参考,只需要很短的时间接近 R 米 年代 D 高原 0 。在每个面板的图5,参数 方法一个几乎恒定的水平大约一半的最大值 o ff 年代 e t (即。,25% of the total simulation time). This might suggest that—regardless of the total simulation length—the fraction of usable (independent) configurations remains the same (final 75% of the trajectory). As the length of the simulation run increases, the criterion derived from the run itself gets increasingly more stringent. From the 200 ns run the first 50 ns should be discarded, which is more than 5times the total length of the 10 ns trajectory. If using the latter as a basis of estimate, however, the last 7.5 ns seem to be trustable.

尽管报告的分析目前工作是特定于TCRpMHC复合体,我们预计落后RMSD分析的方法适用于类似的分子系统,如膜蛋白类似的规模和结构。这里介绍的方法是用来评估的收敛程度MD模拟,因此这种模拟得出的结论的统计质量。

确认

作者感谢托尔。Litov andP。Petkov讨论和建设性的评论。这部分工作是支持保加利亚科学基金在格兰特DCoE-02/1/2009和部分由奥地利科学基金(FWF P22258-B12)。