孟盼陸啟韶趙勇董健衛(wèi)

（1.廣東藥科大學(xué)基礎(chǔ)學(xué)院，廣州 510006）（2.北京航空航天大學(xué)航空科學(xué)與工程學(xué)院，北京 100191）（3.河南理工大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院，焦作 454000）

兩房室神經(jīng)元模型的同步簇放電的動(dòng)力學(xué)分析*

孟盼1?陸啟韶2趙勇3董健衛(wèi)1

（1.廣東藥科大學(xué)基礎(chǔ)學(xué)院，廣州 510006）（2.北京航空航天大學(xué)航空科學(xué)與工程學(xué)院，北京 100191）（3.河南理工大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院，焦作 454000）

研究一個(gè)具有電流反饋的兩房室錐體神經(jīng)元的簇放電動(dòng)力學(xué)問(wèn)題.首先，利用數(shù)值模擬研究了該神經(jīng)元的樹(shù)突和胞體兩個(gè)房室之間的峰相位同步行為.其次，基于快慢動(dòng)力學(xué)分析方法，給出了快子系統(tǒng)的雙參數(shù)分岔曲線圖，結(jié)果表明快子系統(tǒng)產(chǎn)生了3個(gè)余維-2分岔點(diǎn)，它們分別為fold-Hopf分岔點(diǎn)ZH、Cusp分岔點(diǎn)CP和Bogdanov-Takens分岔點(diǎn)BT.最后，根據(jù)這些余維-2分岔點(diǎn)確定系統(tǒng)存在著豐富的簇放電模式，如“fold/fold”型，“Hopf/Hopf”型，“Hopf/homoclinic”型以“fold/homoclinic”型簇振蕩等，并分析研究了簇放電模式的產(chǎn)生以及它們之間相互轉(zhuǎn)遷的動(dòng)力學(xué)機(jī)理.本文的研究結(jié)果進(jìn)一步表明了神經(jīng)元的幾何形態(tài)特性對(duì)于神經(jīng)元放電模式多樣性的重要影響.

簇放電， 同步， 電流反饋， 快慢動(dòng)力學(xué)分析， 分岔

引言

神經(jīng)元作為神經(jīng)系統(tǒng)的基本功能單位，通過(guò)放電活動(dòng)感受外界刺激，從而實(shí)現(xiàn)神經(jīng)信息的傳遞［1］.不同的放電模式是對(duì)不同外界刺激的不同反映［2-3］.研究表明，簇放電（bursting）是神經(jīng)元最重要的放電模式之一，在生物信息編碼中起著關(guān)鍵作用.簇振蕩比單個(gè)峰放電（spiking）可以攜帶更多的信息并且更加可靠，因此理解簇振蕩活動(dòng)的產(chǎn)生和轉(zhuǎn)遷機(jī)制具有重要意義［4］.簇振蕩表現(xiàn)為慢變量引起的靜息態(tài)（silent phase）與快變量引起的重復(fù)放電狀態(tài)（active phase）之間的相互轉(zhuǎn)遷.利用快慢動(dòng)力學(xué)分析來(lái)考察簇放電的靜息態(tài)和放電狀態(tài)之間的分岔，并依此給出簇模式的拓?fù)漕?lèi)型，是理解簇振蕩產(chǎn)生機(jī)理的主要方法［5］.

為了更好的探究神經(jīng)元的各部分在信息傳遞中的作用，神經(jīng)元的形態(tài)結(jié)構(gòu)多利用房室模型來(lái)刻畫(huà)，如錐體神經(jīng)元［6］，蒲肯野神經(jīng)元［7］和中間神經(jīng)元［8］等.房室模型是將空間上連續(xù)的單個(gè)神經(jīng)元分成離散的多個(gè)房室，各個(gè)房室相互連接后整體描述神經(jīng)元的時(shí)空變化.和單房室神經(jīng)元模型相比，多房室模型同樣可以產(chǎn)生豐富的放電模式，并且能夠更好的體現(xiàn)神經(jīng)元的形態(tài)結(jié)構(gòu)特征［6-8］.利用房室模型來(lái)研究神經(jīng)元的電活動(dòng)行為，一直是計(jì)算神經(jīng)領(lǐng)域的熱點(diǎn).本文主要研究具有電流反饋的兩房室錐體神經(jīng)元模型，通過(guò)數(shù)值模擬和理論分析，討論了神經(jīng)元通過(guò)形態(tài)參數(shù)變化時(shí)展現(xiàn)的放電節(jié)律模式以及各種模式之間的相互轉(zhuǎn)遷機(jī)理.

