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.8N):]
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 = yw_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.
, thank you so much for your complete answer
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.
. You wrote the masterpiece! 