# 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…

# 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
```