(* Generated by JWS Online *)

(* This is an experimental feature of JWS Online. Please report any mistakes.*)

(* Note that the following notable SBML entities or features are not supported in notebook outputyet: *)
    (* Events *)
    (* Constraints *)
    (* Units and UnitDefinitions *)
    (* AlgebraicRules *)
    (* conversionFactors *)

variables = {
DHAP[t], DPG[t], F6P[t], FDP[t], G1P[t], G6P[t], GAP[t], NAD[t], NADH[t], P2G[t], P3G[t], PEP[t], PYR[t]
};

initialValues = {
DHAP[0] == 7*^-05, DPG[0] == 6.5*^-05, F6P[0] == 0.0002, FDP[0] == 7*^-05, G1P[0] == 5*^-05, G6P[0] == 0.00075, GAP[0] == 3*^-05, NAD[0] == 0.0005, NADH[0] == 5*^-07, P2G[0] == 5*^-06, P3G[0] == 5*^-05, PEP[0] == 1*^-05, PYR[0] == 0.000144
};

rates = {
v\[LetterSpace]1, v\[LetterSpace]10, v\[LetterSpace]11, v\[LetterSpace]12, v\[LetterSpace]2, v\[LetterSpace]3, v\[LetterSpace]4, v\[LetterSpace]5, v\[LetterSpace]6, v\[LetterSpace]7, v\[LetterSpace]8, v\[LetterSpace]9
};

