Team:DTU-Denmark/Matlab

From 2011.igem.org

(Difference between revisions)
(→Simulation)
(→Simulation)
 
(16 intermediate revisions not shown)
Line 1: Line 1:
-
{{:Team:DTU-Denmark/Templates/Standard_page_begin|Matlab code}}
+
{{:Team:DTU-Denmark/Templates/Standard_page_begin|Matlab}}
 +
==Steady-state==
 +
For the steady-state analysis the following script calls a function that calculates the steady-state solution.
 +
<pre style="height: 200px;">
 +
% Call file
 +
% Model ver. 7
-
== Simulation ==
+
% The model contains the analytical solution in regard to m*max/m for the 3 coupled differential equations, one for ChiXR (r),
-
Temporal simulation is performed using the Systems Biology Toolbox 2 [[http://www.sbtoolbox.org/]] environment with numerical integration using ode45. The used models.
+
% 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.
-
[[File:DTU_Model7.txt]]
+
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!
-
[[File:DTU_Model8.txt]]
+
% 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'])
 +
</pre>
 +
 
 +
The following function calculates the steady-state given a set of parameter inputs.
 +
 
 +
<pre style="height: 200px;">
 +
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
 +
 
 +
 
 +
</pre>
 +
 
 +
The following function plots steady-state values.
 +
 
 +
<pre style="height: 200px;">
 +
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
 +
 
 +
</pre>
 +
 
 +
==Simulation==
 +
Temporal simulation is performed using the Systems Biology Toolbox 2 [http://www.sbtoolbox.org/ www.sbtoolbox.org] environment with numerical integration using ode45. The catalytical model is specified by
 +
 
 +
<pre style="height: 200px;">
 +
********** 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
 +
</pre>
 +
 
 +
The partly stoichiometric model is specified by
 +
 
 +
<pre style="height: 200px;">
 +
********** 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
 +
</pre>
The script running simulation and generating figures.
The script running simulation and generating figures.
 +
 +
<pre style="height: 200px;">
 +
%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)
 +
</pre>
{{:Team:DTU-Denmark/Templates/Standard_page_end}}
{{:Team:DTU-Denmark/Templates/Standard_page_end}}

Latest revision as of 04:08, 22 September 2011

Matlab

Steady-state

For the steady-state analysis the following script calls a function that calculates the steady-state solution. 

% 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'])

The following function calculates the steady-state given a set of parameter inputs. 

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

The following function plots steady-state values. 

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/ 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)

Retrieved from "http://2011.igem.org/Team:DTU-Denmark/Matlab"