"
Team:DTU-Denmark/Matlab
From 2011.igem.org
Matlab
Steady-state
% Call file % Model ver. 7 % The model contains the analytical solution in regard to m*max/m for the 3 coupled differential equations, one for ChiXR (r), % ChiX (s), and ChiP (m) % This solution is only valid for p_s = 0! % In regard to Model ver. 5, now only alfa_s and alfa_r (and not the % relative amount in regard to alfa_m) are used as the axis. clc clear close all %% %% Steady state solution for ChiP %% Model Parameter % TRANSCRIPTION RATES: alfa_m= 2.57; % transcription rate of m (ChiP) From Overgaard fig1; alpha_m = beta_m * m' [nmol/time] alfa_s= 2.57; % transcription rate of m (ChiX) guess! % Degradation and dilution rates: beta_m= 0.0257; % background rate of degradation and dilution of m (ChiP mRNA) From Overgaard [1/min] beta_s= 0.0257; % background rate of degradation and dilution of s (ChiX sRNA) From Overgaard [1/min] beta_r= 0.0257; % background rate of degradation and dilution of r (ChiXR sRNA), guess - no data [1/min] % Constants: k_s= 0.000820; % From Overgaard [1/(nmol*min)] k_r= 0.00820; % Guess - no data [1/(nmol*min)] lambda_s = beta_s*beta_m/k_s; lambda_r = beta_s*beta_r/k_r; m_ss_max = alfa_m/beta_m; % the maximum steady state level of m [nmol] % Probability of co-degradation: p_r= 1; % Probability that r is codegraded with s %% Simulation N_s = 201; % Number of x values (s). N_r = 201; % Number of y values (r). parameter =[lambda_s lambda_r p_r]; [X, Y, SS_cat SS_stoch] = ss_simu_ver7(parameter, N_s, N_r); %% Plotting figure_no =1; % Figure number for easy naming of the figures. figure name=[' catalytic: \lambda_s=', num2str(lambda_s), ' \lambda_r=', num2str(lambda_r), ' p_r=', num2str(p_r)]; % part of the title ss_plot_single_ver7( X, Y, SS_cat, name ); saveas(gcf,['Figure',int2str(figure_no),'_SS_1_subplot.jpg']) figure name=[' stochiometric: \lambda_s=', num2str(lambda_s), ' \lambda_r=', num2str(lambda_r)]; % part of the title ss_plot_single_ver7( X, Y, SS_stoch, name ); saveas(gcf,['Figure',int2str(figure_no),'_SS_1_subplot.jpg'])
function [ X, Y, SS_cat SS_stoch] = ss_simu_ver7( par, N_s, N_r )
%SIMU Simulates model_ver7 based on input parameters and N_s*N_r inputs
% This model was made by the DTU 2011 iGEM team, consult our wiki for info
% regarding deriving the model: https://2011.igem.org/Team:DTU-Denmark
%% Parameters.
lambda_s= par(1);
lambda_r= par(2);
p_r= par(3);
alpha_s_min = 0.01;
alpha_r_min = 0.0001;
alpha_s_max = 10000;
alpha_r_max = 100;
alpha_s_no = N_s;
alpha_r_no = N_r;
SS_cat=zeros(alpha_r_no,alpha_s_no);
SS_stoch=zeros(alpha_r_no,alpha_s_no);
i=1;
j=1;
vector_s=logspace(log10(alpha_s_min), log10(alpha_s_max), alpha_s_no);
vector_r=logspace(log10(alpha_r_min), log10(alpha_r_max), alpha_r_no);
%% Calculate the steady state for the catalytical and stochiometric model
for X_s = vector_s
for X_r=vector_r
alfa_s= X_s;
alfa_r= X_r;
SS_stoch(i,j)= 1/( (alfa_s*p_r - alfa_r - lambda_r + 2*p_r*lambda_s + sqrt((alfa_s*p_r - alfa_r - lambda_r)^2 + 4*p_r*alfa_s*lambda_r ))/(2*p_r*lambda_s) );
SS_cat(i,j)= 1/( 1 + alfa_s/lambda_s/(1 + alfa_r/lambda_r) );
i=i+1;
end;
i=1;
j=j+1;
end;
[X,Y]= meshgrid(vector_s, vector_r);
end
function ss_plot_single_ver7( X, Y, SS,name)
%% Plots the output from call_model_ver7.
% This model was made by the DTU 2011 iGEM team, consult our wiki for info
% regarding deriving the model: https://2011.igem.org/Team:DTU-Denmark
%% Various functions for plotting
%SS_PLOT Plots steady-states as
%surfc(X, Y, SS); % surface plot with contours
%surf(X, Y, SS); % surface plot without contours
%contourf(X,Y,SS,100)
%meshc(X, Y, SS); % surface plot with contours
pcolor(X, Y, SS); %
shading flat %[interp, flat or faceted]
%contour(X, Y, SS, 3) % contour plot, with 3 contours
%% color map
map = 1/255*[225 204 0; % DTU colormap with the 14 colors
225 153 0;
255 102 0 ;
255 0 0;
153 0 0;
255 0 153;
204 51 153;
153 0 102;
102 0 102;
102 0 153;
51 102 204;
51 204 255;
153 204 51;
102 204 0];
map2 = 1/255*[225 204 0; % DTU colormap with in-between colors
225 153 0; %
240 125 0;
255 102 0 ; %
255 50 0;
255 0 0; %
200 0 0;
153 0 0; %
200 0 75;
255 0 153;%
230 25 153
204 51 153;%
175 25 130;
153 0 102;%
130 0 102;
102 0 102;%
102 0 130;
102 0 153;%
75 50 175;
51 102 204;%
51 150 230;
51 204 255;%
102 204 153;
153 204 51;%
130 204 25;
102 204 0];%
%colormap(map2)
caxis ([0 1]); % specify the min and max value for the colorbar
colorbar % add the color bar
%% Add contours
hold on
contour(X, Y, SS, [0.05 0.05], '-.','LineColor', [0.2 0.2 0.2] ,'LineWidth',3) % contour plot, 5% of max value
contour(X, Y, SS, [0.25 0.25], '-.', 'LineColor',[0.5 0.5 0.5] ,'LineWidth',3) % contour plot, 25% of max value
contour(X, Y, SS, [0.5 0.5], '--k' ,'LineWidth',3) % contour plot, 50% maximum value
contour(X, Y, SS, [0.75 0.75], '-.','LineColor', [0.5 0.5 0.5] ,'LineWidth',3) % contour plot, 75% of max value
contour(X, Y, SS, [0.95 0.95], '-.','LineColor', [0.2 0.2 0.2] ,'LineWidth',3) % contour plot, 95% of max value
%% Set axes, legends and titles.
xlabel('\alpha_s','FontSize',30,'FontWeight','bold')
ylabel('\alpha_r','FontSize',30,'FontWeight','bold')
zlabel('1/\phi','FontSize',30,'FontWeight','bold')
title([name],'FontSize',24,'FontWeight','bold')
set(gca,'xscale', 'log','yscale', 'log')
set(gca,'LineWidth',2,'FontSize',24,'FontWeight','bold')
x_min = min(X(1, :));
x_max = max(X(1, :));
y_min = min(Y(:,1));
y_max = max(Y(:,1));
axis([x_min x_max y_min y_max])
saveas(gcf,['Figure_XX''_SS.jpg'])
end
Simulation
Temporal simulation is performed using the Systems Biology Toolbox 2 http://www.sbtoolbox.org/ environment with numerical integration using ode45. The catalytical model is specified by
********** MODEL NAME Dimensionless form. Catalytical. ********** MODEL NOTES Kinetic model of trap-RNA system. Parameters are estimated from literature. ********** MODEL STATES d/dt(m) = 1 - m - k_s*alpha_m*m*s/(beta_m*beta_s) d/dt(s) = (beta_s/beta_m)*(alpha_s/alpha_m - s - k_r*alpha_m * s * r /(beta_s*beta_r)) d/dt(r) = (beta_r/beta_m)*(alpha_r/alpha_m - r) m(0) = 1 s(0) = 0 r(0) = 0 ********** MODEL PARAMETERS alpha_m = 10 alpha_s = 0 alpha_r = 0 beta_m = 0.0257 beta_s = 0.0257 beta_r = 0.0257 k_s = 0.00082 k_r = 0.0082 ********** MODEL VARIABLES ********** MODEL REACTIONS ********** MODEL FUNCTIONS ********** MODEL EVENTS event = gt(time,1), alpha_s, 40 event = gt(time,3), alpha_r, 200 event = gt(time,6), alpha_r, 0 ********** MODEL MATLAB FUNCTIONS
The partly stoichiometric model is specified by
********** MODEL NAME Dimensionless form. Stoichiometric. ********** MODEL NOTES Kinetic model of trap-RNA system. Parameters are estimated from literature. ********** MODEL STATES d/dt(m) = 1 - m - k_s*alpha_m*m*s/(beta_m*beta_s) d/dt(s) = (beta_s/beta_m)*(alpha_s/alpha_m - s - k_r*alpha_m * s * r /(beta_s*beta_r)) d/dt(r) = (beta_r/beta_m)*(alpha_r/alpha_m - r - k_r*alpha_m * s * r /(beta_s*beta_r)) m(0) = 1 s(0) = 0 r(0) = 0 ********** MODEL PARAMETERS alpha_m = 1 alpha_s = 0 alpha_r = 0 beta_m = 0.0257 beta_s = 0.0257 beta_r = 0.0257 k_s = 0.00082 k_r = 0.0082 ********** MODEL VARIABLES ********** MODEL REACTIONS ********** MODEL FUNCTIONS ********** MODEL EVENTS event = gt(time,1), alpha_s, 40 event = gt(time,3), alpha_r, 200 event = gt(time,6), alpha_r, 0 ********** MODEL MATLAB FUNCTIONS
The script running simulation and generating figures.
%ksim runs a dynamic simulation using Systems Biology Toolbox 2 and plots
clear;
model = SBmodel('model7.txt'); %initialize model
%parameters
alpha_m = 1;
alpha_s = 40; %at induced
alpha_r = 200; %at induced
%%Simulation
time = 6; %running time. Glucose event at t = 6
[out] = SBsimulate(model,time);
%%PLot
t = out.time;
m = out.statevalues(:,1); %m
s = out.statevalues(:,2);
r = out.statevalues(:,3);
%scale to max steady_state at induced levels
s = s .* (alpha_m/alpha_s);
r = r .* (alpha_m/alpha_r);
%ss_r = alpha_r/beta_r;
%r = r ./ss_r;
%Binary on off of s and r
%s = gt(t, 1); %Check model for event time
%r = gt(t, 3);
width = 4; %Line width
subplot(3,1,1)
h1 = plot(t,m); %handle
set(h1, 'color', [51/255, 102/255, 204/255], 'LineWidth',width)
set(gca, 'XTickLabel',[])
subplot(3,1,2)
h2 = plot(t,s);
set(h2, 'color', [237/255, 28/255, 36/255],'LineWidth',width)
set(gca, 'XTickLabel',[])
subplot(3,1,3)
h3 = plot(t,r,'g-');
set(h3, 'color', [102/255, 204/255, 0], 'LineWidth',width)
Sponsors
Thanks to:
