KRAS（Kirsten rat sarcoma viral oncogene homolog）蛋白是一种小GTP水解酶，属于RAS蛋白家族，在细胞信号转导过程中发挥重要作用，参与调控细胞的增殖、分化和凋亡等生理过程^［1］。KRAS蛋白在结构上可分为核苷酸结合域（G域）和C末端高变区疏水尾^［2］。核苷酸结合域由第1~165位氨基酸构成，包含P-loop（第10~17位氨基酸）、Switch Ⅰ（第32~38位氨基酸）和Switch Ⅱ（第59~76位氨基酸）等关键区域，负责GTP/GDP的结合和水解，调控下游的信号转导^［2-3］。C末端疏水尾由第166位至最后一位氨基酸构成，负责将KRAS锚定到细胞膜上^［2］。在生理条件下，KRAS在与GTP结合和与GDP结合的状态之间循环，发挥分子开关的作用。KRAS与GTP结合会激活下游的RAF/MEK/ERK和PI3K/AKT/mTOR等信号通路^［4］。在GTP酶激活蛋白的帮助下，KRAS的GTP水解酶活性提高，将GTP水解为GDP，从而失去激活下游信号通路的功能^［5］。鸟嘌呤核苷酸交换因子通过催化KRAS的GDP与GTP交换，使KRAS重新与GTP结合，促进下游信号通路的激活^［6］。这种KRAS-GTP和KRAS-GDP状态之间的循环有序调控着细胞内信号的转导，从而调节细胞的正常生理活动。

然而，KRAS突变打破了正常的GTP与GDP结合循环，造成细胞内信号通路的异常激活，导致肿瘤等疾病的发生^［7］。绝大多数的KRAS突变为单碱基的错义突变，且98%的KRAS突变发生在第12位甘氨酸、第13位甘氨酸和第61位谷氨酰胺等热点突变位点^［8］。这些突变导致KRAS内在的GTP水解能力受损，使KRAS持续激活下游信号通路，促进细胞不断生长和增殖，从而驱动肿瘤的发生发展^［9］。KRAS突变在消化系统肿瘤中最为常见。在胰腺癌中，KRAS突变频率大于90%，其中KRAS G12D（KRAS第12位甘氨酸突变为天冬氨酸）的突变频率最高，约占45%^［10-11］。在胆道系统肿瘤中，KRAS突变频率为20%~30%^［12-13］。在结直肠癌患者中，超过50%的患者携带KRAS突变^［14］。因此，靶向KRAS突变对于肿瘤的治疗具有重要的意义。

MRTX1133是近年来发现的高效、选择性、非共价KRAS^G12D抑制剂，由WANG等^［15］在2021年于美国加利福尼亚州圣地亚哥的Mirati Therapeutics生物制药公司开发，目前在进行用于治疗具有KRAS G12D突变的晚期实体瘤的Ⅰ/Ⅱ期临床试验（ClinicalTrials.gov ID：NCT05737706）。MRTX1133的分子式为C₃₃H₃₁F₃N₆O₂，通过与KRAS蛋白的Switch Ⅱ口袋结合，阻止鸟嘌呤核苷酸交换因子催化的GDP/GTP交换，并且抑制KRAS^G12D/GTP/RAF1复合物的形成，从而抑制依赖突变KRAS的下游信号通路的激活^［15］。细胞实验和动物实验均证明MRTX1133能选择性抑制携带KRAS G12D突变的肿瘤细胞，而对KRAS野生型肿瘤细胞无显著影响^［16］。

