2014年原环境环保部和原国土资源部发布的《全国土壤污染状况调查公报》显示，受快速的工业化和城市化进程影响，我国土壤污染点位超标严重，其中主要污染物为重金属和类金属^[1]。农作物吸收污染土壤中的重金属不仅会对农作物本身造成减产等危害，还会通过食物链传递对人体健康造成不良影响。植物根表质膜是保持农作物内部环境相对平衡稳定的屏障，也是土壤溶液中自由态重金属离子首先接触、结合并进入细胞被农作物吸收乃至产生毒性作用的主要通道^[2]。因此，了解重金属离子与植物根表质膜之间的相互作用规律，是从根本上理解重金属在土壤-植物体系中迁移机制的关键之一。目前，大多数研究是从宏观上的生物配体模型^[3-5]和分子生物学上单纯的金属通道蛋白^[6-8]等角度探讨植物根表质膜表面重金属离子吸附、吸收等过程以及其对植物毒性作用的影响。但重金属离子与细胞质膜相互作用的微观过程，很难用实验手段观测到，因而缺乏分子水平上对其变化过程的描述。

细胞质膜因其磷脂分子头部的磷酸基团、膜蛋白上所带的氨基酸残基等基团^[9-11]的存在，使得质膜表面带负电荷，引起植物根表质膜-土壤溶液界面上电势的变化，影响着重金属离子在质膜表面的吸附及土壤溶液中的分布。Kinraide等^[12-13]将溶液化学中经典的双电层Gouy-Chapman-Stern模型（GCS模型）引入植物根表-土壤溶液体系，建立了表征植物根细胞膜表面离子活度的根表面络合模型。Kinraide等^[13-14]、Wang等^[15-16]和Wang等^[17]的实验和模拟数据表明GCS模型可以用于量化Cu²⁺、Cd²⁺、Ca²⁺、Mg²⁺等阳离子和砷酸根等阴离子对质膜膜表面电势（Ψ₀）的影响和实际膜表面电荷分布情况。经典GCS模型能够成功地揭示膜表面电势影响下不同离子竞争吸附的本质规律，即细胞膜表面分布着负电和电中性等不同结合位点，不同离子在这些位点上有着不同的结合常数，离子在这些结合位点的竞争吸附改变了膜表面电势，并控制离子在扩散层中的分布。但这一概念化数学模型为了简化运算将细胞质膜视为刚性平板电容，而实际的细胞质膜是结构组成复杂的两性双分子层，表面官能团较多，具有流动性等动态特性，因此很难通过实验手段和这一经典模型获取分子水平上重金属离子与质膜相互作用的微观动态和机制。分子动力学（Molecular dynamics，MD）已被证明可以有效模拟细胞质膜与不同分子、离子间的相互作用，可以有效地弥补常规实验和结构分析方法对分子水平上界面动态变化认识的不足^[18-22]。分子动力模拟可以在分子水平上描述重金属离子在细胞质膜表面的分布，以及不同官能团如羰基（C=O）、羟基（-OH）、磷酸基团等对结合构型的影响；模拟得到的运动轨迹可以定量分析质膜表面电势的动态变化，以及其对溶液中离子浓度分布的影响，有助于深入了解分子尺度下重金属离子在实际生物质膜中吸附、穿透和吸收等过程。因此，本文采用分子动力模拟方法：（1）从分子层面解释质膜表面电势对溶液中重金属离子吸附和分布的影响；（2）与经典双电层GCS模型结合，探究不同阳离子或阴、阳离子共存情况下，质膜结构和表面电势的变化规律，从而更好地描述质膜表面电势对溶液中重金属离子吸附和分布影响的分子机制。

1 材料与方法 1.1 模拟体系

模拟体系由磷脂双分子层（图 1）和溶液相构成。溶液分别由约0.05 mol·L^-1的Na⁺、K⁺等一价阳离子，或0.025 mol·L^-1的Ca²⁺、Mg²⁺、Cd²⁺、Cu²⁺等二价阳离子和53 790个水分子组成：磷脂双分子层由114个POPC（1-palmitoyl- 2- oleoyl-sn-glycero-3-phosphocholine）和14个POPG（1 - palmitoyl-2-oleoyl-sn - glycero-3-phosphoglycerol）组成^[23-25]（图 2）。模拟体系设定的周期边界条件PBC（Periodic boundary conditions）：盒子大小为64 Å×64 Å×170 Å，溶液相中另外补充适当的Cl^-，以保持盒子整体的电中性。所有模拟计算由美国Sandia国家实验室开发的开源性大规模原子分子并行模拟器LAMMPS（Large-scale atomic /molecular massively parallel simulator）软件完成，计算体系中原子间相互作用力的力场选用CHARMM 36^[26-27]，平衡后体系中磷脂双分子层表面可变电荷密度（σ₀）约为25.4 mC·m^-2，与实验观测结果相符。

1.2 平衡方法

