the above are from my run of the code presented so far, with some modified plotting.
I added a 40Hz sine to compare to the gamma rhythm mentioned in the original paper:
visually inspecting the rasters, the oscillation seems less pronounced (and less synchronous), but looking at the rates, there does appear to be a 40Hz oscillation in the inhibitory firing rate.
here is my refactored code. I moved some constants around, and added some plotting:
expand for code
from brian2 import *
import numpy as np
start_scope()
# seed(11922)
n = 1000
num_ex = int(0.8*n)
num_inh = n - num_ex
defaultclock.dt = 0.5*ms
refractory_period = 0.1*ms # I don't think the original paper uses a refractory period
duration = 1.0*second
# The model
neuron_eqs = Equations('''
dv/dt = (0.04*v**2+5.*v+140- u + g_exc - g_inh + noise * randn())/ms : 1 (unless refractory)
du/dt = a*(b*v-u)/ms : 1 (unless refractory)
g_exc :1
g_inh :1
noise :1
a :1
b :1
c :1
d :1
''')
reset_eq = 'v = c; u += d; g_exc=0; g_inh=0'
N = NeuronGroup(n, model=neuron_eqs, threshold='v>=30', reset = reset_eq, refractory=refractory_period,
method='euler')
Ne = N[:num_ex]
Ni = N[num_ex:]
re = np.random.random(num_ex) ; ri = np.random.random(num_inh)
Ne.noise = 5.0 ; Ni.noise = 2.0
Ne.a = 0.02 ; Ni.a = 0.02 + 0.08 * ri
Ne.b = 0.2 ; Ni.b = 0.25 - 0.05 * ri
Ne.c = -65.0 + 15.0 * re**2 ; Ni.c = -65.0
Ne.d = 8.0 - 6.0 * re**2 ; Ni.d = 2.0
Ne.v = -65.0 ; Ni.v = -65.0
Ne.u = Ne.v * Ne.b ; Ni.u = Ni.v * Ni.b
Se = Synapses(Ne, N, model= "w:1", on_pre= 'g_exc_post += w')
Si = Synapses(Ni, N, model= "w:1", on_pre= 'g_inh_post += w')
prob_connect_e = 0.1
prob_connect_i = prob_connect_e
w_e = 0.5
w_i = 1.0
Se.connect(p=prob_connect_e)
Si.connect(p=prob_connect_i)
# Initialization
Se.w = 'rand() * w_e'
Si.w = 'rand() * w_i'
M = SpikeMonitor(N)
Ex_rate = PopulationRateMonitor(Ne)
Inh_rate = PopulationRateMonitor(Ni)
#%%
run(duration)
#%%
figure(figsize=(8,7))
plot(M.t/ms, M.i,'k_',ms=2)
plot([0,duration/ms],[num_ex,num_ex],'r--')
xlim(0, duration/ms)
ylim(0, n)
xlabel('Time (ms)')
ylabel('Neuron index');
show()
#%%
t = Ex_rate.t
f = 40
gamma = sin(2*pi*f*t/second)*5
figure(figsize=(12,2))
plot(Ex_rate.t/ms, Ex_rate.smooth_rate(width=2*ms),'r')
plot(Inh_rate.t/ms, Inh_rate.smooth_rate(width=2*ms),'b')
plot(t/ms, gamma,'g')
ylim(-20,70)
legend(['excitatory','inhibitory','40Hz sine'])
xlabel('Time (ms)')
ylabel('Firing Rate (spikes/sec)')
now for some remaining issues / discrepancies:
1.) the paper sets the I→E weights to be uniformly distributed between -1.0 and 0.0:
S=[0.5*rand(Ne+Ni,Ne), -rand(Ne+Ni,Ni)];
the brian code we have seems to use positive weights for the “inhibitory” population:
Se.w = 'rand() * 0.5'
Si.w = 'rand() * 1.0'
maybe you’ve encoded that negative influence somewhere else but I don’t see it.
If I modify that line to Si.w = rand() * -1.0 the raster looks very different, arguably less like the original paper.
2.) the resultant network is NOT scale invariant. Synchrony seems to strongly depend on the number of total neurons
- EDIT: this CAN be corrected by scaling the synaptic weights relative to the population size
- i.e.
w_e = 0.5 * (1000 / n)
3.) I don’t see anywhere in the original Izhikevich paper where they have a probability of connection, I think it’s essentially “all-to-all” whereas the code above uses a 10% probability of connection.
4.) The brian implementation is NOT invariant to changing the size of the timestep
- the timing of the initial transient bump in firing rate depends strongly on what the timestep (dt) is