基于结构的转录因子蛋白沿DNA运动的模拟与采样：从原子尺度步进到粗粒化扩散

本实验方案的目的是以植物转录因子WRKY结构域蛋白为示例体系，揭示蛋白质沿DNA进行一维扩散的结构动态。基于马尔可夫状态模型构建的原子级模拟揭示了蛋白质沿DNA以1个碱基对（1-bp）为步长的原子细节运动过程。而粗粒化模拟则侧重于采样蛋白质在DNA上持续扩散超过数十个碱基对的过程。

首先，使用一个10微秒的全原子分子动力学（MD）轨迹，均匀提取10,000帧，形成单个碱基对步进路径。在VMD中通过点击“File”（文件）并选择“Save Coordinates”（保存坐标）来准备包含10,000帧的过渡路径。然后，在“Selected Atoms”（选定原子）框中输入protein或nucleic。

在“帧”框中选择帧，然后单击“保存”以获取所需的帧。将DNA晶体结构中参考物的长轴对齐至x轴，并通过在VMD中单击“扩展”然后选择“TkConsole”，将完整的34个碱基对DNA的初始质心设置在坐标系原点。随后，在TkConsole命令窗口中输入命令。

然后，点击 VMD，依次进入 Extensions（扩展）、Analysis（分析），选择 RMSD Trajectory Tool（RMSD 轨迹工具），计算蛋白质主链的均方根距离。在原子选择框中，输入 nucleic 以及残基 14 到 23 和 46 到 55。点击 ALIGN（比对），然后点击 RMSD 框。

在 MATLAB 中，为了计算蛋白质围绕 DNA 的旋转角度 theta T（初始角位置定义为 theta 0）在 XY 平面上的旋转度，请执行以下命令。在 MATLAB 中输入指令，使用 K-means 方法将 10,000 个结构分为 25 个聚类。完成后，收集这 25 个聚类中心的结构，用于进一步的分子动力学（MD）模拟。

为了进行第一轮分子动力学（MD）模拟，使用 GROMACS 和构建系统的 sh 脚本为 25 个结构建立原子级体系。在 NPT 系综下，以 2 飞秒的时间步长，通过在 shell 中执行命令，对这 25 个体系各自进行 60 纳秒的 MD 模拟。为了对第一轮 MD 轨迹进行聚类，去除每条模拟轨迹的前 10 纳秒部分，并从 25 条各 50 纳秒的轨迹中收集构象。

对于时间无关的组分分析，在 GROMACS 中输入脚本，然后选择蛋白质与 DNA 之间的距离对作为输入参数进行投影。从 index.ndx 文件导出至一个新的文本文件 index.dat。

为了获取这些原子之间的配对信息，请使用 Python 脚本。在 MSMBuilder 命令窗口中计算每条轨迹的 415 个距离配对。接下来，执行命令进行时间独立成分分析，将数据维度降至前两个时间独立成分或向量上。

使用 MSMBuilder 处理指令后，采用 case center 方法将投影数据集聚类为 100 个簇，并选择每个簇的中心结构。为了进行第二轮分子动力学（MD）模拟，从这 100 个初始结构出发，分别进行 60 纳秒的 MD 模拟。在为所有原子施加随机初始速度后，通过在 MDP 文件中启用速度生成（Velocity generation）来添加随机初始速度。

按照之前所述，去除每次模拟的前10纳秒数据。然后，从100条50纳秒的轨迹中均匀采集2,500,000个快照，用于构建马尔可夫状态模型（MSM）。为了对第二轮分子动力学轨迹进行聚类，请在MSMBuilder中对第二轮轨迹执行时间独立成分分析。

通过执行 Python 脚本计算隐含的时间尺度以验证参数。然后通过更改参数来调整滞后时间 tau 和微观状态数量。执行命令，将构象分类为 500 个簇。

在构建马尔可夫状态模型（MSM）时，将500个微状态聚类为3至6个宏状态。为确定最合适的宏状态数目，可使用MSMBuilder中的PCCAplus算法，并通过Python脚本实现。针对每个微状态，将高维构象映射到蛋白质沿DNA的X轴位置及其旋转角度。

为计算平均首达时间，基于包含500个微状态且滞后时间为10纳秒的马尔可夫状态模型（MSM）的转移概率矩阵，在蒙特卡洛模拟中进行五条各10毫秒的轨迹模拟，并将10纳秒设为蒙特卡洛模拟的时间步长。在Python脚本中计算每对宏观状态之间的平均首达时间，并利用Bash文件计算平均首达时间的标准误差的平均值。在CafeMol 3.0软件中，通过在终端执行命令来运行粗粒化模拟。

在输入文件中指定模块后，使用各个命令设置文件名模块和作业控制模块。接着设置单位与状态模块，然后设置能量函数模块和分子动力学信息模块。所有蛋白质在DNA上的构象均映射为蛋白质沿DNA的纵向移动X和旋转角度，可进一步聚类为三个宏观状态。

S1态的稳定性较差，因为其氢键结构与模拟结构相似，而S3代表一种亚稳态，在一个碱基对步进后所有氢键发生偏移，并呈现出最稳定的状态，占比最高，达63%。中间态S2连接S1和S3，具有中等偏高的占比，为30%。S2向S3的转变过程伴随着氢键的集体断裂与重新形成，耗时约七微秒；而S1到S2的转变可在约0.06微秒内完成。通过计算蛋白质与DNA之间的接触数，共识别出四种状态。在状态1和状态3中，锌指区域朝向Y方向结合。

而在状态2和状态3中，锌指区域则朝向Y方向结合。通过测量不同DNA序列上各保守残基的步进大小，发现这些残基在polyA DNA上的步进大小比在polyAT或随机DNA序列上更为同步。构建马尔可夫状态模型的关键步骤包括：选择蛋白质与DNA数据之间用于修复1-bp步进运动的距离对，以及选择合适的微观状态和宏观状态数量。