用分子动力学模拟和结合自由能计算探索PHD指和H3K4me0相互作用：AIRE-PHD1，一项比较研究

H3K4me0肽结合改变AIRE-PHD1动力学

我们使用MD模拟研究了AIRE-PHD1对肽结合的结构和/或动态变化的可能反应。为了部分补偿单个轨迹的不完整MD采样，我们采用了多副本方法，以允许更详尽的构象采样[31]对于每个系统（即自由和束缚AIRE-PHD1），我们进行了五个分子动力学模拟（每个10 ns），其中每个模拟的最后8 ns被连接成一个轨迹并进行分析。AIRE-PHD1的自由初始结构的Cα根-平方（RMSD）偏差(图1A)和装订形式(图1B)相对稳定，RMSD<2.5º。此外，在整个模拟过程中，H3K4me0在现有的AIRE-PHD1反平行β板上形成的额外β链保持不变，并诱导β1链延伸至Glu307残基，这是通过模拟过程中的二级结构分配进行评估的(图1D、E)（我们分别使用单字母代码和三字母代码来识别肽和蛋白质残基）。

总肽结合限制了AIRE-PHD1探索的构象空间，通过Cα均方根波动（RMSF）的总体降低进行评估(图1C,图S2). 特别是分子间相互作用（总结于表S1)在前9个组蛋白H3残基和结合沟（残基Asn295-Cys310，Ser332-Trp335）之间，AIRE-PHD1在该区域的柔韧性降低。有趣的是，不与肽直接接触的结构域区域，例如连接两条β链（Gly313–Arg316）的环，也降低了它们在复合物中的波动。

接下来，我们想知道H3K4me0肽与AIRE-PHD1的结合是否会影响AIRE-PHD1的协同运动。为此，我们计算了相关矩阵Corr_ij公司所有Cα原子对。该矩阵描述了动力学过程中所有成对Cα原子围绕其平均位置移动时，它们之间的线性相关性，并提供了相关和反相关运动的全面图片。free相关矩阵的比较(图2A)并绑定AIRE-PHD1(图2B)揭示了两个系统相关运动网络的差异。一方面，自由域几乎没有表现出反相关运动(图2A)另一方面，肽的存在抑制了其中一些运动，并在AIRE-PHD1中诱导了一系列新的相关和反相关运动(图2B)减少短程相关性，增加远端残留物的协同运动。我们详细观察到，在结合后，AIRE-PHD1β链（图2A、C)减少，有利于β1和组蛋白H3肽之间的新相关性(图S3). 此外，肽的结合诱导Glu296和Pro315之间形成相关运动（图2B、D)它们通过复合体中稳定的主干相互作用联系在一起(图S4A). 相反，H3K9和Glu298之间的盐桥形成降低了Glu298-Ala300和Arg303-Glu305之间的相关性，这是由自由域中Arg303和Glu288之间的静电相互作用连接的（图2A、C,图S4B). 最后，在存在配体的情况下，我们还观察到残基Cys322-Pro325和Cys337-Cys340之间存在广泛的反相关（图2B，D)这是由Ser324和Ser338之间的两个氢键快速形成和断裂引起的(图S4C).

综合起来，RMSF和相关矩阵分析表明，H3K4me0的结合降低了AIRE-PHD1的灵活性，增加了域的协同运动，包括相关运动和反相关运动。

主成分分析

为了表征整个领域的运动，我们对轨迹产生的协方差矩阵进行了主成分分析（PCA）[32]PCA可以将相关变量的原始空间转换为独立变量的缩减空间（即主成分或特征向量）。主成分分析识别残基组的相关低能位移，并通过将轨迹投影到降维空间来强调主要蛋白质运动的幅度和方向，从而提取轨迹中捕获的慢模式[31]利用这种方法，我们确定了参与最相关的集体构象变化的蛋白质区域，并阐明了配体结合诱导的AIRE-PHD1动态调制。首先，我们观察到，与AIRE-PHD1/H3K4me0复合体相比，自由AIRE-PHD1的前几个特征向量捕获的累积方差始终较低。特别是，在自由约束AIRE-PD1中，总方差的累积比例的70%由前六个和四个特征向量捕获(表S2). 这与肽结合增加AIRE-PHD1相关运动的事实相符(图2B). 通过计算主成分的平方根内积（RMSIP），对模拟结果进行比较[33]揭示了自由和束缚AIRE-PHD1的基本运动特征不同（RMSIP=0.52），并描述了不同的基本子空间。沿着这些组分的MD轨迹投影表明，自由AIRE-PHD1的基本运动主要是由残基Arg328-Thr334组成的环向β1链的运动（图2B，E). 值得注意的是，这个内在领域“呼吸”(图2F)在H3K4me0肽的存在下被强烈还原，该肽在“开放”构象中阻断了结构域(图2G). 在整个模拟过程中，我们始终观察到Glu307和Gly333的Cα原子之间的采样距离分布变窄(图2H).