靶向治疗目前已成为肿瘤治疗的一个重要策略，然而，耐药突变的出现是靶向治疗面临的一大挑战。在临床实践中，KRAS Y96D（KRAS第96位酪氨酸突变为天冬氨酸）继发突变导致了KRAS^G12C共价抑制剂AMG510对非小细胞肺癌的治疗效果下降^［17］。BCR-ABL1 T315I（BCR-ABL1融合蛋白第315位苏氨酸突变为异亮氨酸）突变则会导致慢性粒细胞白血病细胞对受体酪氨酸激酶抑制剂耐药^［18］。KRAS R68G突变为KRAS蛋白中第68位的赖氨酸突变为甘氨酸。CHOI等^［19］使用饱和突变文库筛选的方法发现KRAS R68G继发突变会导致携带KRAS^G12D突变的肿瘤细胞对MRTX1133产生较高的耐药性，然而其具体机制尚未明确。分子动力学模拟是一种通过数值求解牛顿运动方程、模拟原子和分子在一定时间内的运动和相互作用的方法，具有细粒度的时间和空间分辨率优势，可提供原子水平的位置信息，解析分子行为的动态过程等，在蛋白质和核酸研究、药物设计等领域具有广泛的应用。因此，本研究采用分子动力学模拟的方法，探究KRAS R68G继发突变引发KRAS^G12D靶向抑制剂MRTX1133耐药的具体机制，为克服KRAS靶向药的继发耐药提供理论依据和指导。

从结构生物信息学研究合作组织（Research Collaboratory for Structural Bioinformatics，RCSB）蛋白质数据库（Protein Data Bank，PDB）中获取KRAS^G12D与MRTX1133相互作用的复合物的晶体结构数据（PDB ID： 7RPZ）^［15］，作为G12D体系进行分子动力学模拟的初始结构。使用PyMOL软件将G12D体系中KRAS的第68位精氨酸突变为甘氨酸，作为G12D_R68G体系进行分子动力学模拟的初始结构。由于7RPZ晶体结构中缺失了KRAS蛋白第2位氨基酸的坐标，因此对KRAS蛋白的模拟从第3位的谷氨酸开始至第169位的赖氨酸结束。

使用Amber22软件和AmberTools23软件^［20］中的antechamber程序和键电荷校正（bond charge correction，BCC）的方法计算和构建MRTX1133的力场参数。使用LEaP程序构建和描绘KRAS-MRTX1133复合物的拓扑学参数和位置坐标。在LEaP程序中，依次载入KRAS^G12D/KRAS^G12D/R68G蛋白、GDP、MRTX1133和二价镁离子的初始坐标，形成蛋白-GDP-配体-离子的复合物。应用ff19SB力场^［21］计算KRAS中标准氨基酸的力场参数，应用GAFF（general AMBER force field）力场^［22］计算MRTX1133的力场参数，应用TIP3P（intermolecular potential three point）力场计算水分子的力场参数，GDP和二价镁离子的力场参数从Amber参数数据库中获取（http：//amber.manchester.ac.uk/）。为模拟正常的生理盐水环境，将复合物溶解于正八面体的水盒子中，并加入Na⁺和Cl^-对体系中的电荷进行中和，计算整个系统的拓扑学参数，完成G12D和G12D_R68G体系的构建。

首先，对各体系进行2轮能量最小化。在第1轮能量最小化中，对复合物施加约束力以保持其相互之间的位置关系不变，优化溶剂分子的位置坐标以防止原子之间碰撞的产生。在第2轮能量最小化中，解除对复合物施加的约束力，充分优化体系中各个原子的位置坐标。随后，进行300 ps的分子动力学模拟，将体系从0 K逐步升温至300 K。接着，在等温等容的条件下进行700 ps的体系平衡模拟。最后，在等温等压的条件下，进行200 ns的分子动力学模拟。每个体系重复5次，每次重复时体系中每个原子的初速度均为随机产生。在模拟过程中施加周期性边界，使用particle mesh Ewald方法^［23］计算长程静电相互作用。对于短程静电相互作用和范德瓦耳斯力，设置计算的截止距离为10 Å（1 Å=10^-10 m）。使用SHAKE方法对包含氢原子的共价键进行约束。

均方根偏差（root mean square deviation，RMSD）是衡量2个分子结构相似性的指标，在分子动力学模拟中用来衡量模拟过程中的结构与参考结构的偏离程度。RMSD的计算公式如下：

