王 鋒， 焦賢發(fā)

（合肥工業(yè)大學(xué) 數(shù)學(xué)學(xué)院，安徽 合肥 230009）

在神經(jīng)模型中，Hodgkin－Huxley模型（H－H模型）能夠精確地模擬真實神經(jīng)元的放電行為［1］，因此H－H模型對于人們研究神經(jīng)元活動的動力學(xué)機(jī)制具有重要意義。利用H－H模型神經(jīng)元來構(gòu)建神經(jīng)網(wǎng)絡(luò)［2］時，盡管能很好地再現(xiàn)神經(jīng)回路行為，但由于每個神經(jīng)元模型均要用多個非線性微分方程來描述，神經(jīng)回路模型就會非常復(fù)雜，這給模型的數(shù)學(xué)處理和數(shù)值分析帶來很多的不便。在構(gòu)造回路模型時，有時并不需要考慮神經(jīng)元的內(nèi)部機(jī)制，只需考慮神經(jīng)元的輸入－輸出關(guān)系。整合發(fā)放神經(jīng)元模型（integrate－and－fire model，IF模型）是H－H模型的簡化形式，該類模型能夠比較精確地反映單個神經(jīng)元電學(xué)特性和對外部輸入的響應(yīng)。整合發(fā)放神經(jīng)元模型是建立在假定神經(jīng)元的輸入發(fā)生在突觸處以及神經(jīng)元的所有輸入均具有可加性的基礎(chǔ)之上，它把神經(jīng)元的活動分成閾下電位的整合和達(dá)到閾值時的發(fā)放及復(fù)位2個階段。在理論神經(jīng)模型研究方面，如果能夠建立輸入－輸出之間的關(guān)系，就能在給定外部輸入情況下對輸出作出準(zhǔn)確的預(yù)測。

為了建立神經(jīng)元輸入－輸出關(guān)系，確定神經(jīng)模型的內(nèi)部參數(shù)和輸入?yún)?shù)是至關(guān)重要的。另外，模型有效性可以通過比較仿真模擬得到的估計值與真實數(shù)據(jù)來定量檢測。從統(tǒng)計學(xué)角度來看，定量檢測首先是對模型中所涉及的參數(shù)進(jìn)行有效的估計，因此參數(shù)估計是分析和驗證神經(jīng)元模型的一個重要方面。

模型中的參數(shù)［3］可以分成以下2類：

（1）描述神經(jīng)元本身電化學(xué)特性的參數(shù)（如膜時間參數(shù)、離子通道電導(dǎo)），這類參數(shù)通常與細(xì)胞的某些內(nèi)在特性有關(guān)，有具體的生物學(xué)解釋。因此，該類參數(shù)對特定的細(xì)胞來說，值是固定不變的，同時可以利用一些電生理實驗方法測定。

（2）描述神經(jīng)元輸入信號特性的參數(shù)，稱之為輸入?yún)?shù)。輸入?yún)?shù)是可變的，只能通過假定神經(jīng)元膜電位服從給定的模型進(jìn)行估計，因此輸入?yún)?shù)是模型中最重要的參數(shù)。

實驗數(shù)據(jù)從本質(zhì)上分成細(xì)胞內(nèi)記錄和細(xì)胞外記錄［4－6］。對于細(xì)胞內(nèi)記錄數(shù)據(jù)，文獻(xiàn)［7］基于首次放電時間（ISI）分布推導(dǎo)出了帶有抑制性反轉(zhuǎn)電位的OU模型（Feller模型）中參數(shù)估計方法，并提出了在閾下情況和閾上情況下的輸入?yún)?shù)矩估計。首次放電時間分布是經(jīng)驗公式，所以利用ISI數(shù)據(jù)得到的估計方法必然會帶來一定的誤差。文獻(xiàn)［8］提出了利用首次放電時間的密度和轉(zhuǎn)移概率密度來估計輸入?yún)?shù)。該方法計算簡便、運算速度快，并且對小數(shù)據(jù)集十分有效。對于細(xì)胞外記錄數(shù)據(jù)，文獻(xiàn)［9］對未知輸入的OU模型中的輸入?yún)?shù)進(jìn)行了估計，并研究了閾下機(jī)制下的參數(shù)估計，所以在估計過程中不可避免地存在系統(tǒng)誤差。文獻(xiàn)［10－13］對閾下機(jī)制和閾上機(jī)制條件下的3種不同的數(shù)學(xué)模型中輸入?yún)?shù)進(jìn)行估計。為了減小經(jīng)驗公式和神經(jīng)元放電特性帶來的系統(tǒng)誤差，本文盡量避免或減少運用經(jīng)驗公式，并且考慮神經(jīng)元放電閾值特性。

