# Realization the model of network in Brain 2

Hello, everybody!

I am trying to understand Brian simulator and translate my realization of neural network with Izhikevich neurons to Brian simulator.

The code of the model on C++ is look like :

#include <iostream>
#include <fstream>
#include <cstdlib>
#include <cmath>

const float h     = .5f; 	  // integration time step
const int   Tsim  = 1000/.5f; // simulation time in discrete samples
const int   Nexc  = 100; 	  // Number of excitatory neurons
const int   Ninh  = 25;  	  // Number of inhibitory neurons
const int   Nneur = Nexc + Ninh;
const int   Ncon  = Nneur*Nneur*0.1f; //Number of links,
// 0.1 is the probability of communication between 2 random neurons

float Vms[Nneur][Tsim]; // membrane potentials
float Ums[Nneur][Tsim]; // auxiliary variables of the Izhikevich model
float Iex[Nneur]; 		//external direct current applied to the neuron
float Isyn[Nneur]; 		// synaptic current per neuron

int pre_conns[Ncon];    // presynaptic neuron indices
int post_conns[Ncon];   // postsynaptic neuron indices
float weights[Ncon]; 	// bond weights
float y[Ncon][Tsim];    // variable modulating synaptic current depending on the spikes at the presynapse

float psc_excxpire_time = 4.0f;
float minWeight = 50.0f;
float maxWeight = 100.0f;

// Neuron parameters
float Iex_max = 40.0f;
float a		  = 0.02f;
float b		  = 0.5f;
float c		  = -40.0f;
float d 	  = 100.0f;
float k		  = 0.5f;
float Vr	  = -60.0f;
float Vt	  = -45.0f;
float Vpeak	  = 35.0f;
float V0	  = -60.0f;
float U0	  = 0.0f;  /
float Cm      = 50.0f;

float spike_times[Nneur*Tsim]; // spike times
int spike_neurons[Nneur*Tsim]; // corresponding neuron numbers
int spike_num = 0;

using namespace std;

void init_connections(){
for (int con_idx = 0; con_idx < Ncon; ){
// randomly select postsypantic and presynaptic neurons
int pre = rand() % Nneur;
int post = rand() % Nneur;
pre_conns[con_idx] = pre;
post_conns[con_idx] = post;
weights[con_idx] = (rand() % ((int)(maxWeight - minWeight)*10))/10.0f + minWeight;
if (pre >= Nexc){
// if the presynaptic neuron is inhibitory, then the weight of the connection goes with a minus sign
weights[con_idx] = -weights[con_idx];
}
y[con_idx][0] = 0.0f;
con_idx++;
}
}

void init_neurons(){
for (int neur_idx = 0; neur_idx < Nneur; neur_idx++){
// randomly scatter the applied currents
Iex[neur_idx] = ((float) rand() / RAND_MAX) * Iex_max;
Isyn[neur_idx] = 0.0f;
Vms[neur_idx][0] = V0;
Ums[neur_idx][0] = U0;
}
}

float izhik_Vm(int neuron, int time){
return (k*(Vms[neuron][time] - Vr)*(Vms[neuron][time] - Vt) - Ums[neuron][time] + Iex[neuron] + Isyn[neuron])/Cm;
}

float izhik_Um(int neuron, int time){
return a*(b*(Vms[neuron][time] - Vr) - Ums[neuron][time]);
}

void save2file(){
ofstream res_file;
res_file.open("rastr.csv");
for (int k = 0; k < spike_num; k++){
res_file << spike_times[k] << "; " << spike_neurons[k] + 1 << "; " << endl;
}
res_file.close();

}

