刘岩 孙毅 张桂敏 夏健明 张洪明
摘 要：KCa2.3离子通道蛋白是分布在心肌细胞上的生物大分子，其主要作用是控制细胞内钾离子外流速率，前人试验研究发现该分子的表达情况与血压大小有关，但未给出数量级关系。为深入探究KCa2.3离子通道蛋白与血压之间的关系，本文建立了KCa2.3离子通道蛋白模型，运用分子动力学模拟（molecular dynamics， MD）方法对其进行深入分析，将KCa2.3离子通道蛋白置于不同压力环境下进行了扩散性、体系自由能、稳定性和分子聚集特性分析，计算发现：当压力处于120～150 mmHg时，通过KCa2.3离子通道蛋白的钾离子相对其他压力环境下少，说明压力对KCa2.3离子通道蛋白有一定的影响，二者之间的关系呈“U”型，当压力处于120～150 mmHg时，通过KCa2.3离子通道蛋白的钾离子较少，所引起的生物电流较小;不同环境下体系的势能均处在某一个区间内，说明钾离子外流过程中未发生化学反应，仅为各原子之间势能的转移;不同环境下KCa2.3离子通道蛋白均未发生严重变形，说明压力对KCa2.3离子通道蛋白的稳定性没有影响;钾离子在通过KCa2.3离子通道蛋白时未出现聚集现象。
关键词：生物大分子; KCa2.3离子通道;分子动力学模拟;钾离子外流
中图分类号：Q61;Q71
文献标志码：A
KCa2.3离子通道蛋白是分布在人体心肌细胞隶属于小电导钙离子激活的钾离子通道（small-conductance Ca2+ -activated K+ channel，SK Channel），对钾离子有高度选择性[1]。当细胞质内游离的钙离子浓度上升后会被激活，从而开放钾离子通道，产生钾离子外流，钾离子外流是通过各类钾离子通道实现的，钾离子通道是生物体内分布最广泛、最复杂的一类离子通道[2]，SK离子通道就是钾离子通道中的一类。钾离子外流产生的动作电位会引起超极化后电位（afterhyperpolarization potential， AHP）现象， 超极化后电位对动作电位反复激动现象有抑制作用，以保证细胞的各项生命活动[3]。但是KCa2.3离子通道开放程度过大会使心房肌的动作电位时长明显缩短[4]，从而会引发一些心血管疾病。
近幾年的研究表明，通过减小钾离子外流可以有效减少某些心血管疾病尤其是房颤（valvular atrial fibrillation，VAF）的病发率[5]。我国的孙曼青[6]、苑磊[7]、LI[8]等的研究里表述了SK离子通道在人的心肌细胞里广泛分布，且当血压发生变化时，SK离子通道内的电流也有较明显波动。XU等在2003年发现SK通道分布心房优于心室[9]。TAKAI等在2013年发现KCa2.3离子通道受心脏内湍流切应力调控[10]。因此，对KCa2.3离子通道蛋白表达的进一步深入研究有科学应用价值。
本文基于超级计算机和分子动力学模拟软件Gromacs，模拟计算了KCa2.3离子通道蛋白在不同压力下对钾离子外流的调控程度，分析了不同血压下钾离子外流的速度。
1 计算模型及研究方法
1.1 模拟原理
分子在整个环境中所受力场分为动能和势能，势能包括非键结势能（Unb）、键的伸缩项（Ub）、键角弯曲项（Uθ）、二面角扭曲项（Uφ）、离平面振动项（Uχ）、库伦作用项（Uel），总势能公式如下：
U=Unb+Ub+Uθ+Uφ+Uχ+Uel（1）
依照经典力学，系统中任一原子在环境中所受的力为：
Fi=-Uri=-（ixi+jyi+kziU）（2）
式中，r为体系中原子的位置，通过数值求解牛顿方程模拟原子的运动。
ai=d2ridt2=Fimi （3）
最后得出每个原子的速度、位置信息。
1.2 KCa2.3通道工作机理
KCa2.3离子通道是小电导钙激活钾离子通道的一种，由KCNN3编码，α亚单位和β单位组成，如图1所示：
KCa2.3离子通道蛋白有1个单通电导，并受细胞内Ca2+调控。KCa2.3离子通道作用机理如下：
a）当细胞质内钙离子处于高浓度状态时，钙离子及钙调蛋白（CaM）会结合KCa2.3离子通道的α亚单位，见图2（a）;
b）KCa2.3离子通道缓慢激活、开放，见图2（b）;
c）细胞质内钾离子外流，见图2（c）;从而产生动作电位。
1.3 软件及控制变量的数值选取
本文采用Gromacs软件进行模拟，力场选取全原子（OPLS-AA/L all-atom）力场[11]，水分子模型采用SPC/E模型。体系压力通过Parrinello-Rahman[12]控制，其数值分别设置为60、70、80、90、100、110、120、130、140、150、160、170、180 mmHg;温度耦合由Nose-Hoover[13]控制，数值设置为310 K;静电相互作用采用Particle-Mesh-Ewald算法进行处理;范德华力计算采用截断Cut-off方式，计算截断值设为0.12 nm。
1.4 构建分子模型
构建分子模型步骤如下：
a）在蛋白质数据库下载并修改KCa2.3离子通道的分子模型，修改后的模型如图3所示。左图为蛋白质的全原子模型，右图为蛋白质泛素的模型。
b）利用分子建模软件构建钙离子和钾离子模型;正常生理状态下，细胞内、外钙离子和钾离子的浓度如表1所示。
结合该数据，将一定数量的钙离子、钾离子分别固定在KCa2.3离子通道蛋白的两侧，整个模拟体系在6 nm×6 nm×52 nm的盒子中进行，如图4所示。
图中大的离子表示钙离子，小的离子表示钾离子，体系中间为KCa2.3离子通道蛋白。
c）加入正/负离子（Na+/Cl-）以保证体系的电荷守恒。表2是加入Cl-之后，分子动力学模拟体系内各单位的数量。
d）进行能量最小化模拟，优化初始模型。初步进行能量最小化模拟，采用最速下降法，将最大能量设置在1 000 kJ/mol，防止分子模型在初始阶段受分子间势能影响进行其他运动。
1.5 参数设置
本次模拟仅研究不同压力下KCa2.3离子通道对钾离子的调控作用，设置温度耦合为恒温环境，温度场为36.85 ℃（310 K），步数为1×105，步长为1fs（1fs=1×10-15s），耦合总时长为0.1 ns。
压力耦合设置选用60～180 mmHg共 13组压力值，每组数据步数均为1×105，步长为1fs，耦合总时长为0.1 ns。
进行最终分子动力学模拟时，步数选用1.5×106步，步长为1fs，总时间为1.5 ns，每50 ps保存一次轨迹文件。
2 计算结果
2.1 扩散系数
通过观察钾离子的扩散系数，我们可以直观的看到钾离子通过KCa2.3离子通道蛋白流出细胞外的程度，扩散系数越大，钾离子外流速率越快。而Gromacs软件中并不能直接给出钾离子的扩散系数。因此，我们需要根据其均方位移（mean-square displacement，MSD）函数计算出钾离子在不同压力体系下的扩散系数。
均方位移表示的是某一原子相对起始位置的位移量，其公式为：
MSD=〈|r（t）-r（0）|2〉（4）
再结合爱因斯坦的扩散定律得到，
limt→∞〈|r（t）-r（0）|2〉=6Dt（5）
计算得出钾离子的扩散系数为，
D=MSD6t（6）
式中，D为粒子扩散系数（diffusion constant）且为常数，所以模拟时间t与钾离子的均方位移为一次函数关系。图5为同一体系在不同压力下钾离子的MSD函数，其斜率为钾离子的扩散系数。
从表中我们可以看出在压力在90 mmHg时，MSD函数斜率最大;压力在140 mmHg时，MSD函数斜率最小。因此，當血压为90 mmHg时，钾离子外流速率最大，血压为140 mmHg时，钾离子外流速率最低。根据图5的MSD函数曲线，可以得到不同压力对应钾离子扩散系数，如表3所示。
根据表3绘制压力和钾离子扩散系数的样条曲线如图6。
从图中可以看出压力与钾离子扩散系数的关系呈“U”型，当压力在120～150 mmHg之间时，钾离子扩散系数最低，而当压力小于120 mmHg或者大于150 mmHg时，钾离子扩散系数变大，钾离子外流程度增大，较120～150 mmHg时增大了3.19%～6.12%。该数据与邢爱君统计研究[14]发现的当收缩压处于120～140 mmHg时发生心血管疾病死亡风险最低的结果范围大致一致。
2.2 体系的自由能分析
自由能是影响系统变化的重要因素，其中分子非键作用势能占绝对主导地位，其主要包括分子间的范德华力和分子间的静电力两部分。
Unb=∑i<j∑4εijσijrij12-σijrij6+∑i<j∑qiqj4πε0rij（7）
式中，εij是第i个原子和第j个原子之间相互作用势能的最小值;σij是第i个原子和第j个原子之间达到平衡的距离;rij是第i个原子和第j个原子之间的绝对距离;qi，qj是该原子的电荷量;ε0是介电常数，ε0=1/4πk。
图7是模拟体系在进行压力耦合前做的能量最小化数据，在模拟达到6 325步后，系统趋于平衡，其能量最小值（Epot）为-7.215×106 kJ/mol。此时，体系达到最稳定的状态，与模拟过程中体系内相互作用能进行对比（图8），在模拟过程中能量最大为-7 140 348.5 kJ/mol，能量最小为-7 224 514.5 kJ/mol，波动范围在84 116 kJ/mol之间。因此，随时间变化体系内的自由能稳定性能未发生改变，说明在KCa2.3离子通道工作时未发生化学反应，且各原子动能未发生较大变化，各原子的运动速度是一个定值。
2.3 KCa2.3离子通道在不同压力下的稳定性
均方根偏差（root mean square deviation，RMSD）是衡量模拟过程中蛋白质稳定性的重要指标，当模拟后的蛋白质位置与初始状态的位置偏差较大时，RMSD值就会越大，其计算公式为：
RMSD=1N∑Ni=1（rij-ri0）2（8）
式中，rij表示j时刻对应i原子的坐标位置;ri0表示初始时刻原子所在的位置。
图9、图10分别为体系在压力为60 mmHg和180 mmHg时的RMSD变化趋势。从图中我们可以看出，随着模拟时间的增加，蛋白质与其初始位置的距离也变大，且均在模拟0.8 ns后，其结构趋于稳定，其位移大致为0.45 nm和0.49 nm。体系在其他压力状态下KCa2.3离子通道蛋白的位移均在0.45～0.49 nm范围内，对比心肌细胞直径的15 μm可以忽略不计，说明压力对KCa2.3离子通道蛋白的稳定性没有影响或是影响较小。
为了更直观地观察不同压力下KCa2.3离子通道蛋白的变化情况，我们将模拟前后的KCa2.3离子通道蛋白构型放到同一体系内进行观察对比，如图11所示。
图11（a）、11（b）中较薄的构型为体系模拟前的分子构型，图11（c）中较薄的构型为体系在60 mmHg压力下的分子构型。通过对比发现KCa2.3离子通道蛋白在模拟前后，结构上有细微差异;而60 mmHg压力下与180 mmHg压力下KCa2.3离子通道蛋白结构没有发生变化，仅在模拟后的位置上有细微差异，进一步验证了压力对KCa2.3离子通道蛋白的稳定性没有影响。
2.4 分子聚集特性
径向分布函数（radial distribution function，RDF）其物理意义如下：
设中心距离为r到r+dr圆环内的原子数目为dN，RDF函数g（r）为
ρg（r）4πr2=dN （9）
解得，
g（r）=dNρ4πr2dr （10）
当RDF值越接近1时，表示该区域范围内的密度越接近平均密度。我们可以根据径向分布函数查看钾离子之间最短距离，对比不同体系中钾离子之间的径向分布函数，发现其分布函数没有较大差异。图12是在100 mmHg压力下体系以钾离子为中心，距离中心r到dr的圆环内钾离子数目的径向分布函数。
在0.338 nm处RDF函数值开始不为零，即在该系统下两个钾离子之间最短距离为0.338 nm，当距离为0.434 nm时，径向分布函数值达到最大值为1.631。在其他压力情况下，当中心距离越大时，RDF也都是趋近于1。说明钾离子在外流的过程未出现聚集现象，钾离子在整个模拟过程中都是均匀分布在细胞液内。
3 结论
为了研究KCa2.3离子通道蛋白在不同压力环境下各方面性能的差异，本文采用分子动力学模拟的方法，对同一体系在13种不同压力下进行了分子动力学模拟。通过分析和对比KCa2.3离子通道蛋白在表达时不同压力状态下钾离子的扩散系数、体系内的势能与KCa2.3离子通道蛋白的RMSD函数、体系能内势能的变化、钾离子的径向分布函数，得出以下结论：
1）KCa2.3离子通道蛋白受压力影响，二者之间的关系呈“U”型，当压力处于120～150 mmHg时，通过KCa2.3离子通道蛋白的钾离子数量较其他压力状态下少3.19%～6.12%;
2）KCa2.3离子通道蛋白在表达时未发生化学反应，仅为各原子之间势能的转移;
3）压力对KCa2.3离子通道的稳定性没有影响或是影响较小;
4）KCa2.3离子通道蛋白在表达时，钾离子未出现聚集现象，其在整个模拟过程中都是均匀分布在细胞液内。
参考文献：
[1]HIRSCHBERG B， MAYLIE J， ASELMAN J P. Gating of recombinant small-conductance Ca2+-activated K+channel by calcium[J]. Gen Physiol， 1998，111（4）：565-581.
[2] 傅涛，郝选明.钾离子通道的研究进展[J].内蒙古石油化工，2010，36（10）：1-4.
[3] LANCASTER B， NICOLL R A， PERKEL D J. Calcium activates two types of potassium channels in rat hippocampal neurons in culture [J]. Neurosci， 1991，11（1）：23-30.
[4] ZHANG X D， VALERIY T， LI N，et al . Critical roles of a small conductance Ca2+-activated K+ channel （SK3） in the repolarization process of atrial myocytes[J]. Cardiovasc Research， 2014，101（2）：317-325.
[5] DINESS J G， SORENSEN U S， NISSEN J D. Inhibition of small-conductance Ca2+-activated k+ channels terminates and protects against fibrillation [J]. Journal of the European Ceramic Society， 2017，37（5）：2179-2187.
[6] 孙曼青，钱玲玲，王如兴.小电导钙离子激活钾离子通道与心律失常的发生机制[J].中华心血管病杂志，2015，43（2）：187-189.
[7] 苑磊.小电导钙激活钾（SK）通道-心房选择性抗心律失常药物新靶点[C]//2014第十二届国际新药发明科技年会. 苏州：
百奥泰国际会议（大连）有限公司， 2014.
[8] LI M L， LI T， LEI M， et al. Increased small conductance calcium-activated potassium channel （sk2 channel） current in atrial myocytes of patients with persistent atrial fibrillation[J]. 中華心血管病杂志，2011，39（2）：147-151.
[9] XU Y F， TUTEJA D， Zhang Z， et al. Molecular identification and functional roles of a Ca2+-activated K+ channel in human and mouse hearts[J]. Journal of Biological Chemistry， 2003， 278（49）：49085-49094.
[10]TAKAI J， SANTU A， ZHENG H F，et al. Laminar shear stress upregulates endothelial Ca2+-activated K+ channels KCa2.3 and KCa3.1 via a Ca2+/calmodulin-dependent protein kinase kinase/Akt/p300 cascade[J]. AJP Heart and Circulatory Physiology，2013，305（4）：
484-493.
[11]SHOWALTER S A， BRUSCHWEILER R. Validation of molecular dynamics simulations of biomolecules using nmr spin relaxation as benchmarks：
application to AMBER99SB force field [J]. J Chem Theory Com put， 2007，3（3）：
961-975.
[12]CARRETERO G R. A Parrinello-Raham approach to vortex lattices[J]. Physics Letters A， 2005，341（1）：128-134
[13]FAOU E. Nose-Hoover dynamics in a shaker[J]. J Chem Phys，2006，124（18）：184104-184106.
[14]邢愛君，赵权辉，马丽梅，等.收缩压水平对非瓣膜性心房颤动患者死亡的影响[J].中华心血管病杂志，2021，49（3）：236-241.
（责任编辑：于慧梅）
Molecular Dynamics Simulation of the Effect of Blood Pressure
Regulation of KCa2.3 ion Channel on Valvular Atrial Fibrillation
LIU Yan1， SUN Yi2 ， ZHANG Guimin2， XIA Jianming2， ZHANG Hongming*1
（1.Department of Force， Kunming University of Science and Technology， Kunming 650500， China;2.Fuwai Yunnan Cardiovascular Hospital， Kunming 650032， China）
Abstract：
KCa2.3 ion channel protein is a biomolecule distributed in cardiomyocytes， and its main role is to control the intracellular potassium ion outflow rate. Previous experimental studies have found that the expression of this molecule is related to blood pressure， but the order of magnitude relationship has not been given. In order to further explore the relationship between KCa2.3 ion channel protein and blood pressure， the KCa2.3 ion channel protein model was established and analyzed by Molecular Dynamics （MD） method. The diffusion， system free energy， stability and molecular aggregation characteristics of KCa2.3 ion channel protein were analyzed under different pressure environments. The calculation results show that：
When the pressure is 120～150 mmHg， less potassium ions pass through KCa2.3 ion channel protein compared with other pressure environments， indicating that pressure has a certain influence on KCa2.3 ion channel protein， and the relationship between the two is in the shape of "U". When the pressure was 120～150 mmHg， less potassium ions passed through KCa2.3 ion channel protein， resulting in less biological current; The potential energy of the system in different environments is all in a certain interval， indicating that no chemical reaction occurs in the process of potassium ion outflow， only the transfer of potential energy between atoms; KCa2.3 ion channel protein did not undergo serious deformation under different environments， indicating that pressure had no effect on the stability of KCa2.3 ion channel protein; There is no aggregation of potassium ions through KCa2.3 ion channel protein.
Key words：
biomacromolecule; KCa2.3 ion channel; molecular dynamics; potassium outflow
猜你喜欢势能外流原子原子究竟有多小？少儿科学周刊·儿童版(2021年22期)2021-12-11原子可以结合吗？少儿科学周刊·儿童版(2021年22期)2021-12-11带你认识原子少儿科学周刊·儿童版(2021年22期)2021-12-11人力资源流动对农村经济发展的影响农家致富顾问·下半月(2018年5期)2018-10-21不外流领导文萃(2018年13期)2018-08-13势能的正负取值及零势能面选择问题初探高中生学习·高三版(2017年9期)2017-10-26“动能和势能”“机械能及其转化”练习中学生数理化·八年级物理人教版(2016年5期)2016-08-26弹性势能纵横谈新高考·高一物理(2015年3期)2015-08-20关于重力势能和弹性势能理解与运用的几个典型错误新高考·高一物理(2014年4期)2014-09-17
