Brian2开源框架学习笔记----6.不应期

定义不应期

G = NeuronGroup(N, model='...', threshold='...', reset='...', refractory=2*ms)

或者，它也可以是以计算时间为结果的字符串表达式。该表达式将在每次峰值后进行计算，并允许不同的不应期。比如说，以下的定义将允许不应期的时长在1ms到3ms之间。

G = NeuronGroup(N, model='...', threshold='...', reset='...', refractory='(1+2*rand())*ms')

如果说我们想对不同神经元的不应期进行建模，我们需要声明一个状态变量。

G = NeuronGroup(N, model='''...
                    ref : second''',threshold='...',
                    reset='...', refractory='ref')
# set the refractory period for each cell
G.ref =...

这个状态变量本身也可以是动态变量。例如，它可以作为一种适应机制，在每次峰值后增加它，并让它在峰值之间放松到稳态值:

refractory_0 = 2*ms
tau_refractory = 50*ms
G = NeuronGroup(N, model='''...
                            dref/dt = (refractory_0 - ref) / tau_refractory : second''',
                threshold='...', refractory='ref',
                reset='''...
                         ref += 1*ms''')
G.ref = refractory_0

 不太理解这个适应机制的效果，在spyder中运行如下代码：

"""
Spyder Editor

This is a temporary script file.
"""

from brian2 import *
from IPython import get_ipython
get_ipython().run_line_magic('matplotlib', 'inline')


start_scope()

tau = 1*ms
refractory_0 = 1*ms
tau_refractory = 5*ms
eqs = '''
dv/dt = (1-v)/tau : 1 (unless refractory)
dref/dt = (refractory_0 - ref) / tau_refractory : second 
'''

G = NeuronGroup(1, eqs, threshold='v>0.8', 
                reset='''
                v = 0
                ref += 1*ms''', refractory='ref', method='exact')

G.ref = refractory_0

M = StateMonitor(G, 'v', record=0)
run(50*ms)
plot(M.t/ms, M.v[0])
xlabel('Time (ms)')
ylabel('v');



show()

运行显示图像如下图所示，

 因为峰值的输入非常稳定，所以在两个峰值之间并没有发现不应期有什么变化。

创造一个不稳定的输入，为了看清输出的值的图像，在spyder中运行如下代码：

"""
Spyder Editor

This is a temporary script file.
"""

from brian2 import *
from IPython import get_ipython
get_ipython().run_line_magic('matplotlib', 'inline')

start_scope()

N = 100
tau = 10*ms
v0_max = 3.
duration = 100*ms
sigma = 0.2

eqs = '''
dv/dt = (v0-v)/tau+sigma*xi*tau**-0.5 : 1 (unless refractory)
v0 : 1
'''

G = NeuronGroup(N, eqs, threshold='v>1', reset='v=0', refractory=5*ms, method='euler')
M = StateMonitor(G, 'v', record=0)

G.v0 = 'i*v0_max/(N-1)'

run(duration)


plot(M.t/ms, M.v[0])
xlabel('Time (ms)')
ylabel('v');



show()

得到如下所示输入值的图像：

 把适应机制用在此输入上，看看效果：

"""
Spyder Editor

This is a temporary script file.
"""

from brian2 import *
from IPython import get_ipython
get_ipython().run_line_magic('matplotlib', 'inline')

start_scope()

N = 100
tau = 10*ms
v0_max = 3.
duration = 500*ms
sigma = 0.2
refractory_0 = 4*ms
tau_refractory = 50*ms

eqs = '''
dv/dt = (v0-v)/tau+sigma*xi*tau**-0.5 : 1 (unless refractory)
dref/dt = (refractory_0 - ref)/tau_refractory : second
v0 : 1
'''

G = NeuronGroup(N, eqs, threshold='v>0.2', 
                reset='''
                v=0
                ref += 10*ms''', 
                refractory='ref', method='euler')
M = StateMonitor(G, 'v', record=0)

G.v0 = 'i*v0_max/(N-1)'
G.ref = refractory_0

run(duration)


plot(M.t/ms, M.v[0])
xlabel('Time (ms)')
ylabel('v');



show()

这么看，感觉也没有什么变化……

调整参数，让阈值等条件都一样，进行对比……

没有加适应机制

