【Brian2学习教程之三】Introduction to Brian part 3: Simulations

前言：

本次教程主要介绍跟仿真相关的内容。本文内容根据Brian2官网提供的技术文档进行整理总结。

链接：https://brian2.readthedocs.io/en/stable/resources/tutorials/3-intro-to-brian-simulations.html

多次运行
在运行期间改变一些内容
添加输入

正文：

多次运行

让我们从一个非常常见的任务开始：用一些改变的参数来多次运行一个仿真。接下来，我们看一个简单的例子，了解leaky integrate-and-fire(LIF)神经元的firing rate是如何被泊松脉冲神经元驱动随着其膜时间常数tau变化的。执行如下代码。

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallback# remember, this is here for running separate simulations in the same notebook
start_scope()
# Parameters
num_inputs = 100
input_rate = 10*Hz
weight = 0.1
# Range of time constants
tau_range = linspace(1, 10, 30)*ms
# Use this list to store output rates
output_rates = []
# Iterate over range of time constants
for tau in tau_range:# Construct the network each timeP = PoissonGroup(num_inputs, rates=input_rate)eqs = '''dv/dt = -v/tau : 1'''G = NeuronGroup(1, eqs, threshold='v>1', reset='v=0', method='exact')S = Synapses(P, G, on_pre='v += weight')S.connect()M = SpikeMonitor(G)# Run it and store the output firing rate in the listrun(1*second)output_rates.append(M.num_spikes/second)
# And plot it
plot(tau_range/ms, output_rates)
xlabel(r'$\tau$ (ms)')
ylabel('Firing rate (sp/s)');

程序运行结果如下图所示。

在程序执行时，耗时有点久。这是因为在这个程序中对于每一次循环，它都会从头重新创建对象。我们可以通过只设置一次网络来改进它。在循环之前存储一个网络状态的副本，并在每次迭代开始时恢复它。执行如下代码。

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallbackstart_scope()
num_inputs = 100
input_rate = 10*Hz
weight = 0.1
tau_range = linspace(1, 10, 30)*ms
output_rates = []
# Construct the network just once
P = PoissonGroup(num_inputs, rates=input_rate)
eqs = '''
dv/dt = -v/tau : 1
'''
G = NeuronGroup(1, eqs, threshold='v>1', reset='v=0', method='exact')
S = Synapses(P, G, on_pre='v += weight')
S.connect()
M = SpikeMonitor(G)
# Store the current state of the network
store()
for tau in tau_range:# Restore the original state of the networkrestore()# Run it with the new value of taurun(1*second)output_rates.append(M.num_spikes/second)
plot(tau_range/ms, output_rates)
xlabel(r'$\tau$ (ms)')
ylabel('Firing rate (sp/s)');

程序运行结果如下图所示。本例是一个非常简单的使用store和restore的例子，你也可以将它们用在更复杂的仿真中。例如，你想要运行一个较长的训练，之后运行多个测试，就可以在较长的训练后面加上一个store，并在每个测试运行之前加上一个restore。

从上图可以看出，输出曲线存在噪声，没有像我们预期的那样单调递增。这些噪声来源于我们每一次都重新运行一遍泊松神经元组。如果我们只是想看时间常数的影响，我们可以确保每次的峰值都是相同的。通过只运行泊松神经元组一次，然后记录其脉冲，最后创建一个新的SpikeGeneratorGroup输出这些已经被记录的脉冲即可实现。执行如下代码。

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallbackstart_scope()
num_inputs = 100
input_rate = 10*Hz
weight = 0.1
tau_range = linspace(1, 10, 30)*ms
output_rates = []
# Construct the Poisson spikes just once
P = PoissonGroup(num_inputs, rates=input_rate)
MP = SpikeMonitor(P)
# We use a Network object because later on we don't
# want to include these objects
net = Network(P, MP)
net.run(1*second)
# And keep a copy of those spikes
spikes_i = MP.i
spikes_t = MP.t
# Now construct the network that we run each time
# SpikeGeneratorGroup gets the spikes that we created before
SGG = SpikeGeneratorGroup(num_inputs, spikes_i, spikes_t)
eqs = '''
dv/dt = -v/tau : 1
'''
G = NeuronGroup(1, eqs, threshold='v>1', reset='v=0', method='exact')
S = Synapses(SGG, G, on_pre='v += weight')
S.connect()
M = SpikeMonitor(G)
# Store the current state of the network
net = Network(SGG, G, S, M)
net.store()
for tau in tau_range:# Restore the original state of the networknet.restore()# Run it with the new value of taunet.run(1*second)output_rates.append(M.num_spikes/second)
plot(tau_range/ms, output_rates)
xlabel(r'$\tau$ (ms)')
ylabel('Firing rate (sp/s)');