天然和突变AIRE-PHD1/H3K4me0配合物的MM/PBSA计算

接下来，我们研究了驱动AIRE-PHD1和组蛋白H3肽之间天然和突变形式相互作用的能量参数，并将计算结果与可用的实验数据进行了比较。为此，我们利用MD模拟生成的大量构象数据，使用MM/PBSA方法对野生型复合体进行结合自由能计算。该方法表示结合自由能（ΔG_comp公司)作为络合物的自由能与受体加配体的自由能之差，在多个轨迹快照上的平均值[34]接下来，我们在同一组快照上进行丙氨酸突变，并重新计算相关的结合自由能差异（ΔΔG_comp公司).结果然后与ITC测量的可用实验热力学数据进行比较[15],[20]。使用内部脚本进行计算（请参见材料和方法)这提供了程序的快速自动化。

值得注意的是，MM/PBSA方法允许快速估计结合自由能的变化，但通常不会复制绝对结合自由能值。然而，它通常表现出与实验的良好相关性，从而代表了计算和比较结合自由能变化的效率和功效之间的公平折衷。事实上，在我们的计算中，我们观察到计算出的结合自由能比实验得到的大(表1)这是一种在其他系统中已经观察到的高估。这种行为通常被归因于忽略了熵贡献，这是这些计算中的典型近似值[27],[35]–[39]在这里，我们主要关注野生型和突变复合物之间的差异，尽管该方法存在固有的局限性，但我们观察到计算和实验之间具有良好的定性一致性，相关系数为第页0.85(图3).

