"""
Name: hhca_psc_alpha - Hodgkin Huxley neuron model.

Description:

 hhca_psc_alpha is an implementation of a spiking neuron using the Hodkin-Huxley
 formalism.

 (1) Post-syaptic currents
 Incoming spike events induce a post-synaptic change of current modelled
 by an alpha function. The alpha function is normalised such that an event of
 weight 1.0 results in a peak current of 1 pA.


 (2) Spike Detection
 Spike detection is done by a combined threshold-and-local-maximum search: if
 there is a local maximum above a certain threshold of the membrane potential,
 it is considered a spike.

Problems/Todo:

 better spike detection
 initial wavelet/spike at simulation onset

References:

 Spiking Neuron Models:
 Single Neurons, Populations, Plasticity
 Wulfram Gerstner, Werner Kistler,  Cambridge University Press

 Theoretical Neuroscience:
 Computational and Mathematical Modeling of Neural Systems
 Peter Dayan, L. F. Abbott, MIT Press (parameters taken from here)

 Hodgkin, A. L. and Huxley, A. F.,
 A Quantitative Description of Membrane Current
 and Its Application to Conduction and Excitation in Nerve,
 Journal of Physiology, 117, 500-544 (1952)

Sends: SpikeEvent

Receives: SpikeEvent, CurrentEvent, DataLoggingRequest

Authors: Tanguy Fardet
SeeAlso: hh_psc_alpha
"""
neuron hhca_psc_alpha:
  state:
    r integer = 0 # number of steps in the current refractory phase
    V_m mV = -65. mV # Membrane potential
    Ca_i mmol/L = Ca_i_eq

    inline alpha_n_init real = ( 0.01 * ( V_m / mV + 55. ) ) / ( 1. - exp( -( V_m / mV + 55. ) / 10. ) )
    inline beta_n_init  real = 0.125 * exp( -( V_m / mV + 65. ) / 80. )
    inline alpha_m_init real = ( 0.1 * ( V_m / mV + 40. ) ) / ( 1. - exp( -( V_m / mV + 40. ) / 10. ) )
    inline beta_m_init  real = 4. * exp( -( V_m / mV + 65. ) / 18. )
    inline alpha_h_init real = 0.07 * exp( -( V_m / mV + 65. ) / 20. )
    inline beta_h_init  real = 1. / ( 1. + exp( -( V_m / mV + 35. ) / 10. ) )

    Act_m real =  alpha_m_init / ( alpha_m_init + beta_m_init )   # Activation variable m for Na
    Inact_h real = alpha_h_init / ( alpha_h_init + beta_h_init )    # Inactivation variable h for Na
    Act_n real =  alpha_n_init / ( alpha_n_init + beta_n_init ) # Activation variable n for K

    h_Ca  real = y_shape(V_m, V_hhCa, -k_hCa)
    m_Ca  real = y_shape(V_m, V_hmCa, k_mCa)
    m_AHP real = 0.
  end

  function y_shape(V_m mV, Vh mV, k_y mV) real:
    arg_exp real = min(-(V_m-Vh)/k_y, 50.)
    res real     = clip(1./(1. + exp(arg_exp)), 0., 1.)
    return res
  end

  parameters:
    t_ref ms = 2.0 ms      # Refractory period
    g_Na nS = 12000.0 nS    # Sodium peak conductance
    g_K nS = 3600.0 nS      # Potassium peak conductance
    g_L nS = 30 nS          # Leak conductance
    C_m pF = 100.0 pF       # Membrane Capacitance
    E_Na mV = 50 mV         # Sodium reversal potential
    E_K mV = -77. mV        # Potassium reversal potential
    E_L mV = -54.402 mV     # Leak reversal Potential (aka resting potential)
    tau_syn_ex ms = 0.2 ms  # Rise time of the excitatory synaptic alpha function i
    tau_syn_in ms = 2.0 ms  # Rise time of the inhibitory synaptic alpha function

    g_Ca  nS = 0.001 nS
    Ca_env  mmol/L = 1. mmol/L
	Ca_i_eq mmol/L = 5e-5 mmol/L
	tau_Ca      ms = 500. ms

	V_hmCa  mV = -27.5 mV
	k_mCa   mV = 5.7 mV
	tau_mCa ms = 0.5 ms

	V_hhCa  mV = -52.4 mV
	k_hCa   mV = 5.2 mV
	tau_hCa ms = 18. ms

    g_AHP nS = 0.2 nS
	K_AHP   mmol/L = 2e-3 mmol/L
	b_AHP     real = 2.5
	tau_AHP     ms = 30. ms

    E_0     mV = 26.64 mV
    Vcell    um**3 = 4/3.*3.14159*(5*um)**3

    # constant external input current
    I_e pA = 0 pA
  end

  internals:
    RefractoryCounts integer = steps(t_ref) # refractory time in steps
    charge       C = 1.6e-19 C
    avogadro 1/mol = 6.02e23/mol

    c2c mmol/L/fC = 1./(charge*avogadro*Vcell) # convert charge flux to concentration roughly 1e-2 mmol/L/fC
  end

  equations:
    # synapses: alpha functions
    kernel I_syn_in = (e/tau_syn_in) * t * exp(-t/tau_syn_in)
    kernel I_syn_ex = (e/tau_syn_ex) * t * exp(-t/tau_syn_ex)

	inline m_infAHP real = (Ca_i/K_AHP)**2/((Ca_i/K_AHP)**2 + b_AHP)
	inline m_infCa  real = y_shape(V_m, V_hmCa, k_mCa)
	inline h_infCa real = y_shape(V_m, V_hhCa, -k_hCa)

    inline E_Ca mV = 2*E_0*ln(Ca_env/Ca_i)

    inline I_syn_exc pA = convolve(I_syn_ex, spikeExc)
    inline I_syn_inh pA = convolve(I_syn_in, spikeInh)
    inline I_Na  pA = g_Na * Act_m * Act_m * Act_m * Inact_h * ( V_m - E_Na )
    inline I_K   pA  = ( g_K * Act_n*Act_n*Act_n*Act_n + g_AHP*m_AHP**2 ) * ( V_m - E_K )
    inline I_Ca   pA = -g_Ca*m_Ca*h_Ca*(V_m - E_Ca)
    inline I_L   pA = g_L * ( V_m - E_L )

    V_m'  = ( -( I_Na + I_K + I_L ) + I_e + I_stim + I_syn_inh + I_syn_exc ) / C_m
    Ca_i' = c2c*I_Ca - (Ca_i - Ca_i_eq) / tau_Ca

    # Act_n
    inline alpha_n real = ( 0.01 * ( V_m / mV + 55. ) ) / ( 1. - exp( -( V_m / mV + 55. ) / 10. ) )
    inline beta_n  real = 0.125 * exp( -( V_m / mV + 65. ) / 80. )
    Act_n' = ( alpha_n * ( 1 - Act_n ) - beta_n * Act_n ) / ms # n-variable

    # Act_m
    inline alpha_m real = ( 0.1 * ( V_m / mV + 40. ) ) / ( 1. - exp( -( V_m / mV + 40. ) / 10. ) )
    inline beta_m  real = 4. * exp( -( V_m / mV + 65. ) / 18. )
    Act_m' = ( alpha_m * ( 1 - Act_m ) - beta_m * Act_m ) / ms # m-variable

    # Inact_h'
    inline alpha_h real = 0.07 * exp( -( V_m / mV + 65. ) / 20. )
    inline beta_h  real = 1. / ( 1. + exp( -( V_m / mV + 35. ) / 10. ) )
    Inact_h' = ( alpha_h * ( 1 - Inact_h ) - beta_h * Inact_h ) / ms # h-variable

    h_Ca'  = (h_infCa - h_Ca) / tau_hCa
    m_Ca'  = (m_infCa - m_Ca) / tau_mCa
    m_AHP' = (m_infAHP - m_AHP) / tau_AHP
  end

  input:
    spikeInh pA <- inhibitory spike
    spikeExc pA <- excitatory spike
    I_stim pA <- continuous
  end

  output: spike

  update:
    U_old mV = V_m
    integrate_odes()
    # sending spikes: crossing 0 mV, pseudo-refractoriness and local maximum...
    if r > 0: # is refractory?
      r -= 1
    elif V_m > 0 mV and U_old > V_m: # threshold && maximum
      r = RefractoryCounts
      emit_spike()
    end

  end

end