Implementing Izhikevich 2006 in Brian2

Support
1

Description of problem

Hi,

I’m new to Brian (and neuroscience in general) and I have been trying to implement the model from Izhikevich 2006 in Brian2. I have been using this Brian1 code as reference. However, I am struggling to implement the STDP rule and the numerical integration scheme. The STDP rule in Brian1 uses functions that I could not find in Brian2. My state updater also differs from the Brian1 implementation. The main difference in output is that the raster plot from the original paper has delta rhythms which mine lacks.

Minimal code to reproduce problem

start_scope()

defaultclock.dt = 0.5*ms
I_factor = float((1*ms)/dt)

M = 100
D = 20*ms
N_e = 800
a_e = 0.02/ms
d_e = 8*mV/ms
s_e = 6*mV/ms * I_factor

N_i = 200
a_i = 0.1 / ms
d_i = 2 * mV / ms
s_i = -5 * mV / ms * I_factor

b = 0.2/ms
c = -65*mV
sm = 10 * mV / ms * I_factor  

eqs = '''
dv/dt = (0.04/ms/mV)*v**2+(5/ms)*v+140*mV/ms-u+I : volt
du/dt = a*(b*v-u)                                : volt/second
a                                                : 1/second
d                                                : volt/second
I                                                : volt/second
'''

reset = '''
v = c
u += d
'''

new_state_updater = ExplicitStateUpdater('''
    x_support = x + 0.5 * dt *f(x,t)
    x_new = x_support + 0.5 * dt *f(x,t)
    ''')
StateUpdateMethod.register('mymethod', new_state_updater)

G = NeuronGroup(N_e + N_i, model = eqs, threshold='v>30*mV', reset=reset, method='mymethod')
''' Scheme from Brian1
def Izhikevich_scheme(G):
    G.v += 0.5 * dt * ((0.04 / ms / mV) * G.v ** 2 + (5 / ms) * G.v + 140 * mV / ms - G.u + G.I)
    G.v += 0.5 * dt * ((0.04 / ms / mV) * G.v ** 2 + (5 / ms) * G.v + 140 * mV / ms - G.u + G.I)
    G.u += dt * G.a * (b * G.v - G.u)
G._state_updater = Izhikevich_scheme '''


Ge = G[:N_e]
Gi = G[N_e:]
G.v = c
G.u = b * c
Ge.a = a_e
Gi.a = a_i
Ge.d = d_e
Gi.d = d_i

thalamic_input_weight = 20*mV/ms
change_every = int(1*ms/dt)
change_time = 0
thal_ind = 0
@network_operation(dt=1 * ms)
def thalamic_input():
    global change_time, thal_ind
    G.I = 0 * mV / ms
    G.I[thal_ind] = thalamic_input_weight
    if change_time==0:
        change_time = change_every
        thal_ind = randint(N)
        #G.I[randint(N)] = float(thalamic_input_weight*1*ms/dt)
    change_time -= 1


artificial_wmax = 1e10 * mV / ms
taupre = taupost = 20*ms
wmax = artificial_wmax
Apre = (0.1 * mV / ms / artificial_wmax)
Apost = (-0.12 * mV / ms / artificial_wmax)
sm = 10 * mV / ms * I_factor

synapsesE = Synapses(Ge, G, '''
             w_e : volt/second
             dapre/dt = -apre/taupre : volt/second (event-driven)
             dapost/dt = -apost/taupost : volt/second (event-driven)
             ''',
             on_pre='''
             I += w_e
             apre += Apre
             w_e = clip(w_e+apost, 0*mV/ms, wmax)
             ''',
             on_post='''
             apost += Apost
             w_e = clip(w_e+apre, 0*mV/ms, wmax)
             ''')

# Connection matrix based on matlab code from original paper
# Matrix of postsynaptic indices, 100 random connections per neuron
C_e = np.array([np.random.choice(1000, 100, replace=False) for _ in range(800)])
M_e = np.zeros((N_e, N))
#Converting into matrix of 1s and 0s (1 for connection)
for x, y in np.ndenumerate(C_e):
  M_e[x[0], int(y)] = 1

sourcesE, targetsE = M_e.nonzero()
synapsesE.connect(i=sourcesE, j=targetsE)
synapsesE.delay = 'ceil(j/5)*ms' # Delay based on matlab code from original paper
synapsesE.w_e = s_e

@network_operation(dt=1 * ms)
def update_weights_from_derivative():
    synapsesE.w_e = clip(synapsesE.w_e + 0.01 * mV / ms + synapsesE.w_e, 0, sm)


C_i = np.array([np.random.choice(800, 100, replace=False) for _ in range(200)])
M_i = np.zeros((N_i, N))

for x, y in np.ndenumerate(C_i):
  M_i[x[0], int(y)] = 1
