Matlab code - Model B

function p_prime = model_b(time,p)
p_prime = zeros(12,1);

% Initial Reactants
% T + sTa + sTb + sD
%% Units
% [p_dprime] = [mol*dm^-3*s^-1]
% [s_d] = [mol*dm^-3*s^-1]
% [d_d] = [s^-1]
% [p_d] = [mol*dm^-3]

%% Production and Degradation constants
s_t = 2.3808*10^(-10);    %production rate constant of T
d_t = 0.000289;    %degradation rate constant of T

s_sTa = 1.2606*10^(-10);  %production rate constant of sTa
d_sTa = 0.000289;  %degradation rate constant of sTa

s_sTb = 1.2606*10^(-10);  %production rate constant of sTb
d_sTb = 0.000289;  %degradation rate constant of sTb

d_Tsa = 0.000289;  %degradation rate constant od Tsa
%d_TsTa = 1;%degradation rate constant od T-sTa

d_Tsb = 0.000289;  %degradation rate constant od Tsb
%d_TsTb = 1;%degradation rate constant od T-sTb

s_sd = 2.4998*10^(-9);   %production rate constant of sD
d_sd = 0.000289;   %degradation rate constant of sD
%d_TssD = 1;%degradation rate constant od Ts-sD

d_ts = 0.000289;   %degradation rate cpnstant of Ts 

d_d = 0.000289;    %degradation rate constant of D
%% Kinetic rate constants
k011 = 1000; % T + sD -> T-sD 
k012 = 10; % T + sD <- T-sD
k013 = 0.2; % T-sD -> T + D
k051 = 1000; % T + sTa -> T-sTa
k052 = 10; % T + sTa <- T-sTa
k053 = 0.2; % T-sTa -> T + Tsa
k071 = 1000; % T + sTb -> T-sTb
k072 = 10; % T + sTb <- T-sTb
k073 = 0.2; % T-sTb -> T + Tsb
k141 = 0.4*k011; % Ts + sD -> Ts-sD
k142 = 1.4*k012; % Ts + sD <- Ts-sD
k143 = 0.4*k012; % Ts-sD -> Ts + D
k081 = 10; % Tsa + Tsb -> Ts
k082 = 0.1; % Tsa + Tsb <- Ts
%% Abbreviations
%p(1) = [T] = [TEV]
%p(2) = [sD] = [splitDioxygenase]
%p(3) = [T-sD] = [TEV-splitDioxygenase]
%p(4) = [D] = [Dioxygenase]
%p(5) = [sTa] = [splitTEV(a)]
%p(6) = [sTb] = [splitTEV(b)]
%p(7) = [Tsa] = [activated TEV(a)]
%p(8) = [Tsb] = [activated TEV(b)]
%p(9) = [Ts] = [activated TEV]
%p(10) = [T-sTa] = [TEV-splitTEV(a)]
%p(11) = [T-sTb] = [TEV-splitTEV(b)]
%p(12) = [Ts-sD] = [transition state Ts-sD]
%p(13) = [Tsa-Tsb] = [activated TEV(a)-activated TEV(b)]
%% p_prime(1) - TEV production rate
% Terms coming coming from reactions involving: T-sTa T-sD T-sTb
p_prime(1) = -k011*p(1)*p(2) + (k012+k013)*p(3)-k051*p(1)*p(5) + (k052+k053)*p(10) -k071*p(1)*p(6) + (k072+k073)*p(11)+ s_t - d_t*p(1);   %Production rate of TEV

%% p_prime(2) - Production rate of splitDioxygenase
% Terms coming coming from reactions involving: T-sD,Ts-sD
p_prime(2) = -k011*p(1)*p(2) + k012*p(3)-k141*p(9)*p(2) + k142*p(12)+ s_sd - d_sd*p(2);

%% p_prime(3) - Production rate of TEV-splitDioxygenase 
% Terms coming coming from reactions involving: T-sD
p_prime(3) = k011*p(1)*p(2) - (k012+k013)*p(3); %- d_tsd*p(3); - no degradation term for transition state

%% p_prime(4) - Production rate of Dioxygenase
% Terms coming coming from reactions involving: T-sD Ts-sD
p_prime(4) = k013*p(3)+ k143*p(12) - d_d*p(4);

%% p_prime(5) - Production rate of sTa
% Terms coming coming from reactions involving: T-sTa
p_prime(5) = -k051*p(1)*p(5) + k052*p(10) + s_sTa - d_sTa*p(5);

