Brian

1. Simulation of spiking neural networks Brian http://www.briansimulator.org

2. The spirit of “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) Goals: Quick model coding Flexible

3. An example from systems neuroscience Sturzl, W., R. Kempter, and J. L. van Hemmen (2000, June). Theory of arachnid prey localization. Physical Review Letters 84 (24), 5668{71. PMID: 10991021

4. Learning Brian • 1 hour is not long enough for me to say everything • Will include detailed slides, but talk in a more general way • You can download the slides and go through it in more detail later (from Telluride website or Brian website) • Also look at the documentation page on the Brian website • Use the Brian email list for questions

5. Installing Brian • Instructions for installing Brian on our webpage: • http://www.briansimulator.org/docs/installation.html • Recommend using Python(x,y) for Windows users • http://www.pythonxy.com • Includes all packages, various IDEs, etc. • Brian available as a plugin, or can be downloaded separately.

6. Python • No time to introduce Python programming language, but there is an excellent tutorial online: • http://docs.python.org/tutorial/ • For scientific work, use the NumPy (numerical), SciPy (scientific) and Pylab/Matplotlib (plotting) libraries. Documentation (including tutorials) available online: • http://docs.scipy.org/doc/

7. Anatomy of a Brian script Import the Brian package 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

8. Units system • Standard SI names like volt, amp, etc. • Standard SI prefixes, mvolt, namp, etc. • Some short names: mV, ms, nA, etc. • Consistency checking: • 3*mV+2*nA will raise an exception • “dV/dt = -V” will raise an exception (RHS should have units of volt/second) • Useful for catching hard to find bugs • Can be switched off globally: • import brian_no_units before importing Brian

9. Defining a neuron model: Equations • Equations, each of one of the following forms: • dx/dt = -x/tau : voltdifferential equation • y = x*x : volt2equation • z = yalias • w : unitparameter • 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” • Can also get equations from brian.library, but won’t cover this here • Solvers: • Exact for linear DEs (detected automatically) • Euler for nonlinear DEs (method=‘Euler’) • 2nd order Runge-Kutta (method=‘RK’) • Exponential Euler (method=‘exponential_Euler’)

10. NeuronGroup creation group = NeuronGroup(N, eqs, reset=…, threshold=…, refractory=…, method=…) • Resets • Single value • String with Python statement(s), e.g. • ‘V=Vr’ where Vr could be a constant or another state variable • ‘Vt += dVt; V=Vr’ where Vt is another state variable • ‘V=Vr+rand()*sigma’ (rand and randn are known) • Arbitrary Python function f(group, spikes) • Refractory • Single value, note that this only works with single valued resets, otherwise it is ambiguous. In other cases, use CustomRefractoriness (see docs) • Threshold • Single value • String with Python expression, e.g. • ‘V>=Vt’ where Vt could be a constant or another state variable • ‘V+sigma*rand()>=Vt’ (will be in the next release of Brian) • Arbitrary Python function f(var1, var2, …) should return array of spike indices • EmpiricalThreshold(thresh, refrac) fires a spike if V>thresh but won’t fire another spike for period refrac, used for models without an instantaneous reset (e.g. HH)

11. Example: adaptive threshold

12. More on groups • Subgroups • subgp = group.subgroup(N) • subgp = group[i:j] • Can be used whenever a group is used • Use to structure network (e.g. into layers) • More computationally efficient to use one large group with several subgroups than to use several small groups • State variables • Access by name, e.g. group.V, group.ge, group.I • Standard groups • PoissonGroup(N, rates) – N neurons firing with given rates (scalar or vector, can be changed during a run) • SpikeGeneratorGroup(N, spikes) – N neurons that fire spikes specified by the user in the form spikes=[(i1, t1), (i2, t2), …] • PulsePacket(t, N, sigma) – N neurons firing at Gaussian dist. times

13. Connections • C = Connection(source, target, var) • Synapses from source group to target group acting on state variable var. • Weight matrix C.W • When neuron i in source fires: • target.var[j] += C.W[i, j] for each j • Features (next slides): • Building connectivity (full, random, custom) • Delays (homogeneous, heterogeneous) • Weight matrix structures (dense, sparse)

14. Building connectivity • C.connect_full(source, target, weight) • Source and target can be whole group or subgroups • weight can be single value, matrix or function w(i,j) • C.connect_random(source, target, p, weight) • p is the probability of synapse between i and j. • weight can be single value or function w(i, j) • C.connect_one_to_one(source, target, weight) • weight has to be a single value • Multiple connect_* calls allowed for different subgroups, etc. • Can build or modify connectivity by hand, e.g.: • for i in range(N): C.W[i, :] = N-abs(i-arange(N))

15. Delays • Connection(…, delay=…, max_delay=…) • Fixed constant, all synapses have the same (axonal) delay – no need to specify max_delay • delay=True, each synapse has its own delay, have to specify the maximum delay (larger uses more memory) • Specifying delays • Each C.connect_* method has a new delay keyword, which can have different values: • Scalar, all delays equal • Pair (min, max) uniformly distributed • Function f() called for each synapse • Function f(i, j) called for each synapse • Matrix • Specify by hand, C.delay[i,j] = …

16. Matrix structures • C = Connection(…, structure=…) • ‘dense’ • Numpy 2D array • Uses 8NM bytes for matrix of shape (N,M) • ‘sparse’ • Sparse matrix class with fast row access and reasonable column access • Uses 20 bytes per nonzero entry (12 bytes if column access not required) • Cannot be changed at runtime • ‘dynamic’ • Sparse matrix class with reasonable row and column access speeds • Uses 24 bytes per nonzero entry • Can be changed at runtime

17. Example: synfire chain

18. Plasticity • STDP • stdp = ExponentialSTDP(C, taupre, taupost, deltapre, deltapost, …) • stdp = STDP(C, …) • See docs • STP • stp = STP(C, taud, tauf, U) • Tsodyks-Markram model (see docs)

19. Recording activity • Spikes • M = SpikeMonitor(group) • M.spikes = [(i1, t1), (i2, t2), …] • M = SpikeCounter(group) • M.count = [count0, count1, …] • M.nspikes • PopulationSpikeCounter, StateSpikeMonitor • See docs • Others • ISIHistogramMonitor • PopulationRateMonitor • See docs • State variables • M=StateMonitor(group, var, record=…) • Records values of variable var in group • Record= • True, record all neurons (lots of memory) • [i1,i2,…], record given neurons • False (default), only record summary stats • M.plot() • M[i] array of recorded values for neuron i • M.times array of recording times • M.var, M.std summary stats

20. Network operations • Generic mechanism for online control of simulation • @network_operation(when=‘start’)def f(): … • Function f is called every time step, can do anything • Example: stop a network from running if too many spikes being produced: • …M=PopulationSpikeCounter(G)…@network_operationdef check_too_many_spikes(): if M.nspikes>1000: stop()

21. Simulation control • run(duration) creates network of all created objects • net = Network(obj1, obj2, …) for finer control, then call net.run(duration) • stop() and net.stop() • reinit() and net.reinit() • forget(obj) stops run() from using obj • recall(obj) lets run() use obj again • clear() forgets all objects, clear(True) forgets and deletes memory associated to all objects (useful for doing multiple runs)

22. Clocks • By default, Brian uses ‘defaultclock’: • defaultclock.t = 0*ms • defaultclock.dt = 0.1*ms • Can use multiple clocks if some objects need to be simulated with smaller dt than others • Create clock with clk = Clock(dt=1*ms) • Use clock for object with obj=Obj(…, clock=clk) • If using only one clock, don’t need to specify, if using multiple clocks, you have to specify for every object

23. Analysis packages • Brian has some statistical tools for spike trains: CVs, correlelograms, etc. • http://www.briansimulator.org/docs/analysis.html#statistics • NeuroTools is a set of Python packages for analysingneuroscientific data • http://neuralensemble.org/trac/NeuroTools • Disclaimer: I haven’t used it myself • The INCF maintains a list of relevant software • http://software.incf.org/software/ • Disclaimer: likewise

24. Optimisations • Use C code • Build optional C modules • Set the global preference useweave=True for more C code if you have gcc or Visual Studio installed • http://www.briansimulator.org/docs/compiledcode.html • Vectorise your code • http://www.briansimulator.org/docs/efficient.html • Soon: use ‘code generation’ (not yet documented)

25. Model fitting toolbox • Automatically fit spiking neuron models to electrophysiological data • Uses particle swarm optimisation • Adding genetic algorithms and others • Can use GPU for 60x speed improvement Rossant et al., Automatic fitting of spiking neuron models to electrophysiological recordings, Front. Neuroinformatics (2010)

26. “Brian hears” (soon) • Library for auditory modelling • Filtering • High pass, low pass, band pass, etc. • Gammatone, gammachirp, etc. • Head related transfer functions (HRTFs) • Standard neuron models • Meddis, Lyons, etc. • All in Brian framework

27. GPUs • Graphics processing units • Originally made for computer games (and therefore cheap) • Now being used for high-performance scientific computing • Many cores (up to 240 at the moment) but each must do similar computations (like SIMD, but slightly better) • Huge potential speed improvements • Existing papers (early work) show 8-80x improvement • In model fitting library we have 60x improvement • Currently shows a 30x improvement in general case

28. Thanks! http://www.briansimulator.org