其中，N为参与比较的原子对数目，rimodel和riref分别是比较结构和参考结构中第i个原子的坐标。RMSD越小，表示2个结构越相似。本研究使用cpptraj轨迹分析软件^［20］中的“rms”命令计算整个模拟过程中KRAS蛋白中所有氨基酸残基的α碳原子（Cα）与初始构象的RMSD。

均方根波动（root mean square fluctuation，RMSF）用来衡量分子动力学模拟中每个原子或每个氨基酸残基随时间波动的幅度。通过RMSF可了解蛋白质或分子中特定区域的柔性或稳定性。RMSF的计算公式如下：

其中，T是总的模拟时间步数，ri(t)是第i个原子或残基在时间t的位置，ri是第i个原子或残基在整个模拟期间的平均位置。RMSF越小，表示稳定性越高，即原子或残基的波动幅度越低；RMSF越大，表示柔性越大，即原子或残基的波动幅度越高。本研究使用cpptraj轨迹分析软件^［20］中的“atomicfluct”命令计算整个模拟过程中KRAS蛋白中所有氨基酸残基的Cα的RMSF。

在轨迹分析中，主成分分析（principal component analysis，PCA）用于研究分子在运动过程中产生的主要构象变化。通过分析体系的协方差矩阵，PCA可以提取出反映分子大尺度运动的主成分，从而帮助理解分子的动力学行为。PCA流程如下：①轨迹预处理：使用cpptraj软件^［20］加载G12D和G12D_R68G体系的运动轨迹，并对体系进行中心化处理，去除平移和旋转。②使用“matrix covar”命令计算KRAS中Cα原子的协方差矩阵。③使用“analyze matrix”命令对协方差矩阵进行特征值分解，得到特征值和特征向量。④使用“projection”命令将轨迹投影到选定的主成分上。在特征向量中，第1个主成分（PC1）提供了最主要的运动模式。将所有运动轨迹按PC1的投影值进行排序，得到体系中KRAS蛋白在运动中主要的构象变化。

溶剂可及表面积（solvent-accessible surface area，SASA）是指分子表面能够被溶剂分子接触到的总面积，用于评估分子的溶解性、稳定性和相互作用等特性。SASA的计算主要基于球形探针模型，核心为将溶剂探针在分子表面滚动，统计所有可被探针接触到的原子表面积。高SASA值表示该分子更多地暴露于溶剂中，更易与溶剂或其他分子发生相互作用；低SASA值表示该分子更倾向于埋藏在疏水区域或受到遮挡。本研究采用cpptraj软件^［20］中的“surf”命令计算MRTX1133和GDP的SASA。

将各个体系在模拟过程中的每一帧构象根据特征向量投影到构象景观图中。采用hieragglo算法对构象景观图中密度较高的群进行聚类分析，使用average-linkage方法衡量不同类别之间的距离，提取所含成员数量最多的聚类的代表性构象作为该群的代表性构象。为消除平移和旋转对聚类的影响，在分析之前将相关帧中KRAS的所有Cα原子进行叠合。

氢键是一种重要的非共价相互作用，由氢原子（H）作为媒介，在2个电负性较强的原子（如O、N、F等）之间形成的定向静电吸引力，可表示为X—H…Y。氢键在分子间和分子内的结构稳定性及功能上起着关键作用。本研究使用cpptraj软件^［20］中的“hbond”命令来计算区域之间氢键形成的数量，采用的氢键定义的标准为：原子X和原子Y之间的距离≤3.5 Å，X—H…Y的角度≥120°。

动态交叉相关矩阵（dynamic cross-correlation matrix，DCCM）描述了在分子动力学模拟中不同原子或残基的位移之间的相关性，用于揭示分子中不同原子或残基之间的协同运动模式。交叉相关系数（cross-correlation coefficient，CC）的计算公式如下：

