User Tools

Site Tools


biological_systems:examples

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

biological_systems:examples [Thursday, 05 February 2009 : 10:28:12] (current)
Line 1: Line 1:
 +==== Simulation Files ====
 +=== Deterministic model ===
  
 +
 +The hybrid Chi model is:
 +
 +<code chi>
 +model MichaelisMenten()=
 +|[ cont S : real = 0.0000005
 +      , E : real = 0.0000002
 +      , C : real = 0.0
 +      , P : real = 0.0
 + , var  k1: real = 1000000.0
 +     , ​ k2: real = 0.000001
 +     , ​ k3: real = 0.1
 + , alg  A,D: real
 +:: eqn     A = k1 * S * E
 +     , dot S =       ​k2 ​ * C - A
 +     , ​    D = (k2 + k3) * C - A
 +     , dot E =  D
 +     , dot C = -D
 +     , dot P = k3 * C
 +]|
 +</​code>​
 +
 +The MATLAB model is:
 +<code matlab>
 +function MM_plot
 +tspan = [0 50]; 
 +yzero = [5e-7; 2e-7; 0; 0];
 +k1 = 1e6; k2 = 1e-4; k3 = 0.1;
 +
 +options = odeset('​AbsTol',​1e-8);  ​
 +[t,y] = ode15s(@mm,​tspan,​yzero,​options);​
 +
 +plot(t,​y(:,​1),'​g','​LineWidth',​2)
 +hold on
 +plot(t,​y(:,​2),'​b--','​LineWidth',​2)
 +plot(t,​y(:,​3),'​m--','​LineWidth',​2)
 +plot(t,​y(:,​4),'​r','​LineWidth',​2)
 +
 +LEGEND('​[S]','​[E]','​[C]','​[P]'​)
 +xlabel('​Time','​FontSize',​14)
 +ylabel('​Concentration','​FontSize',​14)
 +
 +set(gca,'​FontWeight','​Bold','​FontSize',​12)
 +grid on
 +      %--------------Nested function----------
 +      function yprime = mm(t,y)
 +      % MM    Michaelis-Menten reaction Rate Equation ​
 +      yprime = zeros(4,1);
 +      yprime(1) = -k1*y(1)*y(2) + k2*y(3);
 +      yprime(2) = -k1*y(1)*y(2) + (k2+k3)*y(3);​
 +      yprime(3) =  k1*y(1)*y(2) - (k2+k3)*y(3);​
 +      yprime(4) =  k3*y(3);
 +      end
 +end
 +</​code>​
 +
 +The SBToolbox models are:
 +  * Modeled as reaction rate equations
 +<​code>​
 +********** MODEL NAME 
 + MM_rr
 +
 +********** MODEL NOTES
 +
 +********** MODEL STATES ​
 + ​d/​dt(S) = U 
 + ​d/​dt(E) = V 
 + ​d/​dt(C) = W
 + ​d/​dt(P) = X
 +
 + S(0) = 5e-7 
 + E(0) = 2e-7 
 + C(0) = 0 
 + P(0) = 0
 +
 +********** MODEL PARAMETERS ​
 + k1 = 1e6 
 + k2 = 1e-4 
 + k3 = 0.1
 +
 +********** MODEL VARIABLES
 +
 +********** MODEL REACTIONS
 + U = k2 * C - k1 * S * E 
 + V = (k2 + k3) * C - k1 * S * E
 + W = k1 * S * E - (k2 + k3) * C
 + X = k3 * C
 +
 +********** MODEL FUNCTIONS
 +
 +********** MODEL EVENTS
 +
 +********** MODEL MATLAB FUNCTIONS
 +</​code>​
 +  * Modeled as reaction scheme
 +<​code>​
 +********** MODEL NAME
 + ​MM_reaction_scheme
 +
 +********** MODEL NOTES
 +
 +********** MODEL STATE INFORMATION
 + S(0) = 5e-7 
 + E(0) = 2e-7 
 + C(0) = 0 
 + P(0) = 0
 +
 +********** MODEL PARAMETERS ​
 + k1 = 1e6 
 + k2 = 1e-4 
 + k3 = 0.1
 +
 +********** MODEL VARIABLES
 +
 +********** MODEL REACTIONS
 +
 + S + E <=> C            :R1
 + vf = k1 * S * E
 + vr = k2 * C
 + C => P + E    :R2
 + vf = k3 * C 
 + 
 +********** MODEL FUNCTIONS
 +
 +********** MODEL EVENTS
 +
 +********** MODEL MATLAB FUNCTIONS
 +</​code>​
 +
 +
 +The files needed for the other simulators are listed below:
 +  * COPASI file: {{:​biological_systems:​mm_copasi.zip|:​biological_systems:​mm_copasi.zip}}
 +
 +=== Stochastic model ===
 +The MATLAB model is:
 +<code matlab>
 +%SSA_PLOT.M
 +
 +rand('​state',​100)
 +
 +%stoichiometric matrix
 +V = [-1 1 0; -1 1 1; 1 -1 -1; 0 0 1];
 +
 +%%%%%%%%%% Parameters and Initial Conditions %%%%%%%%%
 +nA = 6.023e23; ​            % Avagadro'​s number
 +vol = 1e-15; ​              % volume of system
 +X = zeros(4,1);
 +X(1) = round(5e-7*nA*vol);​ % molecules of substrate
 +X(2) = round(2e-7*nA*vol);​ % molecules of enzyme ​
 +c(1) = 1e6/​(nA*vol);​ c(2) = 1e-4; c(3) = 0.1;
 +t = 0;
 +tfinal = 50;
 +
 +count = 1;
 +tvals(1) = 0;
 +Xvals(:,1) = X;
 +
 +while t < tfinal
 +     a(1) = c(1)*X(1)*X(2);​
 +     a(2) = c(2)*X(3);
 +     a(3) = c(3)*X(3);
 +     asum = sum(a);
 +     j = min(find(rand<​cumsum(a/​asum)));​
 +     tau = log(1/​rand)/​asum;​
 +     X = X + V(:,j);
 +    ​
 +     count = count + 1;
 +     t = t + tau;
 +     ​tvals(count) = t;
 +     ​Xvals(:,​count) = X;
 +end
 +
 +%%%%%%%%%%% Plots
 +
 +L = length(tvals);​
 +tnew = zeros(1,​2*(L-1));​
 +tnew(1:​2:​end-1) = tvals(2:​end);​
 +tnew(2:​2:​end) = tvals(2:​end);​
 +tnew = [tvals(1),​tnew];​
 +
 +Svals = Xvals(1,:);
 +ynew = zeros(1,​2*L-1);​
 +ynew(1:​2:​end) = Svals;
 +ynew(2:​2:​end-1) = Svals(1:​end-1);​
 +plot(tnew,​ynew,'​go-'​)
 +hold on
 +
 +Pvals = Xvals(4,:);
 +ynew = zeros(1,​2*L-1);​
 +ynew(1:​2:​end) = Pvals;
 +ynew(2:​2:​end-1) = Pvals(1:​end-1);​
 +plot(tnew,​ynew,'​r*-'​)
 +
 +xlabel('​Time','​FontSize',​14)
 +ylabel('​Concentration','​FontSize',​14)
 +Legend('​S',​ '​P'​)
 +
 +set(gca,'​FontWeight','​Bold','​FontSize',​12)
 +grid on
 +</​code>​
biological_systems/examples.txt · Last modified: Thursday, 05 February 2009 : 10:28:12 (external edit)