Can we model the Izhikevich neural network with short-term plasticity?

Support
, ,
1

Description of problem

Can we model the Izhikevich neural network with short-term plasticity?
I tried to model the Izhikevich neural network with the STP synapse model. The program code I took from the HH model reference. The results of the programming are not as I expected, and the synapse model graph does not appear. I wonder if there is another breakthrough in the program code that I have made?.

Minimal code to reproduce problem

eqs = """dv/dt = (0.04*v**2 + 5*v + 140 - u + I + I_noise )/ms : 1
         du/dt = (a*(b*v - u))/ms  : 1
         I : 1
         I_noise : 1
         a : 1
         b : 1
         c : 1
         d : 1
       """

N = NeuronGroup(Ne + Ni, eqs, threshold="v>=30", reset="v=c; u+=d", method="euler")
N.v = -65

N_exc = N[:Ne]
N_inh = N[Ne:]

What you have aready tried

from brian2 import *
import numpy as np
%matplotlib inline

#Izhikevich Model
tfinal = 50 * ms
Ne = 40
Ni = 10

re = np.random.uniform(size=Ne)
ri = np.random.uniform(size=Ni)
weights = np.hstack(
    [
        0.5 * np.random.uniform(size=(Ne + Ni, Ne)),
        -np.random.uniform(size=(Ne + Ni, Ni)),
    ]
).T


start_scope()

#neuron Parameter
duration        = 1*second
sim_dt          = 0.1*ms
defaultclock.dt = sim_dt

V_t             = -30*mV       # Threshold membrane potential
V_r             = -60*mV       # Resting potential
taur            = 3*ms         # Time membrane refractory

#izhikevich model
a = 0.2
b = 0.02
c = -65
d = 8

#HH Model
area            = 20000*umetre**2
g_L             = (5e-5*siemens*cm**-2)*area
g_Ca            = (100*msiemens*cm**-2)*area
g_K             = (30*msiemens*cm**-2)*area
V_L             = -65*mV
C               = (1*ufarad*cm**-2)*area

# Synapse parameters
#from Stimberg et al (2017)
we       = 0.05*nS      # Exc syn conductance
wi       = 1.0*nS       # Inh syn conductance
U_0      = 0.6          # Resting syn release probability
Omegaf   = 3.33/second  # Syn facilitation rate
Omegad   = 2.0/second   # Syn depression rate
taue     = 5*ms         # Exc syn time constant
taui     = 10*ms        # Inh syn time constant
Ee       = 0*mV         # Exc syn reversal potential
Ei       = -80*mV       # Inh syn reversal potential


# izikevich Neuron model
N_eqs= """dv/dt = (g_L*((0.04/mV)*v**2) - g_Ca*(5*v) - g_K*(140*mV - u*ms) + I + 
                  ge*(Ee-v) + gi*(Ei-v))/C                     : volt
                du/dt = (a*(b*v - u))                          : volt/second
                I                                              : amp
                a                                              : Hz
                b                                              : Hz
                c                                              : volt
                d                                              : volt/second
                dge/dt   = -ge/taue                            : siemens
                dgi/dt   = -gi/taui                            : siemens
             """

N = NeuronGroup(Ne + Ni, N_eqs, threshold="v>=V_t", reset="v=c; u+=d", refractory="taur", method="euler")

N.v   = 'V_L + rand()*(V_t-V_L)'
N.ge  = 'rand()*we'
N.gi  = 'rand()*wi'

N_exc = N[:Ne]
N_inh = N[Ne:]

# STP Synapses
syn_eqs = '''
du_S/dt = -Omegaf * u_S     : 1 (event-driven)
dx_S/dt = Omegad * (1-x_S)  : 1 (event-driven)
'''

syn_act = '''
u_S += U_0 * (1-u_S)
r_S  = u_S * x_S
x_S -= r_S
'''

exc_syn = Synapses(N_exc, N, model=syn_eqs,
                   on_pre=syn_act+'ge_post += we*r_S')
