%here we estimate the critical IP3 concentration 
% for propagation

  
function calcium_speed
 global p
 dc = 0.0001;
c = [0:dc:1.5];
 for j = 1:100
     
 p = 0.153+j*0.02/100;
 pj(j) = p;
 flux = current(c);
 ndx = find(flux(1:end-1).*flux(2:end)<=0)
 if length(ndx ==3)
  
 %apply trapezoidal rule
 Flux_int = sum(flux(1:ndx(3)-1)) + flux(ndx(3))/2;
 Flux_int = Flux_int*dc;
 int(j) = Flux_int;
 end
 end

pdx = find(int(1:end-1).*int(2:end)<=0);
p = pj(pdx)
flux = current(c);
ndx = find(flux(1:end-1).*flux(2:end)<=0)
 
figure(1)
plot(c(1:ndx(3)),flux(1:ndx(3)),c(1:ndx(3)),zeros(1,ndx(3)),'--','linewidth',2)
xlabel('c','fontsize',20)
ylabel('Total Flux J_{tot}','fontsize',20)

function c_flux = current(c)
global p
Vp = 0.9;
kf = 1.11;
Jer = 0;
Kp = 0.1;
K1 = 0.13;
K2 = 1.64/20;

ct = 2;
gam = 5.5;

ce = gam*(ct-c);
Jserca = Vp*c.^2./(Kp^2+c.^2);
JIPR = (p*c./((p+K1)*(c+K2))).^3;
Jtot = (kf*JIPR+Jer).*(ce-c)-Jserca;

c_flux = Jtot;


