# -*- coding: utf-8 -*-

from brian2 import *

tau = 0.2 * msecond
tmax = 150 * ms

# Tonic spiking
a = 0.02
b = 0.2
c = -65 * mV
d = 6 * mV
I0 = 14 * mV

# Phasic bursting
# a = 0.02
# b = 0.25
# c = -55 * mV
# d = 0.05 * mV
# I0 = 0.6 * mV

# Tonic bursting
# a = 0.02
# b = 0.2
# c = -50 * mV
# d = 2 * mV
# I0 = 15 * mV


eqs= '''
dv/dt = (0.04*10**3 * volt**-1 * v**2 + 5 * v + 140*mvolt - u + I)/tau : volt
du/dt = a * (b * v - u)/tau : volt
I : volt
'''

reset= '''
v=c
u=u+d
'''

G = NeuronGroup(N=1, model=eqs,
              threshold='v >= 30*mV', reset=reset)

# L'input est de 0 mv pendant 20 ms
# Puis passe à I0
@network_operation(dt=50*ms)
def update_input():
    G.I = I0 - G.I


# On initialise le potentiel au potentiel de reset
G.v = -65*mV
G.u = b * G.v

# On monitore ses spikes et son potentiel
spikemon = SpikeMonitor(G)
vmon = StateMonitor(G, 'v', record=0)
umon = StateMonitor(G, 'u', record=0)
inputmon = StateMonitor(G, 'I', record=0)

run(tmax)

# On trace nos mesures
subplot(311)
plot(spikemon.t/ms, spikemon.i, '.k')
xlim(0*ms, tmax * 10**3) 
subplot(312)
plot(vmon.t, vmon.v.T)
plot(umon.t, umon.u.T)
ylabel('Potential (V)')
xlabel('Time(s)')
legend(['V','U'])
subplot(313)
plot(inputmon.t, inputmon.I.T)
ylim([-I0, 2*I0])
ylabel('Potential (V)')
xlabel('Time(s)')
legend(['I'])

show()