整合發(fā)放神經(jīng)元模型不僅能簡單地描述神經(jīng)元的電生理特性，同時也能夠很好地描述神經(jīng)元動作電位的發(fā)放，并且對計算要求相對較低，因此比其他模型更適合用于數(shù)學(xué)分析處理和數(shù)值模擬。然而，在神經(jīng)系統(tǒng)中噪聲是無處不在的，如熱噪聲、突觸噪聲等［3］，因此本文考慮突觸輸入和噪聲共同作用情況下的整合發(fā)放神經(jīng)元模型突觸輸入?yún)?shù)估計。

1 數(shù)學(xué)模型

整合發(fā)放神經(jīng)元模型把神經(jīng)元的活動分成閾下電位的整合和達(dá)到閾值時的發(fā)放動作電位及膜電位復(fù)位2個階段。模型中假定所有神經(jīng)元均為高爾基Ⅰ型神經(jīng)元，即神經(jīng)元輸入只有突觸輸入。對于第1階段，可以用微分方程來描述［1］，該微分方程如下：

其中，Vi為第i個神經(jīng)元在始段處的跨膜電位；τ為膜時間常數(shù)；Ii（t）為第i個神經(jīng)元接受突觸前神經(jīng)元的突觸輸入總和；R為細(xì)胞膜電導(dǎo)；V0為靜息電位。

然而，在神經(jīng)系統(tǒng)中噪聲是無處不在的，如熱噪聲、突觸噪聲等。因此本文考慮突觸輸入和噪聲共同作用情況下的整合發(fā)放神經(jīng)元模型。令Vt＝Vi－V0，則由（1）式可得：

其中，Vt為神經(jīng)元在t時刻的膜電位；τ為膜時間常數(shù)；μ為突觸輸入總和；σ為噪聲強度；Wt為標(biāo)準(zhǔn)維納過程。

2 突觸輸入?yún)?shù)估計

2.1 不考慮放電閾值的突觸輸入?yún)?shù)估計

不考慮閾值行為，神經(jīng)元的膜電位可以看成是由（2）式所確定的連續(xù)變量。對于一個固定時間t，由（2）式確定的Vt為一個高斯隨機(jī)變量，其均值和方差［11，14］為：

本文所給出的所有估計方法都是利用細(xì)胞內(nèi)記錄數(shù)據(jù)（離散的膜電位數(shù)據(jù)）進(jìn)行估計的，是在離散的時間ti＝ih時刻記錄的膜電位。由于膜電位是一個連續(xù)變量，將離散化的膜電位數(shù)據(jù)與其均值和方差分別做差，并求出最小平方差，即最小二乘法。構(gòu)造的最小函數(shù)為：

對（4）式關(guān)于μ求導(dǎo)，為了使（4）式有最小值，則有：

通過數(shù)學(xué)推導(dǎo)，可得μ的估計為：

同理，為了得到σ2的估計，用（6）式的代替（4）式中的μ，構(gòu)造另一個最小函數(shù)，即

對（7）式關(guān)于σ2求導(dǎo)，并令其值等于0，可得σ2的估計為：

采用四階－五階龍格庫塔算法模擬膜電位演化過程時取S＝20mV。當(dāng)μτ＝S時，稱為閾值；當(dāng)μτ＞S時稱為閾上；當(dāng)μτ＜S時稱為閾下。3種參數(shù)不同取值情況下的膜電位的軌跡如圖1所示。

圖1 不同刺激模式下膜電位演化圖

為了避免偶然性，取3組不同的參數(shù)值。噪聲強度σ取0.3mV／ms時，突觸輸入μ為0.8、1.0、1.4mV／ms，膜時間常數(shù)τ取20ms，初始膜電位V0取10mV，時間步長h取0.01ms。對每組參數(shù)值重復(fù)模擬100次得到100組觀測值，利用上述理論推導(dǎo)和Matlab軟件可以獲得100對），所得估計結(jié)果見表1所列。

由表1可以看出，當(dāng)μ為0.8mV／ms時，突觸輸入的均值為0.806 8，均方誤差為0.041 2，表明在突觸輸入μτ＜S時，閾值情況下突觸輸入的估計效果好；當(dāng)μ為1.4mV／ms時，突觸輸入μ的均值為1.485 8，均方誤差為0.382 9，均方誤差較大，估計效果不佳。

表1 不考慮放電閾值的數(shù)值模擬結(jié)果