sourcesI, targetsI = M_i.nonzero()
synapsesI = Synapses(Gi, Ge, model = 'w_i : volt/second', on_pre='I+=w_i', delay=1*ms)
synapsesI.connect(i=sourcesI, j=targetsI)
synapsesI.w_i = s_i

spikemon = SpikeMonitor(G)


run(1*second)

plot(spikemon.t/ms, spikemon.i, '.k')
xlim([0, 1000*1])
ylim([0, 1000])
xlabel('Time (ms)')
ylabel('Neuron index')

This is parts of the STDP from Brian1


# In order to implement Izhikevich's STDP rule, we need two weight matrices,
# one is the actual weight matrix and the second is a weight matrix derivative.
# To do this with Brian's sparse matrices, we create a second Connection
# Ce_deriv, which we initialise to have the same shape as Ce, but we use the
# forget() function on the Connection so that although the object exists,
# no spikes will be propagated by that Connection. We need the object to exist
# because we create an ExponentialSTDP object that acts on Ce_deriv not
# directly on Ce.

Ce_deriv = Connection(Ge, G, 'I', delay=True, max_delay=D)
forget(Ce_deriv)

# Now we duplicate Ce to Ce_deriv

for i in xrange(Ne):
    Ce_deriv[i, :] = Ce[i, :]
    Ce_deriv.delay[i, :] = Ce.delay[i, :]

# STDP acts directly on Ce_deriv rather than Ce. With this STDP rule alone,
# we would not see any learning, the network operation below propagates changes
# in Ce_deriv to Ce.

artificial_wmax = 1e10 * mV / ms
stdp = ExponentialSTDP(Ce_deriv,
                       taup=20 * ms, taum=20 * ms,
                       Ap=(0.1 * mV / ms / artificial_wmax),
                       Am=(-0.12 * mV / ms / artificial_wmax),
                       wmin= -artificial_wmax,
                       wmax=artificial_wmax,
                       interactions='nearest'
                       )

# Izhikevich's STDP rule has STDP acting on a matrix sd of derivatives, and
# then each second the weight matrix s and derivates sd are updated according
# to the rule:
#   s <- s+0.01+sd
#   sd <- sd*0.9
# Note also that we are using Brian's sparse matrices, but they are have
# exactly the same nonzero entries, and therefore we can do arithmetic
# operations on these matrices using the alldata attribute of Brian sparse
# matrices. The compress() method converts the ConstructionMatrix into a
# ConnectionMatrix, thus freezing the nonzero entries. In the next line, we
# want the actual starting values of Ce_deriv to be zero, but if we had done
# this before compressing to a ConnectionMatrix, all entries would be
# considered zero entries in the sparse matrix, and then Ce and Ce_deriv would
# have a different pattern of non-zero entries.

Ce_deriv.compress()
Ce_deriv.W.alldata[:] = 0
@network_operation(clock=EventClock(dt=1 * second))
def update_weights_from_derivative():
    Ce.W.alldata[:] = clip(Ce.W.alldata + 0.01 * mV / ms + Ce_deriv.W.alldata, 0, sm)
    Ce_deriv.W.alldata[:] *= 0.9
1 Like
2

