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 - 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)

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