1 模型描述

考慮錐體神經(jīng)元，它由樹(shù)突和軸突兩個(gè)房室構(gòu)成.模型描述如下［9］：

其中門(mén)控函數(shù)為：

（1）式和（3）式分別為胞體膜電壓VS和樹(shù)突膜電壓VD的動(dòng)力學(xué)方程，其中胞體房室的電流包含了胞體的突觸輸入電流、樹(shù)突流向胞體的內(nèi)部電流、Na+離子電流、K+離子電流和泄漏電流；樹(shù)突房室的電流包含了樹(shù)突的突觸輸入電流、樹(shù)突流向胞體的內(nèi)部電流和泄漏電流.ENa和EK分別表示Na+離子通道和K+離子通道的平衡電位，ESL和EDL分別為胞體和樹(shù)突膜上泄漏電流的平衡電位；gNa、gK、gSL和gDL分別代表各通道的最大電導(dǎo)；gc表示兩房室連接電導(dǎo)；C表示膜電容；IS和ID分別為胞體和樹(shù)突的突觸輸入電流；p和1-p分別表示胞體和樹(shù)突所占的面積比例.p是錐體神經(jīng)元的一個(gè)重要幾何形態(tài)參數(shù).（2）式描述了K+離子通道門(mén)開(kāi)通概率w的變化規(guī)律.作為外界電流刺激的影響，（1）式中的胞體突觸輸入電流IS由一個(gè)線性反饋的動(dòng)力學(xué)方程（4）支配，這里ε是時(shí)間尺度因子，取為很小的正數(shù)以保證IS的慢變性.方程（1）-（4）稱(chēng)為具有電流反饋的兩房室神經(jīng)元模型.由于ε?1，因此方程（1）-（3）可視為快子系統(tǒng)，其中取IS作為慢變量，p作為控制參數(shù).在這篇文章中，其它參數(shù)取值為：ε＝0.01，gNa＝50，gK＝15，gSL＝2，gDL＝2，gc＝2，ENa＝10，EK＝-90，ESL＝-70，EDL＝-70，Φ＝0.3，C＝2，v0＝-22，v1＝1.2，v2＝18，v3＝10，v4＝6.

為了便于以后的研究，先考慮在上述參數(shù)條件下兩個(gè)房室的同步動(dòng)力學(xué)行為.引入峰相位函數(shù)φS（t）［10］，定義為：

其中ti是第i次峰放電發(fā)放的時(shí)間.因此每經(jīng)過(guò)一次峰放電，峰相位就增加2π.基于相位函數(shù)的定義，引入頻率函數(shù)ω（t）＝˙φ（t）和平均峰頻率函數(shù)Ω（t）＝＜˙φ（t）＞，這里＜·＞表示對(duì)時(shí)間取平均，因此有

如果兩個(gè)房室的平均峰頻率差ΔΩ＝ΩS-ΩD為0，那么兩個(gè)房室的放電頻率一致（在時(shí)間平均的意義上），從而實(shí)現(xiàn)放電相位的同步［11］.

2 神經(jīng)元形態(tài)參數(shù)引起的兩室間簇放電的同步

通過(guò)數(shù)值計(jì)算得到ΔΩ和形態(tài)參數(shù)p的關(guān)系圖（見(jiàn)圖1）.由圖1可以看出在0.18＜p＜1.0的參數(shù)范圍內(nèi)有平均峰頻率差ΔΩ＝0，即兩個(gè)房室總是處于峰相位同步狀態(tài)，下面的圖3～6的數(shù)值結(jié)果也證實(shí)這一點(diǎn).因此在下面的神經(jīng)元簇放電模式研究中，我們可以著重考察其中一個(gè)房室（如胞體）的放電動(dòng)力學(xué)行為，例如只依據(jù)VS的分岔圖就可以進(jìn)行快慢動(dòng)力學(xué)分析.

圖1 平均峰頻率差ΔΩ關(guān)于p的變化圖Fig.1 Relationship of themean spike frequency differenceΔΩ and parameter p