其中，ri和rj分别表示第i个和第j个原子的坐标，ri和rj分别表示第i个和第j个原子在模拟过程中的平均坐标。CCi,j的取值范围为-1~1。其中，CCi,j>0表示i和j之间存在同向协同运动，CCi,j<0表示i和j之间存在反向协同运动，CCi,j≈0表示两者之间没有显著相关性。本研究采用cpptraj软件^［20］中的“matrix correl”命令计算体系中KRAS蛋白的Cα原子之间的CC。在计算之前对轨迹进行拟合，以消除分子的平移和旋转的影响。

应用R 4.4.0统计软件进行数据分析。符合正态分布的定量数据以x±s表示，2组间比较采用Student-t检验。P<0.05表示差异具有统计学意义。

MRTX1133是KRAS^G12D的选择性非共价抑制剂，其结构式如图1A所示。CHOI等^［19］通过饱和突变文库的方法筛选到了8种不同的KRAS^G12D继发耐药突变，其中有7种氨基酸突变位于MRTX1133的结合口袋。通过对这7种氨基酸突变进行分析，发现有5种突变（V9E，第9位的缬氨酸突变为谷氨酸；V9W，第9位的缬氨酸突变为色氨酸；V9Q，第9位的缬氨酸突变为谷氨酰胺；T58Y，第58位的苏氨酸突变为酪氨酸；Y96W，第96位的酪氨酸突变为色氨酸）均是从侧链较小的氨基酸突变为侧链较大的氨基酸，1种突变（Q99L，第99位的谷氨酰胺突变为亮氨酸）为从极性氨基酸突变为非极性疏水氨基酸，而KRAS第68位的精氨酸突变为甘氨酸（R68G）则是从侧链较长的氨基酸突变为侧链最小的氨基酸。配体结合口袋周围的短侧链氨基酸突变为侧链较大的氨基酸会占据更多的空间，对配体的进入和结合形成空间位阻效应，从物理上阻碍配体的结合（图1B）。口袋周围的极性氨基酸突变为疏水氨基酸，可能破坏重要氨基酸之间的氢键和偶极-偶极等相互作用，或产生空间位阻效应以及改变口袋的原有形状等（图1C），导致配体结合困难。然而，从侧链较长的精氨酸突变为侧链最小的甘氨酸，并不会产生明显的空间位阻效应，其具体机制尚未清楚。因此，本研究采用分子动力学模拟的方法探索KRAS R68G继发突变导致肿瘤细胞对MRTX1133耐药的机制。

为了探究KRAS R68G继发突变对MRTX1133耐药的具体机制，本研究构建KRAS^G12D-MRTX1133（简称为G12D）和KRAS^G12D/R68G-MRTX1133（简称为G12D_R68G）2个体系（图2A）。每个体系均进行了5轮200 ns的分子动力学模拟。RMSD结果显示（图2B），G12D体系的RMSD为（1.19±0.34）Å，G12D_R68G体系的RMSD为（2.09±0.49）Å。G12D_R68G体系的RMSD显著高于G12D体系（P<0.001），表明R68G继发突变对KRAS蛋白骨架的运动产生了明显的影响。为了探索R68G突变对KRAS各区域运动的具体影响，计算各体系中KRAS蛋白中每个残基Cα原子的RMSF。如图2C所示，G12D_R68G体系中KRAS的Switch Ⅰ和Switch Ⅱ区域RMSF明显高于G12D体系。为了更直观地展示G12D_R68G体系和G12D体系RMSF的差异，将两体系RMSF的差值投影到KRAS蛋白结构上。如图2D所示，Switch Ⅰ和Switch Ⅱ区域中的柔性部分明显具有较高的RMSF差值，表明G12D_R68G体系中KRAS的Switch Ⅰ和Switch Ⅱ区域的柔性部分摆动幅度显著大于G12D体系，说明R68G继发突变引起了Switch Ⅰ和Switch Ⅱ区域更大的波动。

