Team:DTU-Denmark/Matlab

From 2011.igem.org

Revision as of 02:40, 22 September 2011 by Koplev (Talk | contribs)

Matlab code

Steady-state

% Call file
% Model ver. 7 
% 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

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

Simulation

Temporal simulation is performed using the Systems Biology Toolbox 2 [[1]] 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

File:DTU Model7.txt download

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

File:DTU Model8.txt

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)

File:DTU Ksim.m