3 簇放電模式及轉(zhuǎn)遷

現(xiàn)在考慮前述的兩房室錐體神經(jīng)元模型的復(fù)雜放電模式，圖2給出一些分岔結(jié)果（具體的分析過(guò)程從略）.圖2（a）是全系統(tǒng)（1）-（4）關(guān)于形態(tài)參數(shù)p的分岔圖.可以看出，平衡點(diǎn)的分岔曲線是一條直線，在HB處出現(xiàn)supHopf分岔，系統(tǒng)的靜息態(tài)消失，使得穩(wěn)定平衡點(diǎn)變?yōu)椴环€(wěn)定，同時(shí)產(chǎn)生穩(wěn)定極限環(huán)，系統(tǒng)開(kāi)始振蕩.圖2（b）是雙參數(shù)（IS，p）-平面的分岔曲線圖.其中余維-1分岔曲線包括Hopf分岔曲線（h），平衡點(diǎn)的鞍結(jié)分岔曲線（f1，f2），同宿軌分岔曲線（hc）.余維-2分岔點(diǎn)分別為fold-hopf分岔點(diǎn)ZH（Hopf分岔曲線h和鞍結(jié)分岔曲線f2的交點(diǎn)）、Bogdanov-Takens分岔點(diǎn)BT（fold分岔曲線f1和同宿軌分岔曲線hc相切的點(diǎn)）以及Cusp分岔點(diǎn)CP（兩條鞍結(jié)分岔曲線f1和f2處的交點(diǎn)）.左下角是點(diǎn)BT附近的局部放大圖.右上角是點(diǎn)ZH附近的局部放大圖，右下角是點(diǎn)CP附近的放大圖.

圖2 （a）全系統(tǒng)（1）-（4）關(guān)于形態(tài)參數(shù)p的分岔圖.（b）（IS，p）-參數(shù)平面上快子系統(tǒng)（1）-（3）的分岔曲線.Fig.2 （a）One-parameter bifurcation diagram of the whole system with the variation of parameter p.（b）Bifurcation curves in（IS，p）-parametric plane for the fast subsystem（1）-（3）

為了進(jìn)一步說(shuō)明問(wèn)題，下面結(jié)合圖2來(lái)考慮形態(tài)參數(shù)p的變化對(duì)系統(tǒng)（1）-（4）的影響.通過(guò)取參數(shù)p不同的值來(lái)給出整個(gè)系統(tǒng)的某些典型的放電模式.隨著參數(shù)p的變化，系統(tǒng)呈現(xiàn)不同形式的簇放電模式，如圖3～6所示，其中黑線和紅線分別表示胞體和樹(shù)突的放電時(shí)間序列圖.可以看到這兩個(gè)房室的簇放電行為一直處于峰相位同步狀態(tài).接下來(lái)利用快慢動(dòng)力學(xué)分析，解釋這些簇模式之間轉(zhuǎn)遷的動(dòng)力學(xué)機(jī)理.

首先，當(dāng)形態(tài)參數(shù)p＝0.19時(shí)，系統(tǒng)產(chǎn)生簇放電模式，如圖3（a）所示.快慢動(dòng)力學(xué)分析如圖3（a）所示.快子系統(tǒng)的分岔曲線是一條S形曲線，上支由穩(wěn)定焦點(diǎn)（實(shí)線）和不穩(wěn)定焦點(diǎn)（虛線）組成；中支是由鞍點(diǎn)（虛線）組成；下支是由穩(wěn)定結(jié)點(diǎn)（實(shí)線）組成.由圖2（b）可以看出，現(xiàn)在在慢變參數(shù)區(qū)域內(nèi)的快子系統(tǒng)的分岔曲線從左到右僅有平衡點(diǎn)的鞍結(jié)分岔曲線f1和f2，相應(yīng)的快子系統(tǒng)的分岔點(diǎn)分別記為F1和F2.全系統(tǒng)的軌線也疊加在快子系統(tǒng)的分岔圖上.隨著參數(shù)IS的增加，相應(yīng)于快子系統(tǒng)分岔曲線下支穩(wěn)定結(jié)點(diǎn)的簇放電的靜息態(tài)，經(jīng)鞍結(jié)分岔F2轉(zhuǎn)遷到上支的放電狀態(tài).隨著IS的減小，放電狀態(tài)經(jīng)由鞍結(jié)分岔F1轉(zhuǎn)遷到下支的靜息態(tài).因此簇放電類(lèi)型是“fold/fold”型簇放電［5］.

