%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%            Run Display version 3                           %
%              Anita and Piotr                               %
%              24th August 2010                              %
% Program simulating cleaving of peptide by  TEV protease.   %
% Volume around the cell is calculated using CFU results     %
% from Imperial iGEM 2008. Negelcts transient diffusion      %
% because presence other cells is taken into consideration.  %
% The average resulting diffusive flux is zero.              %
% Assumptions:
% -addition of TEV to solution is done within minutes after  %
% setting the concentration of cells in water to 5e8 CFU/ml  %
% so the cells will not divide a lot by the time TEV is added%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Calculate Diffusion coeffiecient

n=0.0008909;    % viscosity coefficient for H20 in 27C (from H&M notes)
k_b=1.38065e-23;% [J/K] Boltzmann constant
T=310;          % 37C in Kalvin
r=0.000000005;   % [m] average protein radius having 360 residues (E.Coli Statistics)
D=((T*k_b)/(6*pi*n*r)); % Deduced diff. coeff. is twice as big as the smallest in the paper
D=1e-10 % [m^2/s] http://www.life.illinois.edu/crofts/bioph354/diffusion1.html
%D=1.07e-10; % [m^2/s] average foragarose gel from literature for pH 5.6 [http://www.sciencedirect.com/science?_ob=ArticleURL&_udi=B6V5N-4B3MXDC-2&_user=217827&_coverDate=07%2F01%2F2004&_rdoc=1&_fmt=high&_orig=search&_origin=search&_sort=d&_docanchor=&view=c&_searchStrId=1448514265&_rerunOrigin=google&_acct=C000011279&_version=1&_urlVersion=0&_userid=217827&md5=72717f9adb05920bdeb57c11f0d0c6a9&searchtype=a#toc9]

%% Calculate the outer volume using Fick's first law

x= 1.259921e-5; % [m] Diffusion distance determined from CFU/ml value
diff_time=(x^2)/(2*D) %[s] Diffusion time

%% Control volume -> delta_V

V=x^3 %[m^3] - Volume around the cell
v=2.794e-18; % [m^3] - volume of bacterium
delta_V=V-v % [m^3] - control volume

%% Initial conditions and constants
global factor
global s_prot
global deg_c
factor=v/delta_V % unitless
s_prot = 4.13e-8; % production rate of protein to be displayed
deg_c = 0.000289; %common degradation constant due to cell division
%% Run the loading of the cell depending on production and degradation rate
t=0:10:25000;
[t,protein] = ode45(@protein_display,t, 0);
protein_max=max(protein) %The 'protein' variable is already adjusted by the factor
% in the protein_display m file


%% Run the ODE
t0=0;
tf=3000;
c_enz= 4e-6; % concentration of enzyme cleaving protein
options = odeset('NonNegative',[1 1 1 1],'RelTol',1e-12);
[time,p] = ode15s(@display_3, [t0 tf], [c_enz protein_max 0 0], options);


figure
subplot(2,2,1), plot(time,p(:,1))
title('Production of TEV')
xlabel('time [s]')
ylabel('concentration [mol*dm^-3]')

subplot(2,2,2), plot(time,p(:,2))
title('Production of display protein')
xlabel('time [s]')
ylabel('concentration [mol*dm^-3]')

subplot(2,2,3), plot(time,p(:,3))
title('Production of TEV-display protein')
xlabel('time [s]')
ylabel('concentration [mol*dm^-3]')


subplot(2,2,4), plot(time,p(:,4))
title('Production of AIP')
xlabel('time [s]')
ylabel('concentration [mol*dm^-3]')
% Drawing receptor activation threshold line
hold on
thres=zeros(numel(time),1);
thres(:,1)=4.4658e-9;
plot(time,thres,'r') %Plot the threshold value
hold off