Matlab code - Model A

function sd_prime = model_a(time, p)

%% 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
k1 = 10000;
k2 = 620;
k3 = 0.16;           %k3 = kcat;

s_t = 2.3808*1e-10; %production rate constant of TEV
d_t = 0.000289;     %degradation rate constant of TEV

s_sd = 2.4998*1e-0; %production rate constant of splitDioxygenase
d_sd = 0.000289;    %degradation rate constant od splitDioxygenase

d_tsd = 0.000289;   %degradation rate constant of TEV-splitDioxygenase

d_d = 0.000289;     %degradation rate constant of Dioxygenase

sd_prime = zeros(4,1);
%p(1) = [T] = [TEV]
%p(2) = [sD] = [splitDioxygenase]
%p(3) = [T-sD] = [TEV-splitDioxygenase]
%p(4) = [D] = [Dioxygenase]

%% Equations
%Production rate of TEV
sd_prime(1) = -k1*p(1)*p(2) + (k2+k3)*p(3) + s_t - d_t*p(1);     
%Production rate of sD
sd_prime(2) = -k1*p(1)*p(2) + k2*p(3) + s_sd - d_sd*p(2);        
%Production rate of T-sD
sd_prime(3) = k1*p(1)*p(2) - (k2+k3)*p(3) - d_tsd*p(3);      
%Production rate of D
sd_prime(4) = k3*p(3) - d_d*p(4);                

end


time = 0:1:5000;
options = odeset('NonNegative', [1 1 1 1]);
[time,p] = ode15s(@model_a,time, [0 1e-2 0 0], options);

figure
subplot(2,2,1), plot(time,p(:,1))
title('Production of TEV')
xlabel('time [s]')
ylabel('concentration')

subplot(2,2,2), plot(time,p(:,2))
title('Production of splitDioxygenase')
xlabel('time [s]')
ylabel('concentration')

subplot(2,2,3), plot(time,p(:,3))
title('Production of TEV-splitDioxygenase')
xlabel('time [s]')
ylabel('concentration')

subplot(2,2,4), plot(time,p(:,4))
title('Production of Dioxygenase')
xlabel('time [s]')
ylabel('concentration')