# Realization the model of network in Brain 2

**URL:** <https://brian.discourse.group/t/realization-the-model-of-network-in-brain-2/444>\
**Category:** Science, Projects, and Showcases\
**Created:** [5 August 2021 12:30 UTC](https://brian.discourse.group/t/realization-the-model-of-network-in-brain-2/444 "2021-08-05T12:30:57Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![serge25](https://avatars.discourse-cdn.com/v4/letter/s/4af34b/32.png) [@serge25](https://brian.discourse.group/u/serge25)\
**Post date:** [5 August 2021 12:30 UTC](https://brian.discourse.group/t/realization-the-model-of-network-in-brain-2/444/1 "2021-08-05T12:30:58Z")

</div>

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 :

````auto
#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 - 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)

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

---

<div class="post-metadata">

**Author:** ![rth](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/rth/32/16_2.png) [@rth](https://brian.discourse.group/u/rth)\
**Post date:** [9 August 2021 14:18 UTC](https://brian.discourse.group/t/realization-the-model-of-network-in-brain-2/444/2 "2021-08-09T14:18:53Z")

</div>

Hi [@serge25](https://brian.discourse.group/u/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](https://brian2.readthedocs.io/en/stable/user/synapses.html#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:

```python
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

```python
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_

---

<div class="post-metadata">

**Author:** ![serge25](https://avatars.discourse-cdn.com/v4/letter/s/4af34b/32.png) [@serge25](https://brian.discourse.group/u/serge25)\
**Post date:** [9 August 2021 19:47 UTC](https://brian.discourse.group/t/realization-the-model-of-network-in-brain-2/444/3 "2021-08-09T19:47:38Z")

</div>

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

```auto
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

```auto
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:

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

```

---

<div class="post-metadata">

**Author:** ![rth](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/rth/32/16_2.png) [@rth](https://brian.discourse.group/u/rth)\
**Post date:** [9 August 2021 20:01 UTC](https://brian.discourse.group/t/realization-the-model-of-network-in-brain-2/444/4 "2021-08-09T20:01:07Z")

</div>

Here the code which actually works.

```python
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…

 ![brian2-forum](https://global.discourse-cdn.com/free1/uploads/brian/original/1X/74a084a0df077cfab03f48a37489c00a147e86c7.png)

---

<div class="post-metadata">

**Author:** ![rth](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/rth/32/16_2.png) [@rth](https://brian.discourse.group/u/rth)\
**Post date:** [9 August 2021 20:19 UTC](https://brian.discourse.group/t/realization-the-model-of-network-in-brain-2/444/5 "2021-08-09T20:19:29Z")

</div>

> [@serge25](#):
>
> 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)

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](https://brian2.readthedocs.io/en/stable/user/synapses.html), 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 🙃

---

<div class="post-metadata">

**Author:** ![Maryidea](https://avatars.discourse-cdn.com/v4/letter/m/c0e974/32.png) [@Maryidea](https://brian.discourse.group/u/Maryidea)\
**Post date:** [13 November 2021 16:54 UTC](https://brian.discourse.group/t/realization-the-model-of-network-in-brain-2/444/6 "2021-11-13T16:54:35Z")

</div>

Hi  
sorry , i have one problem in the neuron equation , here ‘[Simple Model of Spiking Neurons](https://www.izhikevich.org/publications/spikes.htm)’ 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

---

<div class="post-metadata">

**Author:** ![rth](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/rth/32/16_2.png) [@rth](https://brian.discourse.group/u/rth)\
**Post date:** [14 November 2021 10:25 UTC](https://brian.discourse.group/t/realization-the-model-of-network-in-brain-2/444/7 "2021-11-14T10:25:47Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![rth](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/rth/32/16_2.png) [@rth](https://brian.discourse.group/u/rth)\
**Post date:** [14 November 2021 13:08 UTC](https://brian.discourse.group/t/realization-the-model-of-network-in-brain-2/444/8 "2021-11-14T13:08:08Z")

</div>

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

```python
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

```python
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.

---

<div class="post-metadata">

**Author:** ![Maryidea](https://avatars.discourse-cdn.com/v4/letter/m/c0e974/32.png) [@Maryidea](https://brian.discourse.group/u/Maryidea)\
**Post date:** [14 November 2021 15:40 UTC](https://brian.discourse.group/t/realization-the-model-of-network-in-brain-2/444/9 "2021-11-14T15:40:13Z")

</div>

ohh yes 👍🏻 , thank you so much for your complete answer

---

<div class="post-metadata">

**Author:** ![serge25](https://avatars.discourse-cdn.com/v4/letter/s/4af34b/32.png) [@serge25](https://brian.discourse.group/u/serge25)\
**Post date:** [12 January 2022 19:37 UTC](https://brian.discourse.group/t/realization-the-model-of-network-in-brain-2/444/10 "2022-01-12T19:37:31Z")

</div>

Hi @rth ,

Can you help me with one question. I’m trying to implement Izhikevich network from [Simple Model of Spiking Neurons](https://www.izhikevich.org/publications/spikes.htm). I can’t get the results expected in the article. What’s wrong with my code?

```auto
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()

```

---

<div class="post-metadata">

**Author:** ![rth](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/rth/32/16_2.png) [@rth](https://brian.discourse.group/u/rth)\
**Post date:** [12 January 2022 22:29 UTC](https://brian.discourse.group/t/realization-the-model-of-network-in-brain-2/444/11 "2022-01-12T22:29:21Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![rth](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/rth/32/16_2.png) [@rth](https://brian.discourse.group/u/rth)\
**Post date:** [12 January 2022 23:04 UTC](https://brian.discourse.group/t/realization-the-model-of-network-in-brain-2/444/12 "2022-01-12T23:04:53Z")

</div>

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:

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

```

---

<div class="post-metadata">

**Author:** ![serge25](https://avatars.discourse-cdn.com/v4/letter/s/4af34b/32.png) [@serge25](https://brian.discourse.group/u/serge25)\
**Post date:** [13 January 2022 06:03 UTC](https://brian.discourse.group/t/realization-the-model-of-network-in-brain-2/444/13 "2022-01-13T06:03:27Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![rth](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/rth/32/16_2.png) [@rth](https://brian.discourse.group/u/rth)\
**Post date:** [13 January 2022 15:24 UTC](https://brian.discourse.group/t/realization-the-model-of-network-in-brain-2/444/14 "2022-01-13T15:24:30Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![adam](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/adam/32/197_2.png) [@adam](https://brian.discourse.group/u/adam)\
**Post date:** [13 January 2022 16:38 UTC](https://brian.discourse.group/t/realization-the-model-of-network-in-brain-2/444/15 "2022-01-13T16:38:37Z")

</div>

![Screen Shot 2022-01-13 at 11.01.10 AM](https://global.discourse-cdn.com/free1/uploads/brian/original/1X/0457a33de751f837dbefca62b9c8e824388fbde5.png)  
 ![Screen Shot 2022-01-13 at 11.01.24 AM](https://global.discourse-cdn.com/free1/uploads/brian/original/1X/88d047325ec244f04e12126646f875efc4af5299.png)  
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:  
 ![Screen Shot 2022-01-13 at 11.03.07 AM](https://global.discourse-cdn.com/free1/uploads/brian/original/1X/0e67286275835f9f8fd76a7476a8de9baff72982.png)

## 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**
>
> ```auto
> 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:

```auto
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:

```python
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

---

<div class="post-metadata">

**Author:** ![mstimberg](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mstimberg/32/11_2.png) [@mstimberg](https://brian.discourse.group/u/mstimberg)\
**Post date:** [13 January 2022 16:47 UTC](https://brian.discourse.group/t/realization-the-model-of-network-in-brain-2/444/16 "2022-01-13T16:47:37Z")

</div>

I just had a quick look at the [MATLAB code](https://www.izhikevich.org/publications/net.m) 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](https://nest-simulator.readthedocs.io/en/latest/models/izhikevich.html?#description) and [this paper replicating the polychronization network](https://www.frontiersin.org/articles/10.3389/fninf.2018.00046/full).

---

<div class="post-metadata">

**Author:** ![adam](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/adam/32/197_2.png) [@adam](https://brian.discourse.group/u/adam)\
**Post date:** [13 January 2022 16:55 UTC](https://brian.discourse.group/t/realization-the-model-of-network-in-brain-2/444/17 "2022-01-13T16:55:04Z")

</div>

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)](https://senselab.med.yale.edu/ModelDB/ShowModel?model=39948&file=/izhikevich_simple_nrn/izhi2007Comparison.pdf#tabs-2) hosted through ModelDB )

---

<div class="post-metadata">

**Author:** ![rth](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/rth/32/16_2.png) [@rth](https://brian.discourse.group/u/rth)\
**Post date:** [13 January 2022 20:41 UTC](https://brian.discourse.group/t/realization-the-model-of-network-in-brain-2/444/18 "2022-01-13T20:41:44Z")

</div>

@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! 👏

---

<div class="post-metadata">

**Author:** ![rth](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/rth/32/16_2.png) [@rth](https://brian.discourse.group/u/rth)\
**Post date:** [13 January 2022 20:47 UTC](https://brian.discourse.group/t/realization-the-model-of-network-in-brain-2/444/19 "2022-01-13T20:47:10Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![serge25](https://avatars.discourse-cdn.com/v4/letter/s/4af34b/32.png) [@serge25](https://brian.discourse.group/u/serge25)\
**Post date:** [15 January 2022 01:36 UTC](https://brian.discourse.group/t/realization-the-model-of-network-in-brain-2/444/20 "2022-01-15T01:36:53Z")

</div>

@rth @adam @mstimberg , thank you for your attention.  
With some tricks, I got bursting activity at this model. Also, I have one question about not this model, but Brian2 work. I right understand from the documentation that in the case if I have two variables which one of them dependent on the other through the threshold (not neuron spike), trigger event I can put only in NeuronGroup (custom threshold, like X\>Xthr), not a SynapseGroup. At SynapseGroup I can only update synaptic variables through a new path?

[Next page](https://brian.discourse.group/t/realization-the-model-of-network-in-brain-2/444.md?page=2)