% phase portrait for quorum sensing system
function quorum_sensing_pp
global s g K alp d rfac 

set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7);

 
u = [0:.001:7];
 
s = 0.05; % basal production rate
K = 1; % production rate K_d
g = 1.5; % production rate scale factor
d = 1; % transmembrane permeability rate
alp = 0.0375; % extracellular decay rate



rlist = [0.007;0.02;0.035];
aylist=[0.5;0.8;1.4];

axlist=[2;2;3];
formatSpecF = '%6.3f\n';
f = s + g*u.^2./(K^2+u.^2);

 
for j = 1:3
    r = rlist(j);
  
    rfac = r*d/(1-r);
    sl = alp/rfac+1;
    figure(j)
    plot(u,u-f/d,'--',u,u/sl,'--',  'linewidth',2)
    ax = axlist(j);
    ay = aylist(j);
    %axis([0 ax 0 ay])
    axis([0 3 0 1.4])
    xlabel('u','fontsize',20)
    ylabel('v','fontsize',20)
    hold on 
 title(strcat('\rho = ',sprintf(formatSpecF,r)),'fontsize',20)
switch j
    case 1
  annotation('arrow',[ 0.65,0.55],[.44,0.35])
 annotation('arrow',[ 0.16,0.23],[.29,0.28])
    case 2
         annotation('arrow',[ 0.37,0.47],[.29,0.37])
  annotation('arrow',[ 0.15,0.22],[.29,0.28])
    case 3
         annotation('arrow',[ 0.53,0.6],[.45,0.55])
  annotation('arrow',[ 0.15,0.22],[.29,0.28])
end

 v0list=[0.2;0.15;0.25];
 for jj = 1:length(v0list)
    
 u0=0.;
 w0 = v0list(jj);
s0 = [u0;w0];
tstep = 1;
t_end = 250;

tspan = [ 0:tstep:t_end];
[T,S] = ode23s(@deRHS,tspan,s0); 

 plot( S(:,1),S(:,2),'linewidth',2)
 
  
 end
 
 
end
 

rho = alp*(d*u-f)./(d*f-f*alp+d*alp*u);
v=u-f/d;
figure(4)

plot(rho,u, 'linewidth',2)
axis([0 0.04 0 3])
ylabel('Intracellular Autoinducer','fontsize',20)
xlabel('Cell Density','fontsize',20)
text(0.015,2.5,'c>0','fontsize',20)
text(0.005,2.5,'c<0','fontsize',20)
hold on
 
 rmin = min(rho);
 rmax = max(rho);
 drho = (rmax-rmin)/100;
 rj = [rmin:drho:rmax];
 rk = 0;
 
 for j = 1:length(rj)
     r = rj(j);
geV = gev(u,r);
% identify three zeros
sgc = geV(1:end-1).*geV(2:end);
idx = find(sgc<=0);
 if (length(idx)==3)
     rk = rk+1;
     rhok(rk) = r;
    tst = abs(geV(idx))>abs(geV(idx+1));
    idx = idx+tst;
    Uss = u(idx) % these are the three zeros
    cf =(1-r)*alp/(d*r)+1;
    Vss = Uss/cf;
   
% create the speed sign expressions

Ex1 = -r*d*(Vss-Uss).^2/(2*(1-r))-alp*Vss.^2/2;
DelU = (Uss(3)-Uss(1))/100;
uss = [Uss(1):DelU:Uss(3)];
Fus =  s + g*uss.^2./(K^2+uss.^2);    
Ex2 = r*(sum(Fus)-Fus(1)/2-Fus(end)/2)*DelU/(1-r);
Texp(rk) = Ex1(3)-Ex1(1)+Ex2;
end
 end

 tcrss = Texp(1:end-1).*Texp(2:end);
 tj = find(tcrss<=0)
 Texp(tj)
 Texp(tj+1)
 
 tst = abs(Texp(tj))>abs(Texp(tj+1))
 rhocrss = rhok(tj+tst);
figure(4)
plot(rhocrss*ones(length(u),1),u,'--','linewidth',2)
hold off

 

function s_prime=deRHS(t,y)
global s g K alp d rfac
u = y(1);
v = y(2);
 f = s+g*u^2/(K^2+u^2);
Fu = f+d*(v-u);
Fv = -alp*v+rfac*(u-v);

s_prime = [Fu Fv]';

function evalf=gev(u,r)
global s g K alp d rfac

f = s + g*u.^2./(K^2+u.^2);

cf = 1/((1-r)*alp/(d*r)+1)-1;
evalf = cf*u+f/d;