程序运行结果如下图所示。

从上图中可以看出，现在的噪声比之前要小得多，而且是单调递增的。因为每次输入的脉冲都是一样的，这意味着我们看到的是时间常数的作用，而非随机的脉冲。

注意到在上面的代码中，我们创建了Network对象。因为在循环中，如果我们调用run，它将会尝试模拟所有的对象，包括泊松神经元P，我们只想运行一次。我们使用Network显式地指定了我们想要包含的对象。

到目前为止，我们用到的技术是进行多次运行最简单的方法，但它们并不总是最有效的。由于上面的模型中只有一个输出神经元，所以我们可以简单的复制输出神经元并将时间常数设为这个神经元组的一个参数。执行如下代码。

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallbackstart_scope()
num_inputs = 100
input_rate = 10*Hz
weight = 0.1
tau_range = linspace(1, 10, 30)*ms
num_tau = len(tau_range)
P = PoissonGroup(num_inputs, rates=input_rate)
# We make tau a parameter of the group
eqs = '''
dv/dt = -v/tau : 1
tau : second
'''
# And we have num_tau output neurons, each with a different tau
G = NeuronGroup(num_tau, eqs, threshold='v>1', reset='v=0', method='exact')
G.tau = tau_range
S = Synapses(P, G, on_pre='v += weight')
S.connect()
M = SpikeMonitor(G)
# Now we can just run once with no loop
run(1*second)
output_rates = M.count/second # firing rate is count/duration
plot(tau_range/ms, output_rates)
xlabel(r'$\tau$ (ms)')
ylabel('Firing rate (sp/s)');

程序运行结果如下图所示。可以看出，这次程序运行快了很多。虽然在概念上更加复杂了，但是如果可以实现的话，它会更加高效。

接下来，我们看一下interspike之间的间隔时间的均值和标准差与时间常数之间的关系。执行如下代码。

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallbackstart_scope()
num_inputs = 100
input_rate = 10*Hz
weight = 0.1
tau_range = linspace(1, 10, 30)*ms
num_tau = len(tau_range)
P = PoissonGroup(num_inputs, rates=input_rate)
# We make tau a parameter of the group
eqs = '''
dv/dt = -v/tau : 1
tau : second
'''
# And we have num_tau output neurons, each with a different tau
G = NeuronGroup(num_tau, eqs, threshold='v>1', reset='v=0', method='exact')
G.tau = tau_range
S = Synapses(P, G, on_pre='v += weight')
S.connect()
M = SpikeMonitor(G)
# Now we can just run once with no loop
run(1*second)
output_rates = M.count/second # firing rate is count/durationtrains = M.spike_trains()
isi_mu = full(num_tau, nan)*second
isi_std = full(num_tau, nan)*second
for idx in range(num_tau):train = diff(trains[idx])if len(train)>1:isi_mu[idx] = mean(train)isi_std[idx] = std(train)
errorbar(tau_range/ms, isi_mu/ms, yerr=isi_std/ms)
xlabel(r'$\tau$ (ms)')
ylabel('Interspike interval (ms)');

程序运行结果如下图所示。注意到我们使用SpikeMonitor中的spike_trains()键，它用于指示对应神经元的spike times组成的数组中的值。

在运行期间改变一些内容

