Realization the model of network in Brain 2

[serge25]

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

[rth]

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…

[rth]

@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, 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

[Maryidea]

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

[rth]

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.

[rth]

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.

[Maryidea]

ohh yes , thank you so much for your complete answer

[serge25]

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

[rth]

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?

[rth]

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'

[serge25]

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

[rth]

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.

[adam]

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

[mstimberg]

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.

[adam]

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 )

[rth]

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

[rth]

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.

[serge25]

@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?

[adam]

For a different thread on Izhikevich neurons I worked through an example script with adds units to the (sort-of) dimensionless way the original equations are specified.

In doing so, I added some useful features to the script that may be worth taking a look at, including a dictionary of neuron types:

# a,b,c,d define the type of neuron
# --------------------------------------     a      b       c     d
abcd_neurons = {'regular spiking':       [0.02/ms, 0.2,  -65*mV, 8*mV],
                'intrinsically bursting':[0.02/ms, 0.2,  -55*mV, 4*mV],
                'chattering':            [0.02/ms, 0.2,  -50*mV, 2*mV],
                
                'fast spiking':          [0.10/ms, 0.2,  -65*mV, 2*mV],
                'low-threshold spiking': [0.02/ms, 0.25, -65*mV, 2*mV],
                
                'thalamo-cortical':      [0.02/ms, 0.25, -65*mV, 0.05*mV],
                'resonator':             [0.10/ms, 0.25, -65*mV, 8*mV],
                }
a,b,c,d = abcd_neurons[neuron_type]

of particular interest to @serge25, I was able to get intrinsic bursting:

as well as thalamo-cortical post-inhibitory rebound (bursting)

see the script here
@mstimberg is this worth turning into an example for the documentation?

[mstimberg]

@adam

@mstimberg is this worth turning into an example for the documentation?

Would be great to have this in the example section! I wonder whether it would make sense to have the constants a etc. instead as constant parameters in the equations (i.e. a : 1/second (constant) etc.). This way, one could do a single simulation for all classes and plot them all at once. On the other hand, this would make the code a bit more complex to read, I guess. Having to choose one parametrization in the script would be perfectly fine as well, after all, e.g. Example: Brette_Gerstner_2005 — Brian 2 2.5.0.2 documentation does this as well.