% This file is solving quorum sensing traveling waves with diffusive spread

clear
 %moviename = ['MOVIE.qt']
%MakeQTMovie('start', moviename)
set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7);

 

L = 10; % length of the domain
d = 0.001; % diffusion coefficient
 
N=250;   % number of spatial grid points
M =25;  %number of time steps
kstep = 250;

h = L/N;
dt = 5*h;
 
%parameters
s = 0.05;
K = 1;
g = 1.5;
del = 1.0;
dk = 0.0375;

rlist = [0.007;0.02;0.039];
rho = rlist(2);
%rho=0.0118
formatSpecF = '%6.3f\n';
 
axlist=[2;2;4];
% draw a phase portrait (nullclines)
u = [0:.01:8];
 
f = s + g*u.^2./(K^2+u.^2);
    r = rho;
  
    rfac = r*del/(1-r);
    sl = dk/rfac+1;
    figure(2)
    plot(u,u-f/del,u,u/sl, 'linewidth',2)
   ax = axlist(2);
    axis([0 ax 0 ax/2])
    xlabel('u','fontsize',20)
    ylabel('v','fontsize',20)
    hold on 

    scal = d*dt/h^2;
    X = h*(1:N)';
 
 V = 2*sech(2*(X-L/2)).^2;  %initial autoinducer
 U = zeros(N,1);  % initial extracellular autoinducer
  
t=0;

% uses Crank Nicolson to solve the diffusion equation

% set up matrix

Atm=(1+scal)*ones(N,1);
Atm(1,1) = 1+scal/2;
Atm(N,1) = 1+scal/2;

Am=diag(Atm)-diag(scal/2*ones(N-1,1),1)-diag(scal/2*ones(N-1,1),-1);

% The algorithm
 
for n=2:M
for k = 1:kstep
df = del*(U-V);  

F = s + g*V.^2./(K^2+V.^2)+df;  % intracellular production rate
G = -dk*U-rho*df/(1-rho); %  Extracellular AI rate
U_t =  scal/2 * ([U(2:N );U(N )] - 2*U + [U(1 );U(1:N-1 )]) +dt*G ;

U = Am\(U+U_t);  %update the extracellular AI
V = V +dt*F; % update intracellular AI
t=t + dt;
end
figure(1)
  subplot(2,1,1)
  plot(X,U,'linewidth',2)
   axis([0 L  0 2])
   set(gca,'Ytick',0:0.5:2)
  xlabel('x','fontsize',20)
   ylabel('v(x,t)','fontsize',20)
  hold on
  subplot(2,1,2)
  plot(X, V,'linewidth',2 )
   xlabel('x','fontsize',20)
  ylabel('u(x,t)','fontsize',20)
  axis([0 L  0 2])
  set(gca,'Ytick',0:0.5:2)
 % axis([0 L -.5 1])

hold on
%MakeQTMovie('addframe','moviename')
%MakeQTMovie('quality', 1.0)

  pause(0.05)
  if(n==14)
  figure(2) % make a phase portrait plot of the traveling wave solution
  plot(V,U,'k--','linewidth',2)
   axis([0 3 0 1.4])
hold off
  end
   
end

%MakeQTMovie('finish','moviename')
  
 
TL = t;
 

