# 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 t<duration:
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()
```