int main(){
init_connections();
init_neurons();
float expire_coeff = exp(-h/psc_excxpire_time);
for (int t = 1; t < Tsim; t++){
// go through all neurons
for (int neur = 0; neur < Nneur; neur++){
Vms[neur][t] = Vms[neur][t-1] + h*izhik_Vm(neur, t-1);
Ums[neur][t] = Ums[neur][t-1] + h*izhik_Um(neur, t-1);
Isyn[neur] = 0.0f;

if (Vms[neur][t-1] >Vpeak){
Vms[neur][t] = c;
Ums[neur][t] = Ums[neur][t-1] + d;
spike_times[spike_num] = t*h;
spike_neurons[spike_num] = neur;
spike_num++;
}
}

// go through all connections
for (int con = 0; con < Ncon; con++){
y[con][t] = y[con][t-1]*expire_coeff;

if (Vms[pre_conns[con]][t-1] > Vpeak){
y[con][t] += 1.0f;
}
Isyn[post_conns[con]] += y[con][t]*weights[con];
}
}
save2file();
return 0;
}

My implementation on Brian2:



from brian2 import *

start_scope()

N = 125

seed(11922)

# Parameters of neuron

Iext_max = 40 * mV/ms # maximum current applied to the neuron
a = 0.02 / ms
b = 0.5 /ms
c = -40 * mV # the value of the membrane potential to which it resets after the spike
d = 100 * mV/ms
k = 0.5
Vr = -60 * mV
Vt = -45 * mV
Vpeak = 35 * mV # the maximum value of the membrane potential at which there is a reset to the value with
V0 = -60 * mV # initial value for membrane potentia
U0 = 0 * mV/ms # the initial value for the auxiliary variable
Cm = 50 # electrical capacitance of a neuron, dimension of pcF

minWeight = 50 * mV/ms
maxWeight = 100 * mV/ms
psc_excxpire_time = 4 * ms

# The model

eqs = Equations(’’’
dVm/dt = ((k/ms/mV)(Vm**2 - VmVr - VmVt + VrVt) - Um + Iext + Isyn) / Cm : volt
dUm/dt = a*(b*(Vm - Vr) - Um) : volt/second
Iext : volt/second
Isyn : volt/second
‘’’)

P = NeuronGroup(N, model=eqs, threshold=‘Vm>Vpeak’, reset = ‘Vm = c; Um = Um + d’,
method=‘euler’)
Pe = P[:int(0.8N)]
Pi = P[int(0.8
N):]

eqs_synapse = ‘’’
w_e : volt/second
w_i : volt/second
dy/dt = -y / psc_excxpire_time : 1 (clock-driven)

‘’’
eqs_pre = ‘’’
y += 1

‘’’

Ce = Synapses(Pe, P, model= eqs_synapse, on_pre=eqs_pre+‘Isyn_post = yw_e’,
method=‘euler’)
Ci = Synapses(Pi, P, model= eqs_synapse, on_pre=eqs_pre+'Isyn_post = y
w_i’,
method=‘euler’)

Ce.connect(p=0.2)
Ci.connect(p=0.2)

# Initialization

P.Iext = 'rand() * Iext_max ’

P.Vm = ‘V0’
P.Um = ‘U0’
P.Isyn = 0
Ce.w_e = ‘rand() * (maxWeight - minWeight) + minWeight’
Ci.w_i = ‘(-1)*(rand() * (maxWeight - minWeight) + minWeight)’
Ce.y = 0
Ci.y = 0
seed()

# Record a few traces

#M = SpikeMonitor(P)

M = SpikeMonitor(Pe)
M2 = SpikeMonitor(Pi)

trace = StateMonitor(P, variables =[‘Vm’,‘Iext’], record=True)

trace2 = StateMonitor(Ce, variables =[‘y’], record=True)

run(1*second)
figure(figsize=(16,5))
subplot(121)
#plot(M.t/ms, M.i, ‘.k’)
plot(M.t/ms, M.i, ‘.r’, ms=3)
plot(M2.t/ms, M2.i + int(N * 0.8), ‘.b’, ms=3)
xlim(0, 1000)
ylim(0, N)
xlabel(‘Time (ms)’)
ylabel(‘Neuron index’);
subplot(122)
plot(trace2.t/ms,trace2[9].y)

