首页/文章/ 详情

基于GROMACS与PLUMED的无序蛋白增强采样仿真

9分钟前浏览0

关键词:GROMACS;PLUMED;增强采样;无序蛋白;构象系综

一、文章简要介绍

有一类蛋白质没有固定三维结构,教科书上的折叠图案对它们不适用——内在无序蛋白(IDP)靠动态构象系综行使功能,常规分子动力学在这种体量的构象空间里很难采够样本。这篇哥本哈根大学团队2026年发表在Nature Communications的工作,用GROMACS搭配PLUMED的OPES增强采样,单条模拟副本就覆盖了多温度系综,还在转录共激活因子ACTR里抓到了以往看不见的瞬时三级结构(DOI:10.1038/s41467-026-73067-3)。

图1 OPES多温度系综采样原理与四个无序蛋白体系

OPES单副本覆盖多温度系综,四个无序蛋白采到从紧致螺旋到伸展的构象谱。

二、GROMACS与PLUMED增强采样建模

计算框架是GROMACS 2023.4加PLUMED 2.10,GPU加速运行。增强采样用OPES(在线概率增强采样):把体系势能当作集 合变量,偏置势随模拟在线更新,单条轨迹就能在274到450开尔文的多温度系综里游走,不需要温度副本交换,也省去了模拟回火里估权值的麻烦。体系选了四个:15残基的(AAQAA)3螺旋肽、42残基的ACTR片段、71残基的ACTR和有聚集倾向的HTTex1外显子1,力场用a99SB-disp与a03w,水模型TIP4P/2005,盐浓度150 mM NaCl。

图2 (AAQAA)3肽15微秒采样OPES与REST2对比

(AAQAA)3肽15微秒采样,OPES与REST2的Rg分布高度重合。

工作流分两步:先跑10条短试验模拟让偏置势收敛,小体系约100纳秒、最大的ACTR约600纳秒收敛完成;再把收敛的偏置势当静态偏置,挂到多条独立生产模拟上提吞吐。副本之间不交换,各自在温度空间里按势能涨落扩散,各温度的有效样本量近似均匀。

图3 试验模拟偏置势收敛与统计误差评估

10条试验模拟确认偏置势收敛,小体系约100纳秒、ACTR约600纳秒。

采样效率用(AAQAA)3做了定量对比:同为15微秒总量,OPES的螺旋弛豫时间4.9纳秒比REST2的12.3纳秒短,统计误差收敛快2到3倍;Rg分布与逐残基螺旋倾向两条方法高度重合,说明提速没有牺牲精度。HTTex1上OPES同样追平REST2,部分观测量的收敛快到8.8倍。

图4 HTTex1外显子1的OPES与REST2对比

HTTex1外显子1上OPES与REST2一致,部分观测量收敛快至8.8倍。

三、关键结果与发现

主角是ACTR——与CBP的NCBD结构域结合时,双方各自折叠成三条螺旋。OPES模拟里,游离态ACTR的三条螺旋会协同折叠并形成三级接触,这类稀有序态在轨迹中可逆地出现又消失。

图5 ACTR的OPES与REST2系综对比

ACTR三条螺旋协同折叠形成三级接触,稀有序态可逆出现。

构象靠不靠得住,要拿实验数据说话。团队用贝叶斯/最大熵重加权把化学位移、RDC、PRE和SAXS一起压进系综:化学位移RMSE从1.06降到0.68 ppm,SAXS曲线与PRE标记位点的强度分布都对上,回转半径从1.90纳米修正到2.19纳米。螺旋束状态只占约1%概率,疏水核心由异亮氨酸、亮氨酸、缬氨酸撑起,结构上说得通。

图6 ACTR瞬时三级结构与实验数据重加权验证

贝叶斯重加权后化学位移、RDC、PRE与SAXS全面吻合实验。

自由能面上看得更直观:以天然接触分数Q和Rg为坐标,OPES与无偏模拟、REST2的盆地位置一致,还额外探到高Q的小盆地;沿主要坐标的过渡势垒约13 kJ/mol,折叠-展开在模拟时间尺度上可逆往返。

图7 ACTR的Q-Rg自由能面与方法对比

Q-Rg自由能面:OPES与无偏、REST2盆地一致,势垒约13 kJ/mol。

聚类分析把系综压成可数的构象族:3埃截断下聚出7282个簇,前24个簇就覆盖一半概率,主簇构象与结合态的螺旋1高度相似——这是构象选择机制的采样证据,也给靶向无序蛋白的药物设计留了抓手。

图8 ACTR构象聚类与主簇结构

3埃截断聚出7282个构象簇,前24个簇覆盖一半概率。

来源:320科技工作室
ACT化学UM分子动力学试验METAGROMACS
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2026-11-14
最近编辑:9分钟前
320科技工作室
硕士 | 结构工程师 微❤️ CAE320
获赞 243粉丝 449文章 532课程 0
点赞
收藏
作者推荐

黄永刚晶体塑性UMAT子程序相关教学