想象一个实验，将电流注入一个神经元，然后每10ms随机改变一次振幅。我们来看看是否可以用Hodgkin-Huxley型的神经元模型来建模。执行如下代码。

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallbackstart_scope()
# Parameters
area = 20000*umetre**2
Cm = 1*ufarad*cm**-2 * area
gl = 5e-5*siemens*cm**-2 * area
El = -65*mV
EK = -90*mV
ENa = 50*mV
g_na = 100*msiemens*cm**-2 * area
g_kd = 30*msiemens*cm**-2 * area
VT = -63*mV
# The model
eqs_HH = '''
dv/dt = (gl*(El-v) - g_na*(m*m*m)*h*(v-ENa) - g_kd*(n*n*n*n)*(v-EK) + I)/Cm : volt
dm/dt = 0.32*(mV**-1)*(13.*mV-v+VT)/(exp((13.*mV-v+VT)/(4.*mV))-1.)/ms*(1-m)-0.28*(mV**-1)*(v-VT-40.*mV)/(exp((v-VT-40.*mV)/(5.*mV))-1.)/ms*m : 1
dn/dt = 0.032*(mV**-1)*(15.*mV-v+VT)/(exp((15.*mV-v+VT)/(5.*mV))-1.)/ms*(1.-n)-.5*exp((10.*mV-v+VT)/(40.*mV))/ms*n : 1
dh/dt = 0.128*exp((17.*mV-v+VT)/(18.*mV))/ms*(1.-h)-4./(1+exp((40.*mV-v+VT)/(5.*mV)))/ms*h : 1
I : amp
'''
group = NeuronGroup(1, eqs_HH,threshold='v > -40*mV',refractory='v > -40*mV',method='exponential_euler')
group.v = El
statemon = StateMonitor(group, 'v', record=True)
spikemon = SpikeMonitor(group, variables='v')
figure(figsize=(9, 4))
for l in range(5):group.I = rand()*50*nArun(10*ms)axvline(l*10, ls='--', c='k')
axhline(El/mV, ls='-', c='lightgray', lw=3)
plot(statemon.t/ms, statemon.v[0]/mV, '-b')
plot(spikemon.t/ms, spikemon.v/mV, 'ob')
xlabel('Time (ms)')
ylabel('v (mV)');

程序运行结果如下图所示。

在上面的代码中，我们使用多次运行一个循环来达到这个目的。这很好，但并不是最高效的方法，因为每次运行都要做很多初始化的工作。下面介绍另外一种方法，执行如下代码。

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallbackstart_scope()
# Parameters
area = 20000*umetre**2
Cm = 1*ufarad*cm**-2 * area
gl = 5e-5*siemens*cm**-2 * area
El = -65*mV
EK = -90*mV
ENa = 50*mV
g_na = 100*msiemens*cm**-2 * area
g_kd = 30*msiemens*cm**-2 * area
VT = -63*mV
# The model
eqs_HH = '''
dv/dt = (gl*(El-v) - g_na*(m*m*m)*h*(v-ENa) - g_kd*(n*n*n*n)*(v-EK) + I)/Cm : volt
dm/dt = 0.32*(mV**-1)*(13.*mV-v+VT)/(exp((13.*mV-v+VT)/(4.*mV))-1.)/ms*(1-m)-0.28*(mV**-1)*(v-VT-40.*mV)/(exp((v-VT-40.*mV)/(5.*mV))-1.)/ms*m : 1
dn/dt = 0.032*(mV**-1)*(15.*mV-v+VT)/(exp((15.*mV-v+VT)/(5.*mV))-1.)/ms*(1.-n)-.5*exp((10.*mV-v+VT)/(40.*mV))/ms*n : 1
dh/dt = 0.128*exp((17.*mV-v+VT)/(18.*mV))/ms*(1.-h)-4./(1+exp((40.*mV-v+VT)/(5.*mV)))/ms*h : 1
I : amp
'''
group = NeuronGroup(1, eqs_HH,threshold='v > -40*mV',refractory='v > -40*mV',method='exponential_euler')
group.v = El
statemon = StateMonitor(group, 'v', record=True)
spikemon = SpikeMonitor(group, variables='v')
# we replace the loop with a run_regularly
group.run_regularly('I = rand()*50*nA', dt=10*ms)
run(50*ms)
figure(figsize=(9, 4))
# we keep the loop just to draw the vertical lines
for l in range(5):axvline(l*10, ls='--', c='k')
axhline(El/mV, ls='-', c='lightgray', lw=3)
plot(statemon.t/ms, statemon.v[0]/mV, '-b')
plot(spikemon.t/ms, spikemon.v/mV, 'ob')
xlabel('Time (ms)')
ylabel('v (mV)');

