function [dVdt, dmdt, dhdt] = timeDerivatives(V, m, h, Ve, param)
% function [dVdt, dmdt, dhdt] = timeDerivatives(V, m, h, Ve, param)
% Calculates the time derivatives of transmembrane potential V and gating
% parameters m and h. Ve is the extracellular potential.
% All input (except param) and output arguments are N x 1 vectors, where N
% is the number of Ranvier nodes.
%
% param is a struct containing parameter values. They and their units are:
%  param.gL, S/m2, leakage conductance of active membrane per unit area.
%  param.EL, V, leakage reversal potential
%  param.gNa, S/m2, sodium conductance per unit area when ion channels fully open
%  param.ENa, V, sodium Nernst potential
%  param.sigmai, S/m,  intracellular conductivity
%  param.cm, F/m2, active membrane capacitance per unit area
%  param.d, m, axon diameter
%  param.l, m, node length
%  param.L, m, distance between the nodes of Ranvier
%  param.cm_mye, F/m2, myelin capacitance per unit area
%  param.gL_mye, S/m2, leakage conductance of myelin per unit area
%
% Copyright 2024 Ilkka Laakso (ilkka.laakso@aalto.fi)
% License: The MIT license (https://opensource.org/license/mit)
% doi: 10.5281/zenodo.12652749

GA = pi*(param.d/2)^2 / param.L * param.sigmai; % S, intracellular conductance
Cm = pi * param.d * param.l * param.cm; % F, capacitance of nodal membrane

% Add the effect of myelin capacitance:
Cm_eff = Cm + pi * param.d * (param.L-param.l) * param.cm_mye; % F, capacitance of nodal membrane

% Sodium and leak conductances [S]
GNa = pi * param.d * param.l * param.gNa * h .* m.^2;
GL  = pi * param.d * param.l * param.gL;

% Add the effect of myelin leak conductance:
GL_eff = GL + pi * param.d * (param.L - param.l) * param.gL_mye;

% Time derivative of voltage
dVdt = zeros(size(V));
n = 2 : length(V)-1;
dVdt(n) = 1/Cm_eff * ( ...
  GA * ( V(n-1) - 2*V(n) + V(n+1) + Ve(n-1) - 2*Ve(n) + Ve(n+1) ) ...
  - (GNa(n) .* (V(n) - param.ENa) + GL_eff * (V(n) - param.EL) ) );


% Transmembrane potential in millivolts
V_mV = 1e3 * V;

% Time derivative of gating parameter m
alpham = (126 + 0.363*V_mV) ./ ( 1 + exp( -(V_mV + 49)/5.3 ));
betam = alpham ./ exp( (V_mV + 56.2)/4.17 );
dmdt = 1e3 * (alpham.*(1-m) - betam.*m);

% Time derivative of gating parameter h
betah = 15.6 ./ (1 + exp( -(V_mV + 56)/10 ));
alphah = betah ./ exp( (V_mV + 74.5)/5 );
dhdt = 1e3 * (alphah.*(1-h) - betah.*h);

end