Romain Brette & Dan Goodman Ecole Normale Supérieure Equipe Audition

1. http://www.briansimulator.org Romain Brette & Dan Goodman Ecole Normale Supérieure Equipe Audition romain.brette@ens.fr dan.goodman@ens.fr

2. The spirit of Goals: • Quick model coding • Flexible “A simulator should not only save the time of processors, but also the time of scientists” scientist computer Writing code often takes more time than running it models are defined by equations (rather than pre-defined)

3. briansimulator.org An example Ci P Pi Pe Ce

4. Installation briansimulator.org, click Download Download slides and examples: click Blog and Brian at NeuroComp We have Windows versions on USB sticks (ask Dan)

5. does an adaptive threshold 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, eqs, reset=reset, threshold='v>vt') IF.rest() PG = PoissonGroup(1, 500*Hz) C = Connection(PG, IF, 'v', weight=3*mV) run(100*ms)

6. does synaptic plasticity Song, S., K. D. Miller, and L. F. Abbott (2000). Competitive hebbian learning through spike timing-dependentsynapticplasticity. Nature Neurosci 3, 919-26.

7. does systems neuroscience Sturzl, W., R. Kempter, and J. L. van Hemmen (2000). Theory of arachnid prey localization. Physical Review Letters 84 (24), 5668

8. Work in progress • Parallel computing with graphics processing units (GPUs) • Up to 240 processors • Up to 100x speed improvement • Cheap hardware (a few €100) • Good for vectorised code If you want to contribute: http://groups.google.fr/group/brian-on-gpu

9. Python in 15 minutes seealso: http://docs.python.org/tutorial/

10. The Python console • On Windows, open IDLE. On Linux: type python • interpreted language • dynamic typing • garbage collector • space matters (signals structure) • object-oriented • many libraries

11. Writing a script • If you use IDLE, click on File>New window.Otherwise use anytext editor. Press F5 to execute

12. A simple program comment # Factorialfunction deffactorial(x): if x == 0: return1 else: return x * factorial(x-1) printfactorial(5) function definition untyped argument condition structure by indentation (block = aligned instructions) function call display

13. Numerical objects • Base types: int, long, float, complex • Other numerical types (vectors, matrices) defined in the Numpy library (in a few minutes) x=3+2 x+=1 y=100L z=x*(1+2j) u=2.3/7 x,y,z = 1,2,3 a = b = 123

14. Control structures x = 12 if x < 5or (x > 10and x < 20): print"Good value." if x < 5or10 < x < 20: print‘Good value as well.' for i in [0,1,2,3,4]: print"Number", i for i in range(5): print"Number", i while x >= 0: print"x is not always negative." x = x-1 list the same list

15. Lists mylist = [1,7,8,3] name = ["Jean","Michel"] x = [1,2,3,[7,8],"fin"] heterogeneous list first element = index 0 print name[0] name[1]="Robert" print mylist[1:3] print mylist[:3],mylist[:] print mylist[-1] print x[3][1] « slice »: index 3 not included last element x[3] is a list method (list = object) name.append("Georges") print mylist+name print mylist*2 concatenate Other methods: extend, insert, reverse, sort, remove…

16. List comprehensions carres=[x**2for i in range(10)] pairs=[x for i in range(10) if i % 2 == 0] = list of squares of integers from 0 to 9 = list of even integers between 0 and 9

17. Strings a="Bonjour" b='hello' c=""" Une phrase qui n'en finit pas """ print a+b print a[3:7] print b.capitalize() multiline string ≈ list of characters many methods (find, replace, split…)