In the implementation on Brain2, it is not possible to get network synchronization. What did I do wrong in the implementation on Brian2.

Hi @serge25 ,

I’m not sure that the code you posted above actually works. I’m guessing that you didn’t put the python code in the code python  enclosure. Some multiplications are missing for example on_pre=eqs_pre+'Isyn_post = yw_e' or Pe = P[:int(0.8N)]. So it is quite hard to tell what is going on in your model without the exact code. Could you please update your question and provide a runnable code?

Meanwhile, it seems you attempt to use summed variables to collect synaptic currents.
If this is the case, you need to indicate that this is a ‘summed’ variable. Code should be something like that:

eqs_synapse = '''
w_e : volt/second
w_i : volt/second
dy/dt = -y / psc_excxpire_time : 1 (clock-driven)
Isyn_post = y*w_e : volt/second (summed)
'''


without this (summed), brian2 just puts the last computed synapse in the postsynaptic current, and all other synapses are ignored.

Also, it seems that both, excitatory and inhibitory, synapses have the same time constant: psc_excxpire_time = 4 * ms and they have the first order dynamics:  dy/dt = -y / psc_excxpire_time : 1 (clock-driven). You can move synaptic dynamics into the neuron equations (because they linearly summed) and reduce the computational load ~25 times. The code may look like that

eqs = Equations('''
dVm/dt = ((k/ms/mV)*(Vm**2 - Vm*Vr - Vm*Vt + Vr*Vt) - Um + Iext + y) / Cm : volt
dUm/dt = a*(b*(Vm - Vr) - Um) : volt/second
Iext : volt/second
dy/dt =  -y / psc_excxpire_time : volt/second
''')

Ce = Synapses(Pe, P, model= 'w : volt/second', on_pre='y_post += w ')
Ci = Synapses(Pi, P, model= 'w : volt/second', on_pre='y_post += w ')

Ce.w = w_e
Ci.w = w_i


In this case, you don’t need (summed) variables, because synapses instantaneously update postsynaptic y-variable.

Note that I cannot run this code, so all of this just a suggestion

Thank you for your answer. My apologies for previous edition of my post . The code would be like this

from brian2 import *

start_scope()

N = 125

seed(11922)

# Parameters of neuron

Iext_max = 40 * mV/ms  # maximum current applied to the neuron
a = 0.02 / ms
b = 0.5 /ms
c = -40 * mV  # the value of the membrane potential to which it resets after the spike
d = 100 * mV/ms
k = 0.5
Vr = -60 * mV
Vt = -45 * mV
Vpeak = 35 * mV   # the maximum value of the membrane potential at which there is a reset to the value with
V0 = -60 * mV # initial value for membrane potentia
U0 = 0 * mV/ms    # the initial value for the auxiliary variable
Cm = 50  # electrical capacitance of a neuron, dimension of pcF

minWeight = 50 * mV/ms
maxWeight = 100 * mV/ms
psc_excxpire_time = 4 * ms

# The model
eqs = Equations('''
dVm/dt = ((k/ms/mV)*(Vm**2 - Vm*Vr - Vm*Vt + Vr*Vt) - Um + Iext + Isyn) / Cm : volt
dUm/dt = a*(b*(Vm - Vr) - Um) : volt/second
Iext : volt/second
Isyn : volt/second
''')

P = NeuronGroup(N, model=eqs, threshold='Vm>Vpeak', reset = 'Vm = c; Um = Um + d',
method='euler')
Pe = P[:int(0.8*N)]
Pi = P[int(0.8*N):]

eqs_synapse = '''
w_e : volt/second
w_i : volt/second
dy/dt = -y / psc_excxpire_time : 1 (clock-driven)

'''
eqs_pre = '''
y += 1

'''