%% p_prime(6) - Production rate of sTb
% Terms coming coming from reactions involving: T-sTb
p_prime(6) = -k071*p(1)*p(6) + k072*p(11) + s_sTb - d_sTb*p(6);

%% p_prime(7) - Production rate of Tsa
% Terms coming coming from reactions involving: T-sTa Tsa-Tsb
p_prime(7) = k053*p(10)-k081*p(7)*p(8) + k082*p(9) - d_Tsa*p(7); 

%% p_prime(8) - Production rate of Tsb
% Terms coming coming from reactions involving: T-sTb Tsa-Tsb
p_prime(8) = k073*p(11) -k081*p(7)*p(8) + k082*p(9)- d_Tsb*p(8); 

%% p_prime(9) - Production rate of activated TEV [Ts]
% Terms coming coming from reactions involving: Tsa-Tsb
p_prime(9) = k081*p(7)*p(8) - k082*p(9) -k141*p(9)*p(2) + (k142+k143)*p(12) - d_ts*p(9);

%% p_prime(10) - Production rate of TEV-splitTEV(a) [T-sTa]
% Terms coming coming from reactions involving:  T-sTa
p_prime(10) = k051*p(1)*p(5) - (k052+k053)*p(10); % - d_TsTa*p(10) no degradation term for transition state

%% p_prime(11) - Production rate of TEV-splitTEV(b) [T-sTb]
% Terms coming coming from reactions involving: T-sTb
p_prime(11) = k071*p(1)*p(6) - (k072+k073)*p(11); % - d_TsTb*p(11) no degradation term for transition state

%% p_prime(12) - Production rate of transition state Ts-sD[Ts-sD]
% Terms coming coming from reactions involving: Ts-sD
p_prime(12) = k141*p(9)*p(2) - (k142+k143)*p(12); % - d_TssD*p(12) no degradation term for transition state

%% p_prime(13) - Production rate of activated TEV(a)-activated TEV(b) [Tsa-Tsb]
% This reaction is not ensymatic -> This Ttransition State doers not exist

%% 1) T+sTa+sTb+sD <-> T-sD -> T+D+sTa+sTb
%%Reaction constants
%k011 = 1;
%k012 = 1;
%k013 = 1;
%%Production and degradation rate constants
%s_t = 1;    %production rate constant of TEV
%d_t = 1;    %degradation rate constant of TEV
%s_sd = 1;   %production rate constant of splitDioxygenase
%d_sd = 1;   %degradation rate constant od splitDioxygenase
%d_tsd = 1;  %degradation rate constant of TEV-splitDioxygenase
%d_d = 1;    %degradation rate constant of Dioxygenase
%Equations
%p_prime(1) = -k011*p(1)*p(2) + (k012+k013)*p(3) + s_t - d_t*p(1);   %Production rate of TEV
%p_prime(2) = -k011*p(1)*p(2) + k012*p(3)+ s_sd - d_sd*p(2);         %Production rate of splitDioxygenase   
%p_prime(3) = k011*p(1)*p(2) - (k012+k013)*p(3) - d_tsd*p(3);        %Production rate of TEV-splitDioxygenase 
%p_prime(4) = k013*p(3) - d_d*p(4);                                  %Production rate of Dioxygenase

%% 2) T+D+sTa+sTb <-> T-sTa -> T+Tsa+D+sTb
%Reaction constants
%k021 = 1;
%k022 = 1;
%k023 = 1;
%Production and degradation rate constants
%d_sta = 1;  %degradation rate constant of splitTEV(a)
%d_tsta = 1; %degradation rate constant of TEV-splitTEV(a)
%d_tsa = 1;  %degradation rate constant of activated TEV(a)
%Equations
%p_prime(1) = -k021*p(1)*p(5) + (k022+k023)*p(10) + s_t - d_t*p(1);  %Production rate of TEV
%p_prime(5) = -k021*p(1)*p(5) + k022*p(10) - d_sta*p(5);             %Production rate of splitTEV(a)
%p_prime(10) = k021*p(1)*p(5) - (k022+k023)*p(10) - d_tsta*p(10);     %Production rate of TEV-splitTEV(a)
%p_prime(7) = k023*p(10) - d_tsa*p(7);                               %Production rate of activated TEV(a)

