` Project 6 - "A numerical investigation of Izhikevich's Simple Model of Spiking Neurons"`
` Team Leader: Caroline Hill`
` Team Member(s): Caroline Hill`

read `DMB.txt`
Help6:=proc()  
if nargs=0 then  
("The available functions/procedures are"):

elif nargs=1 and args[1]=SMNS  then
"SMNS(a,b,c,d,u,v): Eugene M. Izhikevich's simple model of neuron spiking with parameters {a,b,c,d} and variables {u,v}"

else
"There is no Help for", args

fi:
end:

`Parameters and Variables:`
`u = membrane recovery of neuron`
`v = membrane potential of neuron`
`a = scales time for u (magnitude of this parameter and recovery speed are inversely proportional)`
`b = the sensitivity of u to the subthreshold fluctuations of v (magnitude of this parameter is proportional to the coupling of the two variables)`
`c = after-spike reset value of  v`
`d = after-spike reset value of u`
`i = current`

SMNS := [0.04*v^2 + 5*v + 140 - u + i, a*(b*v - u)]

convert(EquP(SMNS, [u, v]), radical)
`Equilibrium points given the parameter values typical for mammalian cotrical neurons - observe that the equilibria are complex, meaning the system has oscillatory behavior`
EquP(subs([a = 0.02, b = 0.2, c = -65, d = 2, i = 10], SMNS), [u, v])

`Plotting different neuron spiking dynamics by changing parameters`
`Figure 1: Regular Spiking (RS)`
a := 0.02;
b := 0.2;
c := -65.0;
d := 8.0;
v := -70.0;
u := b*v;
vv := [];
uu := [];
ts := 0.25;
N := round(100/ts);
tspan := [seq(evalf(q*ts), q = 0 .. N)];
T1 := tspan[-1]/10;
for q to N + 1 do
    t := tspan[q];
    if T1 < t then
        i := 14.0;
    else
        i := 0.;
    end if;
    v(x) := v(x) + 0.25*(0.04*diff(v(x), [x $ ~2]) + 5.0*v(x) + 140.0 - u + i);
    u := u + 0.25*a*(b*v(x) - u);
    if 30 < v(x) then
        vv := [op(vv), 30.0];
        v(x) := c;
        u := u + d;
    else
        vv := [op(vv), v(x)];
    end if;
    uu := [op(uu), u];
end do
n := nops(tspan);
P_1 := plots:-listplot([seq([tspan[q], vv[q]], q = 1 .. n)], color = black, style = line, title = "Regular Spiking (RS)", labels = ["Time (ms)", "Membrane Potential (mV)"], labeldirections = ['horizontal', 'vertical']);
plots:-display(P_1);

`Figure 2:  Intrinsically Bursting (IB)`
a := 0.02;
b := 0.2;
c := -55.0;
d := 4.0;
v := -70.0;
u := v*b;
vv := [];
uu := [];
ts := 0.25;
N := round(100/ts);
tspan := [seq(evalf(q*ts), q = 0 .. N)];
T1 := tspan[-1]/10;
for q to N + 1 do
    t := tspan[q];
    if T1 < t then
        i := 14.0;
    else
        i := 0.;
    end if;
    v(x) := v(x) + 0.25*(0.04*diff(v(x), [x $ ~2]) + 5.0*v(x) + 140.0 - u + i);
    u := u + 0.25*a*(b*v(x) - u);
    if 30 < v(x) then
        vv := [op(vv), 30.0];
        v(x) := c;
        u := u + d;
    else
        vv := [op(vv), v(x)];
    end if;
    uu := [op(uu), u];
end do
n := nops(tspan);
P_2 := plots:-listplot([seq([tspan[q], vv[q]], q = 1 .. n)], color = black, style = line, title = "Intrinsically Bursting (IB)", labels = ["Time (ms)", "Membrane Potential (mV)"], labeldirections = ['horizontal', 'vertical']);
plots:-display(P_2);

`Figure 3: Chattering (CH)`
a := 0.02;
b := 0.2;
c := -50.0;
d := 2.0;
v := -70.0;
u := v*b;
vv := [];
uu := [];
ts := 0.25;
N := round(100/ts);
tspan := [seq(evalf(q*ts), q = 0 .. N)];
T1 := tspan[-1]/10;
for q to N + 1 do
    t := tspan[q];
    if T1 < t then
        i := 14.0;
    else
        i := 0.;
    end if;
    v(x) := v(x) + 0.25*(0.04*diff(v(x), [x $ ~2]) + 5.0*v(x) + 140.0 - u + i);
    u := u + 0.25*a*(b*v(x) - u);
    if 30 < v(x) then
        vv := [op(vv), 30.0];
        v(x) := c;
        u := u + d;
    else
        vv := [op(vv), v(x)];
    end if;
    uu := [op(uu), u];
end do
n := nops(tspan);
P_3 := plots:-listplot([seq([tspan[q], vv[q]], q = 1 .. n)], color = black, style = line, title = "Chattering (CH)", labels = ["Time (ms)", "Membrane Potential (mV)"], labeldirections = ['horizontal', 'vertical']);
plots:-display(P_3);