Ce = Synapses(Pe, P, model= eqs_synapse, on_pre=eqs_pre+'Isyn_post = y*w_e',
method='euler')
Ci = Synapses(Pi, P, model= eqs_synapse, on_pre=eqs_pre+'Isyn_post = y*w_i',
method='euler')

Ce.connect(p=0.2)
Ci.connect(p=0.2)

# Initialization
P.Iext = 'rand() * Iext_max '

P.Vm = 'V0'
P.Um = 'U0'
P.Isyn = 0
Ce.w_e = 'rand() * (maxWeight - minWeight) + minWeight'
Ci.w_i = '(-1)*(rand() * (maxWeight - minWeight) + minWeight)'
Ce.y = 0
Ci.y = 0
seed()

# Record a few traces

#M = SpikeMonitor(P)

M = SpikeMonitor(Pe)
M2 = SpikeMonitor(Pi)

trace = StateMonitor(P, variables =['Vm','Iext'], record=True)

trace2 = StateMonitor(Ce, variables =['y'], record=True)

run(1*second)
figure(figsize=(16,5))
subplot(121)
#plot(M.t/ms, M.i, '.k')
plot(M.t/ms, M.i, '.r', ms=3)
plot(M2.t/ms, M2.i + int(N * 0.8), '.b', ms=3)
xlim(0, 1000)
ylim(0, N)
xlabel('Time (ms)')
ylabel('Neuron index');
subplot(122)
plot(trace2.t/ms,trace2[9].y)


This code isn’t working

eqs_synapse = '''
w_e : volt/second
w_i : volt/second
dy/dt = -y / psc_excxpire_time : 1 (clock-driven)
Isyn_post = y*w_e : volt/second (summed)
'''


because of error like multivariable Isyn.

I have and other question: how is right to append at brian simulator variable which dependency from other variable (update with some threshold) and how set dependence Isyn from neurotransmitter concentration (like Tsodyks -Markram model - x,z,u).

For example:

dX/dt = -X / relax_time + beta / (1 + exp(-y + y_thr))


Here the code which actually works.

from brian2 import *

start_scope()

N = 125

seed(11922)

Iext_max = 40 * mV/ms  # maximum current applied to the neuron
a = 0.02 / ms
b = 0.5 /ms
c = -40 * mV  # the value of the membrane potential to which it resets after the spike
d = 100 * mV/ms
k = 0.5
Vr = -60 * mV
Vt = -45 * mV
Vpeak = 35 * mV   # the maximum value of the membrane potential at which there is a reset to the value with
V0 = -60 * mV # initial value for membrane potentia
U0 = 0 * mV/ms    # the initial value for the auxiliary variable
Cm = 50  # electrical capacitance of a neuron, dimension of pcF

minWeight = 50 * mV/ms
maxWeight = 100 * mV/ms
psc_excxpire_time = 4 * ms

# The model
eqs = Equations('''
dVm/dt = ((k/ms/mV)*(Vm**2 - Vm*Vr - Vm*Vt + Vr*Vt) - Um + Iext + y) / Cm : volt
dUm/dt = a*(b*(Vm - Vr) - Um) : volt/second
Iext : volt/second
dy/dt = -y / psc_excxpire_time  : volt/second
''')

P = NeuronGroup(N, model=eqs, threshold='Vm>Vpeak', reset = 'Vm = c; Um = Um + d',
method='euler')
Pe = P[:int(0.8*N)]
Pi = P[int(0.8*N):]

Ce = Synapses(Pe, P, model= "w:volt/second", on_pre= 'y += w')
Ci = Synapses(Pi, P, model= "w:volt/second", on_pre= 'y += w')

Ce.connect(p=0.2)
Ci.connect(p=0.2)

# Initialization
P.Iext = 'rand() * Iext_max '

P.Vm = 'V0'
P.Um = 'U0'
Ce.w = 'rand() * (maxWeight - minWeight) + minWeight'
Ci.w = '(-1)*(rand() * (maxWeight - minWeight) + minWeight)'
seed()