圖3 當(dāng)p＝0.19時(shí)，“fold/fold”型簇放電.（a）兩個(gè)房室的膜電位VS（黑線），VD（紅線）的時(shí)間歷程圖.（b）快慢動(dòng)力學(xué)分岔分析.Fig.3 Bursting of“fold/fold”type when p＝0.19.（a）Time history ofmembrane potentials VSin the black line and VDin the red one.（b）Fast-slow dynamics of the fast subsystem.

圖4 當(dāng)p＝0.35時(shí)，“Hopf/Hopf”型簇放電.（a）兩個(gè)房室的膜電位VS（黑線），VD（紅線）的時(shí)間歷程圖.（b）快慢動(dòng)力學(xué)分岔分析.Fig.4 Bursting of“fold/fold”type when p＝0.35.（a）Time history ofmembrane potentials VSin the black line and VDin the red one.（b）Fast-slow dynamics of the fast subsystem

其次，當(dāng)p＝0.35時(shí)，系統(tǒng)可產(chǎn)生如圖4（a）的簇放電模式.運(yùn)用圖4（b）的快慢動(dòng)力學(xué)分析，來(lái)研究此放電的動(dòng)力學(xué)性質(zhì)及其類(lèi)型.通過(guò)雙參數(shù)（IS，p）-平面分岔圖2（b）可以看出，在兩個(gè)鞍結(jié)分岔點(diǎn)F1和F2之間，產(chǎn)生了兩個(gè)supHopf分岔點(diǎn)，在快子系統(tǒng)的分岔圖上分別記為HB1和HB2.這樣，S形分岔曲線的上支的穩(wěn)定焦點(diǎn)（實(shí)線）經(jīng)由點(diǎn)HB1處的Hopf分岔而失穩(wěn)（虛線），同時(shí)產(chǎn)生穩(wěn)定極限環(huán).隨著參數(shù)IS的增加，該穩(wěn)定極限環(huán)經(jīng)由點(diǎn)HB2處的Hopf分岔消失而轉(zhuǎn)變?yōu)榉€(wěn)定焦點(diǎn).現(xiàn)在結(jié)合全系統(tǒng)的軌線來(lái)討論簇放電的產(chǎn)生過(guò)程.隨著慢變量IS的增加，由于supHopf分岔點(diǎn)HB2處的慢變效應(yīng)［11］，下支的靜息狀態(tài)轉(zhuǎn)遷到上支的放電狀態(tài)；隨著慢變量IS的減小，放電狀態(tài)收斂于經(jīng)由supHopf分岔點(diǎn)HB1產(chǎn)生的穩(wěn)定焦點(diǎn)，再經(jīng)由鞍結(jié)分岔F1轉(zhuǎn)遷到下支的靜息狀態(tài).這種簇振蕩模式稱(chēng)為“Hopf/Hopf”型簇放電［5］.綜上所述，當(dāng)參數(shù)p的取值小于余維-2分岔點(diǎn)BT的縱坐標(biāo)pBT值時(shí)，系統(tǒng)會(huì)產(chǎn)生“fold/fold”型或“Hopf/Hopf”型簇放電.

接下來(lái)，當(dāng)形態(tài)參數(shù)p的取值位于兩個(gè)余維-2分岔點(diǎn)Bogdanov-Takens分岔點(diǎn)BT和fold-hopf分岔點(diǎn)ZH的縱坐標(biāo)值之間，例如p＝0.42時(shí)，兩個(gè)房室的簇放電模式如圖5（a）所示，其相應(yīng)的分岔分析見(jiàn)圖5（b）.和圖4（b）相比，現(xiàn)在快子系統(tǒng)僅有一個(gè)SupHopf分岔點(diǎn)HB.神經(jīng)元的靜息態(tài)經(jīng)由SupHopf分岔而轉(zhuǎn)遷到放電態(tài)，放電態(tài)又經(jīng)極限環(huán)的同宿軌分岔HC轉(zhuǎn)變?yōu)殪o息態(tài)，從而形成了一個(gè)“Hopf/Homoclinic”型簇放電［5］.

