|
|
Demo The following are some examples of Brian in action, with the source code on the left and the output on the right. Note that all these examples are for Brian 1, examples for Brian 2 can be found in the documentation.
Random network
from brian import *
eqs = '''
dv/dt = (ge+gi-(v+49*mV))/(20*ms) : volt
dge/dt = -ge/(5*ms) : volt
dgi/dt = -gi/(10*ms) : volt
'''
P = NeuronGroup(4000, eqs, threshold=-50*mV, reset=-60*mV)
P.v = -60*mV+10*mV*rand(len(P))
Pe = P.subgroup(3200)
Pi = P.subgroup(800)
Ce = Connection(Pe, P, 'ge', weight=1.62*mV, sparseness=0.02)
Ci = Connection(Pi, P, 'gi', weight=-9*mV, sparseness=0.02)
M = SpikeMonitor(P)
run(1*second)
raster_plot(M)
show()
|
 |
Synfire chain
from brian import *
# Neuron model parameters
Vr = -70*mV
Vt = -55*mV
taum = 10*ms
taupsp = 0.325*ms
weight = 4.86 * mV
# Neuron model
eqs=Equations('''
dV/dt=(-(V-Vr)+x)*(1./taum) : volt
dx/dt=(-x+y)*(1./taupsp) : volt
dy/dt=-y*(1./taupsp)+25.27*mV/ms+
(39.24*mV/ms**0.5)*xi : volt
''')
# Neuron groups
P = NeuronGroup(N=1000, model=eqs,
threshold=Vt,reset=Vr,refractory=1*ms)
Pinput = PulsePacket(t=50*ms,n=85,sigma=1*ms)
# The network structure
Pgp = [ P.subgroup(100) for i in range(10)]
C = Connection(P,P,'y')
for i in range(9):
C.connect_full(Pgp[i],Pgp[i+1],weight)
Cinput = Connection(Pinput,Pgp[0],'y')
Cinput.connect_full(weight=weight)
# Record the spikes
Mgp = [SpikeMonitor(p) for p in Pgp]
Minput = SpikeMonitor(Pinput)
monitors = [Minput]+Mgp
# Setup the network, and run it
P.V = Vr + rand(len(P)) * (Vt-Vr)
run(100*ms)
# Plot result
raster_plot(showgrouplines=True,*monitors)
show()
|
 |
Network of sparsely connected inhibitory integrate-and-fire neurons
Dynamics of a network of sparsely connected inhibitory integrate-and-fire neurons. Individual neurons fire irregularly at low rate but the network is in an oscillatory global activity regime where neurons are weakly synchronized.
Reference: Brunel N, Hakim V, Fast Global Oscillations in Networks of Integrate-and-Fire Neurons with Low Firing Rates, Neural Computation 11, 1621–1671 (1999)
from brian import *
# Network parameters
N = 5000
Vr = 10 * mV
theta = 20 * mV
tau = 20 * ms
delta = 2 * ms
taurefr = 2 * ms
duration = .1 * second
C = 1000
sparseness = float(C)/N
J = .1 * mV
muext = 25 * mV
sigmaext = 1 * mV
# Neuron model
eqs = "dV/dt=(-V+muext+sigmaext*sqrt(tau)*xi)/tau : volt"
group = NeuronGroup(N, eqs, threshold=theta,
reset=Vr, refractory=taurefr)
group.V = Vr
# Connections
conn = Connection(group, group, state='V', delay=delta,
weight=-J, sparseness=sparseness)
# Monitors
M = SpikeMonitor(group)
# Run
run(duration)
# Plot
raster_plot(M)
show()
|
 |
Topographically connected network
Topographic map – an example of complicated connections. Two layers of neurons: the first layer is connected randomly to the second one in a
topographical way. The second layer has random lateral connections.
from brian import *
N=100
tau=10*ms
tau_e=2*ms # AMPA synapse
eqs='''
dv/dt=(I-v)/tau : volt
dI/dt=-I/tau_e : volt
'''
rates=zeros(N)*Hz
rates[N/2-10:N/2+10]=ones(20)*30*Hz
layer1=PoissonGroup(N,rates=rates)
layer2=NeuronGroup(N,model=eqs,threshold=10*mV,reset=0*mV)
topomap=lambda i,j:exp(-abs(i-j)*.1)*3*mV
feedforward=Connection(layer1,layer2,sparseness=.5,weight=topomap)
lateralmap=lambda i,j:exp(-abs(i-j)*.05)*0.5*mV
recurrent=Connection(layer2,layer2,sparseness=.5,weight=lateralmap)
spikes=SpikeMonitor(layer2)
run(1*second)
subplot(211)
raster_plot(spikes)
subplot(223)
imshow(feedforward.W.todense(), interpolation='nearest', origin='lower')
title('Feedforward connection strengths')
subplot(224)
imshow(recurrent.W.todense(), interpolation='nearest', origin='lower')
title('Recurrent connection strengths')
show()
|
 |
