%  Keizer-Levine Reduced Model, Fig 5.6 B
%% Carol Lucas

kap=1500; kam=28.8; kbp=1500; kbm=385.9; kcp=1.75; kcm=0.1;

Ka=(kam/kap)^(1/4);Kb=(kbm/kbp)^(1/3);Kc=kcm/kcp;

c=.1;  %% Assume this is the starting C

%%  Then we can find what winf=Po1+Po2+C1 is and assume it's the most
%%  we will have to work with and won't change quickly since C2 is slow

winfss=(1+(Ka/c)^4+(c/Kb)^3)/(1+(1/Kc)+(Ka/c)^4+(c/Kb)^3);

c=[.1:.01:2.5];

peakPo= winfss*(1 + (c/Kb).^3)./(1+(Ka./c).^4+(c/Kb).^3);

winf=(1+(Ka./c).^4+(c/Kb).^3)./(1+(1/Kc)+(Ka./c).^4+(c/Kb).^3);
ssPo=winf.*(1 + (c/Kb).^3)./(1+(Ka./c).^4+(c/Kb).^3);

figure(1);clf;
plot(c,ssPo,c,peakPo);xlabel('[Ca2+]i');ylabel('Po');title('Figure 5.6 (B)');