rateEquations = {
v\[LetterSpace]1 -> fracA*((GLY*Ph*Vfgly)/(KgpAigly*KgpApi*(1 + GLY/KgpAglyb + GLY/KgpAglyf + (GLY*Ph)/(KgpAglyf*KgpAipi) + Ph/KgpApi + G1P[t]/KgpAg1p + (GLY*G1P[t])/(KgpAglyb*KgpAig1p))) - (3.2258064516129035*GLY*Vfgly*G1P[t])/(KgpAigly*KgpApi*(1 + GLY/KgpAglyb + GLY/KgpAglyf + (GLY*Ph)/(KgpAglyf*KgpAipi) + Ph/KgpApi + G1P[t]/KgpAg1p + (GLY*G1P[t])/(KgpAglyb*KgpAig1p)))) + fracB*((50.*amp^nH*GLY*Ph*Vfgly)/(Kgpamp^nH*(1 + (50.*amp^nH)/Kgpamp^nH)*KgpBiglyf*KgpBpi*(1 + GLY/KgpBiglyb + GLY/KgpBipi + Ph/KgpBiglyf + (GLY*Ph)/(KgpBiglyf*KgpBpi) + G1P[t]/KgpBig1p + (GLY*G1P[t])/(KgpBg1p*KgpBiglyb))) - (161.29032258064518*amp^nH*GLY*Vfgly*G1P[t])/(Kgpamp^nH*(1 + (50.*amp^nH)/Kgpamp^nH)*KgpBiglyf*KgpBpi*(1 + GLY/KgpBiglyb + GLY/KgpBipi + Ph/KgpBiglyf + (GLY*Ph)/(KgpBiglyf*KgpBpi) + G1P[t]/KgpBig1p + (GLY*G1P[t])/(KgpBg1p*KgpBiglyb)))), v\[LetterSpace]10 -> ((Vfen*P2G[t])/Ken2pg - (2.0408163265306123*Vfen*PEP[t])/Ken2pg)/(1 + P2G[t]/Ken2pg + PEP[t]/Kenpep), v\[LetterSpace]11 -> ((adp*Vfpk*PEP[t])/(Kpkadp*Kpkpep) - (atp*Vfpk*PYR[t])/(10304*Kpkadp*Kpkpep))/(1 + adp/Kpkadp + atp/Kpkatp + PEP[t]/Kpkpep + (adp*PEP[t])/(Kpkadp*Kpkpep) + PYR[t]/Kpkpyr + (atp*PYR[t])/(Kpkatp*Kpkpyr)), v\[LetterSpace]12 -> (-(LAC*Vfldh*NAD[t])/(16198*Kldhnadh*Kldhpyr) + (Vfldh*NADH[t]*PYR[t])/(Kldhnadh*Kldhpyr))/(1 + LAC/Kldhlac + NAD[t]/Kldhnad + (LAC*NAD[t])/(Kldhlac*Kldhnad) + NADH[t]/Kldhnadh + PYR[t]/Kldhpyr + (NADH[t]*PYR[t])/(Kldhnadh*Kldhpyr)), v\[LetterSpace]2 -> ((Vfpglm*G1P[t])/Kpglmg1p - (0.06016847172081829*Vfpglm*G6P[t])/Kpglmg1p)/(1 + G1P[t]/Kpglmg1p + G6P[t]/Kpglmg6p), v\[LetterSpace]3 -> (-((Vbpgi*F6P[t])/Kpgif6p) + (0.45*Vbpgi*G6P[t])/Kpgif6p)/(1 + F6P[t]/Kpgif6p + G6P[t]/Kpgig6p), v\[LetterSpace]4 -> (atp*Vfpfk*F6P[t]*(1 + ((1 + (amp*en)/Kpfkamp)^4*Kpfkatp*Kpfkf6p*(1 + atp/Kpfkiatp)^4*Lo*((1 + atp/KpfkatpT)*(1 + F6P[t]/Kpfkf6pT) + FDP[t]/KpfkfdpT + (adp*(1 + FDP[t]/KpfkfdpT))/KpfkadpT)^3)/((1 + amp/Kpfkamp)^4*KpfkatpT*Kpfkf6pT*(1 + (atp*dn)/Kpfkiatp)^4*((1 + atp/Kpfkatp)*(1 + F6P[t]/Kpfkf6p) + FDP[t]/Kpfkfdp + (adp*(1 + FDP[t]/Kpfkfdp))/Kpfkadp)^3)))/(Kpfkatp*Kpfkf6p*((1 + atp/Kpfkatp)*(1 + F6P[t]/Kpfkf6p) + FDP[t]/Kpfkfdp + (adp*(1 + FDP[t]/Kpfkfdp))/Kpfkadp)*(1 + ((1 + (amp*en)/Kpfkamp)^4*(1 + atp/Kpfkiatp)^4*Lo*((1 + atp/KpfkatpT)*(1 + F6P[t]/Kpfkf6pT) + FDP[t]/KpfkfdpT + (adp*(1 + FDP[t]/KpfkfdpT))/KpfkadpT)^4)/((1 + amp/Kpfkamp)^4*(1 + (atp*dn)/Kpfkiatp)^4*((1 + atp/Kpfkatp)*(1 + F6P[t]/Kpfkf6p) + FDP[t]/Kpfkfdp + (adp*(1 + FDP[t]/Kpfkfdp))/Kpfkadp)^4))) - (0.00411742907728414*adp*Vfpfk*FDP[t]*(1 + ((1 + (amp*en)/Kpfkamp)^4*Kpfkatp*Kpfkf6p*(1 + atp/Kpfkiatp)^4*Lo*((1 + atp/KpfkatpT)*(1 + F6P[t]/Kpfkf6pT) + FDP[t]/KpfkfdpT + (adp*(1 + FDP[t]/KpfkfdpT))/KpfkadpT)^3)/((1 + amp/Kpfkamp)^4*KpfkatpT*Kpfkf6pT*(1 + (atp*dn)/Kpfkiatp)^4*((1 + atp/Kpfkatp)*(1 + F6P[t]/Kpfkf6p) + FDP[t]/Kpfkfdp + (adp*(1 + FDP[t]/Kpfkfdp))/Kpfkadp)^3)))/(Kpfkatp*Kpfkf6p*((1 + atp/Kpfkatp)*(1 + F6P[t]/Kpfkf6p) + FDP[t]/Kpfkfdp + (adp*(1 + FDP[t]/Kpfkfdp))/Kpfkadp)*(1 + ((1 + (amp*en)/Kpfkamp)^4*(1 + atp/Kpfkiatp)^4*Lo*((1 + atp/KpfkatpT)*(1 + F6P[t]/Kpfkf6pT) + FDP[t]/KpfkfdpT + (adp*(1 + FDP[t]/KpfkfdpT))/KpfkadpT)^4)/((1 + amp/Kpfkamp)^4*(1 + (atp*dn)/Kpfkiatp)^4*((1 + atp/Kpfkatp)*(1 + F6P[t]/Kpfkf6p) + FDP[t]/Kpfkfdp + (adp*(1 + FDP[t]/Kpfkfdp))/Kpfkadp)^4))), v\[LetterSpace]5 -> ((Vfald*FDP[t])/Kaldfdp - (10526.315789473683*Vfald*DHAP[t]*GAP[t])/Kaldfdp)/(1 + DHAP[t]/Kalddhap + FDP[t]/Kaldfdp + GAP[t]/Kaldgap), v\[LetterSpace]6 -> ((-0.052083333333333336*Vftpi*DHAP[t])/Ktpigap + (Vftpi*GAP[t])/Ktpigap)/(1 + DHAP[t]/Ktpidhap + GAP[t]/Ktpigap), v\[LetterSpace]7 -> ((Ph*Vfgad*GAP[t]*NAD[t])/(Kgapdhgap*Kgapdhnad*Kgapdhpi) - (11.2359550561798*Vfgad*DPG[t]*NADH[t])/(Kgapdhgap*Kgapdhnad*Kgapdhpi))/(1 + Ph/Kgapdhpi + DPG[t]/Kgapdh13dpg + GAP[t]/Kgapdhgap + NAD[t]/Kgapdhnad + (GAP[t]*NAD[t])/(Kgapdhgap*Kgapdhnad) + (Ph*GAP[t]*NAD[t])/(Kgapdhgap*Kgapdhnad*Kgapdhpi) + NADH[t]/Kgapdhnadh + (DPG[t]*NADH[t])/(Kgapdh13dpg*Kgapdhnadh)), v\[LetterSpace]8 -> ((57109*adp*Vbpgk*DPG[t])/(Kpgk3pg*Kpgkatp) - (atp*Vbpgk*P3G[t])/(Kpgk3pg*Kpgkatp))/(1 + adp/Kpgkadp + atp/Kpgkatp + DPG[t]/Kpgk13dpg + (adp*DPG[t])/(Kpgk13dpg*Kpgkadp) + P3G[t]/Kpgk3pg + (atp*P3G[t])/(Kpgk3pg*Kpgkatp)), v\[LetterSpace]9 -> ((-5.555555555555555*Vfpgm*P2G[t])/Kpgm3pg + (Vfpgm*P3G[t])/Kpgm3pg)/(1 + P2G[t]/Kpgm2pg + P3G[t]/Kpgm3pg)
};