将结合自由能分解为其组分，包括范德华、静电、极性溶剂化和非极性溶剂化互作用能项，确定了野生型和突变体结合亲和力的主导因素。一方面，库仑（ΔG_库尔)范德瓦尔斯（ΔG_{虚拟数据仓库})和非极性溶剂化项（ΔG_{非营利组织})有利于复合物的形成。另一方面，正极性溶剂化作用（ΔG_秒)由于肽和AIRE-PHD1的极性残基和带电残基的去溶能不利。总的来说，这种高脱溶能量损失不能被库仑相互作用完全补偿，从而导致极性项对结合自由能变化（ΔG_极地的>0). 相反，非极性能量项ΔG_非极性的由范德瓦尔斯项和非极性溶剂化项之和（远小于范德瓦尔项）形成的，促进了络合物的形成（ΔG_非极性的<0). 综上所述，这些数据表明复合物由非极性相互作用稳定，并由极性贡献调节(表1). 总的来说，丙氨酸突变降低了库仑溶剂化和极性溶剂化对野生型复合物的贡献，导致正极性项的总体增加，从而对抗结合。始终ΔG_极地的与实验结果有很好的相关性(第页 = 0.83），而ΔG_非极性的没有显示出实质性的变化。

识别H3K4me0的其他PHD手指的MM/PBSA计算

受AIRE-PHD1/H3K4me0复合物令人鼓舞的结果的启发，我们接下来想知道MM/PBSA方法是否也可以用于计算和排序H3K4me0肽复合物中其他PHD指的结合ΔG，包括CHD4-PHD2、BHC80-PHDD、TRIM24-PHD和BRPF2-PHD1(表2). 五个10 ns长MD模拟的最后一纳秒被连接成一个单一轨迹（5 ns），用于MM/PBSA计算。尽管使用了不同的实验条件（主要差异在于缓冲液类型和离子强度，而pH值和温度是可比较的）和技术（荧光和ITC）来确定离解常数(表2和表S3)，我们观察到结合的实验和计算ΔΔG之间有很好的相关性(第页 = 0.96) (图4A). 重要的是，在所有五个配合物中，H3K4me0的识别通过非极性项（ΔG）稳定_非极性的<0），而ΔG_极地的由负静电和正极性溶剂化组分组成的术语对抗结合。值得注意的是，ΔG_极地的与ΔΔG相关良好_经验(第页 = 0.88），而ΔG_非极性的不同复合体之间没有显著变化(表2).

分子间接触分析

利用对五种不同配合物进行的MD模拟，我们对与肽域分子间接触相关的能量贡献（库仑和范德瓦尔斯能量）进行了比较分析。这些接触是根据它们的分子间距离（<3º）和它们在MD模拟期间的统计发生率（>30%）以及它们的相关能量来定义的(表S1)映射到域表面(图4B). 有趣的是，我们观察到组蛋白肽的前6个残基与等效残基建立了类似的接触，如基于结构的比对所定义的(图4C). 值得注意的是，主要差异涉及残基H3R8和H3K9（H3A7的贡献可以忽略不计），这表明接触次数较高(表S1)和有利的相互作用能(图4B、D;表S1)在那些PHD手指中显示出对H3K4me0的最佳亲和力。相反，缺乏这些相互作用的PHD手指，如BRPF2-PHD1，表现出最低的亲和力。总的来说，这些结果表明对H3R8和H3K9的识别可能会调节相互作用的精细细节，从而赋予不同PHD域的阅读特异性。

AIRE-HD1 MD模拟

自由形式的AIRE-PHD1结构（1xwh）和结合到H3K4me0肽（残基1-10）（2ke1）的结构用于MD模拟。使用Gromacs 4.0.7软件包进行MD模拟[43]液体模拟（OPLS）力场的优化参数[44].通过添加适当数量的钠反离子来中和系统。执行能量最小化程序和200 ps的位置约束阶段，以缓解不利的相互作用；随后在NVT系综中进行300 ps的MD模拟，以平衡系统（T=296 K）。生产运行在NPT装置中进行，P=1 bar，T=296 K，持续10 ns，时间步长为2 fs。用粒子网格Ewald（PME）方法处理长程静电，使用1.1 nm截止。对于Lennard-Jones相互作用，使用1.1纳米截止。每10个MD步骤更新一次配对列表。进行了五个独立的10 ns长MD模拟（50 ns生产运行），以便进行更好的构象取样，并对结合自由能进行统计验证。五个独立模拟的起始结构是从沉积的核磁共振束中提取的。

通过计算平均结构的RMSD值和MM/PBSAΔG的累积平均值来评估结构和能量收敛_绑定值分别作为时间的函数(图S5A、B).

为了进一步验证收敛性，我们还将每个副本的模拟长度扩展到50 ns。50 ns后，所有模拟均具有稳定的RMSD曲线(图S6). 由于10 ns和50 ns模拟所探索的构象空间是相同的（请参见下文），因此对10 ns长的模拟进行了构象分析。

H3K4me0肽复合物中其他PHD手指的MD模拟

如上所述，对CHD4-PHD2（2l75）、BHC80-PHDD（2puy）、TRIM24-PHD（3o37）和BRPF2-PHD1（2I43）与H3K4me0肽的复合物进行MD模拟。对于与H3肽复合物中的CHD4-PHD2，在进行MD模拟之前从H3K9中去除甲基。对于TRIM24-PHD，计算中仅使用了残基Asn825-Pro885，这是基于溴代多巴胺不参与PHD-H3肽相互作用的结构观察[18]对于BRPF2-PHD1，去除残基Gly13-Ser18，生成两条分离链，包括PHD指和12个残基长度的组蛋白尾（结构是通过融合构建物确定的，既包含PHD指又包含组蛋白尾，而结合亲和力是通过用相应的组蛋白肽滴定BRPF2-PHD1来确定的[19]). 对每个综合体进行五次独立的10 ns长MD模拟（50 ns生产运行）。对于NMR测定的配合物，从沉积的结构系综（2l75，2l43）中提取了五个起始结构，而对于晶体结构（2puy，3o37），五个独立的动力学运行使用了不同的种子数。根据RMSD和MM/PBSAΔG的平均结构计算累积平均值，评估结构和能量收敛_绑定值(图S5B，A).

主成分分析

协方差矩阵的主成分分析（PCA）确定沿分子动力学模拟生成的轨迹的主要低频大尺度运动。这种统计方法用于描述最相关的运动，使用一个新的基集直接反映系统所经历的集体运动。该方法可以过滤决定系统运动的主要模式的噪声，从而降低MD轨迹的维数[32]PCA需要构造协方差矩阵覆盖（cov）基于给定原子集相对于系综平均位置的三维位置起伏，通过最小二乘拟合去除整体旋转和平移运动。元素秒_ij公司矩阵的冠状病毒是(方程式1):

(1)

哪里第页_我和第页_j个是向量的Cα_我和Cα_j个轨迹中的位置；<r_我>_t吨和<r_j个>_t吨是第页_我和第页_j个时间平均值(t吨)超过MD轨迹；指数i和j表示氨基酸残基。数据通常由根据相关矩阵计算的相关图表示，其元素由

(2)

该图允许识别具有相关（在同一方向）和反相关（在相反方向）运动的残差对/组。协方差矩阵c的对角化ov型提供了一组正交的特征向量，每个特征向量由一个特征值定义，分别表示运动的方向和幅度，其中第一个特征向量表示对系统总波动的最大贡献，第二个特征向量是第二个最大贡献，依此类推。对AIRE-PHD1（残基Asn295-Thr344）结构系综的Cα坐标进行PCA计算，该结构系综由游离AIRE-PHD1的五个MD的最后8ns的串联产生，并分别与H3K4me0肽复合。

本质子空间与收敛性评估的比较

为了比较自由和束缚AIRE-PHD1模拟中识别的特征向量所描述的基本子空间，我们通过计算两组相应特征向量之间的均方根内积（RMSIP）来计算两个子空间之间的重叠：

(3)

其中v_我和w_j个是我-th和j个-两个集合的第个特征向量。RMSIP的范围从零（正交、非重叠子空间）到一（即相同子空间）。RMSIP通常是在定义基本子空间的十个特征值最大的特征向量上计算的[33].

RMSIP的计算还可以用于比较描述相同仿真的不同时间窗口的特征向量，以验证仿真收敛性[45].

因此，为了评估模拟收敛性，我们将单个模拟分为三个不同的时间窗口（2-5 ns；2-8 ns；2-10 ns），并计算了它们相应的特征向量。接下来，我们计算RMSIP以验证这些特征向量所描述的单个子空间是否相似。重要的是，所有RMSIP值均大于0.85，表明这些特征向量所描述的子空间之间存在较大重叠(表S4). 重要的是，10 ns和50 ns仿真所描述的子空间基本相同，正如两个仿真的前10个特征向量之间的高RMSIP值所示(表S5). 我们的结论是，在10毫微秒时，模拟已经达到收敛，并且2到10毫微斯的生产运行为这些系统的有意义分析提供了足够的采样。

基于MD的结合自由能计算

根据MM/PBSA方法测定结合自由能的方法已在前面描述[34]蛋白质分子与溶液中配体分子的结合自由能定义为：

(4)

进行MD模拟以生成结构的热力学加权集合（在我们的示例中，是时间等距快照的集合）。自由能项计算为所考虑结构的平均值：

(5)

能量项E_{MM（毫米）}定义为：

(6)

其中E_整数表示键能、角能和扭角能，以及E_库尔和E_LJ公司分别表示分子内静电和范德瓦尔斯能量。

溶剂化项G_解决在里面等式（7）被分成极性G_极地的和非极性贡献，G_非极性的 [34]:

(7)

在本工程中G_极地的和G_非极性的使用APBS（自适应泊松-玻耳兹曼解算器程序）进行计算[46].极性贡献G_极地的指将溶质从低介电常数（ε=1）的连续介质转移到介电常数为水的连续介质（ε=80）所需的能量。G公司_极地的使用非线性泊松-玻尔兹曼方程进行计算。栅格间距自动设置为上限0.5º。温度设置为296 K，盐浓度为0.15 M。非极性贡献G_非极性的被认为与溶剂可及表面积（SASA）成正比：

(8)

哪里γ = 0.0227千焦摩尔⁻¹Å⁻²和β = 0千焦摩尔⁻¹ [47]使用半径为1.4º的探针确定介电边界。

基于MM/PBSA方法的结合自由能计算可以根据三轨道法（TTM）或单轨道法（STM）进行。TTM需要对三种系统成分（复合物、游离配体和游离受体）进行三次单独的MD模拟。这是一种需要计算的方法，容易产生结构噪声[28],[30]相反，STM需要复合物的单一轨迹运行，从而直接从复合物结构中提取蛋白质和配体结构[28]从而使E归零_整数术语。在这种情况下，假设蛋白质和配体在结合形式和自由形式中的行为相似。这一假设对于PHD指状物是合理的，因为它们在结合时不会发生结构重排[6]–[8]如果肽结构重排发生在与PHD指结合时。在所研究的系统中，肽与PHD指结合时总是采用β链构象。因此，可以合理假设不同复合体中的熵项非常相似，并且在计算ΔΔG时会合理抵消_comp公司使用AIRE-PHD1/H3K4me0作为参考值（ΔΔG）计算差异_comp公司 = ΔG_AIRE公司-ΔG_博士). 因此，根据STM协议进行了MM/PBSA计算。

