"""
--------------------------------------------------------------------------------
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
'''