建立好的模型体系先用正则系综NVT（N表示原子数，V表示盒子体积，T表示温度，三者不变）平衡，初始温度为300 K，平衡时间1.5 ns，再进行模拟退火0.8 ns，温度由300 K逐步升至450 K再逐步降至300 K，最后采用等温等压系综NPT（P表示压力恒定），压力为0.1 MPa，弛豫时间为1.5 ns。

1.3 模拟过程

平衡后的模型体系继续选用非均向性的NPT系综进行计算并收集数据分析，体系温度设定为300 K，体系压强为0.1 MPa，积分步长设置为2 fs，为保证体系总模拟时长可代表实际情况，模拟过程采用LAMMPS中的SHAKE算法约束53 790个水分子的振动（键长或键角）以达到延长时间步长的目的，体系总模拟时长为60 ns。

1.4 GCS模型计算磷脂表面离子浓度分布

设定质膜双分子层中心为Z=0，GCS理论计算溶液中阳离子Mⁿ⁺和Cl^-在体系中的密度分布公式为^[18]：

(1)

式中：n_b表示溶液中金属离子的浓度；κ_Z=（8πn_be²/ εk_BT）^1/2为Debye-Hückel参数，其中ε表示水的介电常数，k_B为玻尔兹曼常数，T为温度；Φ_s表示Z=0处的膜表面静电势，由式（2）计算得出：

(2)

2 结果与讨论 2.1 质膜表面的离子吸附构型

图 3分别为不同溶液体系中Na⁺、K⁺、Ca²⁺、Mg²⁺、Cd²⁺、Cu²⁺等阳离子及阴离子Cl^-在质膜双分子层表面的吸附模型。模拟结果表明，K⁺在质膜表面吸附较弱（图 3b），而其余阳离子都倾向于与质膜表面位点强烈结合[如羰基（C=O）与磷酸基（-OP）]。这是因为与Na⁺、Ca²⁺、Mg²⁺、Cu²⁺、Cd²⁺相比，K⁺的离子半径更大，离子表面的电荷密度更小，与周围水分子的结合能力更弱，水合半径更小，进而KCl体系的离子强度也比其他体系低，而离子强度的大小又是影响细胞质膜表面电势绝对值的关键。因此K⁺吸附并进入细胞质膜的自由能垒低于Na⁺以及二价金属阳离子^[29]。

除K⁺外，其他阳离子在结合到质膜表面时，倾向于与质膜头部的极性基团羟基、磷酸氧和羰基结合（图 3内插图中内容所示），主要是这些官能团的负电荷与阳离子发生静电络合作用。金属阳离子与质膜头部基团的结合导致了一个与磷脂双分子层自身方向相反偶极矩的生成。质膜头部带有正电荷的胆碱基团，也因为与阳离子之间的排斥作用以及与阴离子的吸引作用发生一定的构型改变，倾向于转向水相或者与磷酸基团呈相反的方向。金属阳离子吸附导致了分子构型变化进而促使磷脂面积减小和质膜厚度增加，促使质膜更加紧密的包裹^{[19, 30-31]}。尽管自然状态下根细胞质膜呈负电性，通过静电排斥作用导致阴离子Cl^-通常分布在溶液中，随着距离质膜表面越远，阳离子的浓度越低，阴离子Cl^-的浓度越高。但图 3c的Ca²⁺体系吸附构型还显示出在Ca²⁺与质膜表面的官能团结合时有Cl^-伴随，证明了溶液相的阴离子可以通过阳离子桥键结合的形式吸附到质膜表面。

2.2 磷脂表面密度分布

不同体系中阳离子、阴离子、水分子和质膜双分子头部的羰基氧、羟基氧、磷酸氧和氮原子在质膜表面密度分布情况如图 4表示。原子的位置对脂质双层结构至关重要。分布峰之间的距离可以直接反映质膜极性基团构型的变化。以质膜中心作为Z 轴的原点（Z=0），图中显示的为质膜双分子层的一半，密度数值显示的是质膜两侧的分布均值。在膜-水界面的周围，K⁺的密度分布没有呈现出明显的峰值（图 4b），这是因为K⁺半径较大，电荷密度低，与负电基团亲和性较弱^{[19, 32]}。Na⁺及二价金属离子体系中（图 4a和图 4c至图 4f深蓝线），阳离子在靠近质膜表面的位置出现一个分布峰，其峰值出现的位置和位于脂质的羰基（Mⁿ⁺-Ocarb）和磷酸氧（Mⁿ⁺-P）最大峰处的分布具有高度的相关性，羟基基团（Mⁿ⁺-OH）密度分布的最大值与金属离子概率密度的峰值重叠部分较少，说明羟基基团对吸附Ca²⁺、Mg²⁺、Cd²⁺、Cu²⁺阳离子的贡献比羰基和磷酸基团少，这一结论也验证了不同金属离子在质膜表面吸附构型的不同。同时我们观察到氮原子（Mⁿ⁺-N）与磷原子（Mⁿ⁺-P）的密度分布相似，并且受阳离子结合的影响方式相似。这种构型和电荷分布的相关变化是由于胆碱氮和阳离子之间的相互排斥作用以及阴离子和胆碱基团之间的吸引力引起的^[19]。