`Figure 4: Fast Spiking (FS)`
a := 0.1;
b := 0.2;
c := -65.0;
d := 2.0;
v := -70.0;
u := v*b;
vv := [];
uu := [];
ts := 0.25;
N := round(100/ts);
tspan := [seq(evalf(q*ts), q = 0 .. N)];
T1 := tspan[-1]/10;
for q to N + 1 do
    t := tspan[q];
    if T1 < t then
        i := 14.0;
    else
        i := 0.;
    end if;
    v(x) := v(x) + 0.25*(0.04*diff(v(x), [x $ ~2]) + 5.0*v(x) + 140.0 - u + i);
    u := u + 0.25*a*(b*v(x) - u);
    if 30 < v(x) then
        vv := [op(vv), 30.0];
        v(x) := c;
        u := u + d;
    else
        vv := [op(vv), v(x)];
    end if;
    uu := [op(uu), u];
end do
n := nops(tspan);
P_4 := plots:-listplot([seq([tspan[q], vv[q]], q = 1 .. n)], color = black, style = line, title = "Fast Spiking (FS)", labels = ["Time (ms)", "Membrane Potential (mV)"], labeldirections = ['horizontal', 'vertical']);
plots:-display(P_4);

`Figure 5: Thalamo-Cortical (TC)`
a := 0.02;
b := 0.25;
c := -65.0;
d := 0.05;
v := -70.0;
u := v*b;
vv := [];
uu := [];
ts := 0.25;
N := round(100/ts);
tspan := [seq(evalf(q*ts), q = 0 .. N)];
T1 := tspan[-1]/10;
for q to N + 1 do
    t := tspan[q];
    if T1 < t then
        i := 14.0;
    else
        i := 0.;
    end if;
    v(x) := v(x) + 0.25*(0.04*diff(v(x), [x $ ~2]) + 5.0*v(x) + 140.0 - u + i);
    u := u + 0.25*a*(b*v(x) - u);
    if 30 < v(x) then
        vv := [op(vv), 30.0];
        v(x) := c;
        u := u + d;
    else
        vv := [op(vv), v(x)];
    end if;
    uu := [op(uu), u];
end do
n := nops(tspan);
P_5 := plots:-listplot([seq([tspan[q], vv[q]], q = 1 .. n)], color = black, style = line, title = "Thalamo-Cortical (TC)", labels = ["Time (ms)", "Membrane Potential (mV)"], labeldirections = ['horizontal', 'vertical']);
plots:-display(P_5);

`Figure 6: Resonator (RZ)`
a := 0.1;
b := 0.26;
c := -65.0;
d := 2.0;
v := -70.0;
u := v*b;
vv := [];
uu := [];
ts := 0.25;
N := round(100/ts);
tspan := [seq(evalf(q*ts), q = 0 .. N)];
T1 := tspan[-1]/10;
for q to N + 1 do
    t := tspan[q];
    if T1 < t then
        i := 14.0;
    else
        i := 0.;
    end if;
    v(x) := v(x) + 0.25*(0.04*diff(v(x), [x $ ~2]) + 5.0*v(x) + 140.0 - u + i);
    u := u + 0.25*a*(b*v(x) - u);
    if 30 < v(x) then
        vv := [op(vv), 30.0];
        v(x) := c;
        u := u + d;
    else
        vv := [op(vv), v(x)];
    end if;
    uu := [op(uu), u];
end do
n := nops(tspan);
P_6 := plots:-listplot([seq([tspan[q], vv[q]], q = 1 .. n)], color = black, style = line, title = "Resonator (RZ)", labels = ["Time (ms)", "Membrane Potential (mV)"], labeldirections = ['horizontal', 'vertical']);
plots:-display(P_6);

`Figure 7: Low-Threshold Spiking (LTS)`
a := 0.02;
b := 0.25;
c := -65.0;
d := 2.0;
v := -70.0;
u := v*b;
vv := [];
uu := [];
ts := 0.25;
N := round(100/ts);
tspan := [seq(evalf(q*ts), q = 0 .. N)];
T1 := tspan[-1]/10;
for q to N + 1 do
    t := tspan[q];
    if T1 < t then
        i := 14.0;
    else
        i := 0.;
    end if;
    v(x) := v(x) + 0.25*(0.04*diff(v(x), [x $ ~2]) + 5.0*v(x) + 140.0 - u + i);
    u := u + 0.25*a*(b*v(x) - u);
    if 30 < v(x) then
        vv := [op(vv), 30.0];
        v(x) := c;
        u := u + d;
    else
        vv := [op(vv), v(x)];
    end if;
    uu := [op(uu), u];
end do
n := nops(tspan);
P_7 := plots:-listplot([seq([tspan[q], vv[q]], q = 1 .. n)], color = black, style = line, title = "Low-Threshold Spiking (LTS)", labels = ["Time (ms)", "Membrane Potential (mV)"], labeldirections = ['horizontal', 'vertical']);
plots:-display(P_7);