function ReactionNMP_Range
%%SIMULATION FOR THE MEDIATOR NMP for NAD/mediator-ranges
clear all
clc

disp('Please enter time period for simulation');
tf = input('>> ');
disp(' ');

tspan = [0 tf];

%% Start concentrations of the network components
           
        y0 = zeros(5,1);  
        y0 (1) = 100; % [µM] NAD      
        y0 (2) = 0;   % [µM] NADH/H+
        y0 (3) = 0 ;  % [µM] Mediator oxidised;
        y0 (4) = 100; % [µM] Mediator reduced;
        y0 (5) = 0;   % [µA] Current output
        
%% kinetic constants
          
k = zeros(5,1);
k(1) = 45;  % [µM]; km for GAP3DH ,  NAD  -> NADH
k(2) = 2.68;  % [µM/s];  vmax, kcat ==0.268 vmax NAD -> NADH , E_0 == 10µM
k(3) = 0.373; % rate constant for mediator reduction: M_ox  ->  M_red , for mediator MB
k(4) = 0.15; % [1/s] rate constant for mediator oxidation at the anode
k(5) = 28.95 ; %0.02895; now: mA

%% ODE-System
            
function dydt = ODESystemFirst(y,k)

NAD = y(1);     
NADH = y(2);
Mox = y(3);
Mred = y(4);
I = y(5);

v1 = k(2)*(NAD/(k(1)+NAD));
v2 = k(3)*(NADH)*(Mox);
v3 = k(4)*(Mred);
v4 = k(5)*(Mred); % current output 
v5 = (I);
dtNAD = - v1 + v2;
dtNADH = + v1 - v2;
dtMox = - v2 + v3;
dtMred = v2 - v3;
dtI = v4-v5;

dydt = [dtNAD;dtNADH;dtMox;dtMred;dtI];

end

%% Solver and Plot characteristics        

options=odeset('NonNegative',[1:5]);

fh=@(t,y)ODESystemFirst(y,k);
[t S]=ode15s(fh,tspan,[y0],options);


% %% Current output dependant on NAD+ start concentration
NAD = [];
NADrange = [10,50,100,500];

for i =  1:length(NADrange);
    
    NAD = NADrange(i);
    y0(1) = NAD; % [µM] NAD

fh=@(t,y)ODESystemFirst(y,k);
[t S]=ode15s(fh,tspan,[y0],options);

plot(t,S(:,5))
   hold on

% plot(t,S(:,4),'r')
%    hold on

   legstr = {};
   
xlabel('Time [s]');
ylabel('Current Output [µA]');
title('');
legend('10','50','100','500');

end


for i =  1:length(NADrange);

     NADb = NADrange(i);

legstr{end+1} = ['M[NAD] =',num2str(NADb)];
legend(legstr); 
%legend off;
end

%% Current output dependant on mediator start concentration

Mox = [];
Moxrange = [10,50,100,500];

for i =  1:length(Moxrange);
    
    Mox = Moxrange(i);
    y0(1) = Mox; % [µM] mediator oxidized

fh=@(t,y)ODESystemFirst(y,k);
[t S]=ode15s(fh,tspan,[y0],options);

plot(t,S(:,5))
   hold on
   
   legstr = {};
   
xlabel('Time [s]');
ylabel('Current Output [µA]');
title('');
legend('10','50','100','500');

end


for i =  1:length(Moxrange);

     Moxb = Moxrange(i);

legstr{end+1} = ['[Mox] =',num2str(Moxb)];
legend(legstr); 
%legend off;
end




end