(1)FCC/BCC改成HCP黄永刚本构模型适用于立方晶体,FCC/BCC,由于三个轴是相互正交的,通过输入三指数来描述,(hkl)用于晶面,[uvw]用于晶向在黄永刚的UMAT中可以根据输入的三指数,自动生成所有的滑移方向和滑移面法向。但是该代码不适用于密排六方的本构晶体。对于密排六方晶体(HCP)结构,滑移面和滑移方向通常采用四指数法来描述:(hkil) 用于晶面, [uvtw] 用于晶向本教学提供基于输入的四指数自动生成所有滑移方向和滑移面法向的代码,适用于HCP的晶体塑性计算,并提供材料参数说明。(2)加JC本构Johnson-Cook本构模型用于描述材料在动态载荷下的本构行为,其核心思想在于材料的流动应力取决于塑性应变、应变率和温度三者耦合作用。可以在黄永刚晶体塑性的UMAT中加入JC本构模型,使该本构模型考虑应变率硬化和温度软化。修改方法在于修改硬化模型,其硬化模量计算公式如下:上述模型中前两项为Bassani硬化模型,黄永刚程序里有,第三项为JC本构中的应变率硬化和温度软化项。修改后的本构计算效果如下图所示。图 { SEQ 图 \* ARABIC |1} 本构模型应力-应变曲线,(a) 单晶压缩模型,(b)不同应变率、温度的应力应变曲线(3)Voronoi多晶结构进行材料随机取向赋值脚本晶体塑性本构用于描述单晶的行为,然而大多数材料是由多晶构成的,采用晶体塑性UMATs进行多晶的仿真,可以借助Voronoi多晶模型,即生成具有多晶结构的代表性体积单元、模型。如下图所示。图 { SEQ 图 \* ARABIC |2} 多晶结构的代表性体积单元、模型当晶粒较多时,采用人工给不同晶粒逐一设置材料取向信息过程繁琐,可以采用脚本进行取向赋值,脚本支持随机取向或者具有织构的取向。(4)在滑移的基础上考虑孪晶系在传统晶体塑性理论中,除了常规滑移变形外,还引入孪晶作为另一种独立的塑性机制。当考虑了孪晶后,模型变为:这类模型还可以拓展为形状记忆合金相变+塑性的晶体塑性本构。(5)修改VUMATUMAT主要应用于隐式分析,而对于大变形接触问题,隐式分析往往计算效率较低。对于接触、碰撞、冲击等问题采用VUMAT往往具有更高的计算效率和收敛速度。可以将UMATs修改为VUMAT,用于隐式分析。本教学提供修改方法、代码及案例。(6)添加损伤本构晶体在滑移的作用下如何产生永久变形。采用黄勇刚晶体塑性本构可以模拟出晶体在恒定载荷下产生的蠕变。由于塑性变形进一步地会出现刚度下降,承载力下降的情况。标准的晶体塑性模型,无法描述这一现象,但如果在此基础上引入连续损伤本构模型,定义损伤因子,并将连续损伤因子与切应力关联,可以实现对材料损伤的仿真。其建模思路如下图所示:图 { SEQ 图 \* ARABIC |3} 晶体塑性蠕变损伤本构建模思路 (7)添加背应力背应力是描述运动硬化和包辛格效应的关键机制,背应力演化公式:其中,是第个滑移系的背应力变化率。添加背应力前后,循环曲线如下图所示:(8)添加蠕变和损伤蠕变模型公式:第一项:纯蠕变贡献,第二项,损伤耦合贡献。添加蠕变前后,计算的位移-时间曲线如下:另外,该模型也可以用于计算晶体塑性低周疲劳,载荷谱及计算结果如下图所示: (不考虑蠕变) (考虑蠕变)(9)晶体塑性损伤+单元删除切削案例(10)位错密度晶体塑性本构模型在黄永刚晶体塑性框架的基础上,进一步开发了位错密度晶体塑性本构子程序(Abaqus UMAT):位错密度以显式状态变量引入,12 个刃型与 6 个螺型几何必需位错(GND)分量逐增量求解,统计存储位错(SSD)走独立演化通道,硬化律采用位错森林投影而非总位错密度。相比只有标量硬化的唯象模型,它能回答三个问题:晶粒尺寸效应从何而来、织构演化多强、晶界处位错如何累积。(1)晶粒尺寸效应。三套 Voronoi 多晶模型仅改变晶粒数(平均尺寸 6.76 / 4.35 / 3.48 μm),单轴拉伸至 10% 应变。细晶的流动应力提升 43%,而总位错密度几乎不变(比值 0.960),位错摩擦应力却达 1.393 倍——细晶强化的载体是 GND 构成的位错森林,不是位错总数;应力与 1/d 的平方根呈线性关系,可据此定量估计 Hall-Petch 斜率。(2)织构演化。拉伸至 30% 应变,取向由每个增量步解出的滑移量积分更新,可直接从 ODB 抽取后绘制极图、计算 ODF。织构指数 J 由 10% 应变时的 +1.8 / +4.8 / +4.6 跃升到 +3.8 / +24.1 / +31.5,76% / 89% / 85% 的晶粒转动超过 5°;初始取向按构造均匀随机,长出来的结构完全是变形诱导的。(3)晶界位错累积。晶界两侧取向不同、滑移系不同,塑性应变梯度的跳变在晶界处造成 GND 累积。18432 个 C3D8 单元、20 个 Voronoi 晶粒拉伸至 5% 应变:内部晶界的 GND 密度达 170.3 μm-2,深部晶内仅 50.0 μm-2(3.4 倍),而自由表面只高 1.36 倍,说明堆积主要归属晶界而非表面效应;SSD 由滑移量决定,晶界处只有晶内的 0.34 倍。同一网格、同一时刻,两类位错给出完全不同的空间形貌。上述(1)-(10)的修改可以根据用户的材料特性,科研需求,组合定制本构子程序,可以提供教学+代码+案例,或者单独提供代码+案例。来源:320科技工作室

未登录
还没有评论
课程
培训
服务
行家
VIP会员 学习计划 福利任务
下载APP
联系我们
帮助与反馈