圖5 當(dāng)p＝0.42時(shí)，“Hopf/Homoclinic”型簇放電.（a）兩個(gè)房室的膜電位VS（黑線）、VD（紅線）的時(shí)間歷程圖.（b）快慢動(dòng)力學(xué)分岔分析.Fig.5 Bursting of“Hopf/Homoclinic”type when p＝0.42.（a）Time history ofmembrane potentials VSin the black line and VDin the red one，（b）Fast-slow dynamics of the fast subsystem.

圖6 當(dāng)p＝0.9時(shí)，“fold/homoclinic”型簇放電.（a）兩個(gè)房室的膜電位VS（黑線）、VD（紅線）的時(shí)間歷程圖.（b）快慢動(dòng)力學(xué)分岔分析.Fig.6 Bursting of“fold/Homoclinic”type when p＝0.9，（a）Time history ofmembrane potentials VSin the black line and VDin the red one.（b）Fast-slow dynamics of the fast subsystem.

最后，當(dāng)形態(tài)參數(shù)p的取值大于余維-2分岔點(diǎn)fold-hopf分岔點(diǎn)ZH的縱坐標(biāo)pZH值，例如p＝0.9時(shí)，系統(tǒng)可產(chǎn)生如圖6（a）的簇放電模式，相應(yīng)的快慢動(dòng)力學(xué)分析如圖6（b）所示.隨著慢變量的增加，快子系統(tǒng)依次經(jīng)過(guò)鞍結(jié)分岔F1、同宿軌分岔HC、鞍結(jié)分岔F2、supHopf分岔HB.值得注意的是，現(xiàn)在鞍結(jié)分岔點(diǎn)F2和supHopf分岔點(diǎn)HB的先后順序和圖5（b）有所不同.因此，系統(tǒng)軌線的靜息態(tài)由平衡點(diǎn)的鞍結(jié)分岔（F2）轉(zhuǎn)變?yōu)榉烹姞顟B(tài)，然后放電狀態(tài)沿著由HB產(chǎn)生的穩(wěn)定極限環(huán)運(yùn)動(dòng)，最后消失于極限環(huán)的同宿軌分岔（HC）.這種簇振蕩模式稱(chēng)為“fold/homoclinic”型簇放電［5］.

4 結(jié)論

通過(guò)數(shù)值模擬和快慢動(dòng)力學(xué)分析，研究了具有電流反饋的兩房室錐體神經(jīng)元的簇放電模式.結(jié)果表明，胞體和樹(shù)突這兩個(gè)房室在一定條件下可以處于峰相位同步狀態(tài)：神經(jīng)元可以表現(xiàn)出豐富的放電模式，如“fold/fold”型、“Hopf/Hopf”型、“Hopf/Homoclinic”型和“fold/homoclinic”型簇放電模式；特別地，神經(jīng)元的幾何形態(tài)參數(shù)對(duì)神經(jīng)元放電節(jié)律的產(chǎn)生和轉(zhuǎn)遷有著重要的影響.還注意到，如果不考慮電流反饋的話，此兩房室神經(jīng)元模型并不存在這些豐富的簇放電模式，因此本文的研究對(duì)進(jìn)一步揭示外界激勵(lì)對(duì)神經(jīng)元放電活動(dòng)的重要作用也有一定意義.

1 汪云九.神經(jīng)信息學(xué)-神經(jīng)系統(tǒng)的理論和模型.北京，高等教育出版社，2006（Wang Y J.Neuro-informationtics：theories and models of neural systems.Beijing：Higher Education Press，2006（in Chinese））

2 Fan Y S，Chay T R.Generation of periodic and chaotic bursting in an excitable cellmodel.Biological Gybernetics，1994，71（5）：417～431

3 Gu H G，Pan B B.Identification of neural firing patterns，frequency and temporal coding mechanisms in individual aortic baroreceptors.Frontiers in Computational Neuroscience，2015，9：108

4 Lisman JE.Bursts as a unit of neural information：making unreliable synapses reliable.Trends in Neuroscience，1997，20（1）：38～43

5 Rinzel J.Bursting oscillations in an excitable membranemodel.In：Ordinary and Partial Differential Equations，Lecture Notes in Mathematics.edited by Sleeman B D，Jarvis R J，New York：Springer，1995：304～316

