(*Values for the physical quantities relevant in the experiment performed by Hodgkin and Huxley with the squid giant axon. \
 Conductances are given in mho/cm^2, voltages in mV*)
 GNa = 120.0*10^-3; (*Sodium Maximum Conductance*)
 GK = 36*10^-3; (*Potassium Maximum Conductance*)
 GL = 0.3*10^-3; (*Leakage Current Maximum Conductance*)
 VNa = 115; (*Sodium Diffusion Potential*)
 VK = -12; (*Potassium Diffusion Potential*)
 VL = 11; (*Leakage Current Diffusion Potential*)
 a = 238*10^-4; (*Radius of the axon, given in cm*)
Cap = 10^-3; (*Capacitance per unit area of the axon membrane, given in mF/cm^2*)
r = 2*10^4; (*Longitudinal resistance per unit length of the axon, given in ohm/cm*)
m0 = 0.053;
h0 = 0.59;
n0 = 0.318;
V0 = 0;

(* Define utility functions *)

(*Voltage-dependent coefficients to be used in the equations for the dynamics of gate-variables m,h and n*)
am[V_] := (0.1 (25 - V))/(E^((25 - V)/10) - 1);
bm[V_] := 4 E^(-(V/18));
ah[V_] := 0.07 E^(-(V/20));
bh[V_] := 1/(E^((30 - V)/10) + 1);
an[V_] := (0.01 (10 - V))/(E^((10 - V)/10) - 1);
bn[V_] := 0.125 E^(-(V/80));
(*Equations for the dynamics of gate-variables m,h and n. dm, dh and dn stand for dm/dt, dh/dt and dn/dt, respectively*)
dm[t_] := am[V[t]] (1 - m[t]) - bm[V[t]] m[t];
dh[t_] := ah[V[t]] (1 - h[t]) - bh[V[t]] h[t];
dn[t_] := an[V[t]] (1 - n[t]) - bn[V[t]] n[t];
(*Aperture of the sodium and potassium channels*)
ChannelsNa[u_] := (m[u])^3 h[u];
ChannelsK[u_] := (n[u])^4;
(*Current density of sodium, potassium and leakage current, together with the resultant ion current density*)
JNa[u_] := GNa ChannelsNa[u] (V[u] - VNa);
JK[u_] := GK ChannelsK[u] (V[u] - VK);
JL[u_] := GL (V[u] - VL);
Jion[u_] := JNa[u] + JK[u] + JL[u];
(*UPS stands for 'Unique Point Solution', since it finds and plots a numerical solution for the behaviour of voltage, currents and gates in a chosen point along the axon. J0 is the stimulus current,dV is the function dV/dt and s is the numerical solution returned by the Mathematica function NDSolve*)

(* Initialize stimulus and plot range values *)

stimulus = 0.02;(* stimulus strength (mA/cm^2) *)
start = 10;(* start of stimulus (ms) *)
duration = 40;(* stimulus duration (ms) *)
domaini = 0;(* initial time *)
domainf = 50;(* final time *)
rangel = -20;(* minimum range *)
rangeh = 120;(* maximum range *)

(* Numerically solve the DE and plot the potential versus time *)

Block[{J0, dV, s},
J0 = Function[t, If[start <= t <= duration, stimulus, 0]];
dV = Function[t, (J0[t] - Jion[t])/Cap];
s = NDSolve[{
    V'[t] == dV[t],
    m'[t] == dm[t],
    h'[t] == dh[t],
    n'[t] == dn[t],
    V[0] == V0,
    m[0] == m0,
    h[0] == h0,
    n[0] == n0},
   {V, m, h, n},
   {t, domainf}
   ];
Plot[
  Evaluate[{V[t]} /. s],
  {t, domaini, domainf},
  PlotRange -> {rangel, rangeh},
  AxesLabel -> {"time (ms)", "membrane voltage (mV)"},
  AxesOrigin -> {domaini, rangel},
  GridLines -> Automatic, GridLinesStyle -> Dashed,
  ImageSize -> 450, PlotStyle -> {Blue}]
]