inh_syn = Synapses(N_inh, N, model=syn_eqs,
                   on_pre=syn_act+'gi_post += wi*r_S')

exc_syn.connect(condition='i!=j', p=0.8)   #pe=0.8
inh_syn.connect(condition='i!=j', p=0.4)   #pi=0.4
exc_syn.x_S = 1
inh_syn.x_S = 1


#Monitor
#statemon
P = StateMonitor(N, ['v','ge','gi'], record=0, when='after_thresholds')
Q = StateMonitor(exc_syn, ['u_S', 'x_S'], record=exc_syn[Ne],
                when='after_synapses')
R = StateMonitor(inh_syn, ['u_S', 'x_S'], record=inh_syn[Ni],
                when='after_synapses')
#spikemon
S = SpikeMonitor(N)
T = SpikeMonitor(N_exc)
U = SpikeMonitor(N_inh)
N_exc.a = 0.02/ms
N_exc.b = 0.2/ms
N_exc.c = (-65 + 15)*volt * re**2
N_exc.d = (8 - 6)*volt/second * re**2

N_inh.a = (0.02 + 0.08)/ms * ri
N_inh.b = (0.25 - 0.05)/ms * ri
N_inh.c = -65*volt
N_inh.d = 2*volt/second

N_exc.u = "b*v"
N_inh.u = "b*v"

# recorded neuron
ne = 0 # excitatory neuron
ni = 0 # inhibitory neuron

# Run Simulation
run (duration, report='text')


# Plotting Neural Synchrony Activity
figure(figsize=(15, 30))
subplot(511)
plot(T.t/ms, T.i, '.r', label='exc_syn')
plot(U.t/ms, U.i + 40, '.k', label='inh_syn')
xlabel('t (ms)')
ylabel('neuron index')
legend()

# Plotting Post-Synaptic Conductance
figure(figsize=(15, 30))
subplot(512)
plot(P.t/ms, P[ne].ge/nS, 'r', label='ge')
plot(P.t/ms, -P[ne].gi/nS, 'k', label='gi')
xlabel('t (ms)')
ylabel('Conductance (nS)')
legend()


# Plotting display for exc_syn use Q,T,ne; inh_syn use R,U,ni
figure(figsize=(15, 30))
subplot(513)
# display for exc_syn use Q,T,ne; inh_syn use R,U,ni
spk_index = np.in1d(Q.t, T.t[T.i == ne])
plot(Q.t[spk_index]/ms, Q.x_S[0][spk_index], 'o',
            ms=4, color='C3')
plot(Q.t[spk_index]/ms, Q.u_S[0][spk_index], 'o',
            ms=4, color='C4')
time = Q.t
tspk = Quantity(Q.t, copy=True)
for ts in T.t[T.i == ni]:
    tspk[time >= ts] = ts
plot(Q.t/ms, 1+(Q.x_S[0] - 1)*exp(-(time-tspk)*Omegad),
    '-', color='C3')
plot(Q.t/ms, Q.u_S[0]*exp(-(time-tspk)*Omegaf),
    '-', color='C4')
xlabel('t (ms)')
ylabel('us, xs')
legend()

subplot(514)
nspikes = np.sum(spk_index)
x_S_spike = Q.x_S[0][spk_index]
u_S_spike = Q.u_S[0][spk_index]
vlines(Q.t[spk_index]/ms, np.zeros(nspikes),
             x_S_spike*u_S_spike/(1-u_S_spike))
xlim(-50, duration/ms + 50)
xlabel('t (ms)')
ylabel('rs')
legend()

subplot(515)
plot(P.t/ms, P[ne].v/mV, 'r', label='potential')
# line -30mV threshold
axhline(V_t/mV, ls='--', c='k')
# spike indeks
#for t in T.t[T.i == ne]:
# axvline(t/ms, ls='--', color='k')
legend()

