%% ODE for full gonadotroph model (Figure 5.16)
%% Carol Lucas

function xdot=lrpmlfun(t,x)
% The Morris-Lecar + Li-Rinzel model

V=x(1,1);
n=x(2,1);
C=x(3,1);
h=x(4,1);
Ct=x(5,1);

% PM parameters
% units: V = mV; Kca = uM; Cm = uF/cm^2; g = uS/cm^2; phi = 1/s 
 Vk=-85;Vca=120;Kca=0.5;Cm=1;
 gk=20;gca=20;gkca=8;phi=12;
 V1=-3;V2=30;V3=-20;V4=30;

% units: V = mV; Kca = uM; Cm = uF/cm^2; g = uS/cm^2; phi = 1/s 
 Vk=-85;Vca=120;Kca=0.5;Cm=1;
 gk=20;gca=20;gkca=8;phi=12;
 V1=-3;V2=30;V3=-20;V4=30;

% ER parameters
fi=0.01;
% Vi = pL
Vi=4;
% L;P = pL/s
L=1.48;P=2960;
% I; C; Ce; Ct;I; Ki; Ka; Ke; Kd; Kp = uM 
I=0.9;Ki=0.0;Ka=0.4;
% Ve; Vp = aMol/s [sic]
Ve=480;Ke=0.2;A=2.0;Kd=0.4;

% sigma;eps; fi = unitless
sigma=0.185;eps=0.01;Vp=400;Kp=0.3;
% alpha = aMol cm^2/nC (approximately Acell/(2*Faraday) in those units)
alpha=.2;



% PM functions
minf=.5*(1+tanh((V-V1)/V2));
ninf=.5*(1+tanh((V-V3)/V4));
tau= 1/(phi*cosh((V-V3)/(2*V4)));
I_ca=gca*minf*(V-Vca);
I_k=gk*n*(V-Vk);
I_kca=gkca*(V-Vk)*C^4/(C^4+Kca^4);
Ce = (Ct - C)/sigma;
jin=-alpha*I_ca;
% PM equations
xdot(1,1)= (-I_k - I_ca - I_kca)/Cm;
xdot(2,1)= (ninf-n)/tau;

% ER equations
xdot(3,1) = fi/Vi*( (L + P*( (I*C*h)/( (I+Ki)*(C+Ka) ) )^3 )*(Ce - C) - Ve*C*C/(Ke*Ke+C*C) + eps*(-alpha*I_ca - Vp*C*C/(Kp*Kp + C*C)) );
xdot(4,1) = A*( Kd - (C + Kd)*h );
xdot(5,1) = fi/Vi*eps*(-alpha*I_ca - Vp*C*C/(Kp*Kp + C*C));