金屬離子-質膜相互作用的分子動力學模擬研究

徐 冰，劉 存，周東美，王玉軍

（1.中國科學院土壤環境與污染修復重點實驗室（南京土壤研究所），南京 210008；2.中國科學院大學，北京 100049）

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模型可以用于量化 Cu2+、Cd2+、Ca2+、Mg2+等陽離子和砷酸根等陰離子對質膜膜表面電勢（Ψ0）的影響和實際膜表面電荷分布情況。經典GCS模型能夠成功地揭示膜表面電勢影響下不同離子競爭吸附的本質規律，即細胞膜表面分布著負電和電中性等不同結合位點，不同離子在這些位點上有著不同的結合常數，離子在這些結合位點的競爭吸附改變了膜表面電勢，并控制離子在擴散層中的分布。但這一概念化數學模型為了簡化運算將細胞質膜視為剛性平板電容，而實際的細胞質膜是結構組成復雜的兩性雙分子層，表面官能團較多，具有流動性等動態特性，因此很難通過實驗手段和這一經典模型獲取分子水平上重金屬離子與質膜相互作用的微觀動態和機制。分子動力學（Mo?lecular 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的 Ca2+、Mg2+、Cd2+、Cu2+等二價陽離子和53 790個水分子組成：磷脂雙分子層由114個POPC（1-palmitoyl-2-oleoyl-sn-glycero-3-phospho?choline）和 14個 POPG（1-palmitoyl-2-oleoyl-snglycero-3-phosphoglycerol）組成[23-25]（圖2）。模擬體系設定的周期邊界條件PBC（Periodic boundary condi?tions）：盒子大小為64 ?×64 ?×170 ?，溶液相中另外補充適當的Cl-，以保持盒子整體的電中性。所有模擬計算由美國Sandia國家實驗室開發的開源性大規模原子分子并行模擬器LAMMPS（Large-scale atomic/molecular massively parallel simulator）軟件完成，計算體系中原子間相互作用力的力場選用CHARMM 36[26-27]，平衡后體系中磷脂雙分子層表面可變電荷密度（σ0）約為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 磷脂雙分子層結構組分POPC和POPG分子式Figure 1 Phospholipid bilayer structure components POPC and POPG formula

1.4 GCS模型計算磷脂表面離子濃度分布

設定質膜雙分子層中心為Z=0，GCS理論計算溶液中陽離子Mn+和Cl-在體系中的密度分布公式為[18]：

式中：nb表示溶液中金屬離子的濃度；κZ=（8πnbe2/εkBT）1/2為Debye-Hückel參數，其中ε表示水的介電常數，kB為玻爾茲曼常數，T為溫度；Φs表示Z=0處的膜表面靜電勢，由式（2）計算得出：

2 結果與討論

2.1 質膜表面的離子吸附構型

圖2 植物根表質膜與土壤溶液界面的GCS宏觀數學模型（a,示意圖修改自文獻[28]）和分子動力學模型（b）Figure 2 GCS model（a,Adapted from[28]）and molecular dynamics（b）for metal ion-lipid bilayer interactions

圖 3分別為不同溶液體系中Na+、K+、Ca2+、Mg2+、Cd2+、Cu2+等陽離子及陰離子Cl-在質膜雙分子層表面的吸附模型。模擬結果表明，K+在質膜表面吸附較弱（圖3b），而其余陽離子都傾向于與質膜表面位點強烈結合[如羰基（C=O）與磷酸基（-OP）]。這是因為與Na+、Ca2+、Mg2+、Cu2+、Cd2+相比，K+的離子半徑更大，離子表面的電荷密度更小，與周圍水分子的結合能力更弱，水合半徑更小，進而KCl體系的離子強度也比其他體系低，而離子強度的大小又是影響細胞質膜表面電勢絕對值的關鍵。因此K+吸附并進入細胞質膜的自由能壘低于Na+以及二價金屬陽離子[29]。

除K+外，其他陽離子在結合到質膜表面時，傾向于與質膜頭部的極性基團羥基、磷酸氧和羰基結合（圖3內插圖中內容所示），主要是這些官能團的負電荷與陽離子發生靜電絡合作用。金屬陽離子與質膜頭部基團的結合導致了一個與磷脂雙分子層自身方向相反偶極矩的生成。質膜頭部帶有正電荷的膽堿基團，也因為與陽離子之間的排斥作用以及與陰離子的吸引作用發生一定的構型改變，傾向于轉向水相或者與磷酸基團呈相反的方向。金屬陽離子吸附導致了分子構型變化進而促使磷脂面積減小和質膜厚度增加，促使質膜更加緊密的包裹[19，30-31]。盡管自然狀態下根細胞質膜呈負電性，通過靜電排斥作用導致陰離子Cl-通常分布在溶液中，隨著距離質膜表面越遠，陽離子的濃度越低，陰離子Cl-的濃度越高。但圖3c的Ca2+體系吸附構型還顯示出在Ca2+與質膜表面的官能團結合時有Cl-伴隨，證明了溶液相的陰離子可以通過陽離子橋鍵結合的形式吸附到質膜表面。

圖3 金屬離子在質膜表面的吸附分布（內插圖放大顯示了不同離子與質膜官能團的典型結合方式）Figure 3 The adsorption configuration of Na+，K+，Ca2+，Mg2+，Cu2+and Cd2+on the surface of the plasma membrane，respectively，and the insets enlarged the coordination environment of cations