%% 3) T+D+Tsa+sTb <-> T-sTb -> T+Tsb+Tsa+D 
%Reaction constants
%k031 = 1;
%k032 = 1;
%k033 = 1;
%%Production and degradation rate constants
%d_stb = 1;  %degradation rate constant of splitTEV(b)
%d_tstb = 1; %degradation rate constant of TEV-splitTEV(b)
%d_tsb = 1;  %degradation rate constant of activated TEV(b)
%%Equations
%p_prime(1) = -k031*p(1)*p(6) + (k032+k033)*p(11) + s_t - d_t*p(1);  %Production rate of TEV
%p_prime(6) = -k031*p(1)*p(6) + k032*p(11) - d_stb*p(6);            %Production rate of splitTEV(b)
%p_prime(11) = k031*p(1)*p(6) - (k032+k033)*p(11) - d_tstb*p(11);    %Production rate of TEV-splitTEV(b)
%p_prime(8) = k033*p(11) - d_tsb*p(8);                              %Production rate of activated TEV(b)%

%% 4) T+Tsa+Tsb+D <-> Tsa-Tsb -> Ts+T+D
% Note that it's actually: Tsa + Tsb <-> Ts
%Reaction constants
%k081 = 1;
%k082 = 1;

%Production rate of Tsa
%p_prime(7) = -k081*p(7)*p(8) + k082*p(9) - d_Tsa*p(7);     
%Production rate of Tsb
%p_prime(8) = -k081*p(7)*p(8) + k082*p(9) - d_Tsa*p(8);
%Production rate of Ts
%p_prime(9) = k081*p(7)*p(8) - k082*p(9) - d_ts*p(9);
%% 5) T + sTa + sTb + sD -- T-sTa --> T + Tsa + sTb + sD
%Reaction constants
%k051 = 1;
%k052 = 1;
%k053 = 1;
%Production rate of sTa
%p_prime(5) = -k051*p(1)*p(5) + k052*p(10) + s_sTa - d_sTa*p(5);     
%Production rate of T
%p_prime(1) = -k051*p(1)*p(5) + (k052+k053)*p(10) + s_t - d_t*p(1);          
%Production rate of T-sTa
%p_prime(10) = k051*p(1)*p(5) - (k052+k053)*p(10) - d_TsTa*p(10);      
%Production rate of Tsa
%p_prime(7) = k053*p(10) - d_Tsa*p(7); 
%% 6) T + Tsa + sTb + sD --  T-sD --> T + Tsa + sTb + D
%Reaction constants
%k061 = 1;
%k062 = 1;
%k063 = 1;
%Production rate of sD
%p_prime(2) = -k061*p(1)*p(2) + k062*p(3) + s_sd - d_sd*p(2);     
%%Production rate of T
%p_prime(1) = -k061*p(1)*p(2) + (k062+k063)*p(3) + s_t - d_t*p(1);          
%Production rate of T-sD
%p_prime(3) = k061*p(1)*p(2) - (k062+k063)*p(3) - d_tsd*p(3);      
%Production rate of D
%p_prime(4) = k063*p(3) - d_d*p(4); 
%% 7) T + Tsa + sTb + D --  T-sTb --> T + Tsa + Tsb + D
%Reaction constants
%k071 = 1;
%k072 = 1;
%k073 = 1;
%Production rate of sTb
%p_prime(6) = -k071*p(1)*p(6) + k072*p(11) + s_sTb - d_sTb*p(6);     
%Production rate of T
%p_prime(1) = -k071*p(1)*p(6) + (k072+k073)*p(11) + s_t - d_t*p(1);          
%Production rate of T-sTb
%p_prime(11) = k071*p(1)*p(6) - (k072+k073)*p(11) - d_TsTb*p(11);      
%Production rate of Tsb
%p_prime(8) = k073*p(11) - d_Tsb*p(8); 

%% 8) T + Tsa + Tsb + D --  Tsa-Tsb --> T + Ts + D
% Note that it's actually: Tsa + Tsb <-> Ts
%Reaction constants
%k081 = 1;
%k082 = 1;