程序运行结果如下图所示。

我们使用run_regularly替换了之前有多个run调用的循环，能够特定地运行单个NeuronGroup，但还有更加灵活的方法。你可以使用network_operation运行任意的Python代码。执行如下代码。

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallbackstart_scope()
# Parameters
area = 20000*umetre**2
Cm = 1*ufarad*cm**-2 * area
gl = 5e-5*siemens*cm**-2 * area
El = -65*mV
EK = -90*mV
ENa = 50*mV
g_na = 100*msiemens*cm**-2 * area
g_kd = 30*msiemens*cm**-2 * area
VT = -63*mV
# The model
eqs_HH = '''
dv/dt = (gl*(El-v) - g_na*(m*m*m)*h*(v-ENa) - g_kd*(n*n*n*n)*(v-EK) + I)/Cm : volt
dm/dt = 0.32*(mV**-1)*(13.*mV-v+VT)/(exp((13.*mV-v+VT)/(4.*mV))-1.)/ms*(1-m)-0.28*(mV**-1)*(v-VT-40.*mV)/(exp((v-VT-40.*mV)/(5.*mV))-1.)/ms*m : 1
dn/dt = 0.032*(mV**-1)*(15.*mV-v+VT)/(exp((15.*mV-v+VT)/(5.*mV))-1.)/ms*(1.-n)-.5*exp((10.*mV-v+VT)/(40.*mV))/ms*n : 1
dh/dt = 0.128*exp((17.*mV-v+VT)/(18.*mV))/ms*(1.-h)-4./(1+exp((40.*mV-v+VT)/(5.*mV)))/ms*h : 1
I : amp
'''
group = NeuronGroup(1, eqs_HH,threshold='v > -40*mV',refractory='v > -40*mV',method='exponential_euler')
group.v = El
statemon = StateMonitor(group, 'v', record=True)
spikemon = SpikeMonitor(group, variables='v')
# we replace the loop with a network_operation
@network_operation(dt=10*ms)
def change_I():group.I = rand()*50*nA
run(50*ms)
figure(figsize=(9, 4))
for l in range(5):axvline(l*10, ls='--', c='k')
axhline(El/mV, ls='-', c='lightgray', lw=3)
plot(statemon.t/ms, statemon.v[0]/mV, '-b')
plot(spikemon.t/ms, spikemon.v/mV, 'ob')
xlabel('Time (ms)')
ylabel('v (mV)');

程序运行结果如下图所示。

现在让我们扩展这个例子，在多个神经元上运行，每个神经元都有不同的电容以观察它是如何影响细胞的行为。执行如下代码。

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallbackstart_scope()
# Parameters
area = 20000*umetre**2
Cm = 1*ufarad*cm**-2 * area
gl = 5e-5*siemens*cm**-2 * area
El = -65*mV
EK = -90*mV
ENa = 50*mV
g_na = 100*msiemens*cm**-2 * area
g_kd = 30*msiemens*cm**-2 * area
VT = -63*mV
N = 3
# The model
eqs_HH_2 = '''
dv/dt = (gl*(El-v) - g_na*(m*m*m)*h*(v-ENa) - g_kd*(n*n*n*n)*(v-EK) + I)/C : volt
dm/dt = 0.32*(mV**-1)*(13.*mV-v+VT)/(exp((13.*mV-v+VT)/(4.*mV))-1.)/ms*(1-m)-0.28*(mV**-1)*(v-VT-40.*mV)/(exp((v-VT-40.*mV)/(5.*mV))-1.)/ms*m : 1
dn/dt = 0.032*(mV**-1)*(15.*mV-v+VT)/(exp((15.*mV-v+VT)/(5.*mV))-1.)/ms*(1.-n)-.5*exp((10.*mV-v+VT)/(40.*mV))/ms*n : 1
dh/dt = 0.128*exp((17.*mV-v+VT)/(18.*mV))/ms*(1.-h)-4./(1+exp((40.*mV-v+VT)/(5.*mV)))/ms*h : 1
I : amp
C : farad
'''
group = NeuronGroup(N, eqs_HH_2,threshold='v > -40*mV',refractory='v > -40*mV',method='exponential_euler')
group.v = El
# initialise with some different capacitances
group.C = array([0.8, 1, 1.2])*ufarad*cm**-2*area
statemon = StateMonitor(group, variables=True, record=True)
# we go back to run_regularly
group.run_regularly('I = rand()*50*nA', dt=10*ms)
run(50*ms)
figure(figsize=(9, 4))
for l in range(5):axvline(l*10, ls='--', c='k')
axhline(El/mV, ls='-', c='lightgray', lw=3)
plot(statemon.t/ms, statemon.v.T/mV, '-')
xlabel('Time (ms)')
ylabel('v (mV)');

