function MnPDesensitization

%% define parameters
p = 5;    % production rate of MnP per unit bacteria per day = 5 unit/min
b0 = 100; % starting concentration of bacteria               = 100 unit
v = 0.879;% volume                                           = 0.879 m^3
j = 0.19; % flow rate                                        = 0.18m^3/min

k1f = 2;  % forward reaction of Enzyme-substrate formation   = unit^-1 * min^-1
k1r = 0.5;% reverse reaction of Enzyme-substrate formation   = min^-1
k2 = 1;   % forward reaction of breaking down of nylon       = min^-1

n = 1000; % Nylon concentration of waster water flowing into chamber = unit

params = {p,v,b0,j,k1f,k1r,k2};

%% Run single simulation
c0 = [0; b0; n; 0];                                                        %starting conc of [MnP; bacteria;nylon;nylon-MnP];
tspan = [0 60];                                                            % 1 hour time period
options = [];
[t,c] = ode15s(@MnPfunc,tspan,c0,options,params);
for tstep=1:length(t),                                                     % loop to evaluate free receptor, given solution (t,y)
    [~,algvars(tstep,:)] = MnPfunc(t(tstep),c(tstep,:),params);
end

plot(t,c(:,1),'b.-',t,c(:,3),'r.-');
xlabel('Time (hr)'); ylabel('Concentration(\unit)'); legend('MnP\unit','Nylon\unit');
% plot(t,c(:,1),'b-.',t,c(:,3),'r-o',t,c(:,2),'g-o');
% xlabel('Time (min)'); ylabel('Concentration(\unit)'); legend('MnP\unit','Nylon\ unit','bacteria\ unit');

%% Run MnP conc response over a range of total Bacteria concentrations, b0
p = 5;    % production rate of MnP per unit bacteria per day = 5 unit/min
b0 = 100; % starting concentration of bacteria               = 100 unit
v = 0.879;% volume                                           = 0.879 m^3
j = 0.19; % flow rate                                        = 0.18m^3/min

k1f = 2;  % forward reaction of Enzyme-substrate formation   = unit^-1 * min^-1
k1r = 0.5;% reverse reaction of Enzyme-substrate formation   = min^-1
k2 = 1;   % forward reaction of breaking down of nylon       = min^-1
n = 100;  % Nylon concentration of waster water flowing into chamber = unit

params = {p,v,b0,j,k1f,k1r,k2};
paramRange = 10.^[0:.1:2];
for i=1:length(paramRange)
    b0 = paramRange(i);
    params = {p,v,b0,j,k1f,k1r,k2};
    c0 = [0; p*b0; n; 0];
    tspan = [0 60];
    options = [];
    [t,c] = ode15s(@MnPfunc,tspan,c0,options,params);
    MnPfinal(i) = c(end,1);
    Nylfinal(i) = c(end,3);
%   BaCfinal(i) = c(end,2);
end

figure(2)
subplot(1,2,1)
semilogx(paramRange,MnPfinal,'b-o');
% plot(paramRange,MnPfinal,'ro',paramRange,BaCfinal,'ro');
xlabel('Starting number of bacteria (\unit)'); ylabel('Concentration after 60 hours(\unit)');
legend('MnP\unit');


%% Run Nyl conc response over a range of total Bacteria concentrations, b0
p = 5;    % production rate of MnP per unit bacteria per day = 5 unit/min
b0 = 100; % starting concentration of bacteria               = 100 unit
v = 0.879;% volume                                           = 0.879 m^3
j = 0.19; % flow rate                                        = 0.18m^3/min

k1f = 2;  % forward reaction of Enzyme-substrate formation   = unit^-1 * min^-1
k1r = 0.5;% reverse reaction of Enzyme-substrate formation   = min^-1
k2 = 1;   % forward reaction of breaking down of nylon       = min^-1
n = 100;  % Nylon concentration of waster water flowing into chamber = unit

params = {p,v,b0,j,k1f,k1r,k2};
paramRange = [0:.1:2];
for i=1:length(paramRange)
    b0 = paramRange(i);
    params = {p,v,b0,j,k1f,k1r,k2};
    c0 = [0; p*b0; n; 0];
    tspan = [0 60];
    options = [];
    [t,c] = ode15s(@MnPfunc,tspan,c0,options,params);
    MnPfinal(i) = c(end,1);
    Nylfinal(i) = c(end,3);
%   BaCfinal(i) = c(end,2);
end
subplot(1,2,2)
plot(paramRange,Nylfinal,'r.-');
xlabel('Starting number of bacteria (\unit)'); ylabel('Concentration after 60 hours(\unit)');
legend('Nyl\unit');


%% sensitivity analysis
p = 5;             % production rate of MnP per unit bacteria per day = 5 unit/min
b0 = 100;          % starting concentration of bacteria               = 100 unit
v = 0.879;         % volume                                           = 0.879 m^3
j = 0.19*60*24;    % flow rate                                        = 0.18m^3/min

k1f = 0.2;         % forward reaction of Enzyme-substrate formation   = unit^-1 * min^-1
k1r = 5;           % reverse reaction of Enzyme-substrate formation   = min^-1
k2 = 10;           % forward reaction of breaking down of nylon       = min^-1

n = 100;   % Nylon concentration of waster water flowing into chamber = unit
delta = 1.01;
params = {p,v,b0,j,k1f,k1r,k2};
c0 = [0; p*b0; n; 0];
tspan = [0 60];
options = [];
[t,c] = ode15s(@MnPfunc,tspan,c0,options,params);
cend0 = c(end,:);

for i=1:length(params)
    newParams = params;
    newParams{i} = delta*params{i};   
    [t,c] = ode15s(@MnPfunc,tspan,c0,options,newParams);
    cfinal(i,:) = c(end,:); 
    S(i,:) = (cfinal(i,:)-cend0)./cend0  .*  params{i}./(newParams{i}-params{i});
end

% plot sensitivity coefficients for MnP
figure(5);
bar(S(:,1));
xlabel('Parameter'); ylabel('MnP Sensitivity Coefficient');
set(gca,'XTickLabel',{'p','v','b0','j','k1f','k1r','k2'},'fontsize',8);

% plot sensitivity matrix with only MnP
figure(6);
imagesc(S');
colormap('gray');
colorbar;
caxis([-1 1]);
xlabel('Perturbation'); ylabel('Output');
set(gca,'XTickLabel',{'p','v','b0','j','k1f','k1r','k2'},'fontsize',8);
set(gca,'YTickLabel',{'','MnP','','Bac','','Nylon','','Nylon-MnP complex',''},'fontsize',8);
