神經元的交流，傳輸與活動，都離不開一個個非常短暫的脈沖-Spike。有各種各樣的模型，可以描述神經元的電位變化，發放，比如HH模型等等。但是如果只考慮比較粗糙的一些性質，比如膜電位的簡單變化和spike的頻率之類，可以直接利用神經元的膜電位性質（與電容類似），以及其動作電位這種比較模式化(stereotypical)，又比較短暫的性質，吧神經元抽象為一個電容，以及超過特定閾值會產生一個大的spike，且有一定的不應期(Refractory Period)的簡單模型。這樣的模型可以描述基本的神經元對于外界輸入電流 (Injected current, 或者建模的突觸輸入)的相應，也可以引入噪聲模型，觀察神經元發放的穩定性等等。在本文內我會主要討論一下神經元的建模，應對不同輸入電流的響應，噪聲模型，以及對發放穩定性(CV, Coefficient of Variation，變異系數)等特性，以及用Python代碼對其實現。

這學期我將會更新一系列計算神經科學相關的內容，主要是整理一些本學期Computational Neuroscience課程上學到的東西分享給大家。也非常歡迎大家一起討論各種計算相關的問題。

LIF 模型的基本原理

一般情況下，神經元的動作電位主要是靠Na+, K+的選擇性 on/off dynamics實現的。經典的HH模型將各種離子的電位，及其門控通道等效為一個電池和其受調控的可變導鈉的電阻，又吧膜電位的被動特性等效為一個電容。再將以上并聯，從而得到一個對神經元電特性的建模。

Hodgkin and Huxley (1952)

LIF模型忽略了比較brief的動作電位具體的Dynamic，只考慮膜電位的被動特性，以及動作電位的觸發，所以LIF模型將HH模型簡化為電容+一個電池&可變電阻的并聯。同時LIF模型想要建模不同輸入電流的影響（用于模擬膜片鉗實驗的current clamp），因此再并聯了一個外部輸入的電流源I_app，如下圖：

Koch (1999)

離子通道的電壓門控特性靠設定一個threshold來實現，除此之外，LIF模型還考慮了動作電位之后的不應期: 在發放了一個動作電位之后: 神經元會位置在一個reset電位幾毫秒。

LIF模型的公式

LIF模型的數學闡述如下: 主要是通過這個公式來描述膜電位的被動特性以及外部輸入電流的影響。

這個公式每個部分都代表上圖的每個并聯的部分，利用基爾霍夫電流結點定律(流入節點的電流=流出節點電流)。在這里

代表細胞膜表面的電容， ， 分別代表Leak的導納，以及被動平衡的電壓 （一般就是靜息電位，-70mV）。如果沒有外界電流輸入，公式還可以改寫為時間常數的版本（純粹的被動特性）。

在存在一個外界的

輸入時，膜電位的平衡電位(穩態電位)會升高，一旦膜電位超過一個給定的閾值 后，會產生一個Spike (比如30mV)，并且會將膜電位重置到一個重置電位 并維持一個給定的不應期時間 (以上給定的參數僅為參考)，然后再利用LIF的被動電位公式進行下一步的計算。

通過LIF模型，我們可以計算出一個理論的神經元活動的發放間歇長度(ISI)，從而通過其倒數推算其理論的發放速度 (Firing Rate)。理論的ISI是這樣計算的：首先有一個確定的不應期時間

，然后存在一個不應期的電位 ，可以通過對微分方程進行求解，得到膜電位從 上升到 所需要的時間，再加上不應期時間 就可以得到每個ISI的時間，再求倒數，得到具體的Firing Rate。可以算出來的 Firing Rate公式推算如下：

其中，

是能觸發Spike所需要的最小電流。

LIF模型的簡單實現

這里給出我自己寫的一個python的函數，利用一階歐拉法迭代計算該微分方程進行求解，并且對閾值、不應期等情況進行一個計算和區分。

import numpy as np import matplotlib.pyplot as pltdef calc_next_step(Vm, I, step_t, remaining_refrac_time):Vl = -70Gl = 0.025C = 0.5if Vm > -50 and Vm < 0: #thresholdVm = 30 #spike potentialelif Vm > 0:Vm = -60 #reset potentialremaining_refrac_time = remaining_refrac_time*0 + 2 #reset everything to 2elif remaining_refrac_time>0:Vm = -60 #reset potentialremaining_refrac_time -= step_telse:Vm = Vm + step_t*(-Gl*(Vm-Vl) + I)/Cif remaining_refrac_time<0:remaining_refrac_time = remaining_refrac_time*0return Vm,remaining_refrac_timestep_t = 0.001 t = np.arange(0,500+step_t,step_t) Vm_out = np.zeros(np.shape(t)[0]) I = 0.9ind = 0 Vm_out[0]=-70 #init state remaining_refrac_time = 0 for tt in t[0:-1]:Vm_out[ind+1],remaining_refrac_time = calc_next_step(Vm_out[ind],I,step_t,remaining_refrac_time)ind += 1plt.plot(t,Vm_out) plt.xlabel("time /ms") plt.ylabel("Memberance Potential /mV") plt.show()

上述代碼模擬外部電流鉗輸入電流

的情況，代碼會產生如下結果：

可以看到，該模擬的細胞會產生一個持續恒定的動作電位。

通過給定不同的輸入電流

，我們可以得到不同Firing Rate的一個序列，如下圖。可以看到再輸入電流小于一定程度（在這個例子里是0.5nA）時，由于穩態的膜電位達不到閾值，不會觸發spike，但是當其高于閾值后，輸入電流越大，Firing Rate也越快。