为了描绘2个体系全局的构象转变，对整个运动过程的轨迹进行PCA。根据每帧构象在前2个主成分（PC1和PC2）上的投影，绘制G12D和G12D_R68G体系的点密度图（图3A）。结果显示，与G12D体系相比，G12D_R68G体系具有明显较高的PC1值［（7.69±8.46） vs （-7.69±5.52），P<0.001，图3A］。将运动轨迹按照PC1值排序，得到KRAS蛋白随着PC1值增加而产生的构象变化（图3B），即Switch Ⅰ区域和Switch Ⅱ区域的柔性部分均远离P-loop区域而向外运动。综上所述，KRAS R68G继发突变破坏了KRAS^G12D的运动模式，促进KRAS向更加开放的构象转变。

为了进一步探索KRAS R68G继发突变对KRAS Switch Ⅰ和Switch Ⅱ区域运动的影响，选择3个特征向量用于描绘Switch Ⅰ、SwitchⅡ和P-loop 3个区域之间的关系（图4A）。KRAS中第12位天冬氨酸和第34位脯氨酸α碳原子之间的距离（d_D12-P34）代表P-loop和Switch Ⅰ区域之间的距离^［24-25］，其一定程度上反映GDP结合口袋的大小。KRAS中第12位天冬氨酸和第60位甘氨酸α碳原子之间的距离（d_D12-G60）代表P-loop和Switch Ⅱ区域柔性部分之间的距离^［24-25］，其一定程度上反映MRTX1133结合口袋的大小。KRAS中第34位脯氨酸和第60位甘氨酸α碳原子之间的距离（d_P34-G60）代表Switch Ⅰ和Switch Ⅱ区域之间的距离。将两体系的全部构象根据d_D12-P34和d_D12-G60特征向量进行投影，绘制构象景观图（图4B）。图中点密度较高的区域为具有相似d_D12-P34和d_D12-G60聚集的区域，依次将这些区域分类为C1（7.0 Å<d_D12-P34<10.0 Å，8.0 Å<d_D12-G60<10.5 Å）、C2（12.5 Å<d_D12-P34<15.0 Å，8.5 Å<d_D12-G60<9.5 Å）、C3（16.0 Å<d_D12-P34<18.5 Å，8.5 Å<d_D12-G60<9.5 Å）和C4（15.5 Å<d_D12-P34<20.5 Å，10.5 Å<d_D12-G60<13.5 Å）。在G12D体系中，C1占47.68%（为G12D体系的主要构象），C2占7.49%，而C3和C4在G12D体系中几乎没有。在G12D_R68G体系中，C2占6.10%，C3占5.20%，C4占34.69%（为G12D_R68G体系的主要构象），而C1在G12D_R68G体系中占比很少。该结果表明R68G继发突变显著影响了KRAS Switch Ⅰ和Switch Ⅱ区域的运动。C1至C4的代表性构象如图4C所示。在C1的代表性构象中，Switch Ⅰ、Switch Ⅱ和P-loop区域之间的距离较近，MRTX1133和GDP的结合口袋均处于关闭状态。在C2和C3的代表性构象中，MRTX1133的结合口袋均处于关闭状态，而GDP结合口袋在C3中比在C2中有更大的开放程度，且开始出现Switch Ⅰ与Switch Ⅱ分离的趋势。在C4的代表性构象中，MRTX1133的结合口袋处于开放状态，Switch Ⅰ向外运动程度更高，GDP结合口袋开放更大，出现Switch Ⅰ与Switch Ⅱ区域的分离。

为了探究KRAS R68G继发突变对小分子溶剂暴露面积的影响，计算两体系整个运动过程中MRTX1133和GDP的SASA。MRTX1133的SASA如图5A所示。在G12D体系中，MRTX1133的SASA为（149.00±37.30）Ų；而在G12D_R68G体系中，MRTX1133的SASA为（181.00±36.60）Ų（图5B）。两体系之间的差异具有统计学意义（P<0.001）。GDP的SASA如图5C所示。G12D_R68G体系中GDP的SASA明显高于G12D体系［（137.00±29.50） Ųvs （98.60±27.70） Ų，P<0.001，图5D）］。综上所述，KRAS R68G继发突变促进了KRAS Switch Ⅰ与Switch Ⅱ区域向外运动，扩大了MRTX1133和GDP的结合口袋，增加了它们的溶剂暴露面积。