在MM/PBSA近似值内<E_{MM（毫米）}>+<G_解决>解释了与络合物形成相关的焓变化。计算确定结合自由能需要计算熵对复杂形成的贡献，包括溶质转动、平移和振动自由度的构象变化。溶质的熵贡献通常通过准手征方法（例如Schlitter方程）或正规模式分析进行估算[48]熵计算需要对自由能源景观进行全面采样，这是一个计算要求极高的步骤，可能会导致不可靠的结果[36]标准误差通常比与其他能量项相关的误差大一个数量级[49]此外，正常模式分析估计通常是非常定性的[50]在短动态时间范围内的构型熵估计可能不重要[51]基于这些考虑，以及我们主要对ΔΔG感兴趣的事实_comp公司（详见下文）我们决定在计算中忽略熵项。MM/PBSA计算考虑了五个MD模拟（即125个等时距帧）最后纳秒的平均值。标准误差（SE）计算如下：

(9)

其中σ是标准偏差，N是计算中使用的结构数量（125）。

丙氨酸突变

在野生型复合物丙氨酸突变后没有重大构象变化的情况下，可以使用MM/PBSA方法对野生型模拟中的快照进行丙氨酸突变复合物的结合自由能计算，而不是对单个突变复合物进行模拟[29],[34]该方案已成功应用于研究各种蛋白质相互作用[28]–[30]简单地说，对MD模拟获得的野生型结构进行后处理，以引入丙氨酸突变。因此，计算结合自由能值可以表示为野生型和突变复合物之间的结合自由能差异（ΔΔG_comp公司). 在这些计算中，忽略了熵，假设相似配体结合自由能的熵贡献在相对比较中抵消了[28],[30],[34],[37],[52].