除K⁺以外，Na⁺和二价金属离子Ca²⁺、Mg²⁺、Cu²⁺、Cd²⁺都在水分子（Mⁿ⁺-H₂O）达到最大分布密度峰附近出现一个倒峰，水分子分布出现最大值的位置为质膜吸附金属阳离子紧密层（双电层中的Stern层）的外侧，二价金属离子在从质膜中心到体相溶液的分布过程中（距质膜双分子层中心距离由0到100 Å），先出现一个最高峰，随后出现一个倒峰，然后趋于和伴随阴离子相似的分布情况，但Na⁺没有倒峰出现。Ca²⁺、Mg²⁺、Cu²⁺、Cd²⁺的倒峰值大小并不相同，这与离子半径和表面电荷密度等特性有关。与阳离子分布不同，阴离子主要分布在质膜的Stern层以外，由于静电斥力的相互作用其分布密度与质膜表面官能团没有相关性，不同阳离子体系中的Cl^-渗透到质膜的能力也不同，表现为Mg²⁺ > Ca²⁺ > Na⁺ > K^+[19]，图 4中显示，在Ca²⁺和Mg²⁺等二价阳离子体系中，Cl^-以与阳离子桥键结合的形式进入质膜区域，而K⁺和Na⁺体系则没有出现阴离子伴随进入质膜的情况。

2.3 磷脂表面离子浓度分布与GCS计算比较

MD模拟可以准确预测不同体系表面电势下的金属离子浓度随距质膜双分子层中心距离变化的分布情况，其结果和经典的GCS理论模型计算结果具有一致性（图 5）。如图 5所示，通过GCS模型中的公式（1）计算得到的膜表面离子浓度分布规律为：与质膜表面带电性相反的金属离子（带正电），在质膜表面的浓度整体呈随距离减小而增大的变化趋势（图 5黑线部分），这是由质膜（负电）与金属离子之间的静电作用引起的；而体系中的阴离子Cl^-在质膜表面的浓度整体呈随着距离的减小而减小的变化趋势（图 5红线部分）。这一变化趋势与MD的模拟结果一致（图 5浅蓝线与黑线重合，深蓝线与红线重合）。MD模拟结果表明，除KCl体系外（图 5a），Stern层内金属离子由于静电吸引力的作用浓度明显比Stern层外大，体系中Cl^-在Stern层内的浓度则相反。此外，CaCl₂、MgCl₂、CuCl₂和CdCl₂等二价金属离子体系中（图 5c至图 5e中浅蓝线），MD计算结果还表明由Stern层外向Stern层内靠近时，出现阳离子浓度先降低后显著增加的过程，这是由于二价金属离子与质膜表面负电荷的官能团吸附强烈，吸附的正电荷超过了膜表面带的负电荷量，使膜表面电势由负变正，所以Stern层外优先静电吸引Cl^-，排斥二价阳离子。通过比较MD模拟和GCS模型计算结果可以发现，MD模拟不但很好地吻合了GCS理论计算得到的Stern层外随着Stern距离增加阴、阳离子浓度分布衰减趋势明显的变化特点，同时能够预测离子在Stern层中的结合和分布情况。

但GCS理论计算离子在膜表面的电荷浓度公式是基于以下假设得到的：（1）体系中的离子为点电荷，没有考虑实际的原子体积和半径^[33]；（2）GCS理论假设背景溶液为稀溶液，界面上水密度均为1 g·cm^-3，不考虑界面相互作用的影响^[34]。MD的模拟结果显示实际平衡后的体系盒子中（图 2b）表面吸附离子浓度较高时，离子会占据相当的体积，且不同离子在体系中占据的体积大小并不相同，不能简单视为点电荷；从图 5中可以看出，体系中的水分子由于与磷脂相互间的作用，在界面表面呈现出随着距双分子层中心距离的增加浓度呈由低到高的变化趋势，而并不是均匀分布的，水分子在界面上的分布情况会影响到膜表面电势的变化和离子的吸附浓度。MD模拟结果直观而有力地弥补了GCS理论因简化假设而造成的微观描述缺陷，为拓展GCS模型的使用范围和验证MD模拟提供了良好的实验数据。

3 结论

（1）MD模拟从分子水平上揭示了不同金属阳离子与质膜结合的构型，发现质膜表面电势强烈影响着离子的分布，二价阳离子主要结合在质膜头部羰基基团和磷酸基团。

（2）MD模拟可以与GCS模型很好地吻合，并且能进一步准确描述双电层中的Stern层中阴、阳金属离子的分布和结合位点，进而更好地为研究重金属的农作物吸收和毒性提供分子层面的机理解释和直观微观图像。