// A model of PKA

p = defn DaveModel

 var Raf, Raftot, RafPP, Mek, Mektot, MekPP, Erk, Erktot, ErkPP,
  C, R, RC,  Apuch_delay,
  mRnaReg, Reg, RegP, cAMP, stim,
  AA, H12pete, Mkk3, Mkk3PP, p38, p38PP, p38tot, Mkk3tot, Pp1,
  Ppka, Perk, Pp38, Apuch;

  // part I - 5HT
  [htInput1]  htX1 -> C;                   0.002 ;  // Constant source
  [htInput2]  htX2 -> R;                   0.002 ;  // Constant source
  [htDeg1]    C -> htWaste1;               0.00048*C;
  [htDegR1]   R -> htWaste2;               0.00048*R;
  [htDeg3]    RC -> htWaste3;              0.00048*RC;
  [htRC]      C + R -> RC;                 (12.0/(1.0+(RegP/0.00064)))*R*C;
  [htPGI]     RC -> C + R;                 105.0*RC*cAMP*cAMP;
  [htapuch]   RC -> C;                      0.007*RC*(Apuch - 0.1);
  [htapuchR]  R -> htWaste4;               0.007*R*(Apuch - 0.1);
  //[htapuch]   RC -> C;                      0.007*RC*(Apuch_delay - 0.1);
  //[htapuchR]  R -> htWaste4;               0.007*R*(Apuch_delay - 0.1);

  // Part II - Erk
  [erkInput1]  erkX1 -> Raf;               0.001*RafP;
  [erkDeg1]    Raf -> erkWaste1;           0.00345*Raf*stim;
  [erkRafPP]   ErkPP -> erkPPwaste;        0.12*ErkPP/(ErkPP + 0.08);
  [erkDeg2]    RafPP -> erkWaste2;         0.012*RafPP;
  [erkInput2]  erkX2 -> Mapkk;             0.12*MapkkP/(MapkkP + 0.08) ;
  [erkDeg3]    Mapkk -> erkWaste3;         0.7*RafP*Mapkk/(Mapkk + 0.08);
  [erkMekPP]   erkX3 -> MapkkPP;           0.7*RafP*MapkkP/(MapkkP + 0.08);
  [erkDeg4]    MapkkPP -> erkWaste4;       0.12*MapkkPP/(MapkkPP + 0.08);
  [erkInput3]  erkX4 -> Erk;               0.12*ErkP/(ErkP + 0.08) ;
  [erkDeg5]    Erk -> erkWaste5;           0.44*MapkkPP*Erk/(Erk + 0.08);
  [erkErkPP]   erkX5 -> ErkPP;             0.44*MapkkPP*ErkP/(ErkP + 0.08);
   RafP   = 0.5 - Raf;
   MapkkP = 0.5 - Mapkk - MapkkPP;
   ErkP   = 0.5 - Erk   - ErkPP;

  // Part III - pka-mid
  [pkaInput1]  pkaX1 -> mRnaReg;           0.00002 ;
  [pkaDeg1]    mRnaReg -> pkaWaste1;       0.00003*mRnaReg;
  [pkaDeg2]    mRnaReg -> pkaWaste2;       0.00225*stim*mRnaReg/(mRnaReg+0.01);
  [pkaInput2]  pkaX2 -> Reg;               4.0 *mRnaReg*ErkAct;
  [pkaDeg3]    Reg -> pkaWaste2;           Reg*0.02;
  [pkaMekPP]   Reg -> RegWaste1;           0.16*Reg/(Reg+0.0015);
  [pkaDeg5]    RegP -> pkaWaste3;          RegP*0.02;
  [pkaDeg6]    RegP -> pkaWaste4;          0.16*RegP/(Reg+0.0015);
  [pkaDeg7]    RegSource -> RegP;          1.0*ErkAct*(Reg-RegP)/(Reg-RegP+1.5);
  [pkaInput3]  pkaX3 -> cAMP;              3.6*stim /(stim+14);
  [pkaDeg4]    cAMP -> pkaWaste4;          cAMP-0.06;
  ErkAct = ErkPP + 0.015;

  // Part V - Ap-uch
  [apInput1]  apX1 -> Ppka;                0.01 * C * (1.0 - Ppka);
  [apDeg1]    Ppka -> apWaste1;            1.5 * 0.1*Ppka;
  [apInput2]  apX2 -> Perk;                0.005 * ErkAct * (1.0 - Perk);
  [apDeg2]    Perk -> apWaste2;            0.5 * 0.1 * Perk;
  [apInput4]  apX4 -> Apuch;               0.0289 * (Ppka^2/(Ppka^2+0.15^2))*(Perk^2/(Perk^2+0.02^2));
  [apDeg4]    Apuch -> apWaste4;           Apuch/100.0;
  [apSource]  SourceA -> Apuch;            0.001;
  
  end;

