calciumregMatlab3

function [xp] = f(t,x)

ca = x(1);
pth = x(2);
dthh = x(3);
ctn = x(4);

IC50 = 46.8;
EC50 = 53.2;

PTHmax = 1.27e-4;
dpth = log(2)/(4/(60*24));
dctn = log(2)/(10/(24*60));
CTNmax = 1e-5;
DHCCmean = 2.1632e-2;
ddhcc = log(2)/(5/24);
pth_mean = 3.58e-5;

K_PTH = PTHmax*dpth;
K_CTN = CTNmax*dctn;
K_DHCC = (DHCCmean*ddhcc)/pth_mean;

K1 = 10000; %pth
K2 = 1; %ctn
K3 = (9.245e-5)*K2 - 0.00165*K1 + 2323.62; %dthh
K4 = 10;
K5 = 0.1;

pthp = K_PTH * (1/(1+(ca/IC50))) - dpth*pth;
dthhp = K_DHCC*pth - ddhcc*dthh;
ctnp = K_CTN * ((ca/EC50)/(1+(ca/EC50))) - dctn * ctn;
cap = K1*pth - K2*ctn + K3*dthh - ca*K4*((ctn/pth)/(1+(ctn/pth))) - K5*ca;

xp = [cap;pthp;dthhp;ctnp];

------- different file: ode_driver.m ---------
[t,y] = ode45(@f,[0,10],[95,35.8e-6,2.163e-2,2e-6]);
subplot(4,1,1);
plot(t,y(:,1));
legend('cap');
subplot(4,1,2);
plot(t,y(:,2));
legend('PTH');
subplot(4,1,3);
plot(t,y(:,3));
legend('dthh');
subplot(4,1,4);
plot(t,y(:,4));
legend('ctn');