// Goldbeter book chapter
p = defn GoldbeterModel
  var n , v, s, a, g, ks, q, L, c, e, ep, e1, e2;

  [a1]  a -> g;  (s*a*e*(1+a*e)*(1+g)^2)/(L*(1+a*c)^2+((1+a*e)^2*(1+g)^2));
  [a2]  e1 -> a;     v;
  [g1]  g -> e2;     ks*g;

  end;
// Initialise all the metabolite concentrations

   p.n = 2; p.v = 0.2; p.s = 1000; p.a = 9; p.g = 6;
   p.ks = 0.1; p.q = 1; p.L = 7500000; p.c = 0.01;
   p.e = 1.0/(1.0+0.1); p.ep = 0; p.e1 = 0; p.e2 = 0;

   t = 0;

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

   println "Computing....";
  
   hours = 1000;
   hStep = 1;
   dataPoints = trunc(hours/hStep);
   m = matrix(dataPoints,3);
   startTime = timer;

   for i = 1 to dataPoints do
      begin
        t = p.sim.OneStep(t,hStep);
        m[i] = {p.time, p.g, p.a};
      end;
   endTime = timer;
   println "Simulation time (sec) :"   , (endTime - startTime)/1000.0;
 
   graph (m);