因此，當(dāng)突觸輸入μτ＞S時，閾上情況下突觸輸入的估計效果不是十分理想。所以最小二乘估計在閾下情況的估計效果要優(yōu)于閾上情況的估計效果。這主要是由于在估計過程中沒有考慮神經(jīng)元放電閾值行為。事實上，神經(jīng)元可看作是存在放電閾值特性的可激發(fā)系統(tǒng)。當(dāng)外部刺激強度在閾值以下時稱為閾下，神經(jīng)元軸突膜上只產(chǎn)生被動的去極化電位，而不產(chǎn)生動作電位［13－14］，膜電位會很快衰減到靜息電位。而當(dāng)刺激強度高于閾值時稱為閾上，神經(jīng)元將會發(fā)放一個動作電位并通過軸突傳導(dǎo)，從而完成從接受信息、處理信息到傳遞信息的過程［15－16］。所以，最小二乘估計依賴于機(jī)制，它僅僅適合閾下情況參數(shù)估計，而對閾上情況估計效果不佳。

2.2 考慮放電閾值的突觸輸入?yún)?shù)估計

如果膜電位的轉(zhuǎn)移概率密度已知并且容易處理，則可以通過極大似然法去估計輸入?yún)?shù)。這樣參數(shù)估計問題就轉(zhuǎn)化為求閾上機(jī)制下膜電位的轉(zhuǎn)移概率密度。為了求解出閾上機(jī)制下的轉(zhuǎn)移概率密度，首先需要求出閾下機(jī)制下的膜電位概率密度［17］。

對應(yīng)于方程（2）的 Fokker－Planck方程［18－19］如下：

其中，f（V，t｜V′，t′）為膜電位由V′變到V時的概率，稱為轉(zhuǎn)移概率密度；D1（V）＝－V／τ；D2（V）＝σ2／2。

為得到微分方程（2）的解析解或數(shù)值解，考慮如下初值條件：

通過數(shù)學(xué)推導(dǎo)可得到閾下機(jī)制下的膜電位轉(zhuǎn)移概率密度為：

其中，N＝ － ［（V－V0e－（t－t′）／τ－μτ＋ （μτ－V′）e－（t－t′）／τ）2］／［σ2τ（1－e－2（t－t′）／τ）］。

為了得到閾上機(jī)制下的膜電位轉(zhuǎn)移概率密度，假定放電閾值為常數(shù)Vth，將閾值條件看成是吸收邊界。考慮帶有吸收邊界條件下的Fokker－Planck方程為：

構(gòu)造似然函數(shù)如下：其中，f（Vi，h｜Vi－1）為轉(zhuǎn)移概率密度函數(shù)；離散膜電位Vi為在離散時間點ih處的膜電位值；h＝t－t′。

將方程（11）式和（12）式代入方程（13）式得：

通過數(shù)學(xué)推導(dǎo)，可得μ和σ2的估計為：

為了研究不同機(jī)制下的估計效果，噪聲強度σ取0.3mV／ms，突觸輸入μ分別取0.8、1.0、1.4mV／ms，膜時間常數(shù)τ取20ms，初始膜電位V0取10mV，時間步長h取0.01ms，閾值Vth取20mV。對每組參數(shù)值重復(fù)模擬100次得到100組觀測值，利用上述理論推導(dǎo)和Matlab軟件可以獲得100對），所得估計結(jié)果見表2所列。

表2 考慮閾值的數(shù)值模擬結(jié)果

由表2可以看出，當(dāng)μ＝0.8mV／ms時，其均值和均方誤差分別為0.790 2和0.040 1；當(dāng)μ＝1.4mV／ms時，其均值和均方誤差分別為1.398 0和0.016 6。由此可以看出極大似然估計的效果非常好。因此極大似然估計使用范圍更廣，而最小二乘估計只使用于閾下活動。為了更好地比較本文給出的最小二乘估計和極大似然估計的優(yōu)劣，依據(jù)膜電位演化方程真實值和估計值得到膜電位演化圖如圖2所示。

圖2 膜電位隨時間演化圖

從圖2a可看出，極大似然估計擬合曲線與真實曲線的匹配程度要高于最小二乘估計擬合曲線與真實曲線的匹配程度，由此可以看出，極大似然估計要優(yōu)于最小二乘估計。由圖2b可看出，極大似然估計擬合曲線與真實曲線相近，基本上將極大似然估計擬合曲線平移一小段距離就可與真實曲線重合。這是因為在估計過程中忽視了不應(yīng)期的存在。事實上，神經(jīng)元發(fā)放一個動作電位后將進(jìn)入一個短暫的絕對不應(yīng)期，在不應(yīng)期內(nèi)無論刺激多么大，神經(jīng)元都不再放電。

3 結(jié)束語

