HH_equations2.py
"""
--------------------------------------------------------------------------------
Date: 29/09/2021
@author: Nikolaos Vardalakis
--------------------------------------------------------------------------------
Implementation Notes
--------------------------------------------------------------------------------
| 1: I removed the {Vm = ...} equations; were not used anywhere in the original model
| 2: Stimulation is weighed by distance {r} from the electrode position
| 3: Re-wrote the equations; simplified the alpham/betam for I_Na in the inhibitory equations set
"""
""" Excitatory Neuron Types """
""" ------------------------------------------------------------------------ """
# Pyramidal CAN
py_CAN_inp_eqs = '''
dv/dt = ( - I_CAN - I_M - I_leak - I_K - I_Na - I_Ca + I_exc + r*I_stim) / ((1 * ufarad * cm ** -2) * (size)) + noise: volt
I_CAN = ((gCAN) * (size)) * mCAN ** 2 * (v - (-20 * mV)) : amp
dmCAN/dt = (mCANInf - mCAN) / mCANTau : 1
mCANInf = alpha2 / (alpha2 + (0.0002 * ms ** -1)) : 1
mCANTau = 1. / (alpha2 + (0.0002 * ms ** -1)) / (3.0 ** ((36. -22) / 10)) : second
alpha2 = (0.0002 * ms ** -1) * (Ca_i / (5e-4 * mole * metre ** -3)) ** 2 : Hz
I_M = ((gM) * (size)) * p * (v - Ek) : amp
dp/dt = (pInf - p) / pTau : 1
pInf = 1. / (1 + exp(- (v + (35 * mV)) / (10 * mV))) : 1
pTau = (1000 * ms) / (3.3 * exp((v + (35 * mV)) / (20 * mV)) + exp(- (v + (35 * mV)) / (20 * mV))) : second
I_leak = ((1e-5 * siemens * cm ** -2) * (size)) * (v - (-70*mV)) : amp
I_K = ((5 * msiemens * cm ** -2) * (size)) * (n ** 4) * (v - Ek) : amp
dn/dt = alphan * (1 - n) - betan * n : 1
alphan = - 0.032 * (mV ** -1) * (v - (-55 * mV) - 15 * mV) / (exp(- (v - (-55 * mV) - 15 * mV) / (5 * mV)) - 1.) / ms : Hz
betan = 0.5 * exp( - (v - (-55 * mV) - 10 * mV) / (40 * mV)) / ms : Hz
I_Na = ((50 * msiemens * cm ** -2) * (size)) * (m ** 3) * h * (v - (50 * mV)) : amp
dm/dt = alpham * (1 - m) - betam * m : 1
dh/dt = alphah * (1 - h) - betah * h : 1
alpham = - 0.32 * (mV ** -1) * (v - (-55 * mV) - 13 * mV) / (exp(- (v - (-55 * mV) - 13 * mV) / (4 * mV)) - 1.) / ms : Hz
betam = 0.28 * (mV ** -1) * (v - (-55 * mV) - 40 * mV) / (exp((v - (-55 * mV) - 40 * mV) / (5 * mV)) - 1.) / ms : Hz
alphah = 0.128 * exp(- (v - (-55 * mV) - 17 * mV) / (18 * mV)) / ms : Hz
betah = 4. / (1 + exp(- (v - (-55 * mV) - 40 * mV) / (5 * mV))) / ms : Hz
I_Ca = ((1e-4 * siemens * cm ** -2) * (size)) * (mCaL ** 2) * hCaL * (v - (120 * mV)) : amp
dmCaL/dt = (alphamCaL * (1 - mCaL)) - (betamCaL * mCaL) : 1
dhCaL/dt = (alphahCaL * (1 - hCaL)) - (betahCaL * hCaL) : 1
alphamCaL = (0.055 * mV ** -1) * ((-27 * mV) - v) / (exp(((-27 * mV) - v) / (3.8 * mV)) - 1.) / ms : Hz
betamCaL = 0.94 * exp(((-75 * mV) - v) / (17 * mV)) / ms : Hz
alphahCaL = 0.000457 * exp(((-13 * mV) - v) / (50 * mV)) / ms : Hz
betahCaL = 0.0065 / (exp(((-15 * mV) - v) / (28 * mV)) + 1.) / ms : Hz
dCa_i/dt = driveChannel + ((2.4e-4 * mole * metre**-3) - Ca_i) / (200 * ms) : mole * meter**-3
driveChannel = (-(1e4) * I_Ca / (cm ** 2)) / (2 * (96489 * coulomb * mole ** -1) * (1 * umetre)) : mole * meter ** -3 * Hz
I_SynE = + ge * (v - (0 * mV)) : amp
dge/dt = (-ge+he) * (1. / (0.3 * ms)) : siemens
dhe/dt=-he/(5*ms) : siemens
I_SynExt = + ge_ext * (v - (0 * mV)) : amp
dge_ext/dt = (-ge_ext+he_ext) * (1. / (0.3 * ms)) : siemens
dhe_ext/dt=-he_ext/(5*ms) : siemens
I_SynHipp = + ge_hipp * (v - (0 * mV)) : amp
dge_hipp/dt = (-ge_hipp+he_hipp) * (1. / (0.3 * ms)) : siemens
dhe_hipp/dt = -he_hipp/(5*ms) : siemens
I_SynI = + gi*(v - 0 * mV)*int(Cl>0.5) + gi*(v - (-80 * mV))*int(Cl<=0.5): amp
dgi/dt = (-gi+hi) * (1. / (1 * ms)) : siemens
dhi/dt=-hi/(10*ms) : siemens
dCl/dt=-Cl/tau_Cl:1
dglu/dt=(1-glu)/(3*second):1
noise = sigma*(2*(0.1e-3 * siemens ) / (1 * ufarad))**.5*randn()/sqrt(tstep) : volt/second (constant over dt)
x_soma:metre
y_soma:metre
z_soma:metre
I_exc = inp_theta(t)*int(z_soma<100*scale_aussel)*int(z_soma>0*mm) : amp
# I_exc = inp_theta(t) : amp
r : 1
I_stim = inputs_stim(t) : amp
sigma : volt
size:metre ** 2
'''
py_CAN_eqs = '''
dv/dt = ( - I_CAN - I_M - I_leak - I_K - I_Na - I_Ca - I_SynE - I_SynExt - I_SynI - I_SynHipp + r*I_stim) / ((1 * ufarad * cm ** -2) * (size)) + noise: volt
I_CAN = ((gCAN) * (size)) * mCAN ** 2 * (v - (-20 * mV)) : amp
dmCAN/dt = (mCANInf - mCAN) / mCANTau : 1
mCANInf = alpha2 / (alpha2 + (0.0002 * ms ** -1)) : 1
mCANTau = 1. / (alpha2 + (0.0002 * ms ** -1)) / (3.0 ** ((36. -22) / 10)) : second
alpha2 = (0.0002 * ms ** -1) * (Ca_i / (5e-4 * mole * metre ** -3)) ** 2 : Hz
I_M = ((gM) * (size)) * p * (v - Ek) : amp
dp/dt = (pInf - p) / pTau : 1
pInf = 1. / (1 + exp(- (v + (35 * mV)) / (10 * mV))) : 1
pTau = (1000 * ms) / (3.3 * exp((v + (35 * mV)) / (20 * mV)) + exp(- (v + (35 * mV)) / (20 * mV))) : second
I_leak = ((1e-5 * siemens * cm ** -2) * (size)) * (v - (-70*mV)) : amp
I_K = ((5 * msiemens * cm ** -2) * (size)) * (n ** 4) * (v - Ek) : amp
dn/dt = alphan * (1 - n) - betan * n : 1
alphan = - 0.032 * (mV ** -1) * (v - (-55 * mV) - 15 * mV) / (exp(- (v - (-55 * mV) - 15 * mV) / (5 * mV)) - 1.) / ms : Hz
betan = 0.5 * exp( - (v - (-55 * mV) - 10 * mV) / (40 * mV)) / ms : Hz
I_Na = ((50 * msiemens * cm ** -2) * (size)) * (m ** 3) * h * (v - (50 * mV)) : amp
dm/dt = alpham * (1 - m) - betam * m : 1
dh/dt = alphah * (1 - h) - betah * h : 1
alpham = - 0.32 * (mV ** -1) * (v - (-55 * mV) - 13 * mV) / (exp(- (v - (-55 * mV) - 13 * mV) / (4 * mV)) - 1.) / ms : Hz
betam = 0.28 * (mV ** -1) * (v - (-55 * mV) - 40 * mV) / (exp((v - (-55 * mV) - 40 * mV) / (5 * mV)) - 1.) / ms : Hz
alphah = 0.128 * exp(- (v - (-55 * mV) - 17 * mV) / (18 * mV)) / ms : Hz
betah = 4. / (1 + exp(- (v - (-55 * mV) - 40 * mV) / (5 * mV))) / ms : Hz
I_Ca = ((1e-4 * siemens * cm ** -2) * (size)) * (mCaL ** 2) * hCaL * (v - (120 * mV)) : amp
dmCaL/dt = (alphamCaL * (1 - mCaL)) - (betamCaL * mCaL) : 1
dhCaL/dt = (alphahCaL * (1 - hCaL)) - (betahCaL * hCaL) : 1
alphamCaL = (0.055 * mV ** -1) * ((-27 * mV) - v) / (exp(((-27 * mV) - v) / (3.8 * mV)) - 1.) / ms : Hz
betamCaL = 0.94 * exp(((-75 * mV) - v) / (17 * mV)) / ms : Hz
alphahCaL = 0.000457 * exp(((-13 * mV) - v) / (50 * mV)) / ms : Hz
betahCaL = 0.0065 / (exp(((-15 * mV) - v) / (28 * mV)) + 1.) / ms : Hz
dCa_i/dt = driveChannel + ((2.4e-4 * mole * metre**-3) - Ca_i) / (200 * ms) : mole * meter**-3
driveChannel = (-(1e4) * I_Ca / (cm ** 2)) / (2 * (96489 * coulomb * mole ** -1) * (1 * umetre)) : mole * meter ** -3 * Hz
I_SynE = + ge * (v - (0 * mV)) : amp
dge/dt = (-ge+he) * (1. / (0.3 * ms)) : siemens
dhe/dt=-he/(5*ms) : siemens
I_SynExt = + ge_ext * (v - (0 * mV)) : amp
dge_ext/dt = (-ge_ext+he_ext) * (1. / (0.3 * ms)) : siemens
dhe_ext/dt=-he_ext/(5*ms) : siemens
I_SynHipp = + ge_hipp * (v - (0 * mV)) : amp
dge_hipp/dt = (-ge_hipp+he_hipp) * (1. / (0.3 * ms)) : siemens
dhe_hipp/dt = -he_hipp/(5*ms) : siemens
I_SynI = + gi*(v - 0 * mV)*int(Cl>0.5) + gi*(v - (-80 * mV))*int(Cl<=0.5): amp
dgi/dt = (-gi+hi) * (1. / (1 * ms)) : siemens
dhi/dt=-hi/(10*ms) : siemens
dCl/dt=-Cl/tau_Cl:1
dglu/dt=(1-glu)/(3*second):1
noise = sigma*(2*(0.1e-3 * siemens ) / (1 * ufarad))**.5*randn()/sqrt(tstep) : volt/second (constant over dt)
x_soma:metre
y_soma:metre
z_soma:metre
r : 1
I_stim = inputs_stim(t) : amp
sigma : volt
size:metre ** 2
'''
#Pyramidal non CAN :
py_eqs = '''
dv/dt = ( - I_M - I_leak - I_K - I_Na - I_Ca - I_SynE - I_SynExt - I_SynI - I_SynHipp + r*I_stim) / ((1 * ufarad * cm ** -2) * (size)) + noise : volt
I_M = ((gM) * (size)) * p * (v - Ek) : amp
dp/dt = (pInf - p) / pTau : 1
pInf = 1. / (1 + exp(- (v + (35 * mV)) / (10 * mV))) : 1
pTau = (1000 * ms) / (3.3 * exp((v + (35 * mV)) / (20 * mV)) + exp(- (v + (35 * mV)) / (20 * mV))) : second
I_leak = ((1e-5 * siemens * cm ** -2) * (size)) * (v - (-70*mV)) : amp
I_K = ((5 * msiemens * cm ** -2) * (size)) * (n ** 4) * (v - Ek) : amp
dn/dt = alphan * (1 - n) - betan * n : 1
alphan = - 0.032 * (mV ** -1) * (v - (-55 * mV) - 15 * mV) / (exp(- (v - (-55 * mV) - 15 * mV) / (5 * mV)) - 1.) / ms : Hz
betan = 0.5 * exp( - (v - (-55 * mV) - 10 * mV) / (40 * mV)) / ms : Hz
I_Na = ((50 * msiemens * cm ** -2) * (size)) * (m ** 3) * h * (v - (50 * mV)) : amp
dm/dt = alpham * (1 - m) - betam * m : 1
dh/dt = alphah * (1 - h) - betah * h : 1
alpham = - 0.32 * (mV ** -1) * (v - (-55 * mV) - 13 * mV) / (exp(- (v - (-55 * mV) - 13 * mV) / (4 * mV)) - 1.) / ms : Hz
betam = 0.28 * (mV ** -1) * (v - (-55 * mV) - 40 * mV) / (exp((v - (-55 * mV) - 40 * mV) / (5 * mV)) - 1.) / ms : Hz
alphah = 0.128 * exp(- (v - (-55 * mV) - 17 * mV) / (18 * mV)) / ms : Hz
betah = 4. / (1 + exp(- (v - (-55 * mV) - 40 * mV) / (5 * mV))) / ms : Hz
I_Ca = ((1e-4 * siemens * cm ** -2) * (size)) * (mCaL ** 2) * hCaL * (v - (120 * mV)) : amp
dmCaL/dt = (alphamCaL * (1 - mCaL)) - (betamCaL * mCaL) : 1
dhCaL/dt = (alphahCaL * (1 - hCaL)) - (betahCaL * hCaL) : 1
alphamCaL = (0.055 * mV ** -1) * ((-27 * mV) - v) / (exp(((-27 * mV) - v) / (3.8 * mV)) - 1.) / ms : Hz
betamCaL = 0.94 * exp(((-75 * mV) - v) / (17 * mV)) / ms : Hz
alphahCaL = 0.000457 * exp(((-13 * mV) - v) / (50 * mV)) / ms : Hz
betahCaL = 0.0065 / (exp(((-15 * mV) - v) / (28 * mV)) + 1.) / ms : Hz
dCa_i/dt = driveChannel + ((2.4e-4 * mole * metre**-3) - Ca_i) / (200 * ms) : mole * meter**-3
driveChannel = (-(1e4) * I_Ca / (cm ** 2)) / (2 * (96489 * coulomb * mole ** -1) * (1 * umetre)) : mole * meter ** -3 * Hz
I_SynE = + ge * (v - (0 * mV)) : amp
dge/dt = (-ge+he) * (1. / (0.3 * ms)) : siemens
dhe/dt=-he/(5*ms) : siemens
I_SynExt = + ge_ext * (v - (0 * mV)) : amp
dge_ext/dt = (-ge_ext+he_ext) * (1. / (0.3 * ms)) : siemens
dhe_ext/dt = -he_ext/(5*ms) : siemens
I_SynHipp = + ge_hipp * (v - (0 * mV)) : amp
dge_hipp/dt = (-ge_hipp+he_hipp) * (1. / (0.3 * ms)) : siemens
dhe_hipp/dt = -he_hipp/(5*ms) : siemens
I_SynI = + gi*(v - 0 * mV)*int(Cl>0.5) + gi*(v - (-80 * mV))*int(Cl<=0.5): amp
dgi/dt = (-gi+hi) * (1. / (1 * ms)) : siemens
dhi/dt = -hi/(10*ms) : siemens
dCl/dt = -Cl/tau_Cl:1
dglu/dt = (1-glu)/(3*second):1
noise = sigma*(2*(0.1e-3 * siemens ) / (1 * ufarad))**.5*randn()/sqrt(tstep) : volt/second (constant over dt)
x_soma:metre
y_soma:metre
z_soma:metre
r : 1
I_stim = inputs_stim(t) : amp
sigma : volt
size:metre ** 2
'''
""" Inhibitory Neuron Types """
""" ------------------------------------------------------------------------ """
inh_inp_eqs = '''
dv/dt = ( - I_leak - I_K - I_Na - I_SynE - I_SynExt - I_SynHipp - I_SynI + I_exc + r*I_stim) / ((1 * ufarad * cm ** -2) * (size)) + noise: volt
I_leak = ((0.1e-3 * siemens * cm ** -2) * (size)) * (v - (-65 * mV)) : amp
I_K = ((9e-3 * siemens * cm ** -2) * (size)) * (n ** 4) * (v - (-90 * mV)) : amp
dn/dt = (n_inf - n) / tau_n : 1
n_inf = alphan / (alphan + betan) : 1
tau_n = 0.2 / (alphan + betan) : second
alphan = 0.01 * (mV ** -1) * (v + 34 * mV) / (1. - exp(- 0.1 * (mV ** -1) * (v + 34 * mV))) / ms : Hz
betan = 0.125 * exp( - (v + 44 * mV) / (80 * mV)) / ms : Hz
I_Na = ((35e-3 * siemens * cm ** -2) * (size)) * (m ** 3) * h * (v - (55 * mV)) : amp
dm/dt = (m_inf - m) / tau_m : 1
dh/dt = (h_inf - h) / tau_h : 1
m_inf = alpham / (alpham + betam) : 1
tau_m = 0.2 / (alpham + betam) : second
h_inf = alphah / (alphah + betah) : 1
tau_h = 0.2 / (alphah + betah) : second
alpham = 0.1 * (mV ** -1) * (v + 35 * mV) / (1. - exp(- (v + 35 * mV) / (10 * mV))) / ms : Hz
betam = 4 * exp(- (v + 60 * mV) / (18 * mV)) / ms : Hz
alphah = 0.07 * exp(- (v + 58 * mV) / (20 * mV)) / ms : Hz
betah = 1. / (exp((- 0.1 * (mV ** -1)) * (v + 28 * mV)) + 1.) / ms : Hz
I_SynE = + ge * (v - (0 * mV)) : amp
dge/dt = (-ge+he) * (1. / (0.3 * ms)) : siemens
dhe/dt = -he/(5*ms) : siemens
I_SynExt = + ge_ext * (v - (0 * mV)) : amp
dge_ext/dt = (-ge_ext+he_ext) * (1. / (0.3 * ms)) : siemens
dhe_ext/dt = -he_ext/(5*ms) : siemens
I_SynHipp = + ge_hipp * (v - (0 * mV)) : amp
dge_hipp/dt = (-ge_hipp+he_hipp) * (1. / (0.3 * ms)) : siemens
dhe_hipp/dt = -he_hipp/(5*ms) : siemens
I_SynI = + gi * (v - (-80 * mV)) : amp
dgi/dt = (-gi+hi) * (1. / (1 * ms)) : siemens
dhi/dt = -hi/(10*ms) : siemens
noise = sigma*(2*(0.1e-3 * siemens ) / (1 * ufarad))**.5*randn()/sqrt(tstep) : volt/second (constant over dt)
x_soma:metre
y_soma:metre
z_soma:metre
I_exc = inp_theta(t)*int(z_soma<100*scale_aussel)*int(z_soma>0*mm) : amp
r : 1
I_stim = inputs_stim(t) : amp
sigma : volt
size:metre ** 2
'''
inh_eqs = '''
dv/dt = ( - I_leak - I_K - I_Na - I_SynE - I_SynExt - I_SynHipp - I_SynI + r*I_stim) / ((1 * ufarad * cm ** -2) * (size)) + noise: volt
I_leak = ((0.1e-3 * siemens * cm ** -2) * (size)) * (v - (-65 * mV)) : amp
I_K = ((9e-3 * siemens * cm ** -2) * (size)) * (n ** 4) * (v - (-90 * mV)) : amp
dn/dt = (n_inf - n) / tau_n : 1
n_inf = alphan / (alphan + betan) : 1
tau_n = 0.2 / (alphan + betan) : second
alphan = 0.01 * (mV ** -1) * (v + 34 * mV) / (1. - exp(- 0.1 * (mV ** -1) * (v + 34 * mV))) / ms : Hz
betan = 0.125 * exp( - (v + 44 * mV) / (80 * mV)) / ms : Hz
I_Na = ((35e-3 * siemens * cm ** -2) * (size)) * (m ** 3) * h * (v - (55 * mV)) : amp
dm/dt = (m_inf - m) / tau_m : 1
dh/dt = (h_inf - h) / tau_h : 1
m_inf = alpham / (alpham + betam) : 1
tau_m = 0.2 / (alpham + betam) : second
h_inf = alphah / (alphah + betah) : 1
tau_h = 0.2 / (alphah + betah) : second
alpham = 0.1 * (mV ** -1) * (v + 35 * mV) / (1. - exp(- (v + 35 * mV) / (10 * mV))) / ms : Hz
betam = 4 * exp(- (v + 60 * mV) / (18 * mV)) / ms : Hz
alphah = 0.07 * exp(- (v + 58 * mV) / (20 * mV)) / ms : Hz
betah = 1. / (exp((- 0.1 * (mV ** -1)) * (v + 28 * mV)) + 1.) / ms : Hz
I_SynE = + ge * (v - (0 * mV)) : amp
dge/dt = (-ge+he) * (1. / (0.3 * ms)) : siemens
dhe/dt = -he/(5*ms) : siemens
I_SynExt = + ge_ext * (v - (0 * mV)) : amp
dge_ext/dt = (-ge_ext+he_ext) * (1. / (0.3 * ms)) : siemens
dhe_ext/dt = -he_ext/(5*ms) : siemens
I_SynHipp = + ge_hipp * (v - (0 * mV)) : amp
dge_hipp/dt = (-ge_hipp+he_hipp) * (1. / (0.3 * ms)) : siemens
dhe_hipp/dt = -he_hipp/(5*ms) : siemens
I_SynI = + gi * (v - (-80 * mV)) : amp
dgi/dt = (-gi+hi) * (1. / (1 * ms)) : siemens
dhi/dt = -hi/(10*ms) : siemens
noise = sigma*(2*(0.1e-3 * siemens ) / (1 * ufarad))**.5*randn()/sqrt(tstep) : volt/second (constant over dt)
x_soma:metre
y_soma:metre
z_soma:metre
r : 1
I_stim = inputs_stim(t) : amp
sigma : volt
size:metre ** 2
'''
# Spike and reset
reset_eqs = '''
glu = glu - 0.
Cl = Cl + 0.2
'''
