import brian2 as br
from brian2 import umetre, ufarad, cm, siemens, mV, msiemens, nS, ms

# Hodgkin Huxley model from Brian2 documentation

# The model

area = 20000 * umetre ** 2
hodgkin_huxley_params = {
    "Cm": 1 * ufarad * cm ** -2 * area,
    "gl": 5e-5 * siemens * cm ** -2 * area,
    "El": -65 * mV,
    "EK": -90 * mV,
    "ENa": 50 * mV,
    "g_na": 100 * msiemens * cm ** -2 * area,
    "g_kd": 30 * msiemens * cm ** -2 * area,
    "VT": -63 * mV
}

hodgkin_huxley_eqs = br.Equations('''
dv/dt = (gl*(El-v) - g_na*(m*m*m)*h*(v-ENa) - g_kd*(n*n*n*n)*(v-EK) + I + ih)/Cm : volt
dm/dt = 0.32*(mV**-1)*(13.*mV-v+VT)/
    (exp((13.*mV-v+VT)/(4.*mV))-1.)/ms*(1-m)-0.28*(mV**-1)*(v-VT-40.*mV)/
    (exp((v-VT-40.*mV)/(5.*mV))-1.)/ms*m : 1
dn/dt = 0.032*(mV**-1)*(15.*mV-v+VT)/
    (exp((15.*mV-v+VT)/(5.*mV))-1.)/ms*(1.-n)-.5*exp((10.*mV-v+VT)/(40.*mV))/ms*n : 1
dh/dt = 0.128*exp((17.*mV-v+VT)/(18.*mV))/ms*(1.-h)-4./(1+exp((40.*mV-v+VT)/(5.*mV)))/ms*h : 1
I : amp
''')

hodgkin_huxley_cond_synapse_weak_noise_eqs = br.Equations('''
dv/dt = (gl*(El-v) - g_na*(m*m*m)*h*(v-ENa) - g_kd*(n*n*n*n)*(v-EK) + I + ih - g_syn*(v-E_i))/Cm + sigma_noise*sqrt(1./tau_noise)*xi : volt
dm/dt = 0.32*(mV**-1)*(13.*mV-v+VT)/
    (exp((13.*mV-v+VT)/(4.*mV))-1.)/ms*(1-m)-0.28*(mV**-1)*(v-VT-40.*mV)/
    (exp((v-VT-40.*mV)/(5.*mV))-1.)/ms*m : 1
dn/dt = 0.032*(mV**-1)*(15.*mV-v+VT)/
    (exp((15.*mV-v+VT)/(5.*mV))-1.)/ms*(1.-n)-.5*exp((10.*mV-v+VT)/(40.*mV))/ms*n : 1
dh/dt = 0.128*exp((17.*mV-v+VT)/(18.*mV))/ms*(1.-h)-4./(1+exp((40.*mV-v+VT)/(5.*mV)))/ms*h : 1
I : amp
g_syn: siemens
''')

hodgkin_huxley_with_external_input_eqs = br.Equations('''
dv/dt = (gl*(El-v) - g_na*(m*m*m)*h*(v-ENa) - g_kd*(n*n*n*n)*(v-EK) + I + ih)/Cm
 : volt
dm/dt = 0.32*(mV**-1)*(13.*mV-v+VT)/
    (exp((13.*mV-v+VT)/(4.*mV))-1.)/ms*(1-m)-0.28*(mV**-1)*(v-VT-40.*mV)/
    (exp((v-VT-40.*mV)/(5.*mV))-1.)/ms*m : 1
dn/dt = 0.032*(mV**-1)*(15.*mV-v+VT)/
    (exp((15.*mV-v+VT)/(5.*mV))-1.)/ms*(1.-n)-.5*exp((10.*mV-v+VT)/(40.*mV))/ms*n : 1
dh/dt = 0.128*exp((17.*mV-v+VT)/(18.*mV))/ms*(1.-h)-4./(1+exp((40.*mV-v+VT)/(5.*mV)))/ms*h : 1
I = I_ext_const + stimulus(t) : amp (constant over dt)
I_ext_const : amp
''')

