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 <http://www.frontiersin.org/neuroinformatics/paper/10.3389/neuro.11/005.2008/>`__,
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 <http://www.briansimulator.org/WordPress/wp-content/uploads/2010/10/algorithms-preprint.pdf>`__,
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 t<duration:
	    # STATE UPDATES
	    S[:]=dot(A,S)
	
	    # Threshold
	    all_spikes=(S[0,:]>Vt).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()