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:
