Run error with modeling izhikevich neuron

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:
izhi_chattering

izhi_thalamo-cortical_burst
izhi_regular spiking

2 Likes