1. 项目概述当计算生物学遇上药物发现新前沿最近几年药物研发领域一个非常有意思的趋势就是计算模拟的力量正以前所未有的深度介入到早期发现环节。我们不再仅仅依赖高通量筛选的“大海捞针”而是尝试用计算机先“算”出有潜力的方向。这次直播要聊的“分子模拟驱动的靶向表观转录组学蛋白的药物发现”就是一个典型的、处于交叉学科前沿的硬核课题。简单来说它要解决的核心问题是如何利用计算机模拟技术去设计和发现能够精准调控RNA表观遗传修饰的关键蛋白或靶向这些蛋白的小分子从而为开发新型疗法尤其是针对癌症、神经退行性疾病等提供全新的武器。这里涉及几个关键概念的交汇。首先是“表观转录组学”你可以把它理解为DNA表观遗传学在RNA层面的延伸。我们都知道DNA的甲基化、乙酰化等修饰能影响基因是否表达而不改变DNA序列本身。同样RNA分子如信使RNA、转运RNA等上也存在着超过170种化学修饰比如常见的m6AN6-甲基腺苷。这些修饰就像贴在RNA上的“小标签”精细地调控着RNA的稳定性、定位、翻译效率等进而深刻影响细胞命运和疾病进程。靶向这些修饰过程的“书写器”writer、“擦除器”eraser和“阅读器”reader蛋白就成了非常热门的药物靶点。然而这些蛋白与RNA、或者与小分子抑制剂之间的相互作用在原子尺度上异常复杂且动态。这就是“分子模拟”大显身手的地方。它不再满足于静态的晶体结构“快照”而是通过分子动力学模拟等方法让我们能像看电影一样观察蛋白质与配体如何相互靠近、结合口袋如何发生构象变化、关键氨基酸如何参与相互作用。这种动态视角对于理解别构调控、设计高选择性抑制剂至关重要。我自己的体会是单纯看一个结合模式图远不如跟踪一段几百纳秒的模拟轨迹来得透彻很多设计灵感就藏在那些瞬间的氢键形成或疏水口袋的坍塌过程中。这场直播适合谁呢如果你是计算化学、生物信息学领域的研究者或学生想了解如何将模拟技术应用于最前沿的生物学问题如果你是药物化学或药物发现领域的从业者希望拓展对RNA靶向药物这一新兴领域的认知并学习如何借助计算工具进行理性设计甚至如果你是相关领域的生物学家想更深入地理解你研究的蛋白靶点背后的物理化学原理那么这次分享都应该能带来不少启发。接下来我会结合这个主题拆解其背后的技术逻辑、实操难点以及我们在这个领域趟过的一些“坑”。2. 核心思路与技术选型为什么是分子模拟表观转录组学2.1 靶点特性决定了方法的选择表观转录组学相关的蛋白尤其是“阅读器”蛋白如YTH结构域家族它们识别RNA修饰的过程具有几个鲜明特点这些特点直接影响了我们为何以及如何运用分子模拟。首先是识别的高度特异性和微弱的结合自由能变化。以m6A阅读器蛋白YTHDF2为例它需要精准区分一个普通的腺苷A和一个甲基化的腺苷m6A。这种区别仅仅是一个甲基基团的有无在结合能上差异可能只有几个kcal/mol。传统的分子对接打分函数对于如此细微的能量差异区分度常常不足容易产生假阳性或假阴性。而分子动力学模拟结合自由能微扰FEP或热力学积分TI等计算方法能够更精确地计算相对结合自由能为理解特异性来源和优化分子结构提供定量依据。其次是结合界面的高度动态性和溶剂效应。RNA-蛋白质复合物界面往往富含带电荷的残基如精氨酸、赖氨酸和极性相互作用水分子的参与和置换在结合过程中扮演关键角色。静态结构无法捕捉水分子的桥接作用或结合过程中“去溶剂化”的能垒。显式溶剂下的分子动力学模拟可以直观展示水分子网络的形成与破坏这是许多基于经验打分函数的快速筛选方法所忽略的。再者存在大量的别构调控和蛋白构象变化。有些“书写器”或“擦除器”蛋白如METTL3/METTL14复合物在催化过程中会发生显著的结构重排。要设计别构抑制剂就必须理解这些构象变化的路径和能量景观。常规的分子对接通常将受体视为刚性体这显然不够。我们需要采用增强采样模拟比如元动力学Metadynamics或高斯加速分子动力学GaMD来主动探索这些罕见的构象变化事件寻找潜在的别构口袋。基于以上特性我们的技术栈选择就清晰了对于初始筛选和结合模式预测采用分子对接如AutoDock Vina, Glide进行快速扫描但必须对结果保持审慎需结合多种对接程序的结果和视觉检查。对于结合稳定性和动态相互作用分析必须进行经典分子动力学模拟使用AMBER, GROMACS, NAMD等软件时长通常需要数百纳秒甚至微秒以观察复合物是否稳定、关键相互作用是否持续。对于结合自由能的精确计算和先导化合物优化引入自由能微扰或热力学积分计算。虽然计算成本高昂但在优化分子结构、预测活性趋势SAR时其准确性远高于经验方法可以大幅减少不必要的化学合成。对于探索构象变化和别构机制应用增强采样方法。例如对怀疑存在别构效应的蛋白可以定义反映其构象变化的集体变量如特定原子间的距离、二面角进行元动力学模拟绘制自由能面从而发现隐藏的构象状态和潜在的抑制剂结合位点。这个组合策略不是简单的堆砌而是一个从“快速初筛”到“动态验证”再到“精准优化”的递进过程。在实际项目中计算资源需要合理分配通常80%的资源会投入到对少数几个最优候选分子的深入模拟和自由能计算上。2.2 数据准备与建模的挑战巧妇难为无米之炊。模拟的起点是一个可靠的初始三维结构。对于表观转录组学靶点我们常面临以下情况每种情况对应不同的建模策略情况一有高分辨率的实验结构晶体或冷冻电镜。这是最理想的情况。但需要注意检查完整性晶体结构可能缺失柔性loop区特别是RNA结合界面附近的区域。需要使用建模软件如MODELER, Rosetta补全这些缺失片段。处理质子化状态实验结构通常不含氢原子。对于精氨酸、组氨酸等可质子化的残基以及RNA上的磷酸基团其质子化状态tautomer对相互作用影响巨大。必须使用诸如H服务器、PROPKA或模拟软件自身的工具如AMBER的LEaP模块在生理pH值通常7.4下预测并添加正确的氢原子。这一步极易出错我曾因为一个组氨酸的质子化状态设错导致整个模拟中关键的盐桥无法形成。添加修饰核苷酸这是本领域的特殊之处。如果结构中含有修饰的RNA如m6A需要确保力场参数正确。对于常见修饰AMBER的OL3和OL15力场有现成参数对于非常见修饰可能需要使用GAFF力场结合antechamber生成参数或者进行量子化学计算来拟合参数这是一个专业性极强的步骤。情况二只有同源模型。许多RNA修饰相关蛋白的独特结构域可能没有直接的结构信息。这时需要同源建模。模板选择优先选择与目标序列相似性高、且含有RNA或类似配体复合物的结构作为模板。序列相似性最好30%。关注结合区域建模后要特别审视预测的RNA结合口袋。可以使用分子对接一个小片段RNA或已知抑制剂来检验口袋的合理性。然后对这个复合物进行短时间的如50ns分子动力学模拟观察模型是否稳定特别是结合界面是否会发生不合理的崩塌或漂移。不稳定的模型需要重新评估或采用更精细的循环建模-模拟优化策略。情况三完全从头预测或复合物结构预测。对于全新的蛋白-RNA复合物AlphaFold2及其衍生工具如AlphaFold-Multimer已经展现了惊人潜力。但需要注意置信度评判重点关注预测结构中对结合界面残基的预测局部距离差异测试pLDDT分数和预测对齐误差PAE。pLDDT低70的区域、PAE图中显示界面区域不确定性高的都需要警惕。作为模拟的起点即使AF2预测了复合物也绝不能将其视为绝对真理。它应该作为一个高质量的初始构型必须经过后续分子动力学模拟的“弛豫”和验证。模拟可以优化侧链构象、溶剂化壳层并检验预测界面在物理力场下的稳定性。注意无论结构来源如何在开始生产级模拟前都必须进行充分的能量最小化和平衡。这个过程是为了消除原子间的冲突并使体系适应溶剂环境。平衡阶段要监控体系的温度、压力、能量是否达到稳定。跳过或缩短平衡步骤是导致模拟崩溃或结果失真的常见原因。3. 核心模拟流程与实操要点3.1 分子动力学模拟的标准化流程一旦我们准备好了可靠的初始结构蛋白-RNA复合物或蛋白-小分子复合物就可以进入核心的分子动力学模拟环节。一个标准的、可重复的流程至关重要。以下是一个基于GROMACS的典型工作流其他软件如AMBER、NAMD在思路上是相通的。第一步体系构建与力场选择确定模拟单元明确你要模拟的对象。是Apo蛋白无配体还是蛋白与RNA片段复合物或是蛋白与小分子抑制剂的复合物将相关分子放入一个模拟盒子中。选择力场这是模拟的“物理定律”。对于蛋白质CHARMM36m和AMBER ff19SB是当前广泛认可的优秀力场。对于RNACHARMM36的核酸力场或AMBER的OL3/OL15是专门优化的。对于小分子抑制剂通常使用GAFF2力场其参数通过antechamberAMBER套件或CGenFFCHARMM生成。务必确保所有组分的力场兼容。我个人的习惯是对于蛋白-RNA体系优先使用CHARMM36m CHARMM36核酸力场因为其参数化的一致性较好对于含有机小分子的体系AMBER ff19SB GAFF2的组合更常见。溶剂化与离子中和将体系放入一个足够大的水盒子如TIP3P水模型中盒子边界距离溶质至少1.2 nm。然后添加离子如Na, Cl-以中和体系净电荷并模拟生理离子浓度如0.15 M NaCl。这一步可以使用GROMACS的gmx solvate和gmx genion完成。第二步能量最小化与平衡能量最小化使用最速下降法或共轭梯度法消除原子间的空间冲突。通常进行5000步左右直到最大力小于某个阈值如1000 kJ/mol/nm。NVT平衡在恒定粒子数、体积和温度下进行平衡通常100 ps。使用V-rescale或Nosé-Hoover热浴将体系温度缓慢升至目标温度如310 K即37℃。这一步是让体系温度稳定。NPT平衡在恒定粒子数、压力和温度下进行平衡通常100-200 ps。使用Parrinello-Rahman等压浴将体系压力调节至1个大气压。这一步是让体系密度达到稳定。关键监控必须实时绘制温度、压力、密度、势能随时间变化的曲线确保它们在平衡阶段末期围绕一个平均值平稳波动没有漂移。第三步生产模拟平衡完成后即可开始长时间的生产模拟。模拟时长的选择取决于科学问题考察结合稳定性对于蛋白-配体复合物通常需要100-200 ns来观察配体是否稳定在结合口袋内还是会解离。对于更柔性的体系可能需要更长时间。观察构象变化如果关注蛋白结构域的打开/关闭或别构运动可能需要微秒µs级别的模拟或者采用增强采样方法。设置参数保存轨迹的频率通常为每10-100 ps一帧。积分步长通常为2 fs使用LINCS算法约束化学键。长程静电相互作用使用粒子网格埃瓦尔德方法处理。第四步轨迹分析模拟产生的是海量的轨迹数据分析是关键。核心分析包括均方根偏差衡量蛋白骨架或整个复合物相对于初始结构的整体漂移。快速上升后达到平台期说明体系已平衡。均方根涨落展示每个氨基酸残基的柔性。结合口袋附近的残基通常RMSF较低刚性而loop区RMSF较高。相互作用分析氢键统计整个模拟过程中蛋白与RNA/小分子之间形成的稳定氢键占有率50%。接触面积计算疏水接触面积的变化。相互作用能使用gmx energy或MM-PBSA/GBSA方法估算结合能注意MM-PBSA/GBSA用于绝对结合能误差较大但用于同系列分子相对排序有一定参考价值。可视化检查定期用VMD或PyMOL打开轨迹动画直观观察结合模式、关键相互作用的动态维持情况这是任何定量分析都无法替代的。3.2 增强采样在探索构象空间中的应用经典MD的时长限制了我们观察稀有事件如蛋白大尺度构象变化、配体结合/解离的能力。对于表观转录组学蛋白其功能往往依赖于这种稀有事件。这时就需要增强采样。元动力学是我最常用的一种方法。它的核心思想是“填坑”在预先定义的、能反映我们感兴趣过程的集体变量CV空间里额外添加一个偏置势能将体系从能量最低点“推”出去从而加速对CV空间的探索。定义合适的集体变量这是成败的关键。例如研究一个RNA结合域的开放与关闭可以定义两个结构域质心间的距离作为CV。研究一个别构口袋的开合可以定义口袋入口处几个关键残基侧链二面角作为CV。CV必须能区分不同的状态且与反应坐标相关。设置偏置参数包括高斯势的高度、宽度和沉积频率。参数设置需要一些经验通常从小高度开始测试观察CV空间是否被有效探索。运行与分析运行足够长的元动力学模拟后我们可以根据添加的偏置势能重构出在CV空间上的自由能面。自由能面上的极小值点对应稳定的构象状态鞍点对应过渡态。这能直接告诉我们从状态A到状态B需要克服多大的能垒以及是否存在中间态。实操心得元动力学模拟非常消耗计算资源且对CV的选择极其敏感。一个实用的技巧是先进行一段较长的经典MD如500 ns观察体系自发涨落从中提取可能的关键运动模式可通过主成分分析再用这些模式作为CV的灵感来源。不要试图用一个CV描述所有复杂运动有时需要2-3个CV联合使用。4. 自由能计算指导先导化合物优化当我们通过模拟筛选出几个有潜力的苗头化合物后下一步就是基于结构进行优化提升其结合亲和力和选择性。自由能微扰/热力学积分计算在这里扮演着“计算显微镜”的角色。假设我们有一个先导化合物L其核心骨架的某个位置是一个甲基-CH3我们想通过FEP计算预测如果把这个甲基换成乙基-CH2CH3、氯原子-Cl或者甲氧基-OCH3结合自由能会如何变化ΔΔG。这个ΔΔG直接对应于结合常数Ki的变化ΔΔG ≈ -RT ln(Ki_new / Ki_old)可以指导化学家优先合成哪个类似物。FEP计算流程简述拓扑文件准备为配体L和它的类似物L‘分别生成拓扑文件。对于小的化学修饰通常采用“双拓扑”或“单拓扑”方法在计算中通过一个耦合参数λ将L逐渐“变形”为L‘。所有力场参数必须一致且准确。设置λ窗口将λ从0对应L到1对应L‘分成许多小窗口如16-24个。在每个λ窗口下进行独立的分子动力学模拟。运行模拟与数据收集在每个λ窗口进行充分的平衡和生产模拟收集每个窗口下体系势能对λ的导数或差值。数据分析与误差评估使用Bennett接受率法或MBAR法整合所有窗口的数据计算出从L到L‘的转化自由能变化ΔG。通过重复模拟或块分析评估计算误差。通常一个预测误差在1 kcal/mol以内的FEP计算被认为是成功的这大约对应着结合常数5-7倍的变化。关键注意事项采样充分性这是FEP计算最大的挑战。配体周围的环境水分子、蛋白侧链必须在每个λ窗口都充分弛豫。如果修饰发生在埋藏较深的疏水口袋采样不足会导致结果不可靠。通常每个λ窗口需要5-20 ns的模拟总计算量巨大。误差分析至关重要必须报告计算的标准误差或置信区间。一个没有误差棒的单点ΔΔG值其参考价值有限。与实验交叉验证在项目初期如果有一组已知活性的化合物先用FEP计算一遍看计算预测的ΔΔG趋势是否与实验测得的活性趋势一致。这是验证你整个计算流程力场、参数、设置是否可靠的最佳方式。我们团队在做一个激酶项目时就是通过这种验证发现对小分子中某个磺酰胺基团的参数需要特别处理修正后预测准确性大幅提升。5. 项目实战中的常见陷阱与解决方案在这个领域摸爬滚打踩坑是不可避免的。下面整理了一些典型问题及其应对策略希望能帮你少走弯路。问题一模拟中配体“飞”出了结合口袋。可能原因1初始结构不合理。对接产生的构象本身就不稳定或者质子化状态错误。解决方案尝试多个对接构象作为起始仔细检查并修正质子化状态先用短时间如10ns模拟多个初始构象选择最稳定的一个进行长时模拟。可能原因2力场参数不准确。特别是对于非标准残基或小分子自动生成的参数可能有偏差。解决方案对小分子进行更高级别的量子化学计算如HF/6-31G*来获得更精确的电荷如RESP电荷和扭转角参数。对比使用不同力场如GAFF vs. CGenFF的结果。可能原因3模拟时间不足或平衡不充分。体系尚未达到平衡剧烈的能量波动将配体“弹”出。解决方案严格监控平衡阶段的各项指标确保完全平衡后再开始生产模拟。对于柔性大的体系适当延长平衡时间。问题二MM-PBSA/GBSA计算出的结合能与实验值趋势相反或量级不符。重要认知MM-PBSA/GBSA尤其是基于单轨迹的计算绝对结合能的误差很大通常5 kcal/mol不应用于预测绝对活性。它更适用于同一系列化合物的相对排序。可能原因1熵的贡献计算不准。通常计算熵贡献-TΔS的准简谐近似或正态模式分析误差很大且计算极其耗时。解决方案很多研究只报告焓变ΔH或有效结合能ΔG without entropy。在比较同系物时有时可以忽略熵变因为其差异可能较小。或者使用更精确但昂贵的计算方法。可能原因2溶剂模型和离子强度影响。PB/GB模型对介电常数、离子强度等参数敏感。解决方案保持同系列化合物计算时所有参数完全一致。尝试不同的溶剂模型如GB-Neck2 vs. GB-OBC和离子浓度看趋势是否稳健。最佳实践不要迷信MM-PBSA/GBSA的绝对值。把它当作一个快速的、定性的筛选工具用于从几十个分子中挑出前5-10个最有希望的然后再对这少数分子进行更精确的FEP计算或直接推进实验验证。问题三增强采样模拟“跑偏了”没有看到预期的转变。可能原因集体变量选择不当。CV没有捕捉到真实的反应坐标。解决方案这是增强采样中最艺术的部分。多尝试不同的CV组合。结合主成分分析的结果。有时一个简单的距离或角度可能不够需要更复杂的CV如路径反应坐标或基于机器学习的方法。可以先做短时间的试探性模拟观察偏置势能是否在有效驱动体系变化。问题四计算资源与时间的矛盾。挑战一个微秒级的常规MD或一个多λ窗口的FEP计算在CPU集群上可能需要数周甚至数月。而药物发现项目周期紧张。解决方案层级化策略对大量化合物用快速对接和短MD10-20 ns初筛。对初筛胜出者进行100 ns级MD深入分析。只对最优的2-3个苗头化合物进行FEP计算。利用GPU加速像AMBER、GROMACS、NAMD的新版本都对GPU有良好支持通常能带来10-50倍的加速。投资或租用GPU计算资源是提高效率的关键。增强采样的高效性对于探索构象变化一段成功的100 ns元动力学模拟其采样效率可能相当于数微秒的常规MD合理使用可以节省总机时。最后我想强调的是分子模拟在药物发现中始终是一个“预测工具”和“解释工具”而不是“判决工具”。它的最大价值在于为实验科学家提供高价值的假设和深层次的机理解释缩小实验筛选的范围降低试错成本。将计算预测与湿实验如SPR结合实验、细胞活性测试紧密闭环不断用实验数据验证和校正计算模型是这个领域工作能够产生实际影响力的唯一路径。每一次计算与实验结果的吻合都会增加我们对这套复杂体系的理解也让下一次的预测更加自信。这个过程既有挑战也充满了探索的乐趣。