GMXAPBS工具

尽管免费提供的软件Gromacs 4.0.7很受欢迎[43]和APBS[46]为了直接使用MD输出作为结合自由能计算的输入，自动组合这两个程序是不可用的。为了方便两个程序之间的接口，我们编写了一系列Bash/Perl脚本，直接对MD模拟生成的结构执行MM/PBSA计算。运行脚本只需要三个MD仿真文件：轨迹文件（TRR或XTC）、拓扑文件（TPR）和索引文件（NDX）。对于自定义的力场，该工具还需要拓扑和参数文件。计算以自动方式组织，可以在PBS队列系统中并行运行。图S7总结了脚本工作流的示意图表示。简而言之，Gromacs MD模拟生成的结构是作为PDB文件从轨迹中提取的。接下来，结构进行能量最小化（长度可由用户确定），在此期间，范德瓦尔斯项以双精度计算。然后，PDB文件被转换为PQR文件编辑confGromacs工具。APBS通过APBS辅助程序计算溶剂化（极性和非极性）和静电项库仑泊松-玻尔兹曼方程需要生成一个具有合适网格的方框。为此，蛋白质在每个维度的极端坐标都会自动从结构文件中提取出来。将20和10Ω添加到每个值中，以分别设置粗网格和细网格的限制。然后，我们的工具自动计算网格小于0.5º时适用于APBS计算的网格点数量。最后，当所有计算完成后，将所有结合能项相加，以获得沿轨道的平均结合自由能。在我们的集群（HP 300核、2.9 GHz、4 Gb RAM、InfiniBand连接、Maui/TORQUE队列系统）中，每天已经进行了3000个结构的计算。GMXAPBS工具也被改编为执行丙氨酸突变的后处理协议：GMXAPBS脚本可以将任何氨基酸突变为丙氨酸，将蛋白质的残基截断到C_β原子，并添加缺失的氢原子以完成C的四面体配位_β原子。随后，进行MM/PBSA计算，以计算丙氨酸突变对结合亲和力的能量影响。GMXAPBS是全自动的，很容易适应任何蛋白质和蛋白质结合系统。脚本被广泛注释以便于定制，输出是报告计算术语的文本或pdf文件。该工具可根据要求提供。

图的准备

所有人物都是用PyMol创建的(http://pymol.org)、VMD[53]。图形是使用Xmgree创建的(http://plasma-gate.weizmann.ac.il/Grace/).