程序运行结果如下图所示。

但是从结果中来看感觉有些不对劲，对所有不同的神经元它们的输入电流都是不同的。让我们检查一下，执行如下代码。

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallbackstart_scope()
# Parameters
area = 20000*umetre**2
Cm = 1*ufarad*cm**-2 * area
gl = 5e-5*siemens*cm**-2 * area
El = -65*mV
EK = -90*mV
ENa = 50*mV
g_na = 100*msiemens*cm**-2 * area
g_kd = 30*msiemens*cm**-2 * area
VT = -63*mV
N = 3
# The model
eqs_HH_2 = '''
dv/dt = (gl*(El-v) - g_na*(m*m*m)*h*(v-ENa) - g_kd*(n*n*n*n)*(v-EK) + I)/C : volt
dm/dt = 0.32*(mV**-1)*(13.*mV-v+VT)/(exp((13.*mV-v+VT)/(4.*mV))-1.)/ms*(1-m)-0.28*(mV**-1)*(v-VT-40.*mV)/(exp((v-VT-40.*mV)/(5.*mV))-1.)/ms*m : 1
dn/dt = 0.032*(mV**-1)*(15.*mV-v+VT)/(exp((15.*mV-v+VT)/(5.*mV))-1.)/ms*(1.-n)-.5*exp((10.*mV-v+VT)/(40.*mV))/ms*n : 1
dh/dt = 0.128*exp((17.*mV-v+VT)/(18.*mV))/ms*(1.-h)-4./(1+exp((40.*mV-v+VT)/(5.*mV)))/ms*h : 1
I : amp
C : farad
'''
group = NeuronGroup(N, eqs_HH_2,threshold='v > -40*mV',refractory='v > -40*mV',method='exponential_euler')
group.v = El
# initialise with some different capacitances
group.C = array([0.8, 1, 1.2])*ufarad*cm**-2*area
statemon = StateMonitor(group, variables=True, record=True)
# we go back to run_regularly
group.run_regularly('I = rand()*50*nA', dt=10*ms)
run(50*ms)
figure(figsize=(9, 4))
for l in range(5):axvline(l*10, ls='--', c='k')
axhline(El/mV, ls='-', c='lightgray', lw=3)
plot(statemon.t/ms, statemon.I.T/nA, '-')
xlabel('Time (ms)')
ylabel('I (nA)');

程序运行结果如下图所示。