xlabel('t (ms)')


# PLOTTING CONNECTIVITY dan SYNAPSE PATTERNS
def visualise_connectivity(exc_syn):
       Ns = len(exc_syn.source)
       Nt = len(exc_syn.target)
       figure(figsize=(20, 8))
       subplot(121)
       plot(zeros(Ns), arange(Ns), '.k', ms=10)
       plot(ones(Nt), arange(Nt), '.k', ms=10)
       for i, j in zip(exc_syn.i, exc_syn.j):
           plot([0, 1], [i, j], '-k')
       xticks([0, 1], ['Source', 'Target'])
       ylabel('Neuron index')
       xlim(-0.1, 1.1)
       ylim(-1, max(Ns, Nt))
       subplot(122)
       plot(exc_syn.i, exc_syn.j, 'or')
       xlim(-1, Ns)
       ylim(-1, Nt)
       xlabel('Source neuron index')
       ylabel('Target neuron index')
visualise_connectivity(exc_syn)

def visualise_connectivity(inh_syn):
       Ns = len(inh_syn.source)
       Nt = len(inh_syn.target)
       figure(figsize=(20, 8))
       subplot(121)
       plot(zeros(Ns), arange(Ns), '.k', ms=10)
       plot(ones(Nt), arange(Nt), '.k', ms=10)
       for i, j in zip(inh_syn.i, inh_syn.j):
           plot([0, 1], [i, j], '-k')
       xticks([0, 1], ['Source', 'Target'])
       ylabel('Neuron index')
       xlim(-0.1, 1.1)
       ylim(-1, max(Ns, Nt))
       subplot(122)
       plot(inh_syn.i, inh_syn.j, 'ok')
       xlim(-1, Ns)
       ylim(-1, Nt)
       xlabel('Source neuron index')
       ylabel('Target neuron index')
visualise_connectivity(inh_syn)

print(T.t[T.i == ne])

show()

Expected output (if relevant)

referring to the reference HH…

download (10)
download (10)940×420 17.2 KB

Actual output (if relevant)

Full traceback of error (if relevant)

Starting simulation at t=0. s for a duration of 1. s
WARNING    (string):33: RuntimeWarning: overflow encountered in square
 [py.warnings]
WARNING    (string):32: RuntimeWarning: invalid value encountered in multiply
 [py.warnings]
WARNING    (string):33: RuntimeWarning: invalid value encountered in add
 [py.warnings]
1. s (100%) simulated in 3s
No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
[0.] s
2

Hi @mufa17. There is no specific issue with STP and the Izhikevich model, i.e. it should work in principle. The issues you are seeing do not seem to be related to STP or even the synapses at all – e.g. if you set all your synaptic weights to zero you still get NaN values in your simulation. It looks like you inserted HH conductances into the Izhikevich model, and also set the usual a, b, etc. parameters but do not use them in your equations ? This cannot work, so you’ll need to fix the neuron model first.

Two minor general Python notes (not directly related to Brian):

  • In your plotting code, if you call figure before calling subplot, each subplot will end up in an independent figure, which is most likely not what you want.
  • You do not need to repeat the visualize_connectivity function (otherwise you wouldn’t need a function). You can define it once with a general syn parameter, and refer to it in the function body. You can then call it twice, once as visualise_connectivity(exc_syn) and once as visualise_connectivity(inh_syn).
1 Like
3

Thank you, Mr Stimberg.
but I am still confused, how to combine the izhikevich neuron model with short term plasticity?

4

Hi @mufa17, as I said, there is no problem with using STP and the Izhikevich equations – you can combine e.g. the STP equations from this example: Example: example_1_COBA — Brian 2 2.6.0 documentation and connect them to a synaptic variable in your Izhikevich neuron model (i.e. to ge or gi). But before you worry about the synapses, you need to get your neuron model to work. Your current model mixes HH and Izhikevich equations in a way that will not work. You need to use either of the two models.