"""
--------------------------------------------------------------------------------
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_SynE - I_SynExt - I_SynI - I_SynHipp + G_sin*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) * (5.*mV) / exprel(-(v + 40.*mV) / (5.*mV)) / ms : Hz
betan = 0.5 * exp(- (v + 45.*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
alpham = 0.32 * (mV**-1) * (4.*mV) / exprel(-(v + 42.*mV) / (4.*mV)) / ms : Hz
betam = 0.28 * (mV**-1) * (5.*mV) / exprel( (v + 15.*mV) / (5.*mV)) / ms : Hz
dh/dt = alphah * (1 - h) - betah * h : 1
alphah = 0.128 * exp(- (v + 38.*mV) / (18.*mV)) / ms : Hz
betah = 4. / (1. + exp(- (v + 15.*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
alphamCaL = 0.055 * (mV**-1) * (3.8*mV) / exprel(-(v + 27.*mV) / (3.8*mV)) / ms : Hz
betamCaL = 0.94 * exp(-(v + 75.*mV) / (17.*mV)) / ms : Hz
dhCaL/dt = (alphahCaL * (1. - hCaL)) - (betahCaL * hCaL) : 1
alphahCaL = 0.000457 * exp(-(v + 13.*mV) / (50.*mV)) / ms : Hz
betahCaL = 0.0065 / (exp(-(v + 15.*mV) / (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
sigma : volt
G_sin = 1.*int(z_soma<15*mm)*int(z_soma>0*mm) : 1 # this is the mask/scaling for which neurons get the sinusoidal input
I_exc : amp (linked) # this is the input theta rhythm from the MS
# I_exc = inp_theta(t)*int(z_soma<15*mm)*int(z_soma>0*mm) : amp
r : 1
I_stim = inputs_stim(t) : amp
size : metre**2 (shared)
'''
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
Vm = (- I_CAN - I_M - I_leak - I_K - I_Na - I_Ca) / ((1.*ufarad*cm**-2) * (size))*tstep : 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) * (5.*mV) / exprel(-(v + 40.*mV) / (5.*mV)) / ms : Hz
betan = 0.5 * exp(- (v + 45.*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
alpham = 0.32 * (mV**-1) * (4.*mV) / exprel(-(v + 42.*mV) / (4.*mV)) / ms : Hz
betam = 0.28 * (mV**-1) * (5.*mV) / exprel( (v + 15.*mV) / (5.*mV)) / ms : Hz
dh/dt = alphah * (1. - h) - betah * h : 1
alphah = 0.128 * exp(- (v + 38.*mV) / (18.*mV)) / ms : Hz
betah = 4. / (1. + exp(- (v + 15.*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
alphamCaL = 0.055 * (mV**-1) * (3.8*mV) / exprel(-(v + 27.*mV) / (3.8*mV)) / ms : Hz
betamCaL = 0.94 * exp(-(v + 75.*mV) / (17.*mV)) / ms : Hz
dhCaL/dt = (alphahCaL * (1. - hCaL)) - (betahCaL * hCaL) : 1
alphahCaL = 0.000457 * exp(-(v + 13.*mV) / (50.*mV)) / ms : Hz
betahCaL = 0.0065 / (exp(-(v + 15.*mV) / (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
sigma : volt
r : 1
I_stim = inputs_stim(t) : amp
size : metre**2 (shared)
'''
#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
Vm = (- I_M - I_leak - I_K - I_Na - I_Ca) / ((1.*ufarad*cm**-2) * (size)) * tstep : 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) * (5.*mV) / exprel(-(v + 40.*mV) / (5.*mV)) / ms : Hz
betan = 0.5 * exp(- (v + 45.*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
alpham = 0.32 * (mV**-1) * (4.*mV) / exprel(-(v + 42.*mV) / (4.*mV)) / ms : Hz
betam = 0.28 * (mV**-1) * (5.*mV) / exprel( (v + 15.*mV) / (5.*mV)) / ms : Hz
dh/dt = alphah * (1. - h) - betah * h : 1
alphah = 0.128 * exp(- (v + 38.*mV) / (18.*mV)) / ms : Hz
betah = 4. / (1. + exp(- (v + 15.*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
alphamCaL = 0.055 * (mV**-1) * (3.8*mV) / exprel(-(v + 27.*mV) / (3.8*mV)) / ms : Hz
betamCaL = 0.94 * exp(-(v + 75.*mV) / (17.*mV)) / ms : Hz
dhCaL/dt = (alphahCaL * (1. - hCaL)) - (betahCaL * hCaL) : 1
alphahCaL = 0.000457 * exp(-(v + 13.*mV) / (50.*mV)) / ms : Hz
betahCaL = 0.0065 / (exp(-(v + 15.*mV) / (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
sigma : volt
r : 1
I_stim = inputs_stim(t) : amp
size : metre**2 (shared)
'''
""" Inhibitory Neuron Types """
""" ------------------------------------------------------------------------ """
inh_inp_eqs = '''
dv/dt = ( - I_leak - I_K - I_Na - I_SynE - I_SynExt - I_SynHipp - I_SynI + G_sin*I_exc + r*I_stim) / ((1.*ufarad*cm**-2) * (size)) + noise: volt
Vm = (- I_leak - I_K - I_Na) / ((1*ufarad*cm**-2) * (size))*tstep : 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 = (alphan * (1 - n) - betan * n) / 0.2: 1
alphan = 0.1 / exprel(-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 = (alpham * (1 - m) - betam * m) / 0.2 : 1
alpham = 1. / exprel(-(v + 35.*mV) / (10.*mV)) / ms : Hz
betam = 4. * exp(- (v + 60.*mV) / (18.*mV)) / ms : Hz
dh/dt = (alphah * (1 - h) - betah * h) / 0.2 : 1
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
G_sin = 1.*int(z_soma<15*mm)*int(z_soma>0*mm) : 1 # this is the mask/scaling for which neurons get the sinusoidal input
I_exc : amp (linked) # same as in the pyCAN group, excitatory input from MS
# I_exc = inp_theta(t)*int(z_soma<15*mm)*int(z_soma>0*mm) : amp
sigma : volt
r : 1
I_stim = inputs_stim(t) : amp
size : metre**2 (shared)
'''
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
Vm = (- I_leak - I_K - I_Na) / ((1.*ufarad*cm**-2) * (size))*tstep : 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 = (alphan * (1 - n) - betan * n) / 0.2 : 1
alphan = 0.1 / exprel(-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 = (alpham * (1 - m) - betam * m) / 0.2 : 1
alpham = 1. / exprel(-(v + 35.*mV) / (10.*mV)) / ms : Hz
betam = 4. * exp(- (v + 60.*mV) / (18.*mV)) / ms : Hz
dh/dt = (alphah * (1 - h) - betah * h) / 0.2 : 1
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
sigma : volt
r : 1
I_stim = inputs_stim(t) : amp
size : metre**2 (shared)
'''
# Spike and reset
reset_eqs = '''
glu = glu - 0.
Cl = Cl + 0.2
'''