从结果中可以看出，每次都是不同的。但是为什么呢？group.run_regularly(I =rand()*50*nA’, dt=10*ms)似乎应该给每个神经元赋予相同的值。但是，实际上每一个神经元运行run_regularly时，都会单独地运行，而I是参数，因此对于每个神经元是不同的。要想解决这个问题，只需要将I变成一个共享变量，也就是说它对每个神经元都有相同的值。执行如下代码。

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallbackstart_scope()
# Parameters
area = 20000*umetre**2
Cm = 1*ufarad*cm**-2 * area
gl = 5e-5*siemens*cm**-2 * area
El = -65*mV
EK = -90*mV
ENa = 50*mV
g_na = 100*msiemens*cm**-2 * area
g_kd = 30*msiemens*cm**-2 * area
VT = -63*mV
N = 3
# The model
eqs_HH_3 = '''
dv/dt = (gl*(El-v) - g_na*(m*m*m)*h*(v-ENa) - g_kd*(n*n*n*n)*(v-EK) + I)/C : volt
dm/dt = 0.32*(mV**-1)*(13.*mV-v+VT)/(exp((13.*mV-v+VT)/(4.*mV))-1.)/ms*(1-m)-0.28*(mV**-1)*(v-VT-40.*mV)/(exp((v-VT-40.*mV)/(5.*mV))-1.)/ms*m : 1
dn/dt = 0.032*(mV**-1)*(15.*mV-v+VT)/(exp((15.*mV-v+VT)/(5.*mV))-1.)/ms*(1.-n)-.5*exp((10.*mV-v+VT)/(40.*mV))/ms*n : 1
dh/dt = 0.128*exp((17.*mV-v+VT)/(18.*mV))/ms*(1.-h)-4./(1+exp((40.*mV-v+VT)/(5.*mV)))/ms*h : 1
I : amp (shared) # everything is the same except we've added this shared
C : farad
'''
group = NeuronGroup(N, eqs_HH_3,threshold='v > -40*mV',refractory='v > -40*mV',method='exponential_euler')
group.v = El
group.C = array([0.8, 1, 1.2])*ufarad*cm**-2*area
statemon = StateMonitor(group, 'v', record=True)
group.run_regularly('I = rand()*50*nA', dt=10*ms)
run(50*ms)
figure(figsize=(9, 4))
for l in range(5):axvline(l*10, ls='--', c='k')
axhline(El/mV, ls='-', c='lightgray', lw=3)
plot(statemon.t/ms, statemon.v.T/mV, '-')
xlabel('Time (ms)')
ylabel('v (mV)');

程序运行结果如下图所示。此时运行结果就没有问题了。

添加输入

现在我们来考虑一个由正弦输入驱动的神经元。神经元模型使用LIF模型以简化方程。执行如下代码。

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallbackstart_scope()
A = 2.5
f = 10*Hz
tau = 5*ms
# Create a TimedArray and set the equations to use it
t_recorded = arange(int(200*ms/defaultclock.dt))*defaultclock.dt
I_recorded = TimedArray(A*sin(2*pi*f*t_recorded), dt=defaultclock.dt)
eqs = '''
dv/dt = (I-v)/tau : 1
I = I_recorded(t) : 1
'''
G = NeuronGroup(1, eqs, threshold='v>1', reset='v=0', method='exact')
M = StateMonitor(G, variables=True, record=True)
run(200*ms)
plot(M.t/ms, M.v[0], label='v')
plot(M.t/ms, M.I[0], label='I')
xlabel('Time (ms)')
ylabel('v')
legend(loc='best');

程序运行结果如下图所示。注意到在这个例子中我们将sin函数直接加入到方程中，我们必须使用参数method=’euler’，因为exact积分器在这里不起作用。然而，TimeArray被认为在时间上保持不变，所以可以使用线性积分。这意味着你讲不能从这两种方法中得到相同的行为，原因有两个。首先，exact和euler方法给出的结果略有不同。第二，sin随时间不是恒定的，然而TimedArray是。

现在，为了证明TimeArray适用于任意电流，让我们创建一个奇怪的“recorded”电流并运行它。执行如下代码。

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallbackstart_scope()
A = 2.5
f = 10*Hz
tau = 5*ms
# Let's create an array that couldn't be
# reproduced with a formula
num_samples = int(200*ms/defaultclock.dt)
I_arr = zeros(num_samples)
for _ in range(100):a = randint(num_samples)I_arr[a:a+100] = rand()
I_recorded = TimedArray(A*I_arr, dt=defaultclock.dt)
eqs = '''
dv/dt = (I-v)/tau : 1
I = I_recorded(t) : 1
'''
G = NeuronGroup(1, eqs, threshold='v>1', reset='v=0', method='exact')
M = StateMonitor(G, variables=True, record=True)
run(200*ms)
plot(M.t/ms, M.v[0], label='v')
plot(M.t/ms, M.I[0], label='I')
xlabel('Time (ms)')
ylabel('v')
legend(loc='best');