Hi @rrm . As I mentioned in a comment in another post (Realization the model of network in Brain 2 - #16 by mstimberg), several of the Izhikevich papers from the early 2000s use somewhat sloppy numerics (like the state update rule in which u is updated based on an already updated v). Unfortunately, it seems some of the results can only be reproduced by using the same numerics, which questions their general validity. For example, reducing the simulation time step should give qualitatively the same results, but this does not seem to be the case here. This paper by Pauli et al. analyzes this model in impressive detail, it is well worth a read if you are trying to reproduce it.

Regarding the Brian 1 code: it uses quite a lot of Brian 1 internals that no longer apply to Brian 2, but from a cursory look I’d say that you translated most of it correctly. Due to the update of u with the already updated value of v, you cannot use the standard StateUpdater mechanism, though, unless you use a complicated workaround with two independent NeuronGroups and linked variables. Instead, I’d suggest a more “manual” approach that exactly reproduces the approach in the paper: instead of defining v and u as differential equations, you define them as simple parameters, and update them manually. I.e., in the equations, they are defined as v : volt and u : volt/second. The actual state update then happens with a run_regularly function:

G.run_regularly('''
v += 0.5 * dt * ((0.04 / ms / mV) * v ** 2 + (5 / ms) * v + 140 * mV / ms - u + I)
v += 0.5 * dt * ((0.04 / ms / mV) * v ** 2 + (5 / ms) * v + 140 * mV / ms - u + I)
u += dt * a * (b * v - u)''', when='groups')

I am not 100% sure I correctly understood the use of the derivative of the weight matrix here, but I don’t think your code is doing this correctly. What I think it should do is something like this:

synapsesE = Synapses(Ge, G, '''
             w_e : volt/second  # the weight
             w_e_change : volt/second  # the "derivative
             dapre/dt = -apre/taupre : volt/second (event-driven)
             dapost/dt = -apost/taupost : volt/second (event-driven)
             ''',
             on_pre='''
             I += w_e
             apre += Apre
             w_e_change = clip(w_e_change+apost, 0*mV/ms, wmax)
             ''',
             on_post='''
             apost += Apost
             w_e_change = clip(w_e_change+apre, 0*mV/ms, wmax)
             ''')
# Every s, apply the accumulated changes to the actual weights and decay the "derivative"
synapsesE.run_regularly('''w_e = clip(w_e + 0.01*mV/ms + w_e_change, 0*mV/ms, sm)
                           w_e_change = w_e_change*0.9''',
                           dt=1*second)

The idea is that w_e is used to propagate spikes (added to I), but the STDP acts on w_e_change. Only once every second (not ms as in your code), the changes are applied to the actual weights w_e.
You shouldn’t need a network_operation here.

It’s a minor point, but since version 2.5, there’s actually a simple way to create 100 random connections per neuron. Instead of manually creating matrices, you can do:

synapsesE.connect(j='index for index in sample(N_post, 100)')

Best,
Marcel

4 Likes
3

Hi @mstimberg ,

Thank you for your help. I definitely would not have been able to figure out the implementation for the integration and STDP on my own. After implementing your suggestions, my output was still different than the plots shown in the paper. I realized that this might be because my implementation for the excitatory delays was not behaving as I expected.

In the excitatory synapse, Izhikevich assigns 5 neurons a 1 ms delay, 5 neurons a 2 ms delay, etc. up until 20 ms. From my understanding, he does this for each pre-synaptic neuron. I tried using a randomly distributed delay between 1ms-20ms, as the Pauli paper did, and that gave a more reasonable output. I also tried implementing Izhikevich’s delays by manually creating the delays array. I am not sure how the array works internally, but here’s the code I wrote:

delays = []
for i in range(800):
  for j in range(20):
     delays += [j+1 for _ in range(5)]

synapsesE.delay = delays*ms

The output in this case is still off. I attached the raster plot Brian generated and the one from the original paper. This weird spiking pattern also persists after running for 100 seconds whereas Izhikevich’s starts showing different activity.
Izhikevich delays
Capture

I am also confused about the use of the I_factor term in the Brian1 code. What purpose does this constant serve? I noticed that it causes the maximal synaptic strength to be 20 instead of 10 as in the paper.

Any guidance is appreciated. Thanks in advance.

4

Hi again. I am certainly not an expert on this model and as you can see from the Pauli et al. paper, reproducing the results is far from trivial.
Regarding the delays: random delays make perfect sense, but if you want to use Izhikevich’s systematic approach, you can use the following:

synapsesE.delay = np.tile(np.arange(20).repeat(5) + 1, N_E)*ms

This means, [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, …, 20, 20, 20, 20, 20]*ms is used for the 100 synapses of each pre-synaptic neurons. This works because synapses are ordered by pre-synaptic index first:

>>> synapsesE.i
<synapses.i: array([  0,   0,   0, ..., 799, 799, 799], dtype=int32)>

You can check that delays are correctly assigned:

>>> synapsesE.delay[0, :]  # delays for outgoing synapses of neuron 0
array([ 1.,  1.,  1.,  1.,  1.,  2.,  2.,  2.,  2.,  2.,  3.,  3.,  3.,
        3.,  3.,  4.,  4.,  4.,  4.,  4.,  5.,  5.,  5.,  5.,  5.,  6.,
        6.,  6.,  6.,  6.,  7.,  7.,  7.,  7.,  7.,  8.,  8.,  8.,  8.,
        8.,  9.,  9.,  9.,  9.,  9., 10., 10., 10., 10., 10., 11., 11.,
       11., 11., 11., 12., 12., 12., 12., 12., 13., 13., 13., 13., 13.,
       14., 14., 14., 14., 14., 15., 15., 15., 15., 15., 16., 16., 16.,
       16., 16., 17., 17., 17., 17., 17., 18., 18., 18., 18., 18., 19.,
       19., 19., 19., 19., 20., 20., 20., 20., 20.]) * msecond

I think this is there because the Brian 1 simulation uses a dt of 0.5, but the original code uses a dt of 1ms. With dt=0.5ms, the current I is applied for half the time, so it needs to be twice as strong. But when you use the run_regularly formulation for the update of v and w, you should use dt=1ms as in the original code (v is updated twice in run_regularly, so as if it had a dt of 0.5ms) – this means that the I_factor becomes irrelevant.

1 Like