 %  DOQCS : http://doqcs.ncbs.res.in/ 
 %  Accession Name = Repressilator 
 %  Accession Number = 68 
 %  Transcriber = Amrita Sekhar and Sharat J. Vayttaden, NCBS 
 %  Developer = Michael B. Elowitz, Stanislas Lelbler 
 %  Species = E. Coli 
 %  Tissue = E.coli Expression, Transposon, Lambda phage 
 %  Cell Compartment = Cytosol 
 %  Notes =  This is a network model termed as Repressilator in which LacI inhibit transcription of tetR which in turn inhibit transcription of Lambda phage. This model simulates the result as per the Fig 1c of the paper : <a href = http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?cmd=Retrieve&db=pubmed&dopt=Abstract&list_uids=10659856&query_hl=1&itool=pubmed_docsum>Elowitz MB, Leibler S., Nature. 2000 Jan 20;403(6767):335-8.</a> 
 
   % Note: All the concentrations are in micro Mole (uM)

function moose_convert
%Matlab simulation for chemical kinetics

% Saved on Mon Feb 13 16:22:40 2006
% Automatically generated from kkit model in MOOSE by conv_matlab

global buff;
global t;
global y;

% Define Tables here
tspan=[0, 20000];
% Molecule names
molnames = {
'1 /kinetics/Lambda-cl_gene/lambda-cl_mRNA'
'2 /kinetics/Lambda-cl_gene/Lambda-cl_operator'
'3 /kinetics/Lambda-cl_gene/Lambda-cl_operator_1st_bound'
'4 /kinetics/Lambda-cl_gene/Lambda-cl_operator_2nd_bound'
'5 /kinetics/Lambda-cl_gene/Lambda-cl'
'6 /kinetics/Lac_gene/lac_mRNA'
'7 /kinetics/Lac_gene/Lac_operator'
'8 /kinetics/Lac_gene/Lac_operator_1st_bound'
'9 /kinetics/Lac_gene/Lac_operator_2nd_bound'
'10 /kinetics/Lac_gene/Lac'
'11 /kinetics/TetR_gene/TetR_operator'
'12 /kinetics/TetR_gene/TetR_operator_1st_bound'
'13 /kinetics/TetR_gene/TetR_operator_2nd_bound'
'14 /kinetics/TetR_gene/tetR_mRNA'
'15 /kinetics/TetR_gene/TetR'
};

bufnames = {
'1 /kinetics/Lambda-cl_gene/Nucleotides'
'2 /kinetics/Lambda-cl_gene/Amino_Acids'
'3 /kinetics/Lambda-cl_gene/Degraded_lambda-cl_mRNA'
'4 /kinetics/Lambda-cl_gene/Degraded_Lambda-cl'
'5 /kinetics/Lac_gene/Nucleotides'
'6 /kinetics/Lac_gene/Degraded_lac_mRNA'
'7 /kinetics/Lac_gene/Degraded_lac'
'8 /kinetics/Lac_gene/Amino_Acids'
'9 /kinetics/TetR_gene/Nucleotides'
'10 /kinetics/TetR_gene/Amino_Acids'
'11 /kinetics/TetR_gene/Degraded_tetR_mRNA'
'12 /kinetics/TetR_gene/Degraded_TetR'
};

sumtotnames = {
};

% Initial concs of 15 molecules
y0(1) = 0;
y0(2) = 0.0222222;
y0(3) = 0;
y0(4) = 0;
y0(5) = 0;
y0(6) = 0;
y0(7) = 0.0222222;
y0(8) = 0;
y0(9) = 0;
y0(10) = 0;
y0(11) = 0.0222222;
y0(12) = 0;
y0(13) = 0;
y0(14) = 0;
y0(15) = 0;
buff(1) = 0.00055555;
buff(2) = 0.00055555;
buff(3) = 0;
buff(4) = 0;
buff(5) = 0.00055555;
buff(6) = 0;
buff(7) = 0;
buff(8) = 0.00055555;
buff(9) = 0.00055555;
buff(10) = 0.00055555;
buff(11) = 0;
buff(12) = 0;

% Call the ODE solver and display results
[t, y] = ode23s(@f, tspan, y0);

hold on;
do_plots;

function do_plots()
global t;
global y;
plot(t,y(:,15));
plot(t,y(:,5));
plot(t,y(:,10));

% Evaluation function
function dydt = f(t, y)
global buff;
dydt = [
 + 0.9 * y(4) * buff(1) / ( 1.00008 + buff(1) ) + 0.9 * y(3) * buff(1) / ( 1.00008 + buff(1) ) + 900 * y(2) * buff(1) / ( 1.00008 + buff(1) ) - 0.005775 * y(1)
 - 1000.01 * y(2) * y(15) + 224 * y(3)
 - 1000.01 * y(3) * y(15) + 9 * y(4) + 1000.01 * y(2) * y(15) - 224 * y(3)
 + 1000.01 * y(3) * y(15) - 9 * y(4)
 + 300.6 * y(1) * buff(2) / ( 1.00007 + buff(2) ) - 1000.01 * y(7) * y(5) + 224 * y(8) - 1000.01 * y(8) * y(5) + 9 * y(9) - 0.016666 * y(5)
 + 900 * y(7) * buff(5) / ( 1.00008 + buff(5) ) + 0.9 * y(9) * buff(5) / ( 1.00008 + buff(5) ) + 0.9 * y(8) * buff(5) / ( 1.00008 + buff(5) ) - 0.005775 * y(6)
 - 1000.01 * y(7) * y(5) + 224 * y(8)
 + 1000.01 * y(7) * y(5) - 224 * y(8) - 1000.01 * y(8) * y(5) + 9 * y(9)
 + 1000.01 * y(8) * y(5) - 9 * y(9)
 - 0.016666 * y(10) + 300.6 * y(6) * buff(8) / ( 1.00007 + buff(8) ) - 1000.01 * y(11) * y(10) + 224 * y(12) - 1000.01 * y(12) * y(10) + 9 * y(13)
 - 1000.01 * y(11) * y(10) + 224 * y(12)
 - 1000.01 * y(12) * y(10) + 9 * y(13) + 1000.01 * y(11) * y(10) - 224 * y(12)
 + 1000.01 * y(12) * y(10) - 9 * y(13)
 + 0.9 * y(13) * buff(9) / ( 1.00008 + buff(9) ) + 0.9 * y(12) * buff(9) / ( 1.00008 + buff(9) ) + 900 * y(11) * buff(9) / ( 1.00008 + buff(9) ) - 0.005775 * y(14)
 + 300.6 * y(14) * buff(10) / ( 1.00007 + buff(10) ) - 0.01666 * y(15) - 1000.01 * y(3) * y(15) + 9 * y(4) - 1000.01 * y(2) * y(15) + 224 * y(3)
];

% end of model