程序运行结果如下图所示。

最后，让我们完成一个从文件中实际读取数据的示例。看看你能不能明白这个例子是如何工作的。执行如下代码。

from brian2 import *
prefs.codegen.target = ' numpy ' #use the Python fallbackstart_scope()
from matplotlib.image import imread
img = (1-imread('brian.png'))[::-1, :, 0].T
num_samples, N = img.shape
ta = TimedArray(img, dt=1*ms) # 228
A = 1.5
tau = 2*ms
eqs = '''
dv/dt = (A*ta(t, i)-v)/tau+0.8*xi*tau**-0.5 : 1
'''
G = NeuronGroup(N, eqs, threshold='v>1', reset='v=0', method='euler')
M = SpikeMonitor(G)
run(num_samples*ms)
plot(M.t/ms, M.i, '.k', ms=3)
xlim(0, num_samples)
ylim(0, N)
xlabel('Time (ms)')
ylabel('Neuron index');

程序运行结果如下。注意，读入的图片文件格式为8bit灰度图像，且图像尺寸不能过大，左图为输入图像，右图为程序运行结果图。

【Brian2学习教程之三】Introduction to Brian part 3: Simulations相关推荐

Brian2学习笔记一 Introduction to Brian part1:Neurons
Brian2学习笔记一 Introduction to Brian part 1 :Neurons 1. 前言 2. 正文 2.1 单位系统(Units system) 2.2 一个简单的模型(A s ...
【Brian2学习教程之一】Introduction to Brian part 1: Neurons
前言: 目前正在学习脉冲神经网络建模,由于研究该领域的人员相对较少,所以很难在网上找到很多系统全面的中文学习教程,只能自己找英文原版资料进行学习.为了能够吸引更多脉冲神经网络建模领域的同学相互交流进步 ...
ZBrush全面入门学习教程 Schoolism – Introduction to ZBrush
ZBrush全面入门学习教程 Schoolism – Introduction to ZBrush ZBrush全面入门学习教程 Schoolism – Introduction to ZBrush ...
【Brian2学习教程之二】Introduction to Brian part 2: Synapses
前言: 本次教程主要介绍跟突触相关的内容.本文内容根据Brian2官网提供的技术文档进行整理总结. 链接:https://brian2.readthedocs.io/en/stable/resourc ...
Brian2学习教程——Intro to Brian part 1: Neurons
参考原英文教程网址:https://brian2.readthedocs.io/en/stable/resources/tutorials/1-intro-to-brian-neurons.html ...
Brian2学习笔记_Introduction to Brian part1:Neurons
本文根据Brian2官方英文教程进行翻译总结. Brian2官方安装教程链接: https://brian2.readthedocs.io/en/stable/introduction/install ...
Blender 3.0基础入门学习教程 Introduction to Blender 3.0
成为Blender通才,通过这个基于项目的循序渐进课程学习所有主题的基础知识. 你会学到什么教程获取:Blender 3.0基础入门学习教程 Introduction to Blender 3.0- ...
ABAP 标准培训教程 BC400 学习笔记之三：ABAP 编程语言的特性和基本构成要素
SAP ABAP 标准培训教程 BC400 中对 ABAP 编程语言特性的总结如下: Is typed - 强类型编程语言,任何 ABAP 变量在其声明时,其数据类型就已经的确定下来了. Enable ...
Brian2学习笔记
Brian2学习笔记前言运行环境写点有用的没用的简介引用安装 python编译安装 pip安装 C++ code generation的安装要求测试使用教程 Tutorial par ...

【Brian2学习教程之三】Introduction to Brian part 3: Simulations

【Brian2学习教程之三】Introduction to Brian part 3: Simulations相关推荐

最新文章

热门文章