function task7CaProbeCameleonnProbeConcOscillation()
%Model based on BJ 77,1244-1256(1999),Thomas Hofer
time=0.001:0.01:500;
s0=zeros(1,4); % Ca conc in plasma membrane, Ca conc in ER, Probe conc in plasma membrane, Ca/Probe complex in plasma membrane
s0(2)=2.2; %ERのカルシウム濃度の初期値を2.2uMに設定
s0(3)=2; %細胞質のプローブ濃度の初期値を2uMに設定 ←順番に2, 20, 200(uM)の値を入れてください

figure;
hold on;
param0=[0.12,0.12,0.02,0.2,4,9,3.6,2.5,2,2,0.1,0.3,0.4,0.2,0.2,0.4,40,0.02,0,9.4*10^2,180];
[t,time_course]=ode15s(@(t,s)ODE(t,s,param0),time,s0);
plot(t,time_course);
legend('PlasmaCa','ERCa','Probe','ProbeComplex');
end

function dsdt=ODE(t,s,param)
x=s(1);%initial Ca conc in plasma memberane
y=s(3);%initial Probe conc in plasma membrane
z=s(2);%initial Ca conc in ER
w=s(4);%initial Ca/Probe complex in plasma membrane

K3=param(1);%Affinity of ER Ca influx pomp to Ca
K4=param(2);%Affinity of ER Ca efflux pomp to Ca
rou=param(3);
new0=param(4);
newc=param(5);
new3=param(6);
new4=param(7);
K0=param(8);
P=param(9);
alpha=param(10);
beta=param(11);
d1=param(12);
d2=param(13);
d3=param(14);
dp=param(15);
da=param(16);
k1=param(17);
k2=param(18);
gamma=param(19);
kon=param(20);
koff=param(21);

Jin=new0+newc*P*(K0+P)^(-1);%Ca influx from Plasma membrane
Jout=new4*x^2*(K4^2+x^2)^(-1);%Ca efflux from Plasma membrane
Jrel=(k1*(d2*P*x*(d1+P)*(d3+P)^(-1))^3*((dp+P)^3*(da+x)^3*(d2*(d1+P)*(d3+P)^(-1)+x)^3)^(-1)+k2)*(z-(1+beta)*x)*beta^(-1);%Ca release from ER
Jserca=new3*x^2*(K3^2+x^2)^(-1); %Ca uptake to ER
Jprobe=-kon*x*y+koff*w; %Complex net reduction rate

v1=rou*(Jin-Jout+Jprobe);%d(Ca)/dt from plasma
v2=rou*alpha*(Jrel-Jserca);%d(Ca)/dt from ER

dsdt(1,:)=v1+v2;
dsdt(2,:)=v1;
dsdt(3,:)=Jprobe;
dsdt(4,:)=-Jprobe;
end
