File: developer-simulationprinciples.txt

package info (click to toggle)
brian 1.4.3-1
 `123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231` ``````Simulation principles --------------------- The following paper outlines the principles of Brian simulation: Goodman, D and Brette R (2008), `Brian: a simulator for spiking neural networks in Python `__, Front. Neuroinform. doi:10.3389/neuro.11.005.2008. This one describes the simulation algorithms, which are based on vectorisation: Brette R and Goodman, DF, `Vectorised algorithms for spiking neural network simulation `__, Neural Computation (in press). Sample script ~~~~~~~~~~~~~ Below we present a Brian script, and a translation into pure Python to illustrate the basic principles of Brian simulations. Original Brian script ..................... A script in Brian:: ''' Very short example program. ''' from brian import * from time import time N=10000 # number of neurons Ne=int(N*0.8) # excitatory neurons Ni=N-Ne # inhibitory neurons p=80./N duration=1000*ms 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(N,model=eqs, threshold=-50*mV,reset=-60*mV) P.v=-60*mV+10*mV*rand(len(P)) Pe=P.subgroup(Ne) Pi=P.subgroup(Ni) Ce=Connection(Pe,P,'ge',weight=1.62*mV,sparseness=p) Ci=Connection(Pi,P,'gi',weight=-9*mV,sparseness=p) M=SpikeMonitor(P) trace=StateMonitor(P,'v',record=0) t1=time() run(1*second) t2=time() print "Simulated in",t2-t1,"s" print len(M.spikes),"spikes" subplot(211) raster_plot(M) subplot(212) plot(trace.times/ms,trace[0]/mV) show() Equivalent in pure Python ......................... The script above translated into pure Python (no Brian):: ''' A pure Python version of the CUBA example, that reproduces basic Brian principles. ''' from pylab import * from time import time from random import sample from scipy import random as scirandom """ Parameters """ N=10000 # number of neurons Ne=int(N*0.8) # excitatory neurons Ni=N-Ne # inhibitory neurons mV=ms=1e-3 # units dt=0.1*ms # timestep taum=20*ms # membrane time constant taue=5*ms taui=10*ms p=80.0/N # connection probability (80 synapses per neuron) Vt=-1*mV # threshold = -50+49 Vr=-11*mV # reset = -60+49 we=60*0.27/10 # excitatory weight wi=-20*4.5/10 # inhibitory weight duration=1000*ms """ Equations --------- eqs=''' dv/dt = (ge+gi-(v+49*mV))/(20*ms) : volt dge/dt = -ge/(5*ms) : volt dgi/dt = -gi/(10*ms) : volt ''' This is a linear system, so each update corresponds to multiplying the state matrix by a (3,3) 'update matrix' """ # Update matrix A=array([[exp(-dt/taum),0,0], [taue/(taum-taue)*(exp(-dt/taum)-exp(-dt/taue)),exp(-dt/taue),0], [taui/(taum-taui)*(exp(-dt/taum)-exp(-dt/taui)),0,exp(-dt/taui)]]).T """ State variables --------------- P=NeuronGroup(4000,model=eqs, threshold=-50*mV,reset=-60*mV) """ S=zeros((3,N)) """ Initialisation -------------- P.v=-60*mV+10*mV*rand(len(P)) """ S[0,:]=rand(N)*(Vt-Vr)+Vr # Potential: uniform between reset and threshold """ Connectivity matrices --------------------- Pe=P.subgroup(3200) # excitatory group Pi=P.subgroup(800) # inhibitory group Ce=Connection(Pe,P,'ge',weight=1.62*mV,sparseness=p) Ci=Connection(Pi,P,'gi',weight=-9*mV,sparseness=p) """ We_target=[] We_weight=[] for _ in range(Ne): k=scirandom.binomial(N,p,1)[0] target=sample(xrange(N),k) target.sort() We_target.append(target) We_weight.append([1.62*mV]*k) Wi_target=[] Wi_weight=[] for _ in range(Ni): k=scirandom.binomial(N,p,1)[0] target=sample(xrange(N),k) target.sort() Wi_target.append(target) Wi_weight.append([-9*mV]*k) """ Spike monitor ------------- M=SpikeMonitor(P) will contain a list of (i,t), where neuron i spiked at time t. """ spike_monitor=[] # Empty list of spikes """ State monitor ------------- trace=StateMonitor(P,'v',record=0) # record only neuron 0 """ trace=[] # Will contain v(t) for each t (for neuron 0) """ Simulation ---------- run(duration) """ t1=time() t=0*ms while tVt).nonzero()[0] # List of neurons that meet threshold condition # PROPAGATION OF SPIKES # Excitatory neurons spikes=(S[0,:Ne]>Vt).nonzero()[0] # In Brian we actually use bisection to speed it up for i in spikes: S[1,We_target[i]]+=We_weight[i] # Inhibitory neurons spikes=(S[0,Ne:N]>Vt).nonzero()[0] for i in spikes: S[2,Wi_target[i]]+=Wi_weight[i] # Reset neurons after spiking S[0,all_spikes]=Vr # Reset membrane potential # Spike monitor spike_monitor+=[(i,t) for i in all_spikes] # State monitor trace.append(S[0,0]) t+=dt t2=time() print "Simulated in",t2-t1,"s" print len(spike_monitor),"spikes" """ Plot ---- subplot(211) raster_plot(M) subplot(212) plot(trace.times/ms,trace[0]/mV) show() Here we cheat a little. """ from brian import raster_plot class M: pass M.spikes=spike_monitor subplot(211) raster_plot(M) subplot(212) plot(arange(len(trace))*dt/ms,array(trace)/mV) show() ``````