Adaptive threshold model
A model with adaptive threshold (increases with each spike).
from brian import *
eqs='''
dv/dt = -v/(10*ms) : volt
dvt/dt = (10*mV-vt)/(15*ms) : volt
'''
reset='''
v=0*mV
vt+=3*mV
'''
IF = NeuronGroup(1, model=eqs,reset=reset,threshold='v>vt')
IF.rest()
PG = PoissonGroup(1, 500*Hz)
C = Connection(PG, IF, 'v',weight=3*mV)
Mv = StateMonitor(IF, 'v', record=True)
Mvt = StateMonitor(IF, 'vt', record=True)
run(100*ms)
plot(Mv.times/ms, Mv[0]/mV)
plot(Mvt.times/ms, Mvt[0]/mV)
show()
|
 |
Hodgkin-Huxley network
from brian import *
# Parameters
area=20000*umetre**2
Cm=(1*ufarad*cm**-2)*area
gl=(5e-5*siemens*cm**-2)*area
El=-60*mV
EK=-90*mV
ENa=50*mV
g_na=(100*msiemens*cm**-2)*area
g_kd=(30*msiemens*cm**-2)*area
VT=-63*mV
# Time constants
taue=5*ms
taui=10*ms
# Reversal potentials
Ee=0*mV
Ei=-80*mV
we=6*nS # excitatory synaptic weight (voltage)
wi=67*nS # inhibitory synaptic weight
# The model
eqs=Equations('''
dv/dt = (gl*(El-v)+ge*(Ee-v)+gi*(Ei-v)-g_na*(m*m*m)*h*(v-ENa)-g_kd*(n*n*n*n)*(v-EK))/Cm : volt
dm/dt = alpham*(1-m)-betam*m : 1
dn/dt = alphan*(1-n)-betan*n : 1
dh/dt = alphah*(1-h)-betah*h : 1
dge/dt = -ge*(1./taue) : siemens
dgi/dt = -gi*(1./taui) : siemens
alpham = 0.32*(mV**-1)*(13*mV-v+VT)/(exp((13*mV-v+VT)/(4*mV))-1.)/ms : Hz
betam = 0.28*(mV**-1)*(v-VT-40*mV)/(exp((v-VT-40*mV)/(5*mV))-1)/ms : Hz
alphah = 0.128*exp((17*mV-v+VT)/(18*mV))/ms : Hz
betah = 4./(1+exp((40*mV-v+VT)/(5*mV)))/ms : Hz
alphan = 0.032*(mV**-1)*(15*mV-v+VT)/(exp((15*mV-v+VT)/(5*mV))-1.)/ms : Hz
betan = .5*exp((10*mV-v+VT)/(40*mV))/ms : Hz
''')
P=NeuronGroup(4000,model=eqs,
threshold=EmpiricalThreshold(threshold=-20*mV,refractory=3*ms),
implicit=True,freeze=True)
Pe=P.subgroup(3200)
Pi=P.subgroup(800)
Ce=Connection(Pe,P,'ge',weight=we,sparseness=0.02)
Ci=Connection(Pi,P,'gi',weight=wi,sparseness=0.02)
# Initialization
P.v=El+(randn(len(P))*5-5)*mV
P.ge=(randn(len(P))*1.5+4)*10.*nS
P.gi=(randn(len(P))*12+20)*10.*nS
# Record a few trace
trace=StateMonitor(P,'v',record=[1,10,100])
run(1000*msecond)
plot(trace.times/ms,trace[1]/mV)
plot(trace.times/ms,trace[10]/mV)
plot(trace.times/ms,trace[100]/mV)
show()
|
 |
|