p = defn tyson6
      [Reaction1] M -> C2 + YP;  k6*M ;
      [Reaction2] C2 -> CP;  C2*k8notP ;
      [Reaction3] CP -> C2;  CP*k9 ;
      [Reaction4] CP + Y -> pM;  CP*k3*Y ;
   //   [Reaction5] M -> pM;  k5notP*M ;
      [Reaction6] $EmptySet -> Y;  k1aa ;
      [Reaction7] Y -> $EmptySet;  k2*Y ;
      [Reaction8] YP -> $EmptySet;  k7*YP ;
     // [Reaction9] pM -> M;  pM*(k4prime+k4*M^2) ;
      [Reaction9] pM -> M;  pM*(k4prime) ;
     // [Reaction9] pM -> M;  pM*(k4*M^2) ;
      [Reaction10] E + M  -> E2;  20.0*E*M ;
      [Reaction20] E2  -> E + M;  2*E2 ;
      [Reaction30] E2 + M  -> E3; 100.0*E2*M ;
      [Reaction40] E3  -> E2 + M;  20.0*E3 ;
      [Reaction50] pM + E3 -> E3M ;  1*k4*E3*pM;
      [Reaction60] E3M   -> E3 + pM;  100.0*E3M ;
      [Reaction70] E3M   -> E3 + M;  1000.0*E3M ;
  end;
  p.E = 0.1; p.E3M=0;
  p.E2 = 0; p.E3 = 0;
  p.C2 = 0;
  p.CP = 1;
  p.M = 0;
  p.pM = 0.3;
  p.Y = 0;
  p.YP = 0;
  p.EmptySet = 0;
  p.k6 = 1;
  p.k8notP = 1000000;
  p.k9 = 1000;
  p.k3 = 200;
//  p.k5notP = 0;
  p.k1aa = 0.015;
  p.k2 = 0;
  p.k7 = 0.6;
p.k4 = 180;
  p.k4prime = 0.005;
  t = 0;
  println "Computing....";
minutes = 200;
hStep = 0.5;
dataPoints = trunc(minutes/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.M, p.pM, p.CP, p.Y};
   end;
 endTime = timer;
 println "Simulation time (sec) :"   , (endTime - startTime)/1000.0;
graph (m);