为了探索KRAS蛋白中各氨基酸残基运动之间的相互关系，计算各体系的DCCM。在G12D体系中（图6A），区域a的CC<0，表明Switch Ⅰ和P-loop区域之间运动的负相关性；区域b的CC接近0且>0，表明Switch Ⅱ区域中的柔性部分与P-loop表现出较弱的同向运动关系；区域c的CC均<0，表明Switch Ⅰ和Switch Ⅱ区域之间运动的方向相反。综合Switch Ⅰ、SwitchⅡ和P-loop区域间运动的相关性以及Switch Ⅰ远离P-loop向外运动的结果，可以得知在G12D体系中Switch Ⅱ区域不发生明显的位移或轻微地向P-loop靠近，即MRTX1133的结合口袋基本不变或略微变小。与G12D体系相比，G12D_R68G体系在区域a、b、c的CC均<0，且均明显低于G12D体系中对应区域的CC（图6B），表明在G12D_R68G体系中Switch Ⅰ、Switch Ⅱ和P-loop区域的运动方向均相反，且运动幅度大于G12D体系。因此，G12D_R68G体系中Switch Ⅱ的柔性部分在运动中远离了P-loop和Switch Ⅰ，导致MRTX1133的结合口袋变大，与之前的分析结果相一致。

为了更好地描绘Switch Ⅰ和Switch Ⅱ的位置关系，选取特征向量d_P34-G60代表Switch Ⅰ和Switch Ⅱ之间的距离（图4A）。在整个运动过程中，G12D_R68G体系Switch Ⅰ和Switch Ⅱ之间的距离显著高于G12D体系［（14.50±4.52）Å vs （8.50±2.55）Å，P<0.001，图6C］。

区域之间氢键的形成一方面表明区域之间存在互作关系，另一方面可以反映两区域之间较近的距离。通过观察两体系在运动过程中氢键的形成和数量，发现G12D体系中在Switch Ⅰ和Switch Ⅱ之间有氢键形成，且形成时间较长，而G12D_R68G体系中Switch Ⅰ和Switch Ⅱ在整个运动过程中几乎无氢键的形成（图6D）。通过对所形成氢键的寿命进行分析，发现在G12D体系中，Switch Ⅰ中的E37（第37位谷氨酸）与Switch Ⅱ中的A59（第59位丙氨酸）之间所形成氢键的平均寿命为31.86 ns，占有率为15.93%；Switch Ⅰ中的I36（第36位的异亮氨酸）与Switch Ⅱ中的A59（第59位丙氨酸）之间所形成氢键的平均寿命为0.85 ns，占有率为0.43%。而在G12D_R68G体系中，Switch Ⅰ中的I36（第36位异亮氨酸）与Switch Ⅱ中的A59（第59位丙氨酸）之间所形成氢键的平均寿命为0.01 ns，占有率<0.01%。通过提取两体系的代表性结构并进行比较，发现G12D体系中Switch Ⅰ和Switch Ⅱ之间的距离较近，为7.9 Å，明显小于G12D_R68G体系的13.9 Å（图6E）。在G12D体系中，Switch Ⅰ区域中的第37位谷氨酸与Switch Ⅱ区域中第59位丙氨酸之间形成1个氢键，与Switch Ⅱ区域中第68位精氨酸形成2个盐桥（图6E）。而在G12D_R68G体系中，由于Switch Ⅰ和Switch Ⅱ区域距离较远，两区域之间无氢键和盐桥等极性相互作用形成（图6E）。综上所述，KRAS R68G继发突变破坏了Switch Ⅰ和Switch Ⅱ区域之间的相互作用网络，导致两区域之间的相互作用减弱，造成两区域均远离P-loop向外运动而出现明显分离现象。