parameters = {
Kadkadp -> 0.00035, Kadkamp -> 0.00032, Kadkatp -> 0.00027, Kalddhap -> 0.002, Kaldfdp -> 0.0005, Kaldgap -> 0.001, KckPCr -> 0.00111, Kckcr -> 0.0038, KckiPCr -> 0.0039, Kckiadp -> 0.000135, Kckiatp -> 0.0035, Ken2pg -> 0.0001, Kenpep -> 0.00037, KeqCK -> 233.0, Kgapdh13dpg -> 8*^-07, Kgapdhgap -> 2.5*^-06, Kgapdhnad -> 9*^-05, Kgapdhnadh -> 3.3*^-06, Kgapdhpi -> 0.00029, KgpAg1p -> 0.0027, KgpAglyb -> 0.00015, KgpAglyf -> 0.0017, KgpAig1p -> 0.0101, KgpAigly -> 0.002, KgpAipi -> 0.0047, KgpApi -> 0.004, KgpBg1p -> 0.0015, KgpBig1p -> 0.0074, KgpBiglyb -> 0.0044, KgpBiglyf -> 0.015, KgpBipi -> 0.0046, KgpBpi -> 0.0002, Kgpamp -> 9.7*^-05, Kldhlac -> 0.01717, Kldhnad -> 0.000849, Kldhnadh -> 2.167*^-06, Kldhpyr -> 0.000335, Kpfkadp -> 0.00271, KpfkadpT -> 0.00271, Kpfkamp -> 6*^-05, Kpfkatp -> 8*^-05, KpfkatpT -> 0.00025, Kpfkf6p -> 0.00018, Kpfkf6pT -> 0.02, Kpfkfdp -> 0.0042, KpfkfdpT -> 0.0042, Kpfkiatp -> 0.00087, Kpgif6p -> 0.000119, Kpgig6p -> 0.00048, Kpgk13dpg -> 1.9*^-06, Kpgk3pg -> 0.0012, Kpgkadp -> 8.3*^-05, Kpgkatp -> 0.00035, Kpglmg1p -> 6.3*^-05, Kpglmg6p -> 3*^-05, Kpgm2pg -> 1.4*^-05, Kpgm3pg -> 0.0002, Kpkadp -> 0.0003, Kpkatp -> 0.00113, Kpkpep -> 8*^-05, Kpkpyr -> 0.00705, Ktpidhap -> 0.00061, Ktpigap -> 0.00032, Lo -> 13.0, Vbpgi -> 0.88, Vbpgk -> 1.12, Vfadk -> 0.88, Vfald -> 0.104, Vfen -> 0.5, Vfgad -> 1.65, Vfgly -> 0.15, Vfldh -> 1.92, Vfpfk -> 0.056, Vfpglm -> 0.48, Vfpgm -> 1.12, Vfpk -> 3.0, Vftpi -> 12.0, VrevCK -> 0.5, amp -> 9.52098310973979*^-06, dn -> 0.01, en -> 0.01, fracA -> 0.4, fracB -> 0.6, k -> 0.075, kout -> 0.2, nH -> 1.75, GLY -> 0.112, LAC -> 0.00205129800531363, Ph -> 0.0307416000593277, adp -> 0.000405079560143532, atp -> 0.00779841945674673, default\[LetterSpace]compartment -> 1.0
};

