function sd_prime = display_3(time, p)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             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.             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% Units
% [p_dprime] = [mol*dm^-3*s^-1]
% [s_d] = [mol*dm^-3*s^-1]
% [d_d] = [s^-1]
% [p_d] = [mol*dm^-3]

%% Reaction constants
global k1 
global k2 
global factor
k1 = 10000;
k2 = 620;
k3 = 0.16;

deg_c = 0.000289; %common degradation constant due to cell division

d_t = deg_c; %degradation rate constant of TEV
d_disp_prot = deg_c;
d_aip = deg_c;


s = 4.13e-8*factor;        %production rate of protein

sd_prime = zeros(4,1);
%p(1) = [T] = [TEV]
%p(2) = [disp_prot] = [display protein]
%p(3) = [T-disp_prot] = [TEV-display protein]
%p(4) = [AIP] = [AIP]

%% Equations
%Production rate of TEV
sd_prime(1) = -k1*p(1)*p(2) + (k2+k3)*p(3) - d_t*p(1);     
%Production rate of disp_prot
sd_prime(2) = -k1*p(1)*p(2) + k2*p(3) +s - d_disp_prot*p(2);        
%Production rate of T-disp_prot
sd_prime(3) = k1*p(1)*p(2) - (k2+k3)*p(3);      
%Production rate of AIP
sd_prime(4) = k3*p(3) - d_aip*p(4);                

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%            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.             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 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=1.7-10; % [m^2/s] average from literature for pH 5.6

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

x=1.259921e-5; % [m] Diffusion distance
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 = 1e-10 %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)

c_p=protein_max*factor % [mol dm^(-3)] this is the maximum cell can produce
% so this will be a initia value for our ODE model


%% Run the ODE
time = 0:1:10000;
global k1 
global k2 
c_enz= 4e-4;
options = odeset('NonNegative',[1 1 1 1]);
[time,p] = ode15s(@display_3, time, [c_enz c_p 0 0], options);

Km=k2/k1 %Report value of Km into command window

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]')
hold on
plot(time,4.4658e-9) %Plot the threshold value
hold off