加了适应机制

 肉眼并不能看出来一些很好的效果，为了更直观的展示效果，运行如下代码：

from brian2 import *
from IPython import get_ipython
get_ipython().run_line_magic('matplotlib', 'inline')

start_scope()

N = 100
tau = 10*ms
v0_max = 3.
duration = 1000*ms
sigma = 0.2
refractory_0 = 4*ms
tau_refractory = 50*ms

eqs = '''
dv/dt = (v0-v)/tau+sigma*xi*tau**-0.5 : 1 (unless refractory)
dref/dt = (refractory_0 - ref)/tau_refractory : second
v0 : 1
'''

G = NeuronGroup(N, eqs, threshold='v>0.2', 
                reset='''
                v=0
                ref += 10*ms''', 
                refractory='ref', method='euler')
M = StateMonitor(G, 'v', record=0)
H = StateMonitor(G, 'ref', record=0)

G.v0 = 'i*v0_max/(N-1)'
G.ref = refractory_0

run(duration)

figure(figsize=(12,4))
subplot(121)

plot(M.t/ms, M.v[0])
xlabel('Time (ms)')
ylabel('v');

subplot(122)
plot(H.t/ms, H.ref[0]/ms)
xlabel('Times(ms)')
ylabel=('ref')



show()

得到如下所示图片：

加了适应机制的

 左边的图片是对输入电压的响应随时间的变化，右边的图片是不应期随时间的变化。

 在Hodgkin-Huxley模型中，不应期可以这样定义：

G = NeuronGroup(N, model='...', threshold='v > -20*mV',
                refractory = 'v >= -20*mV')

refractory=2*ms和以下定义等同：

refractory='(t-lastspike)<=2*ms'

然而，对于不应期是模拟时间步长的倍数的常见情况，该表达式可能给出不一致的结果。由于浮点不精确，它的实际值可能略高于或低于模拟时间步长的倍数;因此，将其直接与不应期进行比较，迟早会导致不应期的一个时间步骤的结束。

为了避免这个情况：

refractory='timestep(t-lastspike,dt) <= timestep(2*ms, dt)'

时间步长函数由Brian提供，负责以安全的方式将时间转换为时间步长。

---

定义模型在不应期期间的行为

在下面的模型中，变量v在不应期继续更新，但如果超过阈值则不会引发峰值:

G = NeuronGroup(N, 'dv/dt = -v / tau : 1',
                threshold='v > 1', reset='v=0',
                refractory=2*ms)

比较一下：

from brian2 import *
from IPython import get_ipython
get_ipython().run_line_magic('matplotlib', 'inline')


start_scope()

tau = 3*ms
eqs = '''
dv/dt = (1-v)/tau : 1 (unless refractory)
'''

G = NeuronGroup(1, eqs, threshold='v>0.5', reset='v = 0',refractory=10*ms, method='exact')

M = StateMonitor(G, 'v', record=0)
run(50*ms)
plot(M.t/ms, M.v[0])
xlabel('Time (ms)')
ylabel('v');

unless refractory

from brian2 import *
from IPython import get_ipython
get_ipython().run_line_magic('matplotlib', 'inline')


start_scope()

tau = 3*ms
eqs = '''
dv/dt = (1-v)/tau : 1 
'''

G = NeuronGroup(1, eqs, threshold='v>0.5', reset='v = 0',refractory=10*ms, method='exact')

M = StateMonitor(G, 'v', record=0)
run(50*ms)
plot(M.t/ms, M.v[0])
xlabel('Time (ms)')
ylabel('v');

without：unless refractory

还有Brian支持的耐火度的第二种实现，一个或几个状态变量可以在不应期被停住。为了模拟这种行为，在耐火期间应该停止更新的变量可以标记为(除非耐火)标志:

G = NeuronGroup(N, '''dv/dt = -(v + w)/ tau_v : 1 (unless refractory)
                      dw/dt = -w / tau_w : 1''',
                threshold='v > 1', reset='v=0; w+=0.1', refractory=2*ms)

 在上面的模型中，v变量在尖峰后的2ms内被箝制在0，但适应变量w在此期间继续更新。此外，处于不应期的神经元的变量是只读的:传入的突触或其他代码将不会对v的值产生影响，直到它离开不应期。