assignments = {

};

events = {

};

speciesAnnotations = {

};

reactionAnnotations = {

};

units = {
{"time" -> "min", "metabolite" -> "mole/L", "extent" -> "M"}
};

(* Time evolution *)
odes = {
DHAP'[t] == 1.0*v\[LetterSpace]6 +1.0*v\[LetterSpace]5 , DPG'[t] == 1.0*v\[LetterSpace]7 -1.0*v\[LetterSpace]8, F6P'[t] == 1.0*v\[LetterSpace]3 -1.0*v\[LetterSpace]4, FDP'[t] == 1.0*v\[LetterSpace]4 -1.0*v\[LetterSpace]5, G1P'[t] == 1.0*v\[LetterSpace]1 -1.0*v\[LetterSpace]2, G6P'[t] == 1.0*v\[LetterSpace]2 -1.0*v\[LetterSpace]3, GAP'[t] == 1.0*v\[LetterSpace]5 -1.0*v\[LetterSpace]6 -1.0*v\[LetterSpace]7, NAD'[t] == 1.0*v\[LetterSpace]12 -1.0*v\[LetterSpace]7, NADH'[t] == 1.0*v\[LetterSpace]7 -1.0*v\[LetterSpace]12, P2G'[t] == 1.0*v\[LetterSpace]9 -1.0*v\[LetterSpace]10, P3G'[t] == 1.0*v\[LetterSpace]8 -1.0*v\[LetterSpace]9, PEP'[t] == 1.0*v\[LetterSpace]10 -1.0*v\[LetterSpace]11, PYR'[t] == 1.0*v\[LetterSpace]11 -1.0*v\[LetterSpace]12
};
timeCourse = NDSolve[Join[odes, initialValues]//.rateEquations//.assignments//.parameters, variables, {t, 0, 100}];

(* Steady-state solution initialized with result of time evolution *)
findRootEquations = odes /.D[_[t],t]->0;
findRootVariables = Partition[Flatten[{#, #/.timeCourse/.t->100} &/@variables],2];
steadyStateVariables = FindRoot[findRootEquations//.rateEquations//.assignments//.parameters, findRootVariables, MaxIterations->100]
fluxes = #//.assignments//.parameters/.steadyStateVariables&/@rateEquations

(* Plot the time evolution of the variables *)
plotTable=Table[Plot[variables[[i]]/.parameters/.timeCourse,{t,0,100},PlotLegends->variables[[i]],PlotRange->Full],{i,Length[variables]}]