% Model name = Snifferomyces % Function: Time-dependant sensitivity analysis function Snifferomyces close all %Initial conditions vector x0=zeros(18,1); x0(1) = 20; %Ligand x0(2) = 0.2; %Ste12 x0(3) = 0; %Ste12a x0(4) = 1.667; %Ste2 x0(5) = 0.0; %D x0(6) = 0.686; %Fus3 x0(7) = 0.0; %Fus3a x0(8) = 1.667; %Gabc x0(9) = 0.0; %GaGDP x0(10) = 0.0; %GaGTP x0(11) = 0.0; %Gbc x0(12) = 0.505; %Sst2 x0(13) = 0.0; %Sst2a x0(14) = 0.0; %Ste2a x0(15) = 0.036; %C x0(16) = 0; %mRNA x0(17) = 0; %mature x0(18) = 0; %nascent k=zeros(22,1); global k; % parameters of rates k(1) = 0.071; %kSte2a k(2) = 2.457; %kSte2ia k(3) = 1.199; %kSte2adg k(4) = 0.724; %kgpa k(5) = 1.259; %kgpia k(6) = 1.124; %kgpiasst2 k(7) = 0.002; %Kgpformation k(8) = 0.021; %kDsyn k(9) = 0.004; %kDdg k(10) = 0.261; %kDdgsst2a k(11) = 0.818; %kfus3a k(12) = 0; %kfus3ia k(13) = 0.017; %kfus3iasst2 k(14) = 0.342; %ksst2a k(15) = 0.643; %ksst2ia k(16) = 300; %kste12a k(17) = 0.17; %kste12ia k(18) = 4.0E-6; %reaction_reaction_18_k k(19) = 0.00214; %reaction_reaction_19_k1 k(20) = 2.0; %reaction_reaction_20_k k(21) = 9.625E-5; %reaction_reaction_21_k1 k(22) = 9.7E-6; %reaction_reaction_22_kmgfpdg [tout,x_std] = odesolver(x0); %------------------------ %sensitivity on initial concentrations co_i = 0.2; % magnitude of perturbation to initial concentrations di = zeros(18,1); % obtain the maximum of each species during time course [tout,x] = odesolver(x0); for z =1:18 di(z) = max(x(:,z)); end di = co_i*di; p_i=zeros(length(tout),18); for i = 1:18 x0(i)=x0(i) + di(i); [tout,x] = odesolver(x0); % p_r(:,i) = abs(x(:,17) - x_std(:,17)); %unnormalized for j = 1:length(tout) %normalized p_i(j,i) = (abs(x(j,17) - x_std(j,17))/x_std(j,17))/co_i; end x0(i)=x0(i) - di(i); end tout = tout/3600; figure semilogy(tout,p_i(:,1),'k*',... tout,p_i(:,2),'k+',... tout,p_i(:,3),'k:',... tout,p_i(:,4),'k--',... tout,p_i(:,5),'c*',... tout,p_i(:,6),'c+',... tout,p_i(:,7),'c:',... tout,p_i(:,8),'c--',... tout,p_i(:,9),'b*',... tout,p_i(:,10),'b+',... tout,p_i(:,11),'b:',... tout,p_i(:,12),'b--',... tout,p_i(:,13),'g*',... tout,p_i(:,14),'g+',... tout,p_i(:,15),'g:',... tout,p_i(:,16),'g--',... tout,p_i(:,17),'m*',... tout,p_i(:,18),'m+'); legend('Ligand','Ste12a','Ste12ia','Ste2','D','Fus3','Fus3a','Gabc',... 'GaGDP','GaGTP','Gbc','Sst2','Sst2a','Ste2a','C','mRNA','mature-GFP',... 'nascent-GFP'); title('Sensitivity analysis of initial concentration(input 20 perturbation 20%)'); xlabel('Time (hours)'); ylabel('Sensitivity of initial concentrations'); %------------------------ %sensitivity on parameters co_p = 0.8; %magnitude of perturbation to parameters dp = co_p*k; p_r=zeros(length(tout),22); % x0(2) = 0.21; for i = 1:22 k(i)=k(i) + dp(i); [tout,x] = odesolver(x0); % p_r(:,i) = abs(x(:,17) - x_std(:,17)); %unnormalized for j = 1:length(tout) %normalized p_r(j,i) = (abs(x(j,17) - x_std(j,17))/x_std(j,17))/co_p; end k(i)=k(i) - dp(i); end tout = tout/3600; figure semilogy(tout,p_r(:,1),'ko',... tout,p_r(:,2),'k+',... tout,p_r(:,3),'k:',... tout,p_r(:,4),'k--',... tout,p_r(:,5),'c*',... tout,p_r(:,6),'c+',... tout,p_r(:,7),'c:',... tout,p_r(:,8),'c--',... tout,p_r(:,9),'b*',... tout,p_r(:,10),'b+',... tout,p_r(:,11),'b:',... tout,p_r(:,12),'b--',... tout,p_r(:,13),'g*',... tout,p_r(:,14),'g+',... tout,p_r(:,15),'g:',... tout,p_r(:,16),'g--',... tout,p_r(:,17),'m*',... tout,p_r(:,18),'m+',... tout,p_r(:,19),'m:',... tout,p_r(:,20),'m--',... tout,p_r(:,21),'r*',... tout,p_r(:,22),'r+'); legend('kSte2a 1','kSte2ia 2','kSte2adg 3','kgpa 4','kgpia 5','kgpiasst2 6',... 'kgpformation 7', 'kDsyn 8','kDdg 9','kDdgsst2a 10','kfus3a 11',... 'kfus3ia 12','kfus3iasst2 13','ksst2a 14','ksst2ia 15','kste12a 16',... 'kste12ia 17','kmRNAsyn 18','kmRNAdg 19',... 'nGFPsyn 20','mGFPsyn 21','mGFPdg 22'); title('Sensitivity analysis of rates(input 2 perturbation 80%)'); xlabel('Time (hours)'); ylabel('Sensitivity of rates'); end function [tout,x] = odesolver(x0) %Creating linespace t=linspace(0,21600,100); %600000 is one week, 1800 is 30 mins, 10800 is 3 hours %remember to change the title of plot as well %Solving equations [tout,x] = ode23s(@f,t,x0); %ploting the results % hold on % plot(tout,x(:,1),'r--',... % tout,x(:,2),'k',... % tout,x(:,3),'k',... % tout,x(:,4),'k',... % tout,x(:,5),'k',... % tout,x(:,6),'k',... % tout,x(:,7),'g--',... % tout,x(:,8),'k',... % tout,x(:,9),'k',... % tout,x(:,10),'k',... % tout,x(:,11),'k',... % tout,x(:,12),'k',... % tout,x(:,13),'k',... % tout,x(:,14),'k',... % tout,x(:,15),'k',... % tout,x(:,16),'k',... % tout,x(:,17),'b--',... % tout,x(:,18),'c--'); % legend('Ligand','Ste12a','Ste12ia','Ste2','D','Fus3','Fus3a','Gabc',... % 'GaGDP','GaGTP','Gbc','Sst2','Sst2a','Ste2a','C','mRNA','mature-GFP',... % 'nascent-GFP'); end function xdot=f(~,x) % Compartment: id = compartment_1, name = Extracellular, constant compartment_Extracellular=1.0; % Compartment: id = compartment_2, name = Cell, constant compartment_Cell=1.0; % Reaction: id = reaction_1, name = Receptor Activation global k; % Local Parameter: id = kSte2a, name = kSte2a reaction_reaction_1=compartment_Cell*function_function_1(x(1), k(1), x(4)); % Reaction: id = reaction_2, name = Receptor Inactivation % Local Parameter: id = kSte2ia, name = kSte2ia reaction_reaction_2=compartment_Cell*function_function_2(x(14), k(2)); % Reaction: id = reaction_4, name = Receptor-Ligand Degradation % Local Parameter: id = kSte2adg, name = kSte2adg reaction_reaction_4=compartment_Cell*function_function_4(x(14), k(3)); % Reaction: id = reaction_5, name = G-Protein Activation % Local Parameter: id = kgpa, name = kgpa reaction_reaction_5=compartment_Cell*function_function_5(x(8), x(14), k(4)); % Reaction: id = reaction_6, name = G-Protein Inactivation % Local Parameter: id = kgpia, name = kgpia reaction_reaction_6=compartment_Cell*function_function_7(x(10), k(5)); % Reaction: id = reaction_7, name = G-Protein Inactivation Under phosphorylated Sst2 % Local Parameter: id = kgpiasst2, name = kgpiasst2 reaction_reaction_7=compartment_Cell*function_function_8(x(10), x(13), k(6)); % Reaction: id = reaction_8, name = Heterotrimeric G-Protein Formation % Local Parameter: id = Kgpformation, name = Kgpformation reaction_reaction_8=compartment_Cell*function_function_9(x(11), x(9), k(7)); % Reaction: id = reaction_9, name = complex D synthesis % Local Parameter: id = kDsyn, name = kDsyn reaction_reaction_9=compartment_Cell*function_function_10(x(15), x(11), k(8)); % Reaction: id = reaction_10, name = complex D degradation % Local Parameter: id = kDdg, name = kDdg reaction_reaction_10=compartment_Cell*function_function_11(x(5), k(9)); % Reaction: id = reaction_11, name = complex D degradation under phosphorylated Sst2 % Local Parameter: id = kDdgsst2a, name = kDdgsst2a reaction_reaction_11=compartment_Cell*function_function_12(x(5), x(13), k(10)); % Reaction: id = reaction_12, name = Fus3 Activation % Local Parameter: id = kfus3a, name = kfus3a reaction_reaction_12=compartment_Cell*function_function_13(x(6), k(11), x(5)); % Reaction: id = reaction_13, name = Fus3 Inactivation % Local Parameter: id = kfus3ia, name = kfus3ia reaction_reaction_13=compartment_Cell*function_function_14(x(7), k(12)); % Reaction: id = reaction_14, name = Fus3 Inactivation by Sst2 % Local Parameter: id = kfus3iasst2, name = kfus3iasst2 reaction_reaction_14=compartment_Cell*function_function_15(x(7), k(13)); % Reaction: id = reaction_15, name = Activation of Sst2 % Local Parameter: id = ksst2a, name = ksst2a reaction_reaction_15=compartment_Cell*function_function_16(x(12), x(7), k(14)); % Reaction: id = reaction_16, name = Inactivation of Sst2 % Local Parameter: id = ksst2ia, name = ksst2ia reaction_reaction_16=compartment_Cell*function_function_17(x(13), k(15)); % Reaction: id = reaction_3, name = Activation of Ste12 % Local Parameter: id = kste12a, name = kste12a reaction_reaction_3=compartment_Extracellular*function_function_3(x(2), x(7), k(16)); % Reaction: id = reaction_17, name = Inactivation of Ste12 % Local Parameter: id = kste12ia, name = kste12ia reaction_reaction_17=compartment_Extracellular*function_function_6(x(3), k(17)); % Reaction: id = reaction_18, name = GFP-mRNA Synthesis % Local Parameter: id = k, name = k reaction_reaction_18=compartment_Cell*function_function_18(x(3), k(18)); % Reaction: id = reaction_19, name = GFP-mRNA Degradation % Local Parameter: id = k1, name = k1 reaction_reaction_19=compartment_Cell*k(19)*x(16); % Reaction: id = reaction_20, name = Nascent GFP Synthesis % Local Parameter: id = k, name = k reaction_reaction_20=compartment_Cell*function_function_18(x(16), k(20)); % Reaction: id = reaction_21, name = Mature GFP Synthesis % Local Parameter: id = k1, name = k1 reaction_reaction_21=compartment_Cell*k(21)*x(18); % Reaction: id = reaction_22, name = Mature GFP degradation % Local Parameter: id = kmgfpdg, name = kmgfpdg reaction_reaction_22=compartment_Cell*function_function_19(x(17), k(22)); xdot=zeros(18,1); % Species: id = species_3, name = L % Warning species is not changed by either rules or reactions xdot(1) = 0; % Species: id = species_14, name = Ste12, affected by kineticLaw xdot(2) = (1/(compartment_Extracellular))*((-1.0 * reaction_reaction_3) + ( 1.0 * reaction_reaction_17)); % Species: id = species_15, name = Ste12a, affected by kineticLaw xdot(3) = (1/(compartment_Extracellular))*(( 1.0 * reaction_reaction_3) + (-1.0 * reaction_reaction_17)); % Species: id = species_1, name = Ste2, affected by kineticLaw xdot(4) = (1/(compartment_Cell))*((-1.0 * reaction_reaction_1) + ( 1.0 * reaction_reaction_2)); % Species: id = species_10, name = D, affected by kineticLaw xdot(5) = (1/(compartment_Cell))*(( 1.0 * reaction_reaction_9) + (-1.0 * reaction_reaction_10) + (-1.0 * reaction_reaction_11)); % Species: id = species_11, name = Fus3, affected by kineticLaw xdot(6) = (1/(compartment_Cell))*((-1.0 * reaction_reaction_12) + ( 1.0 * reaction_reaction_13) + ( 1.0 * reaction_reaction_14)); % Species: id = species_12, name = Fus3a, affected by kineticLaw xdot(7) = (1/(compartment_Cell))*(( 1.0 * reaction_reaction_12) + (-1.0 * reaction_reaction_13) + (-1.0 * reaction_reaction_14)); % Species: id = species_5, name = Gabc, affected by kineticLaw xdot(8) = (1/(compartment_Cell))*((-1.0 * reaction_reaction_5) + ( 1.0 * reaction_reaction_8)); % Species: id = species_7, name = GaGDP, affected by kineticLaw xdot(9) = (1/(compartment_Cell))*(( 1.0 * reaction_reaction_6) + ( 1.0 * reaction_reaction_7) + (-1.0 * reaction_reaction_8)); % Species: id = species_4, name = GaGTP, affected by kineticLaw xdot(10) = (1/(compartment_Cell))*(( 1.0 * reaction_reaction_5) + (-1.0 * reaction_reaction_6) + (-1.0 * reaction_reaction_7)); % Species: id = species_6, name = Gbc, affected by kineticLaw xdot(11) = (1/(compartment_Cell))*(( 1.0 * reaction_reaction_5) + (-1.0 * reaction_reaction_8) + (-1.0 * reaction_reaction_9) + ( 1.0 * reaction_reaction_10) + ( 1.0 * reaction_reaction_11)); % Species: id = species_8, name = Sst2, affected by kineticLaw xdot(12) = (1/(compartment_Cell))*((-1.0 * reaction_reaction_15) + ( 1.0 * reaction_reaction_16)); % Species: id = species_13, name = Sst2a, affected by kineticLaw xdot(13) = (1/(compartment_Cell))*(( 1.0 * reaction_reaction_15) + (-1.0 * reaction_reaction_16)); % Species: id = species_2, name = Ste2a, affected by kineticLaw xdot(14) = (1/(compartment_Cell))*(( 1.0 * reaction_reaction_1) + (-1.0 * reaction_reaction_2) + (-1.0 * reaction_reaction_4)); % Species: id = species_9, name = C, affected by kineticLaw xdot(15) = (1/(compartment_Cell))*((-1.0 * reaction_reaction_9) + ( 1.0 * reaction_reaction_10) + ( 1.0 * reaction_reaction_11)); % Species: id = species_16, name = GFP-mRNA, affected by kineticLaw xdot(16) = (1/(compartment_Cell))*(( 1.0 * reaction_reaction_18) + (-1.0 * reaction_reaction_19)); % Species: id = species_17, name = mature-GFP, affected by kineticLaw xdot(17) = (1/(compartment_Cell))*(( 1.0 * reaction_reaction_21) + (-1.0 * reaction_reaction_22)); % Species: id = species_18, name = nascent-GFP, affected by kineticLaw xdot(18) = (1/(compartment_Cell))*(( 1.0 * reaction_reaction_20) + (-1.0 * reaction_reaction_21)); end function z=function_function_1(L,kSte2a,Ste2), z=(Ste2*L*kSte2a);end function z=function_function_2(Ste2a,kSte2ia), z=(Ste2a*kSte2ia);end function z=function_function_4(Ste2a,kSte2adg), z=(Ste2a*kSte2adg);end function z=function_function_5(Gabc,Ste2a,kgpa), z=(Gabc*Ste2a*kgpa);end function z=function_function_7(GaGTP,kgpia), z=(GaGTP*kgpia);end function z=function_function_8(GaGTP,Sst2a,kgpiasst2), z=(GaGTP*Sst2a*kgpiasst2);end function z=function_function_9(GaGDP,Gbc,Kgpformation), z=(GaGDP*Gbc*Kgpformation);end function z=function_function_10(C,Gbc,kDsyn), z=(C*Gbc*kDsyn);end function z=function_function_11(D,kDdg), z=(D*kDdg);end function z=function_function_12(D,Sst2a,kDdgsst2a), z=(D*Sst2a*kDdgsst2a);end function z=function_function_13(Fus3,kfus3a,D), z=(Fus3*D*kfus3a);end function z=function_function_14(Fus3a,kfus3ia), z=(Fus3a*kfus3ia);end function z=function_function_15(Fus3a,kfus3iasst2), z=(Fus3a*kfus3iasst2);end function z=function_function_16(Sst2,Fus3a,ksst2a), z=(Sst2*Fus3a*ksst2a);end function z=function_function_17(Sst2a,ksst2ia), z=(Sst2a*ksst2ia);end function z=function_function_19(GFP,kmgfpdg), z=(GFP*kmgfpdg);end function z=function_function_3(Ste12,Fus3a,kste12a), z=(Ste12*Fus3a*kste12a);end function z=function_function_6(Ste12a,kste12ia), z=(Ste12a*kste12ia);end function z=function_function_18(Activator,k), z=(k*Activator);end