6 馬丹，劉深泉，汪雷.兩房室神經(jīng)元模型的分岔現(xiàn)象.北京生物醫(yī)學(xué)工程，2011，30（6）：567～573（Ma d，Liu S Q，Wang L.Bifurcation phenomenon in a two-compartment neuron model.Beijing Biomedical Engineering，2011，30（6）：567～573（in Chinese））

7 張煦晨，劉深泉，任會(huì)霞.浦肯野神經(jīng)元回路中的同步性研究.動(dòng)力學(xué)與控制學(xué)報(bào)，2015，13（5）：377～381（Zhang X C，Liu SQ，Ren H X.The research on synchronization in purinje neural circuit.Journal of Dynamics and Control，2015，13（5）：377～381（in Chinese））

8 江小芳，劉深泉，張煦晨.中等多棘神經(jīng)元的多房室模型分析.生物物理學(xué)報(bào)，2013，29（7）：527～542（Jiang X F，Liu SQ，Zhang X C.Multi-compartmentmodelanalysis of the medium spiny neuron.Biphysics Reports，2013，29（7）：527～542（in Chinese））

9 伊國(guó)勝，王江，魏熙樂(lè)，鄧斌.電場(chǎng)調(diào)制神經(jīng)元放電動(dòng)力學(xué)的生物物理機(jī)制.中國(guó)科學(xué)：物理學(xué)·力學(xué)·天文學(xué)，2013，42（2）：210～216.（Yin G S，Wang J，Wei X L，Deng B.The mechanism of electrical field modulates the neuron firing pattern dynamics.Science China Physics，Mechanics＆Astronomy，2013，42（2）：210～216（in Chinese））

10王海俠，陸啟韶.對(duì)稱(chēng)電耦合連接方式下的多尺度同步行為和轉(zhuǎn)遷過(guò)程.動(dòng)力學(xué)與控制學(xué)報(bào)，2011，9（3）：214～218（Wang H X，Lu Q S.Multi-time-scale synchronization and transitions with symmetrical and electrical coupling schemes.Journal of Dynamics and Control，2011，9（3）：214～218（in Chinese））

11 Baer S M，Erneux T，Rinzel J.Slow passage through a Hopf bifurcation：Delay，memory effects，and resonance. SIAM Journal on Applied Mathematics，1989，49：55～71

DYNAM ICS ANALYSISOF SYNCHRONOUS BURSTING IN TWO-COMPARTMENT NEURON MODEL*

Meng Pan1?Lu Qishao2Zhao Yong3Dong Jianwei1
（1.School of Basic Courses，Guangdong Pharmaceutical University，Guangzhou 510006，China）（2.School ofAeronautical Science and Engineering，Beihang University，Beijing 100083，China）（3.School of Mathematics and Information Science，Henan Polytechnic University，Jiaozuo 454000，China）

A two-compartment neuronmodelwith outside current feedback is considered in order to study the dynamics of bursting patterns.Firstly，on the base of numerical simulation，the behavior of phase synchronization for the spikes in the soma and dendrite compartments is investigated.Secondly，by using fast-slow dynamics analysis，two-parameter bifurcation analysis of the fast subsystem is performed in the two-parameter plane.It is shown that three condimension-2 bifurcation points are generated in the fast subsystem，i.e.，fold-Hopf bifurcation ZH，cusp bifurcation CP and Bogdanov-Takens bifurcation BT.Finally，these condimension-2 bifurcation points determine that different bursting oscillations occurs in the system，such as“fold/fold”type，“Hopf/Hopf”type，“Hopf/homoclinic”type and“fold/homoclinic”type.Furthermore，the generation and transition mechanisms of different bursting patterns are discussed.The results offer significant data for understanding the important role of the geometric parameter in the variety of neuronal firing activities

bursting， synchronization， current-feedback， fast-slow analysis， bifurcation

10.6052/1672-6553-2016-023

2016-01-31收到第1稿，2016-03-26收到修改稿.

*國(guó)家自然科學(xué)基金資助項(xiàng)目（11402057，11502073，11272034，11371046）

?通訊作者E-mail：mengpan200e＠163.com

Received 31 January 2016，revised 26 March 2016.

*The project supported by the National Natural Science Foundation of China（11402057，11502073，11272024 and 11371046）

?Corresponding author E-mail：mengpan200e＠163.com