Summary
该协议的目标是揭示蛋白质沿DNA的一维扩散的结构动力学，使用植物转录因子WRKY结构域蛋白作为示例系统。为此，已经实施了原子和粗粒度分子动力学模拟以及广泛的计算采样。
Abstract
转录因子（TF）蛋白沿DNA的一维（1-D）滑动对于促进TF的扩散以定位遗传调控的靶DNA位点至关重要。检测TF滑动或踩踏DNA的碱基对（bp）分辨率在实验上仍然具有挑战性。我们最近进行了全原子分子动力学（MD）模拟，捕获了小的WRKY结构域TF蛋白沿DNA的自发1-bp步进。基于从此类仿真中获得的10 μs WRKY步进路径，此处的协议显示了如何通过构建用于1-bp蛋白质步进的马尔可夫状态模型（MSM）对TF-DNA系统进行更广泛的构象采样，并测试各种数量的微观和宏观状态用于MSM构建。为了检查TF蛋白与结构基础的DNA的加工一维扩散搜索，该方案进一步展示了如何进行粗粒度（CG）MD模拟以采样系统的长期尺度动力学。与全原子模拟揭示的亚微秒到微秒的蛋白质步进运动相比，这种CG建模和模拟对于揭示蛋白质-DNA静电对TF蛋白质在几十微秒以上的过程扩散运动的影响特别有用。
Introduction
转录因子（TF）寻找靶DNA以结合和调节基因转录及相关活性1.除了三维（3D）扩散之外，TF的促进扩散被认为是靶DNA搜索所必需的，其中蛋白质还可以沿着一维（1D）DNA滑动或跳跃，或者在DNA2，3，4，5，6，7上进行节间转移。
在最近的一项研究中，我们对植物TF进行了数十微秒（μs）全原子平衡分子动力学（MD）模拟 - DNA8上的WRKY结构域蛋白。在微秒内捕获了WRKY在多晶A DNA上的完整1-bp步进。已经观察到蛋白质沿DNA槽的运动和氢键（HBs）断裂重整动力学。虽然这样的轨迹代表了一条采样路径，但整体蛋白质步进景观仍然缺乏。在这里，我们展示了如何使用构建的马尔可夫状态模型（MSM）围绕最初捕获的蛋白质步进路径扩展计算采样，该模型已广泛实施，用于模拟涉及实质性构象变化和时间尺度分离的各种生物分子系统9，10，11，12，13，14，15，16，17，18，19.目的是揭示TF蛋白沿DNA扩散的构象集合和亚稳态，用于一个循环步骤。
虽然上述MD模拟揭示了DNA上1 bp的蛋白质运动的原子分辨率，但TF以相同的高分辨率沿DNA的长期过程扩散的结构动力学几乎不容易获得。然而，在残留物水平上进行粗粒度（CG）MD仿真在技术上是可行的。CG模拟时间尺度可以有效地扩展到比原子模拟长20、21、22、23、24、25、26、27、28、29的几十倍或几百倍。在这里，我们展示了通过实施Takada lab30开发的CafeMol软件进行的CG模拟。
在目前的方案中，我们首先介绍了沿多A DNA和MSM构建的WRKY结构域蛋白的原子模拟，其重点是对沿DNA仅1 bp的蛋白质步进运动进行采样。然后，我们提出了同一蛋白质 - DNA系统的CG建模和模拟，该系统将计算采样扩展到蛋白质沿着DNA的数十个bps上的过程扩散。
在这里，我们使用GROMACS31，32，33 软件进行MD模拟，并使用MSMbuilder34 构建用于采样构象快照的MSM，以及使用VMD35 可视化生物分子。该协议要求用户能够安装和实现上述软件。然后，CafeMol30 软件的安装和实施对于进行CG MD模拟是必要的。在VMD中还对轨迹和可视化进行了进一步的分析。
Subscription Required. Please recommend JoVE to your librarian.
Protocol
1. 从原子MD模拟构建马尔可夫状态模型（MSM）
- 自发蛋白质步进途径和初始结构收集
- 使用先前获得的10μs全原子MD轨迹8 从“向前”的1-bp步进路径（即，每纳秒一帧）中均匀地提取10000帧。帧的总数需要足够大，以包括所有代表性构象。
- 在 VMD 中准备具有 10000 帧的过渡路径，方法是单击 “文件”>“保存坐标”，在“所选原子”框中键入蛋白质或细胞核，然后在“帧”框中选择帧，单击“ 保存 ”以获取所需的帧。
注意：先前获得的10 μs全原子MD模拟轨迹（此处称为“正向步进轨迹”），用于在34 bp均相多晶A DNA8上进行WRKY步进1-bp距离，用作启动进一步构象采样的初始路径。然而，请注意，在大多数实践中，通过执行定向或有针对性的MD模拟，或实现一般的路径生成方法等，可以构造初始路径。 - 将参考DNA的长轴（从晶体结构）对齐到x轴，并将完整的34-bp DNA的初始质心（COM）设置在坐标空间的原点，以方便进一步数据分析。为此，请单击 VMD 中 Tk 控制台>扩展 ，然后在 Tk 控制台命令窗口中键入：
源旋转.tcl
tcl 脚本可以在 补充文件 3 中找到。 - 然后通过将中心10 bp DNA（A 14至23和T 14'至23'）与晶体结构40的DNA对齐来计算蛋白质主链的均方根距离（RMSD），RMSD表示系统的几何测量值（见 图1A）。为此，请单击 “VMD >扩展>分析> RMSD 轨迹工具 ，然后在原子选择框中键入核和残基 14 到 23 和 46 到 55，单击 对齐 ，然后单击 RMSD 框以计算 RMSD 值。
- 在 MATLAB 中，通过键入命令，计算蛋白质围绕 y-z 平面上 DNA Θ（t） 的旋转度
rad2deg（atan（z/y））
初始角度定位定义为 Θ（0）=0，如之前执行的8 所示。 - 在 MATLAB41 中键入以下命令以使用 K-means 方法42，43，44 ，并通过键入将 10000 个结构分类为 25 个聚类：
[idx， C]=kmeans（ X， 25）
这里的 X 是RMSD的2D矩阵和WRKY在DNA上的旋转角度。收集这 25 个集群中心的结构，以进行进一步的 MD 模拟。
注意：由于相对于DNA采样的蛋白质RMSD覆盖约25 Å的范围，因此我们选择25个簇以每埃具有一个簇。
- 进行第一 轮MD模拟和模拟设置
- 在 parmbsc1 力场45 下使用 GROMACS 5.1.2 软件32 并使用 shell 中补充文件 2 中的 buildsystem.sh 文件，为 25 个结构构建原子系统。
- 通过在shell中键入以下命令，在NPT集成下对这25个系统进行60-ns MD仿真，时间步长为2 fs：
gmx_mpi grompp -f md.mdp -c npt.gro -p topol.top -o md.tpr
gmx_mpi mdrun -deffnm md
- 聚类 1圣 圆形 MD 轨迹
- 通过在 shell 中键入来删除每个仿真轨迹的前 10 ns：
gmx_mpi trjcat -f md.xtc -b 10000 -e 600000 -o newtraj.xtc
并从 25 × 50 ns 轨迹中收集构象以进行聚类，以便为后续更广泛的采样（第 2 轮 MD 模拟）准备输入结构。
注意：为了减少初始路径的影响并允许局部平衡，删除了10-ns的初始仿真周期。 - 选择蛋白质和DNA之间的距离对作为时间无关成分分析（tICA）46，47，48 投影的输入参数。使用GROMACS中的 make_ndx 命令来执行此操作：
gmx_mpi make_ndx -f 输入.pdb -o index.ndx
注意：这里，选择了残基Y119，K122，K125，R131，Y133，Q146，K144，R135，W116，R117，Y134，K118，Q121的蛋白质CA原子和重原子（NH1，NH2，OH，NZ，NE2，ND2），它们可以与DNA核苷酸形成氢键（HBs），它们与DNA核苷酸的O1P O2P和N6原子配对（A14-20， T19-23）。选定的氨基酸可以形成稳定的HBs或与DNA的盐桥。 - 将上述选定的原子索引从 index.ndx 文件复制到新的文本文件（索引.dat）。通过 python 脚本从 补充文件 1 中获取这些原子之间的对信息，generate_atom_indices.py并键入：
python2.6 generate_atom_indices.py索引.dat >原子索引.txt
这在蛋白质和DNA之间生成415个距离对。 - 通过在 MSMbuilder 命令窗口中键入以下命令，计算每个轨迹的 415 个距离对：
msmb AtomPairsFeaturizer -out pair_features --pair_indices AtomIndices.txt --top references.pdb --trjs “trajectories/*.xtc” --transformed pair_features --stride 5 - 执行 tICA，通过键入以下内容，将数据维度减小到前 2 个与时间无关的分量 （tIC） 或向量上：
msmb tICA -i ../tica_rc_a/tmp/ -o tica_results --n_components 2 --lag_time 10 --伽马 0.05 -t tica_results.h5
注意：tICA是一种降维方法，它计算时间滞后相关矩阵
的特征值，通过以下公式确定仿真系统最慢的松弛自由度：
其中 Xi（t） 是时间 t 处第 i 个反应坐标的值，Xj（t+Δt） 是时间 t+Δt 时第 j 个反应坐标的值。
是 Xi（t） 和 X j（t + Δ t） 整体仿真轨迹的乘积的期望值。沿最慢松弛自由度的方向对应于上述时间滞后相关矩阵的最大特征值。在这里，2个tIC似乎是在我们的MSM构造上区分三个宏观状态的最小集合（稍后讨论）。例如，还可以计算广义矩阵瑞利商（GMRQ）得分49，以探索要使用的最佳分量集。 - 使用 MSMbuilder 中的命令，通过 K-center43，44 方法将投影数据集聚类到 100 个聚类中（参见图 1B）：
msmb KCenters -i ./tica_results.h5 -o kcenters_output -t kcenters_output --n_clusters 100.
选择每个聚类的中心结构作为第 二轮 MD 模拟的初始结构。维护模拟的100个结构的模拟信息，包括位置、温度、压力等，速度除外。
注意：在第一轮25次模拟之后，初始路径的内存已经减少，因此我们在第二轮中生成了更多的聚类，例如100个聚类，以大幅扩展构象采样。
- 通过在 shell 中键入来删除每个仿真轨迹的前 10 ns：
- 进行第 二轮广泛的MD模拟
- 从这100个初始结构开始，在所有原子上施加随机初始速度后，进行60-ns MD仿真。通过在 mdp 文件中打开速度生成来添加随机初始速度，即将 md.mdp 文件gen_vel = 否 更改为 gen_vel = yes。
- 按照步骤 1.3.1 中所述，删除每个仿真的前 10 个 ns，从 100 个× 50 ns 轨迹中平均收集 2，500，000 个快照以构建 MSM。
注意：请注意，在后来的宏观状态构造中，发现了少量具有特别低总体（X-Θ平面底部约0.2%）的偏离路径状态。当宏态的总数设置为 3 到 6 时，这些偏离路径状态被分类为一个宏状态（图 2B）。由于如此低的总体宏观状态仅包括3个轨迹，这些轨迹最终被移除，因此该协议中显示的结果确实是从97×50 ns轨迹中获得的，总共有2，425，000帧或快照。
- 聚类第 2 轮 MD 轨迹
- 像以前一样，对第 2轮轨迹进行tICA。在 MSMbuilder 中键入：
msmb tICA -i ../tica_rc_a/tmp/ -o tica_results --n_components 2 --lag_time 10 --伽马 0.05 -t tica_results.h5 - 计算隐含的时间尺度，以验证相关延迟时间Δt和微观状态数的参数（见 图1C），
其中 τ 表示用于构建转移概率矩阵 （TPM） 的滞后时间;μk（τ） 表示滞后时间 τ 下 TPM 的第 k 个特征值。将 补充文件 1 中的 python 脚本用于此 python BuildMSMsAsVaryLagTime.py -d .。/ -f ../trajlist_num -i 50 -m 1000 -t 10 -n 20 -s 500. - 通过更改上面使用的参数来改变滞后时间 τ 和微状态数：
python BuildMSMsAsVaryLagTime.py -d ../ -f ../trajlist_num -i 50 -m 1000 -t 5 10 20 30 40 -n 20 -s 20 200 400 500 800 2000
注意：当隐含的时间尺度曲线开始随着时间尺度分离而趋于平稳时，系统被视为马尔可夫系统。然后，选择 Dt 作为相关延迟时间，选择 τ 作为隐含时间刻度开始趋于平稳以构建 MSM 的延迟时间。 - 因此，选择一个相对较大（但不太大）的状态数，N = 500，以及相对较短的相关延迟时间Δt = 10 ns。发现构建 MSM 的滞后时间为 τ =10 ns。
- 使用以下命令将构象分类为 500 个簇（参见 图 1D）：
msmb KCenters -i ./tica_results.h5 -o kcenters_output -t kcenters_output --n_clusters 500
- 像以前一样，对第 2轮轨迹进行tICA。在 MSMbuilder 中键入：
- 中微管理建筑
- 将 500 个微观状态集中为 3–6 个宏观状态，以找出最适合根据 MSMbuilder 中的 PCCA+ 算法50 的宏观状态的数量，方法是使用补充文件 1 python msm_lumping_usingPCCAplus.py中的 python 脚本。通过构建少量的宏观状态，即在动力学上聚集数百个微观状态，如下所述，确定生物分子最基本的构象变化的模型的简化动力学网络。
- 如步骤1.1.3和1.1.4所述，将高维构象映射到X（蛋白质沿DNA长轴的运动）和蛋白质沿DNA的旋转角度（例如，没有状态具有过低的种群<1%;见 图2C）。然后找到最能代表系统的3个宏状态（图1E）。参见 图2D ，了解蛋白质沿DNA的运动和蛋白质围绕DNA的旋转角度的快照。
注意：在之前生成10 μs自发蛋白质正向步进路径的工作中，我们还进行了5 x 4 μs平衡MD模拟，以适度扩展采样。我们显示了原始正向路径的映射（见 图2A 左图）以及先前进行的正向路径上进一步的4μs采样轨迹（见 图2A 右）8。图中显示了本工作中使用的原始100×50 ns（见 图2B 左）8 和97×50 ns轨迹的映射（见 右图2B ）。
- 计算平均首次通过时间 （MFPT）
- 根据 500 微状态 MSM 的 TPM 执行五个 10 ms 蒙特卡罗 （MC） 轨迹，并将 10 ns 的延迟时间设置为 MC 的时间步长。通过补充文件 1 python python mfpt_msm3.py中的 python 脚本计算每对宏状态之间的 MFPT52（图 3）。
- 使用 补充文件 2 中的 bash 文件计算 MFPT 的平均误差和标准误差，键入：
sh mfpt_analysis.bash
的特征值，通过以下公式确定仿真系统最慢的松弛自由度：
是 Xi（t） 和 X j（t + Δ t） 整体仿真轨迹的乘积的期望值。沿最慢松弛自由度的方向对应于上述时间滞后相关矩阵
2. 进行粗粒度（CG）模拟，对长时间动力学进行采样
- 使用CafeMol 3.0软件30进行CG模拟。查看在扩展名为 .inp 的输入配置文件中指定的 CG 模拟设置，包括输入结构、模拟参数、输出文件等。在终端上键入以下命令以运行 CG 模拟：
咖啡酚 XXX.inp - 在输入文件中指定以下块，每个块都以标签 < and ending with >>>>开头。
- 设置文件名块（必需）以指定工作目录和输入/输出文件存储路径。为这些模拟的文件名块键入以下内容：
<<<<文件名
路径 = XXXXX（工作路径）
文件名 = wrky（输出文件名）
输出 psf pdb 电影 dcd rst
path_pdb = XXXXX（输入本机结构路径）
path_ini = XXXXX（输入初始结构路径）
path_natinfo = XXXXX（本机信息文件路径）
path_para = XXXXX（参数文件路径）
>>>>
注意：由于Go-model53 在CG建模中使用，即蛋白质将偏向于天然构象，因此需要将建模结构设置为天然构象。在这里，输入晶体结构被设置为原生构象。 - 设置作业控制块（必需）以定义模拟的运行模式。键入以下命令：
<<<< job_cntl
i_run_mode = 2（= 2 恒温模拟）
i_simulate_type = 1 （=1 朗格文动力学）
i_initial_state = 2（=2 表示初始配置为本机配置）
>>>>
选择恒温朗格文动力学仿真。 - 设置单位和状态块（必需）以定义输入结构的信息。键入以下命令：
<<<< unit_and_state
i_seq_read_style = 1（=1 表示从 PDB 文件读取序列）
i_go_native_read_style = 1（=1 表示本机结构来自 PDB 文件）
1蛋白质蛋白质.pdb（单位和状态molecular_type native_structure）
2-3 DNA DNA.pdb（单位和状态molecular_type native_structure）
>>>>
注意：需要初始输入结构文件（蛋白质.pdb和DNA.pdb此处）。这些结构以 pdb 格式编写。这里需要两个pdb文件：一个是包含WRKY重原子坐标（单元1）的蛋白质结构文件，另一个是200-bp双链（ds）DNA（单元2-3）的坐标。蛋白质最初放置在距离DNA15埃的位置。 - 设置energy_function块中定义的能量功能块（必需）。键入以下命令：
<<<< energy_function
本地（1） L_GO
本地（2-3） L_DNA2
NLOCAL（1/1） GO EXV ELE
NLOCAL（2-3/2-3） ELE DNA
NLOCAL（1/2-3） EXV ELE
i_use_atom_protein = 0
i_use_atom_dna = 0
i_para_from_ninfo = 1
i_triple_angle_term = 2
>>>>
注意：在CG模拟中，蛋白质由Go-model53 粗粒度表示，每个氨基酸由放置在其Cα位置的CG颗粒表示。蛋白质构象将在Go电位下偏向天然结构或晶体结构（图4A 左）。DNA由3SPN.2模型54描述，其中每个核苷酸由3个CG颗粒S，P，N表示，分别对应于糖，磷酸盐和含氮碱基（图4A 右）。考虑不同链之间的静电和vdW相互作用。在CG模拟中，蛋白质和DNA之间的静电相互作用近似于Debye-Hückel电位55。vdW斥力能量采用与围棋模型中相同的形式。 - 设置md_information块（必需）以定义仿真信息。键入以下命令：
<<<< md_information
n_step_sim = 1
n_tstep（1） = 500000000
tstep_size = 0.1
n_step_save = 1000
n_step_neighbor = 100
i_com_zeroing = 0
i_no_trans_rot = 0
温度 = 300.0
n_seed = -1
>>>>
n_tstep是模拟步骤。将tstep_size设置为每个MD步骤的时间长度，每个CG Cafemol时间步骤约为200 fs30，因此这里的每个MD步骤原则上为200×0.1 fs。每 100 MD 步更新一次邻居列表（n_step_neighbor = 100）。将模拟温度设置为300 K.通过使用Berendsen恒温器56更新蛋白质结构的速度型Verlet算法来控制温度。
注意：n_step_sim是基于Go模型势的盆地数，或能量曲线的局部最小数。多盆电位允许蛋白质构象偏向于不同的构象，以便蛋白质构象可以从一个局部最小值变化到另一个局部最小值。这里只使用单个盆地Go模型，这意味着在模拟中蛋白质只有一个偏置构象（晶体结构）。同时，由于在CG背景下没有蛋白质 - DNA氢键相互作用等，因此分子运动的采样速度甚至更快，即比原子模拟中的>10倍。 - 设置静电块（仅在使用静电相互作用时才需要）作为不同链之间的静电相互作用的考虑因素，因此使用此块通过键入来定义静电相互作用的参数：
<<<<静电
cutoff_ele = 10.0
ionic_strength = 0.15
>>>>
将静电相互作用中的 Debye 长度设置为 10 Å，对应于溶液条件。将离子强度设置为0.15 M，如生理状况一样。
- 设置文件名块（必需）以指定工作目录和输入/输出文件存储路径。为这些模拟的文件名块键入以下内容：
Subscription Required. Please recommend JoVE to your librarian.
Representative Results
旋转耦合滑动或 1 bp 步进的 WRKY 从 MSM 结构
DNA上的所有蛋白质构象都映射到蛋白质COM沿DNA的纵向运动X和旋转角度（见 图3A）。这两个度的线性偶联表示WRKY结构域蛋白在DNA上的旋转耦合步进。构象可以在 MSM 中进一步聚类为 3 个宏状态（S1、S2 和 S3）。然后，WRKY 的正向步进遵循宏观状态转换 S1->S2->S3。S1是指由建模结构（基于WRKY-DNA复合物40的晶体结构）引发的亚稳态，其总体约为6%。请注意，在当前的建模中，初始蛋白质构象是从蛋白质与特定的W-box DNA序列40结合的晶体结构中采用的。因此，这种建模的蛋白质 - 聚A-DNA复合物导致比阶梯或最终松弛结构（S3）更不利的初始结构（S1）。然而，人们可以发现蛋白质 - DNA界面处的氢键（HBs）在S3中心附近恢复，就像在S1中心附近一样（见 图3B）。处于 S1 状态的 HB 维护良好：K125 与 A15、R131、Q146 和 Y133 与 A16、K144 和 Y119 与 A17、R135 与 A18（图 3B 左上角）。S3是指1-bp蛋白步进后的亚稳态，几乎所有的HBs都移位了1-bp距离（图3B 底部），并且结构看起来稳定，群体最高（63%）。中间状态 S2 连接 S1 和 S3，具有中高人口（~30%）。我们发现R135和K144在这种中间状态下非常灵活，通常可以用当前的核苷酸破坏HBs，并用下一个核苷酸进行重整（右上图3B ）。总体而言，WRKY蛋白COM移动了约2.9 Å并旋转约55°以达到1 bp。WRKY步进的限速步骤是S2->S3，它基本上允许对HB进行集体分解和重整，平均需要约7 μs。相比之下，S1至S2可以在~0.06 μs或60-ns的时间内非常快速地转载（图3B），主要涉及蛋白质COM波动（例如，由于DNA上的蛋白质取向变化）。
CG模型中加工扩散过程中WRKY的单链偏置
在我们最近的研究中，我们发现WRKY结构域蛋白优先与dsDNA的一条链结合，无论在1-bp步进或静态结合期间;单链偏倚变得非常突出，特别是在特定的DNA序列结合8上。同时，目前尚不清楚这种趋势在蛋白质沿着DNA的工艺扩散过程中是否仍然存在。在这里，我们试图通过CG模拟来检查潜在的 链偏置 。有趣的是，在加工扩散过程中，在WRKY的CG模拟中已经确定了重要的单链DNA结合构型。为了看到这一点，在各自的DNA链上计算蛋白质和DNA之间的接触号（见 图4B）。当蛋白质CG颗粒和DNA CG P（磷酸基团）颗粒之间的距离小于7 Å时，考虑接触。蛋白质确实对其中一条DNA链表现出偏倚（例如，~4个接触一条链，~1个接触到另一条链），也就是说，即使没有模拟蛋白质 - DNA界面处的详细相互作用，例如HBs。
然而，优选的DNA链可以在DNA的两条链之间不时切换，这取决于DNA上蛋白质的结合方向或配置。特别是，根据蛋白质和DNA各链之间形成的接触号，这里主要有4种状态（如图 4B，C中标有1，2，3和4）。在状态1和3中，锌指区域与-Y方向结合，首选链是蓝色链。在状态2和3中，锌指区域向+Y方向结合，首选链成为红色链。还发现锌 - 无花果区与DNA主要相互作用（见 图4D）。因此，与锌指区域紧密结合的DNA链确实是首选的DNA链。根据上述采样，似乎链偏倚仍然存在，但在加工蛋白扩散的CG模型中在两条DNA链之间切换。
CG模拟中的蛋白质个体残留步进
之前从我们的CG模拟中注意到，WRKY的步进尺寸可能因不同的DNA序列8而异。蛋白质COM倾向于在均质的多A DNA上执行步骤1 bp。在具有2 bp周期性的多效AT DNA上，2-bp步进的比例似乎增加。
此外，在这里，我们检查了单个蛋白质残基是否在蛋白质 - DNA界面处同步移动。我们计算了WRKY图案（WRKYGQK）中每1000次步进中每个高度保守的残基的步进尺寸（图5A）。因此，每个保守残留物的残余步进尺寸可以从CG模拟中测量。结果确实表明，这些单个残基的步进尺寸在poly-A DNA上比在poly-AT或随机DNA序列上更同步（图5B）。
图1：构象生成和微观状态/宏观状态构建。 （A）在蛋白质-DNA RMSD上绘制的初始正向步进路径和蛋白质围绕DNA的旋转角度。最初选择的25个结构用红色圆圈标记。（B）将来自第 1轮25 x 50 ns MD模拟轨迹的100个构象聚类中心映射到两个最高特征值tIC方向上。（C） 隐含时间尺度的绘图，作为使用所选距离对作为输入的 TICA 构建 MSM 构造的滞后时间函数。对于每个集合，通过将构象投影到前2个tIC上，然后进行K中心聚类以产生20至2000个微状态（从左到右列），tICA的相关延迟时间为5至40 ns（从上到下行）来构建MSM。（D）构建的500个微观状态和（E）进一步构建的3个宏观状态，相应的微观状态中心沿最高的两个tIC方向映射。 请点击此处查看此图的大图。
图2：宏观状态的构建（A）初始前向步进路径轨迹（左）的映射和少量额外的微秒轨迹采样（右）在蛋白质质心（COM）上沿DNA长轴（X）和围绕DNA的旋转角度的运动（先前获得8）。（B）当前MSM施工中使用的原始100×50 ns轨迹和97×50 ns轨迹的映射。（C）在广泛的采样图上标注了3-6个宏观状态的构造及其来自构建的MSM的种群。（D）分别显示了蛋白质运动X和围绕DNA的旋转角度。采样的构象最终被归为3个宏观状态，其中红色，蓝色和灰色分别对应于宏观状态1，2和3。请点击此处查看此图的大图。
图3：WRKY结构域蛋白的MSM踩在多A DNA上。 （A）MD构象快照的投影到蛋白质COM运动X的坐标上，并相对于DNA的旋转角度。3 种宏态 S1、S2 和 S3 分别以红色、蓝色和灰色着色。（B）构造的3个宏观态的代表性构象和过渡均值-首次通过时间（MFPT）。图中显示了蛋白质和DNA之间的关键氢键。 请点击此处查看此图的大图。
图4：CG模型中的粗粒（CG）模型和蛋白质和DNA链之间形成的接触。 （A）蛋白质（左）和DNA的粗粒化（右）。（B）WRKY与模拟中每条DNA链的接触号。（C）4种接触模式的分子观点。锌指附近的蛋白质区域以灰色着色，另一个区域以绿色显示。（D）每种蛋白质氨基酸与DNA的接触概率。当氨基酸的CG颗粒与任何DNA CG颗粒之间的距离小于7 Å时，氨基酸被认为与DNA接触。 请点击此处查看此图的大图。
图5：单个蛋白质氨基酸在WRKY基序中的扩散步长，如WRKY沿DNA移动。 （A）原子结构（左）和粗粒化后的高度保守残基（WRKYGQK）（右）。（B）不同DNA序列（poly-A;poly-AT;随机序列）上每个保守残基的步进尺寸 请点击这里查看该图的放大版本。
补充文件 1： 该协议中使用的python代码和软件。MSM主要是通过使用MSMbuilder构建的，附加了必要的python代码。 请点击此处下载此文件。
补充文件 2： 原子分子动力学模拟由GROMACS进行，还附有构建全原子模拟的命令和必要文件。粗粒度模拟由CafeMol软件进行。通过VMD和MATLAB对仿真结果进行分析。 请点击此处下载此文件。
补充文件 3： 在VMD中旋转和移动蛋白质的tcl脚本。 请点击此处下载此文件。
Subscription Required. Please recommend JoVE to your librarian.
Discussion
这项工作解决了如何进行基于结构的计算模拟和采样，以揭示转录因子或TF蛋白沿着DNA移动，不仅在步进的原子细节下，而且在过程扩散中，这对于TF在DNA靶标搜索中的促进扩散至关重要。为此，首先构建了小TF结构域蛋白WRKY沿着均匀的poly-A DNA步进1-bp的马尔可夫状态模型或MSM，以便可以揭示DNA上的蛋白质构象集合以及蛋白质 - DNA界面处的集体氢键或HB动力学。为了获得MSM，我们沿着自发蛋白质步进路径（从先前的10μs模拟中获得）进行了两轮广泛的全原子MD模拟，电流采样聚合为7.5μs（125 x 60 ns）。如此广泛的采样为我们提供了构象聚类成数百个微观状态的快照，利用蛋白质 - DNA界面对距离作为聚类的几何测量。通过检测时间尺度与针对各个MD模拟的各种长度或滞后时间计算的隐含时间尺度的分离，部分验证了MSM构造的马尔可夫性质。然后测试并比较20-2000个微态的时间尺度分离特性，并为MSM构建选择500个微态。此外，将500个微观状态动态地归为少量的宏观状态，为此我们测试了各种数量的状态，发现三个宏观状态足以满足当前系统的需求。三态模型简单地表明，状态S1向S2的过渡相对较快（在几十ns内），以DNA上的蛋白质质心（COM）波动为主，而状态S2向S3的转移速度较慢并且是限速的（平均约7 μs），主要由集体HB动力学作为步进。请注意，将微观状态的动力学归为少数动力学上不同的宏观状态仍然受到方法学发展的影响，测试了不同的算法和改进的机器学习技术57，58，59，60，61，62，63.构建 MSM 的关键步骤包括选择 tICA 中使用的距离对，并确定用于构建微状态的参数。距离对的选择是基于知识的，选择最基本的交互对很重要。构造微状态的参数，如相关延迟时间、滞后时间、微状态的多倍数，需要正确设置，以确保系统为马尔可夫。
通过这样的努力，可以系统地揭示蛋白质沿DNA步进1-bp的亚微到微秒的蛋白质结构动力学以及原子细节。原则上，通过从MSM构造获得的转移概率矩阵，系统可以进化到超过微秒的长时间尺度，或者说，接近毫秒及以上的13，17，64。然而，MSM采样和构造存在内在局限性，它依赖于围绕某个初始路径的亚微秒级个体模拟，并且马尔可夫性质可能无法很好地保证65，66。在大多数实践中，初始路径是在强迫或加速度下构建的，尽管在当前系统中，我们利用从10ms平衡模拟8获得的自发蛋白质步进路径（无强迫或加速度）。由于原子模拟的高计算成本，总体构象采样仍然受到数十微秒的限制。蛋白质步进的这种微秒采样不太可能提供足够的构象以出现在长时间尺度的工艺TF扩散上。如果实现当前获得的转移概率矩阵超过一定的时间尺度，并且不能保证马尔可夫性质确保正确使用当前的MSM14，52，66，则内存问题将变得严重。因此，为了对TF沿DNA的长期尺度过程扩散进行采样，改为进行残基水平粗粒度或CG建模和模拟，以平衡结构基础和降低计算成本。
在CG建模和模拟中，蛋白质残基和DNA核苷酸由磁珠表示（即，一个氨基酸的一个珠子和一个核苷酸的三个珠子），通过Go模型维持蛋白质构象，朝向天然或预平衡构型30，53。虽然HB相互作用的原子水平在CG模型中变得不存在，但蛋白质 - DNA静电相互作用得到了很好的维持，这似乎能够捕获蛋白质沿DNA 67，68，69，70的过程扩散中的主要动力学特征。这里给出了详细的实现协议，用于对WRKY-DNA系统进行建模和仿真。具有代表性的结果表明，首先，在先前的WRKY-DNA系统原子模拟中呈现的单链DNA偏倚在CG模型中持续存在，而在过程扩散过程中采样的各种蛋白质取向/构型导致两条链之间的偏倚不时切换。因此，这种DNA链偏倚不一定与HB关联有关，但似乎主要依赖于蛋白质 - DNA静电相互作用，这些相互作用因DNA上的各种蛋白质配置或取向而异。接下来，蛋白质-DNA界面处或附近的单个氨基酸，例如高度保守的WRKQGQK基序，显示出不同DNA序列的不同步进大小或同步模式。在我们之前的研究中，仅显示了蛋白质COM的步进尺寸变化，因为蛋白质被建模为沿着不同的DNA序列扩散。请注意，DNA的当前CG模型支持具有不同参数化的DNA序列变异54，71，72，尽管缺少原子细节。因此，在基于结构的蛋白质-DNA系统建模中，适当的DNA序列依赖性参数化对于揭示多个时间和长度尺度上的蛋白质-DNA搜索和识别机制至关重要。
Subscription Required. Please recommend JoVE to your librarian.
Disclosures
作者没有利益冲突。
Acknowledgments
这项工作得到了NSFC Grant #11775016和#11635002的支持。JY得到了UCI的CMCF通过NSF DMS 1763272和Uci的西蒙斯基金会赠款#594598和启动基金的支持。LTD得到了上海市自然科学基金#20ZR1425400和#21JC1403100。我们也感谢北京计算科学研究中心（CSRC）的计算支持。
Materials
| Name | Company | Catalog Number | Comments |
| CafeMol | Kyoto University | coarse-grained (CG) simulations | |
| GROMACS | University of Groningen Royal Institute of Technology Uppsala University | molecular dynamics simulations software | |
| Matlab | MathWorks | Numerical calculation software | |
| MSMbuilder | Stanford University | build MSM | |
| VMD | UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN | molecular visualization program |
References
- Latchman, D. S.
- Berg, O. G., von Hippel, P. H. Selection of DNA binding sites by regulatory proteins. Statistical-mechanical theory and application to operators and promoters. Journal of Molecular Biology. 193 (4), 723-750 (1987).
- von Hippel, P. H., Berg, O. G. Facilitated target location in biological systems. The Journal of Biological Chemistry. 264 (2), 675-678 (1989).
- Halford, S. E., Marko, J. F. How do site-specific DNA-binding proteins find their targets. Nucleic Acids Research. 32 (10), 3040-3052 (2004).
- Slusky, M., Mirny, L. A. Kinetics of protein-DNA interaction: facilitated target location in sequence-dependent potential. Biophysical Journal. 87 (6), 4021-4035 (2004).
- Bauer, M., Metzler, R. Generalized facilitated diffusion model for DNA-binding proteins with search and recognition states. Biophysical Journal. 102 (10), 2321-2330 (2012).
- Shvets, A. A., Kochugaeva, M. P., Kolomeisky, A. B. Mechanisms of Protein Search for Targets on DNA: Theoretical Insights. Molecules. 23 (9), Basel, Switzerland. 2106 (2018).
- Dai, L., Xu, Y., Du, Z., Su, X. D., Yu, J. Revealing atomic-scale molecular diffusion of a plant-transcription factor WRKY domain protein along DNA. Proceedings of the National Academy of Sciences of the United States of America. 118 (23), 2102621118 (2021).
- Chodera, J. D., Singhal, N., Pande, V. S., Dill, K. A., Swope, W. C. Automatic discovery of metastable states for the construction of Markov models of macromolecular conformational dynamics. The Journal of Chemical Physics. 126 (15), 155101 (2007).
- Pan, A. C., Roux, B. Building Markov state models along pathways to determine free energies and rates of transitions. The Journal of Chemical Physics. 129 (6), 064107 (2008).
- Bowman, G. R., Huang, X., Pande, V. S. Using generalized ensemble simulations and Markov state models to identify conformational states. Methods. 49 (2), San Diego, California. 197-201 (2009).
- Prinz, J. H., et al. Markov models of molecular kinetics: Generation and validation. The Journal of chemical physics. 134 (17), 174105 (2011).
- Chodera, J. D., Noé, F. Markov state models of biomolecular conformational dynamics. Current Opinion in Structural Biology. 25, 135-144 (2014).
- Malmstrom, R. D., Lee, C. T., Van Wart, A. T., Amaro, R. E. On the Application of Molecular-Dynamics Based Markov State Models to Functional Proteins. Journal of Chemical Theory and Computation. 10 (7), 2648-2657 (2014).
- Husic, B. E., Pande, V. S. Markov State Models: From an Art to a Science. Journal of the American Chemical Society. 140 (7), 2386-2396 (2018).
- Sittel, F., Stock, G. Perspective: Identification of collective variables and metastable states of protein dynamics. The Journal of chemical physics. 149 (15), 150901 (2018).
- Wang, W., Cao, S., Zhu, L., Huang, X. Constructing Markov State Models to elucidate the functional conformational changes of complex biomolecules. WIREs Computational Molecular Science. 8, 1343 (2018).
- Peng, S., et al. Target search and recognition mechanisms of glycosylase AlkD revealed by scanning FRET-FCS and Markov state models. Proceedings of the National Academy of Sciences of the United States of America. 117 (36), 21889-21895 (2020).
- Tian, J., Wang, L., Da, L. T. Atomic resolution of short-range sliding dynamics of thymine DNA glycosylase along DNA minor-groove for lesion recognition. Nucleic Acids Research. 49 (3), 1278-1293 (2021).
- Chu, J. -W., Izveko, S., Voth, G. The multiscale challenge for biomolecular systems: coarse-grained modeling. Molecular Simulation. 32 (3-4), 211-218 (2006).
- Marrink, S. J., Risselada, H. J., Yefimov, S., Tieleman, D. P., De Vries, A. H. The MARTINI force field: coarse grained model for biomolecular simulations. The Journal of Physical Chemistry B. 111 (27), 7812-7824 (2007).
- Givaty, O., Levy, Y. Protein sliding along DNA: dynamics and structural characterization. Journal of Molecular Biology. 385 (4), 1087-1097 (2009).
- Khazanov, N., Levy, Y. Sliding of p53 along DNA can be modulated by its oligomeric state and by cross-talks between its constituent domains. Journal of Molecular Biology. 408 (2), 335-355 (2011).
- Riniker, S., Allison, J. R., van Gunsteren, W. F. On developing coarse-grained models for biomolecular simulation: a review. Physical Chemistry Chemical Physics : PCCP. 14 (36), 12423-12430 (2012).
- Kmiecik, S., et al. Coarse-Grained Protein Models and Their Applications. Chemical Reviews. 116 (14), 7898-7936 (2006).
- Bhattacherjee, A., Krepel, D., Levy, Y. Coarse-grained models for studying protein diffusion along DNA. WIREs Computational Molecular Science. 6, 515-531 (2016).
- Wang, J., et al. Machine Learning of Coarse-Grained Molecular Dynamics Force Fields. ACS Central Science. 5 (5), 755-767 (2019).
- Joshi, S. Y., Deshmukh, S. A. A review of advancements in coarse-grained molecular dynamics simulations. Molecular Simulation. 47 (10-11), 786-803 (2021).
- Bigman, L. S., Greenblatt, H. M., Levy, Y. What Are the Molecular Requirements for Protein Sliding along DNA. The Journal of Physical Chemistry B. 125 (12), 3119-3131 (2021).
- Kenzaki, H., et al. CafeMol: A Coarse-Grained Biomolecular Simulator for Simulating Proteins at Work. Journal of Chemical Theory and Computation. 7 (6), 1979-1989 (2011).
- Berendsen, H. J. C., vander Spoel, D., van Drunen, R. GROMACS: a message-passing parallel molecular dynamics implementation. Computer Physics Communications. 91 (1-3), 43-56 (1995).
- vander Spoel, D., et al. GROMACS: fast, flexible, and free. Journal of Computational Chemistry. 26 (16), 1701-1718 (2005).
- Abraham, M. J., et al. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX. 1-2, 19-25 (2015).
- Harrigan, M. P., et al. MSMBuilder: Statistical Models for Biomolecular Dynamics. Biophysical journal. 112 (1), 10-15 (2017).
- Humphrey, W., Dalke, A., Schulten, K.
- Izrailev, S., et al. Steered Molecular Dynamics. Computational Molecular Dynamics: Challenges, Methods, Ideas. 4, Springer. Berlin, Heidelberg. 39-65 (1999).
- Schlitter, J., Engels, M., Krüger, P. Targeted molecular dynamics: a new approach for searching pathways of conformational transitions. Journal of Molecular Graphics. 12 (2), 84-89 (1994).
- Maragliano, L., Fischer, A., Vanden-Eijnden, E., Ciccotti, G. String method in collective variables: minimum free energy paths and isocommittor surfaces. The Journal of Chemical Physics. 125 (2), 24106 (2006).
- Weiss, D. R., Levitt, M. Can morphing methods predict intermediate structures. Journal of Molecular Biology. 385 (2), 665-674 (2009).
- Xu, Y. P., Xu, H., Wang, B., Su, X. D. Crystal structures of N-terminal WRKY transcription factors and DNA complexes. Protein. 11 (3), 208-213 (2020).
- Higham, D. J., Higham, N. J.
- Hartigan, J. A., Wong, M. A. Algorithm AS 136: A K-Means Clustering Algorithm. Journal of the Royal Statistical Society. Series C (Applied Statistics). 28 (1), 100-108 (1979).
- Gonzalez, T. F. Clustering to minimize the maximum intercluster distance. Theoretical Computer Science. 38, 293-306 (1985).
- Zhao, Y., Sheong, F. K., Sun, J., Sander, P., Huang, X. A fast parallel clustering algorithm for molecular simulation trajectories. Journal of Computational Chemistry. 34 (2), 95-104 (2013).
- Ivani, I., et al. Parmbsc1: a refined force field for DNA simulations. Nature Methods. 13 (1), 55-58 (2016).
- Naritomi, Y., Fuchigami, S. Slow dynamics of a protein backbone in molecular dynamics simulation revealed by time-structure based independent component analysis. The Journal of Chemical Physics. 139 (21), 215102 (2013).
- Naritomi, Y., Fuchigami, S. Slow dynamics in protein fluctuations revealed by time-structure based independent component analysis: the case of domain motions. The Journal of Chemical Physics. 134 (6), 065101 (2011).
- Pérez-Hernández, G., Paul, F., Giorgino, T., De Fabritiis, G., Noé, F. Identification of slow molecular order parameters for Markov model construction. The Journal of Chemical Physics. 139 (1), 015102 (2013).
- McGibbon, R. T., Pande, V. S. Variational cross-validation of slow dynamical modes in molecular kinetics. The Journal of Chemical Physics. 142 (12), 124105 (2015).
- Deuflhard, P., Weber, M. Robust Perron cluster analysis in conformation dynamics. Linear Algebra and its Applications. 398, 161-184 (2005).
- Silva, D. A., et al. Millisecond dynamics of RNA polymerase II translocation at atomic resolution. Proceedings of the National Academy of Sciences of the United States of America. 111 (21), 7665-7670 (2014).
- Swope, W. C., Pitera, J. W., Suits, F. Describing Protein Folding Kinetics by Molecular Dynamics Simulations. 1. Theory. The Journal of Physical Chemistry B. 108 (21), 6571-6581 (2004).
- Clementi, C., Nymeyer, H., Onuchic, J. N. Topological and energetic factors: what determines the structural details of the transition state ensemble and "en-route" intermediates for protein folding? An investigation for small globular proteins. Journal of molecular biology. 298 (5), 937-953 (2000).
- Hinckley, D. M., Freeman, G. S., Whitmer, J. K., De Pablo, J. J. An experimentally-informed coarse-grained 3-Site-Per-Nucleotide model of DNA: structure, thermodynamics, and dynamics of hybridization. The Journal of chemical physics. 139 (14), 144903 (2013).
- Debye, P., Huckel, E. The theory of the electrolyte II-The border law for electrical conductivity. Physikalische Zeitschrift. 24, 305-325 (1923).
- Berendsen, H. J., Postma, J. V., van Gunsteren, W. F., DiNola, A., Haak, J. R. Molecular dynamics with coupling to an external bath. The Journal of Chemical Physics. 81, 3684-3690 (1984).
- Bowman, G. R. Improved coarse-graining of Markov state models via explicit consideration of statistical uncertainty. The Journal of Chemical Physics. 137 (13), 134111 (2012).
- Jain, A., Stock, G. Identifying metastable states of folding proteins. Journal of Chemical Theory and Computation. 8 (10), 3810-3819 (2012).
- Röblitz, S., Weber, M. Fuzzy spectral clustering by PCCA+: application to Markov state models and data classification. Advances in Data Analysis and Classification. 7, 147-179 (2013).
- Mardt, A., Pasquali, L., Wu, H., Noé, F. VAMPnets for deep learning of molecular kinetics. Nature Communications. 9 (1), 5 (2018).
- Wang, W., Liang, T., Sheong, F. K., Fan, X., Huang, X. An efficient Bayesian kinetic lumping algorithm to identify metastable conformational states via Gibbs sampling. The Journal of Chemical Physics. 149 (7), 072337 (2018).
- Chen, W., Sidky, H., Ferguson, A. L. Nonlinear discovery of slow molecular modes using state-free reversible VAMPnets. The Journal of Chemical Physics. 150 (21), 214114 (2019).
- Gu, H., et al. RPnet: a reverse-projection-based neural network for coarse-graining metastable conformational states for protein dynamics. Physical Chemistry Chemical Physics :PCCP. 24 (3), 1462-1474 (2022).
- Lane, T. J., Bowman, G. R., Beauchamp, K., Voelz, V. A., Pande, V. S. Markov state model reveals folding and functional dynamics in ultra-long MD trajectories. Journal of the American Chemical Society. 133 (45), 18413-18419 (2011).
- Konovalov, K. A., Unarta, I. C., Cao, S., Goonetilleke, E. C., Huang, X. Markov State Models to Study the Functional Dynamics of Proteins in the Wake of Machine Learning. JACS Au. 1 (9), 1330-1341 (2021).
- Cao, S., Montoya-Castillo, A., Wang, W., Markland, T. E., Huang, X. On the advantages of exploiting memory in Markov state models for biomolecular dynamics. The Journal of Chemical Physics. 153 (1), 014105 (2020).
- Brandani, G. B., Takada, S. Chromatin remodelers couple inchworm motion with twist-defect formation to slide nucleosomal DNA. PLoS Computational Biology. 14 (11), 1006512 (2018).
- Tan, C., Terakawa, T., Takada, S. Dynamic Coupling among Protein Binding, Sliding, and DNA Bending Revealed by Molecular Dynamics. Journal of the American Chemical Society. 138 (27), 8512-8522 (2016).
- Terakawa, T., Takada, S. p53 dynamics upon response element recognition explored by molecular simulations. Scientific reports. 5, 17107 (2015).
- Brandani, G. B., Niina, T., Tan, C., Takada, S. DNA sliding in nucleosomes via twist defect propagation revealed by molecular simulations. Nucleic Acids Research. 46 (6), 2788-2801 (2018).
- Knotts, T. A., Rathore, N., Schwartz, D. C., de Pablo, J. J. A coarse grain model for DNA. The Journal of Chemical Physics. 126 (8), 084901 (2007).
- Freeman, G. S., Hinckley, D. M., Lequieu, J. P., Whitmer, J. K., de Pablo, J. J.
