Introduction to Neural Networks

Programming Lab 5: Spiking neural networks


Spiking neural networks (SNNs)



Today, we will scale up our simulations from single neurons to networks. This will allow us to simulate the dynamics of the networks of neurons in a small patch of cortex, using networks of LIF neurons. We will study the spiking network simulation below:

from brian2 import *

# network parameters
N = 1000; Ne = int( 0.8*N ); Ni = int( 0.2*N );

# membrane parameters
C = 200 * pF
taum = 20 * msecond
gL = C / taum
El = -70 * mV
Vth = -50 * mV

# synapse parameters
Ee = 0 * mvolt
Ei = -80 * mvolt
taue = 5 * msecond
taui = 5 * msecond

# excitatory and inhibitory weights
g = 4 							# balance between excitation and inhibition
we = 0.5 * nS 						# excitatory synaptic weight
wi = ( 7 * g * (we/nS) ) * nS 				# inhibitory synaptic weight

# input parameters
input_rate = 2*Hz					# Poisson input rate

# membrane equation
eqs = Equations('''
dv/dt = ( gL*(El-v) + ge*(Ee-v) + gi*(Ei-v) ) * (1./C) : volt
dge/dt = -ge*(1./taue) : siemens
dgi/dt = -gi*(1./taui) : siemens''')

# setup network
P = NeuronGroup( N, model=eqs, threshold='v>Vth', reset='v=El', refractory=5*ms )
Pe = P[0:Ne]; Pi = P[Ne:]

# setup excitatory and inhibitory connections
Ce = Synapses( Pe, P, on_pre='ge+=we' )
Ce.connect( True, p=0.01 )
Ci = Synapses( Pi, P, on_pre='gi+=wi' )
Ci.connect( True, p=0.01 )

# setup Poisson input
Ie = PoissonInput( P, 'ge', Ne, input_rate, weight=we )

# initialize network
P.v = randn( len(P) ) * 5 * mV - 70 * mV

# record spikes + rates (excitatory population)
M = SpikeMonitor( Pe )
R = PopulationRateMonitor( Pe )

# run simulation
run( 1 * second, report='text' )

# Make plot
figure( figsize=(15,5) )
xlabel( 'time (ms)' ); ylabel( 'cell id' )
plot( M.t/ms, M.i, '.k' )
show()

This parameter regime creates a low-rate, irregular spiking regime called the asynchronous-irregular state. These highly disordered states are an important model in computational neuroscience, as they recreate the firing patterns of neurons in cortex. Specifically, during spontaneous waking activity, neurons fire with low rate (< 10 Hz) and high irregularity (\(C_V \approx 1\)). We will study how this state emerges as a function of the balance between excitation and inhibition.

The most important parameters determining the state of the network are g (the balance of excitation and inhibition) and input_rate (the firing rate of the external inputs). We can study the mean firing rate in the network using the PopulationRateMonitor object in Brian (R in the code above). Try plotting the firing rate as a function of time for one run of the simulation and study how that matches up with the spike train raster. Also try running the network for different numbers of neurons, to understand how this affects the spiking state.

Finally, we can study how states of regular, highly synchronous spiking emerge in the network at different values of g and input_rate. To quantify the irregularity in the network, we can calculate \(C_V\) using the code below:

def coefficient_of_variation( spike_monitor, N ):
	# calculate Cv of inter-spike intervals
	cv = np.zeros( N )
	for ii in range(N):
		tau = diff( spike_monitor.spike_trains()[ii] )
		cv[ii] = std( tau ) / mean( tau )

If you need to run multiple simulations of the network, you can put the network setup and running code into a function

Experiments: (1) Simulate the network with increasing values of input_rate and observe how the network state changes. (2) Simulate the network for increasing values of g and again observe the change in network state. (3) Quantify the dependence on both control parameters, either in a one- or two-dimensional scan across values.