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

Description:
The dynamics are given by:

   C_m dV_m/dt = g_L*(V-E_L) + I + I_syn_ex + I_syn_in
   tau_e depsilon/dt = (1-epsilon/(alpha*epsilon_0))**3 + (V-E_f)/(E_d-E_f)

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

   if v >= V_th:
     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 elif_psc_alpha_fast:

  state:
    r integer    = 0                # number of steps for refractory phase
    V_m mV       = E_0              # Membrane potential
    epsilon real = alpha*epsilon_0  # Energy
  end

  parameters:
    C_m pF = 100. pF        # Membrane capacitance
    g_L nS = 9. nS          # leak conductance
    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    # standard resting energy level
    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 detection threshold (reset condition)
    V_reset mV = -62. mV    # Reset potential
    tau_syn_ex ms = 0.2 ms  # Synaptic Time Constant Excitatory Synapse
    tau_syn_in ms = 2.0 ms  # Synaptic Time Constant for Inhibitory Synapse
    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
  end

  equations:
    inline eps_bound real = max(epsilon, 0.)  # non-negative energy
    # 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 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)

    V_m' = (g_L*(E_L - V_m) + I_e + I_in + I_ex + currents) * invCm

    epsilon' = ((1 - eps_bound*invae)*(1 - eps_bound*invae)*(1 - eps_bound*invae) - (V_m-E_f)*invEdEf) * 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 >= V_th and epsilon > epsilon_c:
      V_m  = V_reset
      epsilon -= delta
      emit_spike()
      r = RefractoryCounts
    end
  end

end