KRAS蛋白作为一种分子开关，通过不断地在失活构象（GDP结合状态）和激活构象（GTP结合状态）之间循环，精细地调控细胞内的信号转导，参与细胞的生长、增殖、分化、凋亡等生理过程^［1］。而KRAS特定位点的突变则会打破这种循环，破坏KRAS的GTP水解能力，导致下游信号通路的持续激活，造成细胞的增殖和分裂失控，最终导致肿瘤的发生^［7］。因此，特异性靶向KRAS基因突变能够有效抑制其持续激活的下游信号通路，在肿瘤治疗中具有重要意义。长期以来因KRAS蛋白表面光滑且缺乏小分子结合的深口袋，以及其与GDP/GTP较高的亲和力，认为其不可成药^［26］。直至发现KRAS^G12C中新的变构口袋，并开发出具有较高特异性的抑制剂后，才证明可以在结构上靶向KRAS来开发药物^［26-27］。KRAS^G12C的共价小分子抑制剂AMG510（sotorasib）是第一个被美国FDA批准的KRAS靶向药，其在治疗携带KRAS G12C突变的非小细胞肺癌患者中表现出显著的效果^［28-29］。在临床实践中，TANAKA等^［17］从获得KRAS^G12C抑制剂耐药的非小细胞肺癌患者的游离DNA中发现了KRAS Y96D突变，并通过细胞实验证明了KRAS Y96D继发突变导致细胞对包括AMG510在内的KRAS G12C抑制剂的耐药。ZHUANG等^［25］使用分子动力学模拟的方法进一步解释了Y96D继发突变导致携带KRAS G12C突变的细胞对AMG510耐药的具体机制，为开发新的克服KRAS耐药突变的靶向药提供了理论依据。

MRTX1133是首个成功开发的针对KRAS G12D突变的小分子抑制剂，在胰腺癌细胞和动物模型都取得了良好的治疗效果^［15-16］。目前，MRTX1133用于治疗具有KRAS G12D突变的晚期实体瘤的Ⅰ/Ⅱ期临床试验（ClinicalTrials.gov ID： NCT05737706）正在进行中。除了KRAS特定点突变的抑制剂外，近些年来也成功开发了泛KRAS抑制剂（pan-KRAS inhibitors）^［30-31］。泛KRAS抑制剂可以抑制一系列的KRAS突变体，包括G12A/C/D/F/V/S、G13C/D、V14I、L19F、Q22K、D33E、Q61H、K117N和A146V/T^［30］。然而，泛KRAS抑制剂对于野生型KRAS的抑制作用与突变体相当，可能会非选择性地杀伤人体的正常细胞，导致明显的不良反应^［30］。此外，与KRAS G12D的特异性抑制剂MRTX1133相比，肿瘤细胞对于泛KRAS抑制剂BI-2865的敏感性是对MRTX1133的敏感性的1%~10%^［19］。因此，对于明确的KRAS G12D突变的肿瘤患者，KRAS G12D特异性抑制剂MRTX1133的药物敏感性和不良反应均优于泛KRAS抑制剂。

值得注意的是，CHOI等^［19］使用饱和突变文库的方法筛选到了8种不同的对于MRTX1133耐药的KRAS^G12D继发突变，为克服MRTX1133耐药提供了指导。在这8种继发耐药突变中，绝大多数突变分布在MRTX1133结合的口袋周围，并产生了明显的空间位阻效应。而R68G突变不产生明显的空间位阻效应，其机制尚未明确。因此，本研究采用分子动力学模拟的方法，从原子层面探索和解释了R68G继发突变导致携带KRAS G12D突变的蛋白对MRTX1133耐药的具体机制。