18. Dictionaries dico={‘one':1,‘two':2,‘three':‘several'} print dico[‘three'] dico[‘four']=‘many‘ del dico[‘one'] key value for key in dico: print key,'->',dico[key] iterate all keys

19. Functions def power(x,exponent=2): return x**exponent print power(3,2) print power(7) print power(exponent=3,x=2) default value call with named arguments inline definition carre=lambda x:x**2

20. Modules loads the file ‘math.py’ or ‘math.pyc’ (compiled) import math print math.exp(1) from math import exp print exp(1) from math import * print exp(1) import only the exp object import everything You can work with several files (= modules), each one can define any number of objects (= variables, functions, classes)

21. Scipy & Pylab

22. Scipy & Numpy • Scientific libraries • Syntax ≈ Matlab • Many mathematical functions from scipy import * x=array([1,2,3]) M=array([[1,2,3],[4,5,6]]) M=ones((3,3)) z=2*x y=dot(M,x) from scipy.optimize import * print fsolve(lambda x:(x-1)*(x-3),2) vector matrix matrix product

23. Vectors et matrices • Base type in SciPy: array(= vector or matrix) from scipy import * x=array([1,2,3]) M=array([[1,2,3],[4,5,6]]) M=ones((3,2)) z=2*x+1 y=dot(M,x) vector (1,2,3) 1 2 3 4 5 6 matrix 1 1 1 1 1 1 matrix matrix product

24. Operations x+y x-y x*y x/y x**2 exp(x) sqrt(x) dot(x,y) dot(M,x) M.T M.max() M.sum() size(x) M.shape element-wise x² dot product matrix product transpose total number of elements

25. Indexing Vector indexing  lists (first element= 0) x[i] M[i,j] x[i:j] M[i,:] M[:,i] x[[1,3]] x[1:3]=[0,1] M[1,:]+=x (i+1)th element slice from x[i] to x[j-1] (i+1)th row (i+1)th column elements x[1] and x[3] x[1]=0, x[3]=1 add vector x to the 2nd row of M is a « view » on matrix M  copy ( reference) M[i,:] y=M[0,:] y[2]=5 x=z x[1]=3 M[0,2] is 5 copy: z[1] is 3 x=z.copy()

26. Construction x=array([1,2,3]) M=array([[1,2,3],[4,5,6]]) x=ones(5) M=zeros((3,2)) M=eye(3) M=diag([1,3,7]) x=rand(5) x=randn(5) x=arange(10) x=linspace(0,1,100) from lists vector of 1s zero matrix identity matrix diagonal matrix random vector in (0,1) random gaussian vector 0,1,2,...,9 100 numbers between 0 and 1

27. SciPy example: optimisation • Simple example, least squares polynomial fit

28. Vectorisation • How to write efficient programs? • Replace loops by vector operations for i in range(1000000): X[i]=1 X=ones(1000000) for i in range(1000000): X[i]=X[i]*2 X=X*2 for i in range(999999): Y[i]=X[i+1]-X[i] Y=X[1:]-X[:-1] for i in range(1000000): if X[i]>0.5: Y[i]=1 Y[X>.5]=1

29. Pylab

30. Pylab • Plotting library • Syntax ≈ Matlab • Many plotting functions plot from pylab import * plot([1,2,3],[4,5,6]) show() x y last instruction of script polar more examples: http://matplotlib.sourceforge.net/gallery.html

31. Plotting with PyLab hist(randn(1000)) contour(gaussian_filter(randn(100, 100), 5)) specgram(sin(10000*2*pi*linspace(0,1,44100)), Fs=44100) imshow(gaussian_filter(randn(100, 100), 5))

32. Brian At last!

36. Online help

37. Anatomy of a Brian script Import the Brian package + Scipy and Pylab Define some parameters, you can use physical units Define neuron equations, the “volt” term says that V has units of volts, and Brian checks the consistency of your equations. Create synapses – here recurrent random connectivity with probability .1 for each pair of neurons and given synaptic weight, spikes cause an instantaneous increase of amount weight in variable V of the target neuron Create neurons with given equations and fixed threshold and reset value. Record spiking activity Initialise state variables to random values. Run simulation Analyse results – do a raster plot

38. Model equations tau=10*ms eqs=''' dv/dt = (v0-v)/tau : volt u=v*v : volt**2 x=u v0 : volt ''' units differential equation equation alias parameter • Non-autonomous equations, use the reserved symbol “t” for time • Stochastic DEs • Use the term “xi” for noise with: • Has units s-1/2. Use typically as “sigma*xi*(2/tau)**0.5” • Solvers: • Exact for linear DEs (detected automatically, method=‘linear’) • Euler for nonlinear DEs (default, method=‘Euler’) • 2nd order Runge-Kutta (default, method=‘RK’) • Exponential Euler (implicit=True, method=‘exponential_Euler’)

39. Neuron groups units of v model equations group=NeuronGroup(N, model=’dv/dt = -v/tau : volt’, threshold=’v>15*mV’, reset=’v=0*mV’, refractory=2*ms, method=’implicit’) threshold condition what happens after spiking optional refractory period integration method is a vector with values of v for all neurons group.v group.v[5] is the value of v for neuron 5 P=group[:100] P=group.subgroup(100) defines a subgroup of 100 neurons Special groups: N Poisson processes PoissonGroup(N,rates=10*Hz) fires at specific times SpikeGeneratorGroup(N,spikes) spikes=[(3,2*ms),(1,10*ms),(2,15*ms)]

40. Connections modified variable synapses=Connection(P,Q,’v’) synapses[5,7]=2*mV P Q when neuron 5 spikes: Q.v[7]+=2*mV Building connections synapses.connect_full(P,Q,weight=2*mV) synapses.connect_random(P,Q,sparseness=.02,weight=2*mV) synapses.connect_full(P[50:100],Q[20:40], weight=lambda i,j:w0*exp(-(i-j)**2)) synapses=Connection(P,Q,’v’,sparseness=.02,weight=2*mV)

41. Delays Homogeneous delays synapses=Connection(P,Q,’v’,delay=2*ms) Heterogeneous delays synapses=Connection(P,Q,’v’,delay=True,max_delay=10*ms) synapses.delay[5,7]=2*ms synapses.connect_full(P,Q,weight=2*mV,delay=(1*ms,2*ms)) synapses.connect_full(P,Q,weight=2*mV, delay=lambda i,j:w0*exp(-(i-j)**2))

42. Monitors • Recording spike trains • Recording variables M=SpikeMonitor(P) run(100*ms) raster_plot(M) show() M=StateMonitor(P,'v',record=True) run(100*ms) plot(M.times/ms,M[0]/mV) show()

43. Network operations …M=SpikeMonitor(G)…@network_operationdef check_too_many_spikes(): if M.nspikes>1000: stop() … run(10*second) this function is called at every time step

44. Examples

45. Integrate-and-fire model • Stimulate an integrate-and-fire model with regularly spaced spikes • Tip: use SpikeGeneratorGroup

46. Integrate-and-fire model

47. Current-frequency curve • Plot the firing rate of an IF model as a function of input current • Tip: counter=SpikeCounter(group) counts spikes for every neuron

48. Current-frequency curve

49. Cortical column • Ring of noisy IF neurons with distance-dependent connections dv/dt=(v0-v)/ + / (t) w(i,j)=cos(2(i-j)/N)w0

50. Cortical column