// Initialise all the metabolite concentrations
   // 5ht
p.R = 0.1; p.C = 0.1; p.RC = 3.2;     p.Apuch_delay = 0.0;
p.htX1 = 0; p.htX2 = 0;
p.htWaste1 = 0; p.htWaste2 = 0; p.htWaste3 = 0;  p.htWaste4 = 0;
//erk
p.Raf = 0.5;    p.RafPP = 0.0; p.RafP = 0.0; p.Mapkk = 0.5;   p.MapkkPP = 0.0;
p.Erk = 0.05;   p.ErkPP = 0.45;
p.erkX1 = 0; p.erkX2 = 0;   p.erkX3 = 0; p.erkX4 = 0; p.erkX5 = 0;
p.erkWaste1 = 0; p.erkWaste2 = 0; p.erkWaste3 = 0; p.erkWaste4 = 0; p.erkWaste5 = 0;
//pka
p.mRnaReg = 0.67;    p.RegP = 0.0; p.Reg = 0.0;   p.cAMP = 0.0;
p.pkaX1 = 0; p.pkaX2 = 0;   p.pkaX3 = 0;
p.pkaWaste1 = 0; p.pkaWaste2 = 0; p.pkaWaste3 = 0; p.pkaWaste4 = 0;
//apuch
p.Ppka= 0.0;    p.Perk = 0.0;  p.Apuch = 0.1;
p.apX1 = 0; p.apX2 = 0;   p.apX4 = 0;
p.apWaste1 = 0; p.apWaste2 = 0;  p.apWaste4 = 0;

t = 0;
p.stim = 0;

println "Computing....";
minutes = 800;
hStep = 0.1;

dataPoints = trunc(minutes/hStep);

m = matrix(dataPoints,5);
startTime = timer;
 for i = 1 to dataPoints do
    begin
     t = p.sim.OneStep(t,hstep);

       m[i] = {p.time, p.Perk, p.Ppka, p.Apuch, p.C}

        if((t > 50+5.0)  and  (t < 50+10.0)) then p.stim = 10.0;
	if((t > 50+10.0) and  (t < 50+25.0)) then p.stim = 0.0;
        if((t > 50+25.0) and  (t < 50+30.0)) then p.stim = 10.0;
	if((t > 50+30.0) and  (t < 50+45.0)) then p.stim = 0.0;
	if((t > 50+45.0) and  (t < 50+50.0)) then p.stim = 10.0;
	if((t > 50+50.0) and  (t < 50+65.0)) then p.stim = 0.0;
	if((t > 50+65.0) and  (t < 50+70.0)) then p.stim = 10.0;
	if((t > 50+70.0) and  (t < 50+85.0)) then p.stim = 0.0;
	if((t > 50+85.0) and  (t < 50+90.0)) then p.stim = 10.0;
	if (t > 50+90.0)  then  p.stim = 0.0;

    end;
 endTime = timer;
 println "Simulation time (sec) :"   , (endTime - startTime)/1000.0;
 
graph (m);