hodgkin_huxley_eqs_with_synaptic_conductance = br.Equations('''
dv/dt = (gl*(El-v) - g_na*(m*m*m)*h*(v-ENa) - g_kd*(n*n*n*n)*(v-EK) + I + ih - g_syn*(v-E_i))/Cm : volt
dm/dt = 0.32*(mV**-1)*(13.*mV-v+VT)/
    (exp((13.*mV-v+VT)/(4.*mV))-1.)/ms*(1-m)-0.28*(mV**-1)*(v-VT-40.*mV)/
    (exp((v-VT-40.*mV)/(5.*mV))-1.)/ms*m : 1
dn/dt = 0.032*(mV**-1)*(15.*mV-v+VT)/
    (exp((15.*mV-v+VT)/(5.*mV))-1.)/ms*(1.-n)-.5*exp((10.*mV-v+VT)/(40.*mV))/ms*n : 1
dh/dt = 0.128*exp((17.*mV-v+VT)/(18.*mV))/ms*(1.-h)-4./(1+exp((40.*mV-v+VT)/(5.*mV)))/ms*h : 1
I : amp
g_syn: siemens
''')

# # First H-current from Izhikevich p. 48
# ghbar = 40. * nS #Other values would be 0.5, 2, 3.5, 20 depending on neuron type (Rothman, Manis)
# Eh = -43*mV
# V_half_h = -75.0
# k_h = -5.5
# V_max_h = -75.
# sigma_h = 15.
# C_amp_h = 1000.
# C_base_h = 100.
# eqs_ih = """
# ih = ghbar*r*(Eh-v) : amp
# dr/dt= (rinf-r)/rtau : 1
# rinf = 1. / (1+exp((V_half_h - v/mV) / k_h)) : 1
# rtau = (C_base_h + C_amp_h * exp(-(V_max_h-v/mV)**2./sigma_h**2)) * ms : second
# """
# Second H-current from Izhikevich p. 48


# eqs_ih = """
# ih = ghbar*r*(Eh-v) : amp
# dr/dt=(rinf-r)/rtau : 1
# rinf = 1. / (1+exp((v/mV + 90.) / 7.)) : 1
# # rtau = ((100000. / (237.*exp((v/mV+60.) / 12.) + 17.*exp(-(v/mV+60.) / 14.))) + 25.)*ms : second
# rtau = 0.01*(100. + 1000. * exp(-(-76.-v/mV)**2./15.**2)) * ms : second #From Izhikevich
# """


ih_params = {
    "ghbar": 40. * nS,  # Other values would be 0.5, 2, 3.5, 20 depending on neuron type (Rothman, Manis)
    "Eh": -1. * mV,
    "V_half_h": -95.0,
    "k_h": -5.5,
    "V_max_h": -75.,
    "sigma_h": 20.,
    "C_amp_h": 50.,
    "C_base_h": 10.,
}
eqs_ih = """
ih = ghbar*r*(Eh-v) : amp
dr/dt= (rinf-r)/rtau : 1
rinf = 1. / (1+exp((V_half_h - v/mV) / k_h)) : 1
rtau = 0.1*(C_base_h + C_amp_h * exp(-(V_max_h-v/mV)**2./sigma_h**2)) * ms : second
"""

lif_eqs = """
dv/dt =1.0/tau* (-v + u_ext) :volt (unless refractory)
u_ext = u_ext_const : volt
"""

lif_params = {
    "tau": 7 * ms,
    "v_threshold": -40 * mV,
    "v_reset": -60 * mV,
    "tau_refractory": 0.0 * ms,
    "u_ext_const": -50 * mV
}

lif_options = {
    "threshold": "v>v_threshold",
    "reset": "v=v_reset",
    "refractory": "tau_refractory",
    'method': 'euler'
}

delta_synapse_model = 'synaptic_strength: volt'
delta_synapse_on_pre = 'v+=synaptic_strength'
delta_synapse_param = {}
exponential_synapse = """
dg/dt = -g/tau_syn : siemens (clock-driven)
g_syn_post = g :siemens (summed)
synaptic_strength : siemens
"""
exponential_synapse_on_pre = "g+=synaptic_strength"
exponential_synapse_params = {
    "tau_syn": 1 * ms
}

noise_params = {
    'sigma_noise': 1.0 * mV,
    'tau_noise': 10. * ms
}