本文綜合考慮突出輸入和噪聲共同作用下的整合發(fā)放神經(jīng)元模型，該研究工作是基于模型中內(nèi)在參數(shù)已知，尤其是膜時間常數(shù)，并且輸入?yún)?shù)是由給定的微分方程（2）所確定。數(shù)值分析結(jié)果表明，對于2個輸入?yún)?shù)μ和σ均可以通過本文所給出的最小二乘估計和極大似然估計進(jìn)行精確地估計。由本文分析可知，在閾下機(jī)制下對于輸入?yún)?shù)μ的估計效果，極大似然估計要優(yōu)于最小二乘估計；而對于輸入?yún)?shù)σ的估計效果，極大似然估計要遠(yuǎn)遠(yuǎn)優(yōu)于最小二乘估計。在其他機(jī)制下，對于2個參數(shù)的估計，極大似然估計都要遠(yuǎn)遠(yuǎn)優(yōu)于最小二乘估計。

本文針對突觸輸入和噪聲共同作用下的整合發(fā)放神經(jīng)元模型，在不考慮放電閾值前提下，采用最小二乘法估計突觸輸入?yún)?shù)；當(dāng)考慮神經(jīng)元放電閾值特性，將放電閾值看成一個吸收邊界，導(dǎo)出膜電位轉(zhuǎn)移概率密度函數(shù)，然后利用極大似然法估計突觸輸入?yún)?shù)。最小二乘估計依賴于機(jī)制，它僅適合閾下的參數(shù)估計，而對閾上無效。極大似然估計適用于所有情況。無論是從適用范圍還是估計精度來說，極大似然估計都要優(yōu)于最小二乘估計。

［1］ 顧凡及，梁培基.神經(jīng)信息處理［M］.北京：北京工業(yè)大學(xué)出版社，2009：129－138.

［2］ 丁亞明，王樹忠，張志紅，等.基于改進(jìn)神經(jīng)網(wǎng)絡(luò)的模糊聚類算法 ［J］.合 肥 工 業(yè) 大 學(xué) 學(xué) 報：自 然 科 學(xué) 版，2007，30（8）：934－938.

［3］ Lansky P，Sanda P，He J.The parameters of the stochastic leaky integrate－and－fire neuronal model［J］.J Comput Neurosci，2006，21：211－223.

［4］ Paninski L，Pillow J，Simoncelli E.Comparing integrate－andfire－like models given intracellular and extracellular data［J］.Neurocomputing，2005，65：379－385.

［5］ H¨opfner R.On a set of data for the membrane potential in a neuron［J］.Math Biosci，2007，207：275－301.

［6］ Mullowney P，Iyengar S.Parameter estimation for a leaky integrate－and－fire neuronal model from ISI data［J］.J Comput Neurosci，2008，24：179－194.

［7］ Ditlevsen S，Lansky P.Estimation of the input parameters in the Feller neuronal model［J］.Phys Rev E，2006，

73：061910.

［8］ Ditlevsen S，Lansky P.Parameters of stochastic processes estimated from observation of first－h(huán)itting times：application to the leaky integrate－fire neuronal model［J］.Phys Rev E，2007，76：041906.

［9］ Ditlevsen S，Ditlevsen O.Parameter estimation from observations of first－passage times of the Ornstein－Uhlenbeck process and the Feller process［J］.Probabilistic Engineering Mechanics，2008，23：170－179.

［10］ Ditlevsen S，Lansky P.Estimation of the input parameters in the Ornstein－Uhlenbeck neuronal model［J］.Phys Rev E，2005，71：011907.

［11］ Lansky P，Ditlevsen S.A review of the methods for signal estimation in stochastic diffusion leaky integrate－and－fire neuronal models［J］.Biol Cybern，2008，99：253－262.

［12］ Bibbona E，Sacerdote L.Errors in estimation of the input signal for integrate－and－fire neuronal models［J］.Phys Rev E，2008，78：011918.

［13］ Joseph H，Guckenheimer J.Parameter estimation for bursting neural models［J］.J Comput Neurosci，2008，24：358－373.

［14］ Bibbona E，Lansky P，Sirovich P.Estimating input parameters from intracellular recordings in the feller neuronal model［J］.Phys Rev E，2010，81：031916.

［15］ Pawlas Z，Lansky P.Distribution of interspike intervals estimated from multiple spike trains observed in a short time window［J］.Phys Rev E，2011，83：011910.

［16］ Ditlevsen S，Lansky P.Only through perturbation can relaxation times be estimated［J］.Phys Rev E，2012，86：050102.

［17］ Kleinhans D.Estimation of drift and diffusion functions from time series data：a maximum likelihood framework［J］.Phys Rev E，2012，85：026705.

［18］ Haken H.The fokker－planck equation［M］.Berlin：Springer－Verlag，1989：100－105.

［19］ 朱位秋.非線性隨機(jī)動力學(xué)與控制［M］.北京：科學(xué)出版

社，2003：56－65.