% Quorum sensing for a bounded cell culture
function quorum_sensing_ss
global s0 K gamma del alp D rho

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

 

% set parameter values
s0 = 0.05;  % baasal production rate
gamma = 1.5; % max production rate
alp = 0.0375; % decay rate;
% without loss of generality:
K = 1; % production rate half saturation level
del = 1; % membrane porosity
D = 1; % diffusion coeffision

% set rho
rho = 0.04;
% for this value of rho, determine the largest u0 can be
U = [0.:0.005:3];  
%F = s0 +gamma*U.^2./(K^2+U.^2);
F = prod(U);% production rate
 
Fp = prodp(U); % derivative of production rate
Up=-sqrt(alp*D)*(U-F/del)./(D*(1-Fp/del)); % This is the boundary condition at L for u;
Vp = -sqrt(alp/D)*U; % This is the boundary condition  at L for v
V = U-F/del; % this is v as a function of u
rmdx = min(find(V>0));  %v must be positive to make physical sense
 
rhot = alp*(del*U-F)./(del*F-F*alp+del*alp*U); % this is the rho vs. U curve
 
rdx = max(find(rhot<=rho)); % This establishes the upper limit for U
 % find the u nullcline
 % add the nullcline
 
Un = [0:.001:K/sqrt(3)];
F = prod(Un);
Fpp = 2*gamma*K^2*(K^2-3*Un.^2)./((K^2+Un.^2).^3);
UPsq = del*(-alp*Un + (alp/del+rho/(1-rho))*F)./(D*Fpp);
UPpos = find(UPsq>0);
UPnull = -sqrt(UPsq(UPpos));

for j = 1:rdx-rmdx+1  
 ju=j+rmdx-1;  
tstep = 0.01;
t_end = 30;
tspan = [0:tstep:t_end];
% set initial data
    u0 = U(ju);
   up0 = 0;
s = [u0;up0]; % solve the differential equations with these initial data

[T,S] = ode23s(@deRHS,tspan,s,odeset('Events',@events)); 
 u = S(:,1);
 up = S(:,2);
 L(j) = T(end); %record the length of the interval for this initial value
uL(j) = u(end);

% plot the u-up trajectory
if (fix(j/20) ==j/20)  %plot only every twentieth trajectory
 figure(1)
  
 plot(u,up,'linewidth',2)
  
 xlabel('u','fontsize',20)
 ylabel('u_x','fontsize',20)
 
  annotation('arrow',[ 0.87,0.75],[0.87,0.8])
   

hold on
 % and the v- vp trajectory
 F=prod(u);
Fp = prodp(u);

 v = u-F/del;
 vp = (1-Fp/del).*up;
 figure(2)
 plot(v,vp,'linewidth',2)
xlabel('v','fontsize',20)
 ylabel('v_x','fontsize',20)
  annotation('arrow',[ 0.75,0.68],[0.65,0.55])
 hold on
 figure(3)
 plot(T,u, 'linewidth',2)
  xlabel('x','fontsize',20)
 ylabel('u(x)','fontsize',20)
 hold on
end
end

% add the boundary condition curve
figure(1)
plot(U,Up,'--', 'linewidth',2)
axis([0 U(rdx) -1.2 0])
hold off
 
   
figure(2)
plot(U,Vp,'--','linewidth',2)
axis([0 V(rdx) -0.2 0])
 

figure(3)
plot(L,uL,'--','linewidth',2)
axis([0 20 0 3])
hold off
 

figure(5)
plot(L,U(1:rdx-rmdx+1),L,uL,'--','linewidth',2)
xlabel('L','fontsize',20)
 
axis([0 20 0 3])
legend('u(0)','u(L)','fontsize',20)
 
function s_prime=deRHS(t,y)
global s0 K gamma del alp D rho
u = y(1);
up = y(2);
 
Fu = up;

% evaluate the production function and its derivatives
F=prod(u);
Fp = prodp(u);

Fpp = 2*gamma*K^2*(K^2-3*u^2)/((K^2+u^2)^3); % second derivative wrt u

Fup = (D*Fpp*up^2/del + alp*u - (alp/del+rho/(1-rho))*F)/(D*(1-Fp/del));

s_prime = [Fu Fup]';


% event function
function [value,isterminal,direction]=events(t,y)
global s0 K gamma del alp D rho

u = y(1);
up = y(2);

F=prod(u);
Fp = prodp(u);

value = D*(1-Fp/del)*up+sqrt(alp*D)*(u-F/del);
% this is the equation for the ending condition
isterminal = 1;
direction = 0;

function F=prod(u)
global s0 K gamma del alp D rho

F = s0 + gamma*u.^2./(K^2+u.^2); % production rate

function Fp=prodp(u)
global s0 K gamma del alp D rho

Fp = 2*gamma*K^2*u./((K^2+u.^2).^2); % derivative wrt u % production rate