// Smolen et al book chapter

p = defn CircadianModel

 var mRNA,  ImRNA,   Ptot,
 P0, P1, P2, P3, P4, P5, P6, P7, P8, P9, P10, P11, P12, WP12;

  [ImRNA]    ImRNA -> mRNA;          7.0*kr/(P12+kr);
  [WmRNA]    mRNA -> WmRNA;          0.2*mRNA;
  [IP0]      IP0 -> P0;              0.2*mRNA;
  [IP1]      P0 -> P1;               kph*P0 ;
  [WP1]      P1 -> P2;               kph*P1;
  [IP3]      P2 -> P3;               kph*P2 ;
  [WP3]      P3 -> P4;               kph*P3;
  [IP5]      P4 -> P5;               kph*P4 ;
  [WP5]      P5 -> P6;               kph*P5;
  [IP6]      P6 -> P7;               kph*P6 ;
  [WP6]      P7 -> P8;               kph*P7;
  [IP9]      P8 -> P9;               kph*P8 ;
  [WP9]      P9 -> P10;              kph*P9;
  [IP11]     P10 -> P11;             kph*P10 ;
  [WP11]     P11 -> P12;             kph*P11;
  [WP12]     P12 -> WP12;            3.0*P12/(0.0001+P12);

  Ptot = P0 + P1 + P2 + P3 + P4 + P5 + P6 + P7 + P8 + P9 + P10 + P11 + P12;
  end;

// Initialise all the metabolite concentrations

p.ImRNA = 0;   p.kr = 0.0005;   p.kph = 2.0;  p.mRNA = 0.0;     p.Ptot = 0;
p.P0 = 0; p.P1 = 0; p.P2 = 0; p.P3 = 0; p.P4 = 0; p.P5 = 0; p.P6 = 0; p.P7 = 0;
p.P8 = 0; p.P9 = 0; p.P10 = 0; p.P11 = 0; p.P12 = 0; p.WP12 = 0;

t = 0;

//p.xml("PaulModelSBML.xml");

println "Computing....";
  
hours = 60;
hStep = 0.5;
dataPoints = trunc(hours/hStep);
m = matrix(dataPoints,5);
startTime = timer;

 for i = 1 to dataPoints do
    begin
     t = p.sim.OneStep(t,hStep);
  if (mod(i,trunc(1.0/hstep))==0) then print"+";
       m[i] = {p.time, p.mRNA, p.P0, p.P12, p.Ptot};
   end;
 endTime = timer;
 println "Simulation time (sec) :"   , (endTime - startTime)/1000.0;
 
graph (m);