在這里 輸入電流0.5nA相當于一個bifurcation point，在這里分叉出來兩種截然不同的情況：firing or not firing。

LIF模型的性能分析，Firing Rate，以及噪聲

接下來，我試圖量化LIF模型，對不同輸入電流的情況會產生不同Firing Rate的關系：我通過分別計算1）直接模擬神經元，計算實際Firing Rate，以及2）LIF模型理論可以產生的Firing Rate來進行對照。得到的結果如下圖：可以看到兩種方式計算出來的Firing Rate，和輸入電流的關系基本是一致的

理論的LIF模型僅模擬了理想情況，無噪聲時會產生的情形。于是接下來我對LIF模型加入噪聲進行模擬，此時LIF模型的微分方程如下：

這種對于噪音的微分方程迭代求解方法如下：需要注意跟號

這里給出一個微分方程求解的函數：

def calc_next_step_noise(Vm, I, step_t, remaining_refrac_time,sigma):Vl = -70Gl = 0.025C = 0.5wn = np.random.normal()if Vm > -50 and Vm < 0: #thresholdVm = 30 #spike potentialelif Vm > 0:Vm = -60 #reset potentialremaining_refrac_time = remaining_refrac_time*0 + 2 #reset everything to 2elif remaining_refrac_time>0:Vm = -60 #reset potentialremaining_refrac_time -= step_telse:Vm = Vm + step_t*(-Gl*(Vm-Vl) + I)/C + sigma*np.sqrt(step_t)*wnif remaining_refrac_time<0:remaining_refrac_time = remaining_refrac_time*0return Vm,remaining_refrac_time

帶著有噪聲的模型，再跑一遍各種輸入電流的情況，得到的結果如下圖：

我們可以看到，存在噪聲的情況下，即使原先不會達到閾值電位的比較低的輸入電流的神經元也可能會達到閾值并觸發動作電位。對于輸入電流在閾值附近的情況，噪聲會產生比較大的影響，也會對ISI產生很大的波動(紅色的線)。而對于紫色的，比較大的輸入電流的情況，噪聲對于ISI影響有限。這里畫出引入噪聲和不引入噪聲，Firing Rate和輸入電流關系的對比，如下圖：

可以看到，再閾值附近，容易產生比較大的波動，不過隨著輸入電壓的提高，噪聲帶來的相對誤差會不明顯很多。（其實我應該畫一下相對的誤差，用Firing Rate做個normalization，畢竟對于80Hz/82Hz，1Hz/3Hz差別還是很大的）

噪音也會讓bifurcation point變得不明顯。

在這里引入一個統計量：變異系數(coefficient of variation, CV)，來衡量Firing Rate的隨機程度。類似的統計量還有Fano Factor，都是衡量隨機程度的。

這里提一句題外話：對于一個泊松隨機的序列，CV=1，但是其實在大腦，尤其是前額葉的實際測量可以發現CV其實是大于1的，如果以泊松隨機來建模的話，我們可以認為泊松隨機的參數

也是隨時間而變化的，所以可以看到大腦內部加載的信息量所帶來的混亂程度是很高的（Goris et al. 2014）。

CV計算方法如下：

我這里還有一個Python函數的例子，通過給定電壓序列，先計算ISI（通過求Spike差分），再計算CV：

def extract_CV_from_voltage(voltage,step_t):spike_cv = np.zeros(np.shape(voltage)[0]) for ii in np.arange(0,np.shape(voltage)[0],1):idx_spike = np.transpose(np.argwhere(voltage[ii,:]==30))idx_spike_diff = np.diff(idx_spike)*step_tif len(idx_spike_diff[0]):mean_isi = np.mean(idx_spike_diff) #mean ISIstd_isi = np.std(idx_spike_diff)spike_cv[ii] = std_isi/mean_isielse:spike_cv[ii] = 1return spike_cv

然后我畫出Firing Rate和CV的關系，如下圖：每個點代表一次模擬，對于不同的輸入電流

我們可以看到，對于LIF模型，隨著Firing Rate的提高，神經元發放的穩定性也會逐漸提高，CV會逐漸下降。對于數學好的朋友，這個CV也是可以推導出來的(Ricciardi 1977; Renart et al. 2003)公式如下：

總而言之，LIF模型是一種建模神經元spike的比較簡單基本的模型，忽略了動作電位自身的快速特性：基本只考慮了細胞的被動特性，以及動作電位threshold的特性和refreactory period的特性。但是他可以很好的展示一些神經元動作電位序列的一些特征：比如CV，對噪聲的影響之類。

我不會講HH模型的具體內容，這個很多人說。下一篇文章我開始談談對于神經元之間突觸傳遞的建模。有兩種方法：Current Based or Conducatance Based. 同時還有三種模型: Kick Model, Filter Model and Kinetic Model. 以及一些對于LTP，LTD的建模。歡迎大家持續關注。

References

主要來自于不知道怎么引用的上課課件和講義。。。還有

Hodgkin, A. L., & Huxley, A. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve.The Journal of physiology,117(4), 500.

Koch, C., & Schutter, E. D. (1999). Biophysics of computation: Information processing in single neurons.Nature,398(6729), 678-678.

Ricciardi, L. M. (1977). Diffusion Processes and Related Topics in Biology. Springer, Berlin.

Renart, A., Brunel, N., & Wang, X. J. (2003). Mean-field theory of recurrent cortical networks: working memory circuits with irregularly spiking neurons.Computational neuroscience: A comprehensive approach, 432-490.

總結

以上是生活随笔為你收集整理的求给定精度的简单交错序列部分和_单个神经元的简单模型：Leaky integrate and fire (LIF) model...的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。