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.