%%  Program for Figure 5.13, Li and Rinzel Model (closed-cell gonadotroph model
%%  Carol Lucas
%%  calls:  function xdot=lrclosefun(t,x,flag,I,L,Ve)


fi=0.01;Vi=4;P=26640;
Ki=1.0;Ka=0.4;
Ke=0.2;A=0.5;Kd=0.4;
sigma=0.185;L=.37;I=.6;Ct=2;Ve=400;

C=[.015:.01:1];Ce=(Ct-C)/sigma;
hnull=Kd./(C+Kd);

Iray=[.4 .8 2];
cnull=zeros(3,length(C));

mess=' Lots of pauses, hit return to keep going'
for k=1:3
    I=Iray(k);
    vars=['[IP3]= ' num2str(I) ' L= ' num2str(L) '  Vserca= ' num2str(Ve)];
    [t,x]=ode45('lrclosefun',[0 250],[.2;.8],[],I,L,Ve);
    figure(1);clf;plot(t,x);
    xlabel('t');ylabel('C and h'); title(vars);pause;
    figure(2);clf;plot(x(:,1),x(:,2));
    xlabel('[C]');ylabel('h'); title(['phase plane' vars]);pause;
    cnull(k,:)=(((Ve*C.*C./((Ke*Ke+C.*C).*(Ce - C))-L)/P).^(1/3)).*(I+Ki).*(C+Ka)./(I*C);
    figure(3);clf;plot(C,hnull,C,real(cnull(k,:)),x(:,1),x(:,2));axis([0 1 .2 1.4]);
    xlabel('[C]');ylabel('h'); title(['phase plane  ' vars]);pause;
end
figure(4);clf;plot(C,hnull,C,real(cnull));axis([0 1 .2 1.4]);
    xlabel('[C]');ylabel('h'); title(['F5.13A: [IP3] Variations = '  num2str(Iray)]);pause;

Veray=[400 120];I=.4;L=.37;
cnull=zeros(2,length(C));
for k=1:2
    Ve=Veray(k);
    vars=['[IP3]= ' num2str(I) ' L= ' num2str(L) '  Vserca= ' num2str(Ve)];
    [t,x]=ode45('lrclosefun',[0 250],[.2;.8],[],I,L,Ve);
    figure(1);clf;plot(t,x);
    xlabel('t');ylabel('C and h'); title(vars);pause(1);
    figure(2);clf;plot(x(:,1),x(:,2));
    xlabel('[C]');ylabel('h'); title(['phase plane' vars]);pause(1);
    cnull(k,:)=(((Ve*C.*C./((Ke*Ke+C.*C).*(Ce - C))-L)/P).^(1/3)).*(I+Ki).*(C+Ka)./(I*C);
    figure(3);clf;plot(C,hnull,C,real(cnull(k,:)),x(:,1),x(:,2));axis([0 1 .2 1.4]);
    xlabel('[C]');ylabel('h'); title(['phase plane  ' vars]);pause;
end
figure(4);clf;plot(C,hnull,C,real(cnull));axis([0 1 .2 1.4]);
    xlabel('[C]');ylabel('h'); title(['F5.13B: Vserca Variations = ' num2str(Veray)]);pause;

Lray=[3.7 .37];Ve=400;I=.6;
cnull=zeros(2,length(C));
for k=1:2
    L=Lray(k);
    vars=['[IP3]= ' num2str(I) ' L= ' num2str(L) '  Vserca= ' num2str(Ve)];
    [t,x]=ode45('lrclosefun',[0 250],[.2;.8],[],I,L,Ve);
    figure(1);clf;plot(t,x);
    xlabel('t');ylabel('C and h'); title(vars);pause;
    figure(2);clf;plot(x(:,1),x(:,2));
    xlabel('[C]');ylabel('h'); title(['phase plane' vars]);pause;
    cnull(k,:)=(((Ve*C.*C./((Ke*Ke+C.*C).*(Ce - C))-L)/P).^(1/3)).*(I+Ki).*(C+Ka)./(I*C);
    figure(3);clf;plot(C,hnull,C,real(cnull(k,:)),x(:,1),x(:,2));axis([0 1 .4 1.4]);
    xlabel('[C]');ylabel('h'); title(['phase plane  ' vars]);pause;
end
figure(4);clf;plot(C,hnull,C,real(cnull));axis([0 1 .4 1.4]);
    xlabel('[C]');ylabel('h'); title(['F5.13C: L Variations = ' num2str(Lray)]);pause;