%Production rate of Tsa
%p_prime(7) = -k081*p(7)*p(8) + k082*p(9) - d_Tsa*p(7);     
%Production rate of Tsb
%p_prime(8) = -k081*p(7)*p(8) + k082*p(9) - d_Tsb*p(8);
%Production rate of Ts
%p_prime(9) = k081*p(7)*p(8) - k082*p(9) - d_Ts*p(9);
%% 9) T + sTa + sTb + sD -- T-sTb --> T + Tsb + sTa + sD
%Reaction constants
%k091 = 1;
%k092 = 1;
%k093 = 1;
%Production rate of sTb
%p_prime(6) = -k091*p(1)*p(6) + k092*p(11) + s_sTb - d_sTb*p(6);     
%Production rate of T
%p_prime(1) = -k091*p(1)*p(6) + (k092+k093)*p(11) + s_t - d_t*p(1);          
%Production rate of T-sTb
%p_prime(11) = k091*p(1)*p(6) - (k092+k093)*p(11) - d_TsTb*p(11);      
%Production rate of Tsb
%p_prime(8) = k093*p(11) - d_Tsb*p(8); 
%% 10) T + Tsb + sTa + sD --  T-sD -->  T + Tsb + sTa + D
%Reaction constants
%k101 = 1;
%k102 = 1;
%k103 = 1;
%Production rate of sD
%p_prime(2) = -k101*p(1)*p(2) + k102*p(3) + s_sd - d_sd*p(2);     
%%Production rate of T
%p_prime(1) = -k101*p(1)*p(2) + (k102+k103)*p(3) + s_t - d_t*p(1);          
%Production rate of T-sD
%p_prime(3) = k101*p(1)*p(2) - (k102+k103)*p(3) - d_tsd*p(3);      
%Production rate of D
%p_prime(4) = k103*p(3) - d_d*p(4); 
%% 11) T + Tsb + sTa + D --  T-sTa --> T + Tsb + Tsa + D
%Reaction constants
%k051 = 1;
%k052 = 1;
%k053 = 1;
%Production rate of sTa
%p_prime(5) = -k051*p(1)*p(5) + k052*p(10) + s_sTa - d_sTa*p(5);     
%Production rate of T
%p_prime(1) = -k051*p(1)*p(5) + (k052+k053)*p(10) + s_t - d_t*p(1);          
%Production rate of T-sTa
%p_prime(10) = k051*p(1)*p(5) - (k052+k053)*p(10) - d_TsTa*p(10);      
%Production rate of Tsa
%p_prime(7) = k053*p(10) - d_Tsa*p(7); 

%% 12) T + Tsb + Tsa + D --  Tsa-Tsb --> T + Ts + D
% Note that it's actually: Tsa + Tsb <-> Ts
%Reaction constants
%k081 = 1;
%k082 = 1;

%Production rate of Tsa
%p_prime(7) = -k081*p(7)*p(8) + k082*p(9) - d_Tsa*p(7);     
%Production rate of Tsb
%p_prime(8) = -k081*p(7)*p(8) + k082*p(9) - d_Tsa*p(8);
%Production rate of Ts
%p_prime(9) = k081*p(7)*p(8) - k082*p(9) - d_ts*p(9);


%% 13) T + sTa + sTb + sD -- T-sTa --> T + Tsa + sTb + sD
%Reaction constants
%k051 = 1;
%k052 = 1;
%k053 = 1;
%Production rate of sTa
%p_prime(5) = -k051*p(1)*p(5) + k052*p(10) + s_sTa - d_sTa*p(5);     
%Production rate of T
%p_prime(1) = -k051*p(1)*p(5) + (k052+k053)*p(10) + s_t - d_t*p(1);          
%Production rate of T-sTa
%p_prime(10) = k051*p(1)*p(5) - (k052+k053)*p(10) - d_TsTa*p(10);      
%Production rate of Tsa
%p_prime(7) = k053*p(10) - d_Tsa*p(7); 
%% 14) Ts + T + sD -- Ts-sD --> Ts + T + D
%Reaction constants
%k141 = 1;
%k142 = 1;
%k143 = 1;
%Production rate of sD
%p_prime(2) = -k141*p(9)*p(2) + k142*p(12) + s_sd - d_sd*p(2);     
%Production rate of Ts
%p_prime(9) = -k141*p(9)*p(2) + (k142+k143)*p(12) - d_ts*p(9);          
%Production rate of Ts-sD
%p_prime(12) = k141*p(9)*p(2) - (k142+k143)*p(12) - d_TsTb*p(12);      
%Production rate of D
%p_prime(4) = k143*p(12) - d_d*p(4);

end


time = 0:0.1:20;
[time,p] = ode45(@model_b,time,[0 1 0 0 0 0 0 0 0 0 0 0]);

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

subplot(3,2,2), plot(time,p(:,5))
title('Production of splitTEV(a)')
xlabel('time [s]')
ylabel('concentration')

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

subplot(3,2,4), plot(time,p(:,9))
title('Production of activated TEV (TEVs)')
xlabel('time [s]')
ylabel('concentration')

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

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