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 :upside_down_face: