(************** Content-type: application/mathematica **************
                    CreatedBy='SBW JDesigner Export                 

 This notebook can be used with any Mathematica-compatible          
 application, such as Mathematica, MathReader or Publicon. The data 
 for the notebook starts with the line containing stars above.      

synopsis:
      xdot = Goldbeter2008_Somite_Segmentation_Clock_Notch_Wnt_FGF (time, x)
      x0 = Goldbeter2008_Somite_Segmentation_Clock_Notch_Wnt_FGF
 
   Goldbeter2008_Somite_Segmentation_Clock_Notch_Wnt_FGF can be used with odeN functions as follows:
 
   x0 = untitled;
   [t,x] = ode23s(@Goldbeter2008_Somite_Segmentation_Clock_Notch_Wnt_FGF, [0 100], Goldbeter2008_Somite_Segmentation_Clock_Notch_Wnt_FGF);
   plot (x);
 
   where 100 is the end time
 
   When Goldbeter2008_Somite_Segmentation_Clock_Notch_Wnt_FGF is used without any arguments it returns a vector of
   the initial concentrations of the 22 floating species.
   Otherwise Goldbeter2008_Somite_Segmentation_Clock_Notch_Wnt_FGF should be called with two arguments:
   time and x.  time is the current time. x is the vector of the
   concentrations of the 22 floating species.
   When these parameters are supplied Goldbeter2008_Somite_Segmentation_Clock_Notch_Wnt_FGF returns a vector of 
   the rates of change of the concentrations of the 22 floating species.
 
   the following table shows the mapping between the vector
   index of a floating species and the species name.

   index    Name
   x[1][t] := n           
   x[2][t] := na          
   x[3][t] := nan         
   x[4][t] := mf          
   x[5][t] := f           
   x[6][t] := bp          
   x[7][t] := bn          
   x[8][t] := a           
   x[9][t] := k           
   x[10][t] := b           
   x[11][t] := max         
   x[12][t] := rasa        
   x[13][t] := erka        
   x[14][t] := xa          
   x[15][t] := mdusp       
   x[16][t] := dusp        
   x[17][t] := rast        
   x[18][t] := erkt        
   x[19][t] := xt          
   x[20][t] := d           
   x[21][t] := kt          
   x[22][t] := fgf         
*)


initdata = {
   x'[1][0.] == 0.5,                   (* n *)
   x'[2][0.] == 0.2,                   (* na *)
   x'[3][0.] == 0,                     (* nan *)
   x'[4][0.] == 0.1,                   (* mf *)
   x'[5][0.] == 0.001,                 (* f *)
   x'[6][0.] == 0.1,                   (* bp *)
   x'[7][0.] == 0.001,                 (* bn *)
   x'[8][0.] == 0.1,                   (* a *)
   x'[9][0.] == 3,                     (* k *)
   x'[10][0.] == 0.1,                   (* b *)
   x'[11][0.] == 0.1,                  (* max *)
   x'[12][0.] == 0.5,                  (* rasa *)
   x'[13][0.] == 0.2,                  (* erka *)
   x'[14][0.] == 0.1,                  (* xa *)
   x'[15][0.] == 0.1,                  (* mdusp *)
   x'[16][0.] == 0.1,                  (* dusp *)
   x'[17][0.] == 2,                    (* rast *)
   x'[18][0.] == 2,                    (* erkt *)
   x'[19][0.] == 2,                    (* xt *)
   x'[20][0.] == 2,                    (* d *)
   x'[21][0.] == 3,                    (* kt *)
   x'[22][0.] == 1,                    (* fgf *)
}

   vsFK = vsF * (KIG1 / (KIG1 + K));

   AK = Kt - K;

   Rasi = Rast - Rasa;

   ERKi = ERKt - ERKa;

   Xi = Xt - Xa;

Array[p, 77];
   p[1] = 0.5;                      (* __src__*)
   p[2] = 0.5;                      (* __waste__*)
   p[3] = 1.8;                      (* a1*)
   p[4] = 0.5;                      (* ak*)
   p[5] = 0.1;                      (* d1*)
   p[6] = 0.3;                      (* epsilon*)
   p[7] = 0.5;                      (* erki*)
   p[8] = 0.3;                      (* eta*)
   p[9] = 2;                        (* j*)
   p[10] = 0.28;                    (* k1*)
   p[11] = 0.03;                    (* k2*)
   p[12] = 0.05;                    (* ka*)
   p[13] = 0.7;                     (* kab*)
   p[14] = 0.05;                    (* kaerk*)
   p[15] = 0.5;                     (* kafgf*)
   p[16] = 0.5;                     (* kamdusp*)
   p[17] = 0.103;                   (* karas*)
   p[18] = 0.05;                    (* kax*)
   p[19] = 0.05;                    (* kaxa*)
   p[20] = 3.45;                    (* kc*)
   p[21] = 1.35;                    (* kcdusp*)
   p[22] = 0;                       (* kd1*)
   p[23] = 7.062;                   (* kd2*)
   p[24] = 0.63;                    (* kdax*)
   p[25] = 0.5;                     (* kddusp*)
   p[26] = 0.05;                    (* kderk*)
   p[27] = 0.37;                    (* kdf*)
   p[28] = 0.5;                     (* kdmdusp*)
   p[29] = 0.768;                   (* kdmf*)
   p[30] = 1.4;                     (* kdn*)
   p[31] = 0.001;                   (* kdna*)
   p[32] = 0.001;                   (* kdnan*)
   p[33] = 0.1;                     (* kdras*)
   p[34] = 0.05;                    (* kdx*)
   p[35] = 0.5;                     (* kid*)
   p[36] = 0.5;                     (* kif*)
   p[37] = 2.5;                     (* kig1*)
   p[38] = 0.48;                    (* kmd*)
   p[39] = 0.02;                    (* ksax*)
   p[40] = 0.5;                     (* ksdusp*)
   p[41] = 0.3;                     (* ksf*)
   p[42] = 0.1;                     (* kt1*)
   p[43] = 0.1;                     (* kt2*)
   p[44] = 0.7;                     (* kt3*)
   p[45] = 1.5;                     (* kt4*)
   p[46] = 2;                       (* m*)
   p[47] = 2;                       (* n*)
   p[48] = 2;                       (* p*)
   p[49] = 2;                       (* q*)
   p[50] = 2;                       (* r*)
   p[51] = 0.5;                     (* rasi*)
   p[52] = 1.5;                     (* theta*)
   p[53] = 0.06;                    (* v0*)
   p[54] = 0.6;                     (* vdax*)
   p[55] = 2;                       (* vddusp*)
   p[56] = 0.39;                    (* vdf*)
   p[57] = 2.82;                    (* vdn*)
   p[58] = 0.01;                    (* vdna*)
   p[59] = 0.1;                     (* vdnan*)
   p[60] = 3.3;                     (* vmaerk*)
   p[61] = 4.968;                   (* vmaras*)
   p[62] = 1.6;                     (* vmax*)
   p[63] = 1.64;                    (* vmb*)
   p[64] = 0.8;                     (* vmd*)
   p[65] = 0.5;                     (* vmdmdusp*)
   p[66] = 0.41;                    (* vmdras*)
   p[67] = 0.5;                     (* vmdx*)
   p[68] = 1.92;                    (* vmf*)
   p[69] = 5.08;                    (* vmk*)
   p[70] = 1;                       (* vmp*)
   p[71] = 0.9;                     (* vmsmdusp*)
   p[72] = 0.5;                     (* vmxa*)
   p[73] = 0.087;                   (* vsb*)
   p[74] = 3;                       (* vsf*)
   p[75] = 0;                       (* vsfk*)
   p[76] = 0.23;                    (* vsn*)
   p[77] = 0.5;                     (* xi*)

Array[rates, 33];
   rates[1] := cytosol*p[5]*p[75];
   rates[2] := p[5]*cytosol*p[56]*x[0][t]/(p[29]+x[0][t]);
   rates[3] := p[5]*cytosol*p[19]*x[0][t]*Power(p[35],p[8])/(Power(p[35],p[8])+Power(x[4][t],p[8]));
   rates[4] := p[5]*cytosol*p[57]*x[1][t]/(p[30]+x[1][t]);
   rates[5] := p[5]*cytosol*(p[41]*x[1][t]-p[42]*x[2][t]);
   rates[6] := p[5]*cytosol*p[58]*x[2][t]/(p[31]+x[2][t]);
   rates[7] := p[5]*cytosol*p[74]*Power(x[2][t],p[47])/(Power(p[11],p[47])+Power(x[2][t],p[47]));
   rates[8] := p[5]*cytosol*p[67]*x[3][t]/(p[28]+x[3][t]);
   rates[9] := p[5]*cytosol*p[40]*x[3][t];
   rates[10] := p[5]*cytosol*p[55]*x[4][t]/(p[26]+x[4][t]);
   rates[11] := p[51]*cytosol*(p[4]*p[3]-p[2]*x[7][t]*x[8][t]);
   rates[12] := p[51]*cytosol*p[72];
   rates[13] := p[51]*cytosol*p[21]*x[9][t];
   rates[14] := p[51]*cytosol*p[68]*p[34]/(p[34]+x[19][t])*x[9][t]/(p[9]+x[9][t])*p[3]/x[20][t];
   rates[15] := p[51]*cytosol*p[69]*x[5][t]/(p[10]+x[5][t]);
   rates[16] := p[51]*cytosol*(p[44]*x[6][t]-p[43]*x[9][t]);
   rates[17] := p[51]*cytosol*p[22]*x[5][t];
   rates[18] := p[51]*cytosol*p[52];
   rates[19] := p[51]*cytosol*(p[62]*Power(x[6][t],x[0][t])/(Power(p[12],x[0][t])+Power(x[6][t],x[0][t])));
   rates[20] := p[51]*cytosol*(p[71]*Power(x[13][t],p[45])/(Power(p[18],p[45])+Power(x[13][t],p[45])));
   rates[21] := p[51]*cytosol*p[63]*x[10][t]/(p[37]+x[10][t]);
   rates[22] := p[51]*cytosol*p[38]*x[10][t];
   rates[23] := p[51]*cytosol*p[53]*x[7][t]/(p[23]+x[7][t]);
   rates[24] := p[7]*cytosol*p[60]*Power(x[21][t],p[49])/(Power(p[14],p[49])+Power(x[21][t],p[49]))*p[50]/(p[16]+p[50]);
   rates[25] := p[7]*cytosol*p[65]*x[11][t]/(p[32]+x[11][t]);
   rates[26] := p[7]*cytosol*p[59]*x[11][t]/x[16][t]*p[6]/(p[13]+p[6]);
   rates[27] := p[7]*cytosol*p[20]*x[15][t]*x[12][t]/(p[25]+x[12][t]);
   rates[28] := p[7]*cytosol*p[61]*x[12][t]/x[17][t]*p[76]/(p[17]+p[76]);
   rates[29] := p[7]*cytosol*p[66]*x[13][t]/(p[33]+x[13][t]);
   rates[30] := p[7]*cytosol*p[70]*Power(x[13][t],p[48])/(Power(p[15],p[48])+Power(x[13][t],p[48]));
   rates[31] := p[7]*cytosol*p[64]*x[14][t]/(p[27]+x[14][t]);
   rates[32] := p[7]*cytosol*p[39]*x[14][t];
   rates[33] := p[7]*cytosol*p[54]*x[15][t]/(p[24]+x[15][t]);

   x'[1][t] == {  + rates[1] - rates[2] - rates[3]}/.rates;
   x'[2][t] == {  + rates[3] - rates[4] - rates[5]}/.rates;
   x'[3][t] == {  + rates[5] - rates[6]}/.rates;
   x'[4][t] == {  + rates[7] - rates[8]}/.rates;
   x'[5][t] == {  + rates[9] - rates[10]}/.rates;
   x'[6][t] == {  + rates[14] - rates[15] - rates[17]}/.rates;
   x'[7][t] == {  - rates[16]}/.rates;
   x'[8][t] == {  + rates[11] + rates[22] - rates[23]}/.rates;
   x'[9][t] == {  + rates[11]}/.rates;
   x'[10][t] == {  + rates[12] - rates[13] - rates[14] + rates[15] + rates[16]}/.rates;
   x'[11][t] == {  + rates[18] + rates[19] + rates[20] - rates[21]}/.rates;
   x'[12][t] == {  + rates[24] - rates[25]}/.rates;
   x'[13][t] == {  + rates[26] - rates[27]}/.rates;
   x'[14][t] == {  + rates[28] - rates[29]}/.rates;
   x'[15][t] == {  + rates[30] - rates[31]}/.rates;
   x'[16][t] == {  + rates[32] - rates[33]}/.rates;
   x'[17][t] == { }/.rates;
   x'[18][t] == { }/.rates;
   x'[19][t] == { }/.rates;
   x'[20][t] == { }/.rates;
   x'[21][t] == { }/.rates;
   x'[22][t] == { }/.rates;
 

sol = NDSolve[Evaluate[Join[Table[x'[i][t], {i,22}], Table[x'[i][0], {i, 22}]]/.p, Table[x[i], {i, 22}], {t, 0, tmax}]];

Plot[Evaluate[Table[x[i][t], {i, 22}]/.sol], {t, 0, tmax}];