本研究构建了KRAS^G12D和KRAS^G12D/R68G分别与MRTX1133相互作用的体系，并在模拟生理盐水的环境中进行了5轮共1 000 ns的动力学模拟。结果显示，G12D_R68G体系的运动轨迹相对于初始构象的RMSD值显著高于G12D体系，证明R68G继发突变使体系骨架在运动中变得更加松散。通过比较两体系的RMSF，发现KRAS的Switch Ⅰ和Switch Ⅱ区域的波动幅度最大。Switch Ⅰ和Switch Ⅱ为KRAS中负责鸟嘌呤核苷酸交换和GTP水解的关键区域，其构象会不断发生改变以识别和结合下游分子^［3］，这与本研究中其较高的RMSF相一致。PCA证明G12D和G12D_R68G体系的运动轨迹有明显区别，且G12D_R68G体系中KRAS的Switch Ⅰ和Switch Ⅱ的柔性部分均出现远离P-loop的大幅度外向运动，与G12D_R68G体系中Switch Ⅰ和Switch Ⅱ的RMSF明显大于G12D体系中对应的部分这一结果相吻合。DCCM分析进一步证实了G12D_R68G体系中Switch Ⅰ和Switch Ⅱ的这种运动模式。通过测量和比较两体系中Switch Ⅰ和Switch Ⅱ之间的距离和氢键形成数量，发现在G12D体系中，Switch Ⅰ和Switch Ⅱ之间的距离一直处于较小的范围，且在运动的大多数时间内有氢键的形成。在G12D_R68G体系中，Switch Ⅰ和Switch Ⅱ之间的距离波动范围很大，明显高于G12D体系，而Switch Ⅰ和Switch Ⅱ之间始终无氢键形成的结果证明了Switch Ⅰ和Switch Ⅱ之间相互作用减弱，且在运动中逐渐分离。Switch Ⅰ和Switch Ⅱ的外向运动使其大幅度远离P-loop，扩大了GDP和MRTX1133的结合口袋，因此G12D_R68G体系中GDP和MRTX1133的溶剂暴露面积与G12D体系相比明显变大。

值得注意的是，根据CHOI等^［19］的结果可知，KRAS^G12D/R68G突变除了对KRAS^G12D特异性抑制剂MRTX1133耐药外，对泛KRAS抑制剂BI-2865也同样有耐药作用，这可能和泛KRAS抑制剂BI-2865的结合口袋与MRTX1133的结合口袋相同有关。R68G突变打破了Switch Ⅰ与Switch Ⅱ之间的相互作用网络，导致药物结合口袋开放。与R68G突变引起耐药的机制不同，Y96D突变通过破坏KRAS的第96位氨基酸残基与KRAS^G12C抑制剂AMG510之间的范德瓦耳斯力，且引起Switch Ⅱ与α3-helix区域的分离运动，导致AMG510的结合口袋开放和AMG510的结合能力减弱，最终导致AMG510的耐药^［25］。

本研究存在一定的局限性。目前仅KRAS-MRTX1133-GDP复合物的结构被解析，而KRAS-MRTX1133-GTP的结构尚未知，这可能与MRTX1133选择性抑制KRAS^G12D的作用机制有关。MRTX1133通过阻断KRAS与鸟嘌呤核苷酸交换因子的相互作用，从而将KRAS锁定在与GDP结合的状态，阻断其对下游信号通路的激活^［15］。本研究只描绘了R68G继发突变对KRAS失活态中MRTX1133和GDP的影响，无法探究其对KRAS激活态中MRTX1133和GTP的影响。此外，从模拟中所到的KRAS R68G继发突变导致MRTX1133解离速度加快的结论尚需要表面等离子体共振等实验进一步验证。

综上所述，本研究阐明了R68G继发突变导致携带KRAS G12D突变的蛋白对MRTX1133耐药的机制，即R68G突变破坏了Switch Ⅰ和Switch Ⅱ之间的相互作用网络，减弱了Switch Ⅰ和Switch Ⅱ之间的相互作用，导致两区域在远离P-loop向外运动的过程中逐渐分离。Switch Ⅰ和Switch Ⅱ的大幅度向外运动导致GDP和MRTX1133的结合口袋扩大，GDP和MRTX1133的溶剂暴露面积增加，从而加速其与溶剂或其他物质的反应，进而促进其从KRAS结合口袋中解离，最终导致肿瘤细胞对MRTX1133的耐药。这一结果可为克服KRAS^G12D 耐药突变、开发下一代KRAS靶向药以及设计联合用药方案（如靶向药和稳定Switch区域相互作用的分子联用）提供理论指导。