Run error with modeling izhikevich neuron

[Maryidea]

Hi everybody , i try to model izhikevich neuron with brian2 but i get this error and don’t know how fix it

from brian2 import*

%matplotlib inline

duration = 8*second

#neurons parameter

a=0.02/msecond

b = 0.2/msecond

c = -65 * mV

d = 6

v_th = 30 * mV

v_0 = -60 * mV

u_0 = 0 * mV/ms

I_app = 12*uamp

eqs = (’’’ dv/dt = ( 0.04/ms/mV) * v**2 +(5/ms) * v + 140*mV/ms - u + I_app : volt

                 du/dt = a * ( b * v - u ) : volt/second

                 I_app : volt/second

                

            ''')

start_scope()

G = NeuronGroup(1, model = eqs , threshold = ’ v > v_th ’ , reset = ’ v = c ; u = u + d ’ , method = ‘euler’ )

start_scope()

G = NeuronGroup(1, eqs, method=‘exact’)

M = StateMonitor(G, variables=[‘v’,‘I_app’], record=True,dt=1*ms)

run(10*second)

plot(M.t/ms, M.v[0])

xlabel(‘Time (ms)’)

ylabel(‘v’)

##############################################################################

---

UnsupportedEquationsException Traceback (most recent call last)
/usr/local/lib/python3.7/dist-packages/brian2/core/network.py in before_run(self, run_namespace)
891 try:
→ 892 obj.before_run(run_namespace)
893 except Exception as ex:

16 frames
UnsupportedEquationsException: The expression ‘I_app + 140mV/ms - u + 5v/ms + 0.04v**2/(mVms)’, defining the variable ‘v’, could not be separated into linear components.

The above exception was the direct cause of the following exception:

BrianObjectException Traceback (most recent call last)
/usr/local/lib/python3.7/dist-packages/brian2/core/network.py in before_run(self, run_namespace)
892 obj.before_run(run_namespace)
893 except Exception as ex:
→ 894 raise BrianObjectException(“An error occurred when preparing an object.”, obj) from ex
895
896 # Check that no object has been run as part of another network before

BrianObjectException: Error encountered with object named ‘neurongroup_4_stateupdater’.
Object was created here (most recent call only, full details in debug log):
File ‘’, line 27, in
G = NeuronGroup(1, eqs, method=‘exact’)

An error occurred when preparing an object. (See above for original error message and traceback.)

[adam]

Hi Mary,
I think there a couple of things going on here.
The error message you’re getting seems to be directly related to specifying the group G to use exact numerical integration - which it can’t do because the equations have a nonlinear term:

\frac{dv}{dt} = 0.04 v^2 ...

This is a relatively quick fix, you can just replace that line with

G = NeuronGroup(1, eqs, method='euler')

but you’ve also defined G twice and started the scope twice which means maybe you meant to combine those:

start_scope()
G = NeuronGroup(1, model = eqs , threshold = 'v > v_th' , reset = 'v = c ; u = u + d' , method = 'euler' )

and just delete the 2 lines after.

---

After fixing those, the simulation runs, but appears to go unstable still. I suspect this has to do with the units of the constants. I’m looking into that now

[Maryidea]

@Maryidea

M = StateMonitor(G, variables=[‘v’,‘I_app’], record=True,dt=1*ms)

run(10*second)

plot(M.t/ms, M.v[0])

xlabel(‘Time (ms)’)

ylabel(‘v’)

thank you adam , i changed all you said but i still get this error
“MagicError: The magic network contains a mix of objects that has been run before and new objects, Brian does not know whether you want to start a new simulation or continue an old one. Consider explicitly creating a Network object. Also note that you can find out which objects will be included in a magic network with the collect() function.”

[adam]

So I worked on getting the dimensionality of the constants correct, and here’s what I came up with (partially based on some notes by Naureen Ghani )

"dimensionalized" izhikevich equations coded in Brian (click to expand script)

from brian2 import*
%matplotlib inline

import random
import matplotlib.pyplot as plt
#%%

start_scope()
# SIMULATION PARAMETERS
duration = 0.25*second
# defaultclock.dt = 0.001*ms

# Currents should be specified in amps...
# - if you do that, change C5 to be a resistance, in ohm
# - also change > I_app : volt → I_app : amp  
I_off = 0*mV
I_on = 10*mV #*uamp
I_app_start = 20*ms #150*ms
neuron_type = 'chattering'

'''
This script attempts to "dimensionalize" or add units to the basic Izhikevich model neuron.
In the original model "Simple Model of Spiking Neurons" (Izhikevich, 2003) 
 equations are presented in a (sort-of) dimesnionless form.
https://www.izhikevich.org/publications/spikes.pdf

Some lecture notes by Naureen Ghani add dimensions to the relevant constants, 
although with some inconsistencies. I used this as a starting point 
http://www.columbia.edu/cu/appliedneuroshp/Spring2018/Spring18SHPAppliedNeuroLec5.pdf

- by Adam Willats, Jan. 2021
'''
#%%

# 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],
                }
if neuron_type == 'random':
    neuron_type = random.choice(list(abcd_neurons.keys())) #uncomment to get a random neuron type
a,b,c,d = abcd_neurons[neuron_type]

# NEURON CONSTANTS
v_th = 30 * mV

C1 = 0.04 * 1/mV
C2 = 5.0
C3 = 140 * mV
C5 = 1.0 #* kohm # resistance to injected current

eqs = '''
 dv/dt =   (C1*(v**2) + C2*v + C3 - u + C5*I_app) / ms: volt
 du/dt = a * ( b * v - u )  : volt 
 I_app = I_app_fn(t) : volt
'''

v_0 = c
u_0 = 0 * mV
# %%
# SET UP NETWORK, RUN SIMULATION
# create a function to simulate a step-change in current 
I_app_fn = TimedArray( [I_off, I_on], dt=I_app_start)

G = NeuronGroup(1, model = eqs , threshold = 'v > v_th' , reset = 'v = c ; u = u + d' , method = 'euler' )
G.v = v_0
G.u = u_0
M = StateMonitor(G, variables=['v','I_app','u'], record=True)
SpikeMon = SpikeMonitor(G)
run(duration)
#%%
# PLOTTING THE RESULTS
pair = np.ones((2,1)) 
fake_spike = [-45,60]
f,axs = subplots(2,1,figsize=(8,6), gridspec_kw={'height_ratios': [3, 1]})

axs[0].plot(M.t/ms, M.v[0]/mV,'k') #plot voltage

# for st, si in zip(SpikeMon.t,SpikeMon.t):
    # plot vertical lines for spikes
    # axs[0].plot(st/ms*pair, fake_spike, 'k-')
    # pass
# axs[0].plot([0,duration/ms],[v_th/mV,v_th/mV],'r:') #plot voltage threshold
axs[1].plot(M.t/ms, M[0].I_app/mV,'r')

# cleaning up plot formatting
axs[0].set_ylabel('voltage [mV]')
axs[0].set_title(neuron_type,fontsize=20)
axs[0].set_xticklabels([])
axs[1].set_ylabel('external stim. (mV) ')
axs[1].set_xlabel('Time (ms)')
axs[0].spines['right'].set_visible(False)
axs[0].spines['top'].set_visible(False)
axs[1].spines['right'].set_visible(False)
axs[1].spines['top'].set_visible(False)
plt.savefig(f'izhi_{neuron_type}.png')

some highlights:

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]

C1 = 0.04 * 1/mV
C2 = 5.0
C3 = 140 * mV
C5 = 1.0

eqs = '''
 dv/dt =   (C1*(v**2) + C2*v + C3 - u + C5*I_app) / ms: volt
 du/dt = a * ( b * v - u )  : volt 
 I_app = I_app_fn(t) : volt
'''

some example plots:

[Maryidea]

thank you adam ,it’s so usefull , thank you so much for your help

[Maryidea]

Hi adam , sorry for my many questions i ask , according to the model you siad and post here , when i want to set a synapse_eqs between three neurons and then plot it , again i face with the run error !!
how can i solve it ?

[adam]

Hi Mary,
which run error did you get?

[Maryidea]

thank you for your attention .

i seprated three izhikevich neurons in two groups :

from brian2 import*

import numpy

import numpy as np

%matplotlib inline

import random

import matplotlib.pyplot as plt

#%%

start_scope()

SIMULATION PARAMETERS

duration = 0.25*second

defaultclock.dt = 0.001*ms

Currents should be specified in amps…

- if you do that, change C5 to be a resistance, in ohm

- also change > I_app : volt → I_app : amp

I_off = 0*mV

I_on = 10*mV #*uamp

I_app_start = 20ms #150ms

neuron_type = ‘chattering’

#%%

a,b,c,d define the type of neuron

-------------------------------------- a b c d

abcd_neurons = {‘regular spiking’: [0.02/ms, 0.2, -65mV, 8mV],

            '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],

            }

if neuron_type == ‘random’:

neuron_type = random.choice(list(abcd_neurons.keys())) #uncomment to get a random neuron type

a,b,c,d = abcd_neurons[neuron_type]

NEURON CONSTANTS

v_th = 30 * mV

C1 = 0.04 * 1/mV

C2 = 5.0

C3 = 140 * mV

C5 = 1.0 #* kohm # resistance to injected current

eqs = ‘’’

dv/dt = (C1*(v**2) + C2v + C3 - u + C5I_app) / ms: volt

du/dt = a * ( b * v - u ) : volt

I_app = I_app_fn(t) : volt

‘’’

v_0 = c

u_0 = 0 * mV

%%

SET UP NETWORK, RUN SIMULATION

create a function to simulate a step-change in current

I_app_fn = TimedArray( [I_off, I_on], dt=I_app_start)

G_pre = NeuronGroup( 2, model = eqs , threshold = ‘v > v_th’ , reset = ‘v = c ; u = u + d’ , method = ‘euler’ )

G_pre.v = [v_0 , v_0 ]

G_pre.u = [u_0 , u_0 ]

G_post = NeuronGroup( 1, model = eqs , threshold = ‘v > v_th’ , reset = ‘v = c ; u = u + d’ , method = ‘euler’ )

G_post.v = v_0

G_post.u = u_0

#synpses parameter

tau_c = 25 * ms

tau_f = 30 * ms # Facilitation time constant

tau_d = 50 * ms

tau_p = 120 * second

O_p = 0.9 * umolar / second

synapses_eqs = ‘’’

# Usage of releasable neurotransmitter per single action potential:

  du_S/dt = -u_S / tau_f  : 1 (event-driven)

# Fraction of synaptic neurotransmitter resources available for release:

        dx_S/dt = (1 - x_S)/tau_d : 1 (event-driven)

     '''

synapses_action = ‘’’

u_S += U_0 * ( 1 - u_S )

r_S = u_S * x_S

x_S -= r_S

'''

s= Synapses( G_pre, G_post , model= synapses_eqs , on_pre=synapses_action, dt=I_app_start , method=‘exact’)

M = StateMonitor(G_pre, variables=[‘v’,‘I_app’,‘u’], record=True)

SpikeMon = SpikeMonitor(G_pre)

run(duration)

###########################################
error :

---

KeyError Traceback (most recent call last)
/usr/local/lib/python3.7/dist-packages/brian2/core/network.py in before_run(self, run_namespace)
891 try:
→ 892 obj.before_run(run_namespace)
893 except Exception as ex:

17 frames
KeyError: ‘The identifier “U_0” could not be resolved.’

The above exception was the direct cause of the following exception:

BrianObjectException Traceback (most recent call last)
/usr/local/lib/python3.7/dist-packages/brian2/core/network.py in before_run(self, run_namespace)
892 obj.before_run(run_namespace)
893 except Exception as ex:
→ 894 raise BrianObjectException(“An error occurred when preparing an object.”, obj) from ex
895
896 # Check that no object has been run as part of another network before

BrianObjectException: Error encountered with object named ‘synapses_2_pre’.
Object was created here (most recent call only, full details in debug log):
File ‘’, line 86, in
s= Synapses( G_pre, G_post , model= synapses_eqs , on_pre=synapses_action, dt=I_app_start , method=‘exact’)

An error occurred when preparing an object. (See above for original error message and traceback.)

[adam]

Hey Mary,
it will be easier to compare code back and forth if, on this forum you use 3 “backticks” ( ` `` but without the space) before and after your code

print('code here')

this makes it easier to see the formatting in your code, but also if you don’t do that, the forum thinks * are italics for instance.

[adam]

the problem seems to be this line:
u_S += U_0 * ( 1 - u_S )
U_0 doesn’t exist. u_0, the initial value for the recovery value of the izhikevich neuron exists, but I don’t think this is what you want to use for u_S.

So I added a new variable u_S_0 instead. But I’m not sure what you want that to be. For now it’s 0, but that means u_S += u_S_0 * ( 1 - u_S ) doesn’t do anything yet.

Are you using a set of equations as a reference for the synapse? if so we can look at those and see what to do.

---

from brian2 import*
import numpy
import numpy as np
%matplotlib inline

import random
import matplotlib.pyplot as plt
'''
SUMMARY OF CHANGES: 
+ u_S_0 = 0.6
see: https://brian2.readthedocs.io/en/stable/examples/frompapers.Stimberg_et_al_2018.example_1_COBA.html
- u_S += U_0 * ( 1 - u_S ) 
+ u_S += u_S_0 * ( 1 - u_S )
'''
#%%

start_scope()

#SIMULATION PARAMETERS
duration = 0.25*second

defaultclock.dt = 0.001*ms
#Currents should be specified in amps…
#- if you do that, change C5 to be a resistance, in ohm
#- also change > I_app : volt → I_app : amp
I_off = 0*mV

I_on = 10*mV #*uamp

I_app_start = 20*ms #150ms

neuron_type = 'chattering'

#%%

#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],
            }
if neuron_type == 'random':
    neuron_type = random.choice(list(abcd_neurons.keys())) #uncomment to get a random neuron type
a,b,c,d = abcd_neurons[neuron_type]

# NEURON CONSTANTS
v_th = 30 * mV
C1 = 0.04 * 1/mV
C2 = 5.0
C3 = 140 * mV
C5 = 1.0 #* kohm # resistance to injected current

eqs = '''
dv/dt = (C1*(v**2) + C2*v + C3 - u + C5*I_app) / ms: volt
du/dt = a * ( b * v - u ) : volt
I_app = I_app_fn(t) : volt
'''

v_0 = c
u_0 = 0 * mV

u_S_0 = 0.6 #synaptic release probability at rest

#%%
# SET UP NETWORK, RUN SIMULATION
# create a function to simulate a step-change in current
I_app_fn = TimedArray( [I_off, I_on], dt=I_app_start)

G_pre = NeuronGroup( 2, model = eqs , threshold = 'v > v_th' , reset = 'v = c ; u = u + d' , method = 'euler' )
G_pre.v = [v_0 , v_0 ]
G_pre.u = [u_0 , u_0 ]

G_post = NeuronGroup( 1, model = eqs , threshold = 'v > v_th' , reset = 'v = c ; u = u + d' , method = 'euler' )
G_post.v = v_0
G_post.u = u_0

#synpses parameter

tau_c = 25 * ms
tau_f = 30 * ms # Facilitation time constant
tau_d = 50 * ms
tau_p = 120 * second
O_p = 0.9 * umolar / second


#CHANGED
synapses_eqs = '''
# Usage of releasable neurotransmitter per single action potential:
  du_S/dt = -u_S / tau_f  : 1 (event-driven)

# Fraction of synaptic neurotransmitter resources available for release:
dx_S/dt = (1 - x_S)/tau_d : 1 (event-driven)
'''

synapses_action = '''
u_S += u_S_0 * ( 1 - u_S )
r_S = u_S * x_S
x_S -= r_S
'''

s = Synapses( G_pre, G_post , model= synapses_eqs , on_pre=synapses_action, dt=I_app_start , method='exact')
s.connect()
M = StateMonitor(G_pre, variables=['v','I_app','u'], record=True)

SpikeMon = SpikeMonitor(G_pre)

run(duration)
#%%
print(SpikeMon.t)

[mstimberg]

@adam

Are you using a set of equations as a reference for the synapse? if so we can look at those and see what to do.

It looks like the Tsodyks&Markram model for short-term synaptic plasticity as e.g. implemented here: Example: example_1_COBA — Brian 2 2.5.4 documentation. The U_0 parameter is the “synaptic release probability at rest” and should be a value between 0 and 1 (in our example 0.6).

[Maryidea]

yes , sorry i didn’t know this

[Maryidea]

oh , yes , i didn’t pay attention to this . thank you so much , u_0 in synapses_action is about “the average synaptic release” with value of 0.8 .
i want to use the equations of this paper :[https://urldefense.com/v3/Modulation of Synaptic Plasticity by Glutamatergic Gliotransmission: A Modeling Study;!!CjcC7IQ!e8g52zrANCsN9p-DQ-TgnyyEQ5eYt7fs7y1L-cikVeZiFk1SbFsPoxZFy_hmbfyIOViEvaKbq2k$]
but instead of HH neuron i want to use three izhikevich neurons

[adam]

great, so setting u_S_0 to 0.8 in my most recent code snippet should work!

[Maryidea]

yes , i done it , thank you for your attention.