2.2 磷脂表面密度分布

不同體系中陽離子、陰離子、水分子和質膜雙分子頭部的羰基氧、羥基氧、磷酸氧和氮原子在質膜表面密度分布情況如圖4表示。原子的位置對脂質雙層結構至關重要。分布峰之間的距離可以直接反映質膜極性基團構型的變化。以質膜中心作為Z軸的原點（Z=0），圖中顯示的為質膜雙分子層的一半，密度數值顯示的是質膜兩側的分布均值。在膜-水界面的周圍，K+的密度分布沒有呈現出明顯的峰值（圖4b），這是因為K+半徑較大，電荷密度低，與負電基團親和性較弱[19，32]。Na+及二價金屬離子體系中（圖4a和圖4c至圖4f深藍線），陽離子在靠近質膜表面的位置出現一個分布峰，其峰值出現的位置和位于脂質的羰基（Mn+-Ocarb）和磷酸氧（Mn+-P）最大峰處的分布具有高度的相關性，羥基基團（Mn+-OH）密度分布的最大值與金屬離子概率密度的峰值重疊部分較少，說明羥基基團對吸附 Ca2+、Mg2+、Cd2+、Cu2+陽離子的貢獻比羰基和磷酸基團少，這一結論也驗證了不同金屬離子在質膜表面吸附構型的不同。同時我們觀察到氮原子（Mn+-N）與磷原子（Mn+-P）的密度分布相似，并且受陽離子結合的影響方式相似。這種構型和電荷分布的相關變化是由于膽堿氮和陽離子之間的相互排斥作用以及陰離子和膽堿基團之間的吸引力引起的[19]。

圖4 不同體系中的金屬離子K+、Na+、Ca2+、Mg2+、Cu2+、Cd2+和Cl-（-Cl）、水分子（-H2O）、羰基氧（-Ocarb）、磷酸氧（-P）、質膜頭部氮原子（-N）、羥基氧（-OH）等在質膜表面的離子濃度分布Figure4 TheionconcentrationdistributiononthesurfaceoftheplasmamembraneofK+，Na+、Ca2+、Mg2+、Cu2+、Cd2+andCl-（-Cl），water molecule（-H2O），carbonyloxygen（-Ocarb），phosphoricacidoxygen（-P），plasmamembraneheadnitrogenatom（-N）andhydroxyoxygen（-OH）

除 K+以外，Na+和二價金屬離子 Ca2+、Mg2+、Cu2+、Cd2+都在水分子（Mn+-H2O）達到最大分布密度峰附近出現一個倒峰，水分子分布出現最大值的位置為質膜吸附金屬陽離子緊密層（雙電層中的Stern層）的外側，二價金屬離子在從質膜中心到體相溶液的分布過程中（距質膜雙分子層中心距離由0到100 ?），先出現一個最高峰，隨后出現一個倒峰，然后趨于和伴隨陰離子相似的分布情況，但Na+沒有倒峰出現。Ca2+、Mg2+、Cu2+、Cd2+的倒峰值大小并不相同，這與離子半徑和表面電荷密度等特性有關。與陽離子分布不同，陰離子主要分布在質膜的Stern層以外，由于靜電斥力的相互作用其分布密度與質膜表面官能團沒有相關性，不同陽離子體系中的Cl-滲透到質膜的能力也不同，表現為Mg2+＞Ca2+＞Na+＞K+[19]，圖4中顯示，在Ca2+和Mg2+等二價陽離子體系中，Cl-以與陽離子橋鍵結合的形式進入質膜區域，而K+和Na+體系則沒有出現陰離子伴隨進入質膜的情況。

2.3 磷脂表面離子濃度分布與GCS計算比較

圖5 不同體系中金屬陽離子和陰離子在質膜表面的離子濃度分布和GCS計算比較Figure 5 The surface ion concentration distribution of metal cations and anion plasma membranes in different systems and compared with GCS model calculation

MD模擬可以準確預測不同體系表面電勢下的金屬離子濃度隨距質膜雙分子層中心距離變化的分布情況，其結果和經典的GCS理論模型計算結果具有一致性（圖5）。如圖5所示，通過GCS模型中的公式（1）計算得到的膜表面離子濃度分布規律為：與質膜表面帶電性相反的金屬離子（帶正電），在質膜表面的濃度整體呈隨距離減小而增大的變化趨勢（圖5黑線部分），這是由質膜（負電）與金屬離子之間的靜電作用引起的；而體系中的陰離子Cl-在質膜表面的濃度整體呈隨著距離的減小而減小的變化趨勢（圖5紅線部分）。這一變化趨勢與MD的模擬結果一致（圖5淺藍線與黑線重合，深藍線與紅線重合）。MD模擬結果表明，除KCl體系外（圖5a），Stern層內金屬離子由于靜電吸引力的作用濃度明顯比Stern層外大，體系中Cl-在Stern層內的濃度則相反。此外，CaCl2、Mg?Cl2、CuCl2和CdCl2等二價金屬離子體系中（圖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層中陰、陽金屬離子的分布和結合位點，進而更好地為研究重金屬的農作物吸收和毒性提供分子層面的機理解釋和直觀微觀圖像。