"""
Name: madexp_psc_alpha - Energy-based leaky integrate-and-fire neuron.

Description:
The dynamics are given by:
   C_m dV_m/dt   = g_L*(V-E_L) - w + I_e + I_syn_ex + I_syn_in
   tau_w dw/dt   = a(V-E_L) - epsilon/epsilon_0 w + I_KATP*epsilon_0/(epsilon_0 + epsilon)
   tau_e depsilon/dt = (1-epsilon/(alpha*epsilon_0))**3 - (V-E_f)/(E_d-E_f) - gamma*w

   E_L = E_0 + (E_u - E_0)(1-epsilon/epsilon_0)

   if V_m >= V_th and epsilon > epsilon_c:
     V_m is set to V_reset

   On each spike arrival, the membrane potential feels an alpha-shaped current
   of the form:
     I_syn = I_0 * t * exp(-t/tau_syn) / tau_syn.

Sends: SpikeEvent

Receives: SpikeEvent, CurrentEvent, DataLoggingRequest
FirstVersion: 2019
Author: Tanguy Fardet
"""

neuron madexp_psc_alpha:

  state:
    r integer    = 0                                          # number of steps for refractory phase
    epsilon real = alpha*epsilon_0                            # Energy
    V_m   mV     = E_0 + (E_u - E_0)*(1 - epsilon/epsilon_0)  # Membrane potential
    w     pA     = 0*pA                                       # Adaptation current
  end

  parameters:
    C_m pF = 100. pF           # Membrane capacitance
    g_L nS = 9. nS             # leak conductance
    Delta_T mV = 1.0 mV        # Slope factor
    V_peak mV = 0 mV           # Spike detection threshold
    E_0 mV = -65. mV           # resting potential
    E_u mV = -58. mV           # upper potential
    E_d mV = -50. mV           # energy depletion potential
    E_f mV = -60. mV           # energy inflexion potential
    epsilon_0 real = 0.5       # standard resting energy level
    epsilon_c real = 0.2       # energy threshold for spike generation
    alpha real = 1.            # energetic health
    delta real = 0.1           # energy consumption per spike
    tau_e ms = 1000. ms        # time constant for energy production
    I_e pA = 0. pA             # Constant input current
    V_th mV = -55. mV          # Spike generation threshold
    a nS = 1. nS               # subthreshold adaptation
    b pA = 2. pA               # spike-triggered adaptation
    gamma pA = 1000. pA        # normalization of adaptation energy
    I_KATP pA = 10. pA         # peak ATP-gated potassium current
    tau_w ms = 300. ms         # timescale of the adaptation current
    V_reset mV = -62. mV       # Reset potential
    tau_syn_ex ms = 0.2 ms     # Synaptic Time Constant (excitatory)
    tau_syn_in ms = 2.0 ms     # Synaptic Time Constant (inhibitory)
    t_ref ms = 2.0 ms          # Refractory period
  end

  internals:
    RefractoryCounts integer = steps(t_ref) # refractory time in steps
    invae real = 1./(alpha*epsilon_0)
    inveps real = 1./epsilon_0
    invEdEf 1/mV = 1./(E_d-E_f)
    invte 1/ms = 1./tau_e
    invCm 1/pF = 1./C_m
    invDT 1/mV = 1./Delta_T
    invtw 1/ms = 1./tau_w
    invgamma 1/pA = 1./gamma
    invtau_syn_in 1/ms = 1./tau_syn_in
    invtau_syn_ex 1/ms = 1./tau_syn_ex
    Vspike real = Delta_T == 0 ?  min(V_th, V_peak) : V_peak
  end

  function I_spike(epsilon real, V_m mV) pA:
    Ispk pA = 0. pA
    arg_exp real = 0.

    if Delta_T > 0. mV:
      arg_exp = clip((V_m - V_th) * invDT, -20., 20.)
      Ispk = (epsilon-epsilon_c)*g_L*Delta_T*exp(arg_exp) * inveps
    end

    return Ispk
  end

  equations:
    inline eps_bound real = max(epsilon, 0.)  # non-negative energy
    inline V_bound mV = min(V_m, V_peak)  # prevent exponential divergence

    # synapses: alpha functions
    kernel I_syn_in = e * invtau_syn_in * t * exp(-t*invtau_syn_in)
    kernel I_syn_ex = e * invtau_syn_ex * t * exp(-t*invtau_syn_ex)

    inline E_L mV = E_0 + (E_u - E_0)*(1 - eps_bound*inveps)

    inline I_in pA = convolve(I_syn_in, spikesInh)
    inline I_ex pA = convolve(I_syn_ex, spikesExc)

    inline Ispike pA = I_spike(eps_bound, V_bound)

    V_m' = (g_L*(E_L - V_bound) + Ispike - w + I_e + I_in + I_ex + currents) * invCm
    w'   = (a*(V_bound-E_L) - w + I_KATP*epsilon_c/(2*eps_bound + epsilon_c)) * invtw

    epsilon' = ((1 - eps_bound*invae)*(1 - eps_bound*invae)*(1 - eps_bound*invae) - (V_bound-E_f)*invEdEf - w * invgamma) * invte
  end

  input:
    spikesInh pA <- inhibitory spike
    spikesExc pA <- excitatory spike
    currents pA <- continuous
  end

  output: spike

  update:
    integrate_odes()

    if epsilon < 0.:
      epsilon = 0.
    end

    # refractoriness and threshold crossing
    if r > 0: # is refractory?
      r -= 1
      V_m  = V_reset
    elif V_m >= Vspike and epsilon > epsilon_c:
      # this test is necessary to support Delta_T = 0
      V_m  = V_reset
      epsilon -= delta
      w += b
      emit_spike()
      r = RefractoryCounts
    end
  end

end