# Record a few traces

#M = SpikeMonitor(P)

M = SpikeMonitor(Pe)
M2 = SpikeMonitor(Pi)

trace = StateMonitor(P, variables =['Vm','Iext'], record=True)

trace2 = StateMonitor(Ce, variables =['y'], record=True)

run(1*second)
figure(figsize=(16,5))
subplot(121)
#plot(M.t/ms, M.i, '.k')
plot(M.t/ms, M.i, '.r', ms=3)
plot(M2.t/ms, M2.i + int(N * 0.8), '.b', ms=3)
xlim(0, 1000)
ylim(0, N)
xlabel('Time (ms)')
ylabel('Neuron index');
subplot(122)
plot(trace2.t/ms,trace2[9].y)
show()


Note that I moved variable y into neuron’s equations, and that made everything much simpler. Model synchronizes pretty well…

Well, it isn’t clear to me what exactly you are trying to achieve. I would strongly recommend, to consult with documentation first. Everything, including Synapses, is pretty well documented. In a very general picture, I want to move all possible equations to the neuron’s side. It is faster. If I need any local synaptic variables (for example STDP, or facilitation/depression), I forested to move equations to the synapse’s side. Please check with documentation and try to implement what you want using examples from there.

Please note, I’m NOT Brian’s father, just a Brian user who is trying to help a colleague. Pretty sure that @mstimberg and @dan have much more efficient and better-suited code

Hi
sorry , i have one problem in the neuron equation , here ‘Simple Model of Spiking Neurons’ is something different with what you wrote here and may you tell me what is ‘y’ parameter here ? because i use your code to model three izhikevich neurouns with astrocyte and i want to know this equation can help me or not

Hi @Maryidea,

Sorry for the corrections. I didn’t get what was the question. My apology for previous version of the answer.

y here is a dynamic variable, which is synaptic conductance. Because, and only if, synaptic time constants are the same, linear summation of derivatives is equal to the sum of right-hand-sides of differential equations. So, you can move differential equation into neuron equation and just sum presynaptic spikes.

Here how it works.
Say, for example, you have a network of Izhikevich’s models connected through inhibitory synapses. Equations are

\frac{dv}{dt} = 0.04v^2 + 5v + 140 -u -\sum_i w_i s_i \left(v - E_{syn} \right)
\frac{du}{dt} = a(bv-u)

and for each synapse, you have

\frac{ds_i}{dt} = -\frac{s_i}{\tau_s} + \sum_j\delta(t-t'_{ji})

where t'_{ji} is time of a j-th spike of i-th presinaptic neuron.

You can code this model directly like that

equ="""
dv/dt = (0.04*v**2+5.*v+140- u-gs*(u-Es) + I)/ms : 1 (unless refractory)
du/dt = a*(b*v-u)/ms                             : 1 (unless refractory)
gs                                               : 1
"""
n = NeuronGroup(100, equ, threshold="v>30", reset="v=c; u=d", refractory=0.1*ms)
s = Synapses(n,n,model="""
w                   : 1
ds/dt = -s/tau_s/ms : 1 (clock-driven)
gs = w*s            : 1(summed)
""",
on_pre = "s+1"


All above will work just fine, but slow. It is because for each connection, Brian will solve a differential equation. It will be 2 differential equations for neuron dynamics and 100 - differential equations for synapses. 102*100 = 10,200 equations in total!

Now if and only if all tua_s are the same, \sum_i w_i s_i(t), ds_i/dt = \sum_j \delta(t'_{ji}-t) - s_i/\tau_s is equal to ds/dt = \sum_{i,j} w_i \delta(t'_{ji}-t) - s/\tau_s. So … you can use one differential question instead of 100!
The code looks like that

equ="""
dv/dt = (0.04*v**2+5.*v+140- u-s*(u-Es) + I)/ms : 1 (unless refractory)
du/dt = a*(b*v-u)/ms                             : 1 (unless refractory)
ds/dt = -s/tau_s                                 : 1
"""
n = NeuronGroup(100, equ, threshold="v>30", reset="v=c; u=d", refractory=0.1*ms)
s = Synapses(n,n,model="w : 1", on_pre = "s_post+w"


These two implementations produce the same results, but the second one uses 3*100 = 300 equation instead of 10200! Obviously, this is much faster.

2 Likes

Hi @rth ,

Can you help me with one question. I’m trying to implement Izhikevich network from Simple Model of Spiking Neurons. I can’t get the results expected in the article. What’s wrong with my code?

from brian2 import *
import numpy as np

start_scope()
n = 1000
seed(11922)

defaultclock.dt = 0.5*ms

# 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

''')

N = NeuronGroup(n, model=neuron_eqs, threshold='v>=30', reset = 'v = c; u += d', refractory=0.1*ms,
method='euler')

Ne = N[:int(0.8*n)]
Ni = N[int(0.8*n):]

re = np.random.random(int(0.8*n))      ; ri = np.random.random(int(0.2*n))
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')

Se.connect(p=0.1)
Si.connect(p=0.1)

# Initialization
Se.w = 'rand() * 0.5'
Si.w = 'rand() * 1'

M = SpikeMonitor(N)

run(1*second)

figure(figsize=(16,5))
plot(M.t/ms, M.i,'.', ms=3)
xlim(0, 1000)
ylim(0, n)
xlabel('Time (ms)')
ylabel('Neuron index');

show()


Hi @serge25

Well, at this moment, I cannot run your code to check what exactly isn’t working. Your code looks pretty innocent and should work well. Could you please specify what exactly you want to achieve and what this code does not do?

The default time step in the code defaultclock.dt = 0.5*ms but refractory is refractory=0.1*ms, so neurons do not have refractory at all just one simulation step.

For the Izhikevich model, the time-step must be 0.02 ms or even less. The best performance/accuracy is when dt=0.01*ms

Another catch may be in reset. The code resets only v and u variables, but g_x continues to sum inputs. So you may need to reset them too, just to avoid overflow. Something like this:

reset = 'v = c; u += d; g_exc=0; g_inh=0'


As a result of the simulation, I want to reproduce the synchronization of neurons in alpha and gamma rhythms.

I’ve played a bit with the code. And it seems, with corrected reset, it produces a similar rasterplot to what is in the paper. However, all oscillations are gone, just a chaotic network activity. I’m wondering whether @mstimberg or @adam can spot something I have overseen.

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
1 Like

I just had a quick look at the MATLAB code on the linked page, and the model has some peculiarities. It uses a dt of 1ms, but with two 0.5ms update steps for v, and the synaptic input gets reset to 0 at every time step (so in a way it is a delta synapse implementation, but since the time step is very long one should probably rather think of it as a 1ms rectangular step current). There is no refractory period in that model.

Just in case you are not aware of it: some of Izhikevhich’s results have been a bit hard to replicate, and in particular the observed oscillations seem to depend on the details of the integration method and the integration time step. See e.g. the comment in the NEST documentation and this paper replicating the polychronization network.

2 Likes

yeah … I think my conclusion from comparing the brian code and original 2003 paper is that the original model has some odd features that might not be worth replicating.

( just as another source, if replicating Izhikevich / enumerating discrepancies is the way to go, here are some matlab files (and python wrappers) hosted through ModelDB )

1 Like

@mstimberg Thank you for this very polite note. I had tried to put it nicely, but after years of frustration with Izhikevich’s model, I decided to . You wrote the masterpiece!

1 Like

I ran original code in octave this morning, and it doesn’t always oscillate. When I reduced time-step twice, with appropriate scaling of noise amplitudes, oscillation were gone. I think these oscillations maybe a result of numerical instability, but I’m not sure.