(* 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 = {
ACE[t], BPG[t], F16P[t], F6P[t], G6P[t], GLCi[t], NADH[t], P2G[t], P3G[t], PEP[t], PHOS[t], PYR[t], Prb[t], TRIO[t]
};

initialValues = {
ACE[0] == 0.0474837, BPG[0] == 0.111, F16P[0] == 0.836, F6P[0] == 0.75247, G6P[0] == 3.085, GLCi[0] == 0.087, NADH[0] == 0.044, P2G[0] == 0.13771, P3G[0] == 0.825, PEP[0] == 0.1404, PHOS[0] == 10.0, PYR[0] == 0.884031, Prb[0] == 5.0, TRIO[0] == 0.5177
};

rates = {
v1, v10, v11, v12, v13, v14, v15, v16, v17, v18, v2, v3, v4, v5, v6, v7, v8, v9
};

rateEquations = {
v1 -> (VMAXHK*((GLCi[t]*(-SUMAXP + (1 - 4*KEQAK)*Prb[t] + (4*(1 - 4*KEQAK)*KEQAK*Prb[t]^2 + (SUMAXP - (1 - 4*KEQAK)*Prb[t])^2)^0.5))/(2*(1 - 4*KEQAK)*KMHKATP*KMHKGLCi) - (G6P[t]*(Prb[t] - (-SUMAXP + (1 - 4*KEQAK)*Prb[t] + (4*(1 - 4*KEQAK)*KEQAK*Prb[t]^2 + (SUMAXP - (1 - 4*KEQAK)*Prb[t])^2)^0.5)/(1 - 4*KEQAK)))/(KEQHK*KMHKATP*KMHKGLCi)))/(compartment*(1 + G6P[t]/KMHKG6P + (wt*G6P[t])/KiHKG6P + GLCi[t]/KMHKGLCi)*(1 + (-SUMAXP + (1 - 4*KEQAK)*Prb[t] + (4*(1 - 4*KEQAK)*KEQAK*Prb[t]^2 + (SUMAXP - (1 - 4*KEQAK)*Prb[t])^2)^0.5)/(2*(1 - 4*KEQAK)*KMHKATP) + (Prb[t] - (-SUMAXP + (1 - 4*KEQAK)*Prb[t] + (4*(1 - 4*KEQAK)*KEQAK*Prb[t]^2 + (SUMAXP - (1 - 4*KEQAK)*Prb[t])^2)^0.5)/(1 - 4*KEQAK))/KMHKADP)), v10 -> (VmENO*(P2G[t] - PEP[t]/KeqENO))/(compartment*KmENOP2G*(1 + P2G[t]/KmENOP2G + PEP[t]/KmENOPEP)), v11 -> (VmPYK*PEP[t]*(1 + PEP[t]/KmPYKPEP)^(-1 + npyk)*(Prb[t] - (-SUMAXP + (1 - 4*KeqAK)*Prb[t] + (4*(1 - 4*KeqAK)*KeqAK*Prb[t]^2 + (SUMAXP - (1 - 4*KeqAK)*Prb[t])^2)^0.5)/(1 - 4*KeqAK)))/(compartment*KmPYKPEP*(KmPYKADP + Prb[t] - (-SUMAXP + (1 - 4*KeqAK)*Prb[t] + (4*(1 - 4*KeqAK)*KeqAK*Prb[t]^2 + (SUMAXP - (1 - 4*KeqAK)*Prb[t])^2)^0.5)/(1 - 4*KeqAK))*((1 + PEP[t]/KmPYKPEP)^npyk + L0PYK*((1 + (-SUMAXP + (1 - 4*KeqAK)*Prb[t] + (4*(1 - 4*KeqAK)*KeqAK*Prb[t]^2 + (SUMAXP - (1 - 4*KeqAK)*Prb[t])^2)^0.5)/(2*(1 - 4*KeqAK)*KmPYKATP))/(1 + F16P[t]/KmPYKF16P))^npyk)), v12 -> (VmPDC*PYR[t]^nPDC)/(compartment*KmPDCPYR^nPDC*(1 + PYR[t]^nPDC/KmPDCPYR^nPDC)), v13 -> (KSUCC*ACE[t])/compartment, v14 -> (VmGLT*(GLCo - GLCi[t]/KeqGLT))/(compartment*KmGLTGLCo*(1 + GLCo/KmGLTGLCo + GLCi[t]/KmGLTGLCi + (0.91*GLCo*GLCi[t])/(KmGLTGLCi*KmGLTGLCo))), v15 -> -((VmADH*(ETOH*(NADt - NADH[t]) - (ACE[t]*NADH[t])/KeqADH))/(compartment*KiADHNAD*KmADHETOH*(1 + (ETOH*KmADHNAD)/(KiADHNAD*KmADHETOH) + (KmADHNADH*ACE[t])/(KiADHNADH*KmADHACE) + (NADt - NADH[t])/KiADHNAD + (ETOH*(NADt - NADH[t]))/(KiADHNAD*KmADHETOH) + (ETOH*ACE[t]*(NADt - NADH[t]))/(KiADHACE*KiADHNAD*KmADHETOH) + (KmADHNADH*ACE[t]*(NADt - NADH[t]))/(KiADHNAD*KiADHNADH*KmADHACE) + NADH[t]/KiADHNADH + (ETOH*KmADHNAD*NADH[t])/(KiADHNAD*KiADHNADH*KmADHETOH) + (ACE[t]*NADH[t])/(KiADHNADH*KmADHACE) + (ETOH*ACE[t]*NADH[t])/(KiADHETOH*KiADHNADH*KmADHACE)))), v16 -> (VmG3PDH*(-((GLY*(NADt - NADH[t]))/KeqG3PDH) + (NADH[t]*TRIO[t])/(1 + KeqTPI)))/(compartment*KmG3PDHDHAP*KmG3PDHNADH*(1 + (NADt - NADH[t])/KmG3PDHNAD + NADH[t]/KmG3PDHNADH)*(1 + GLY/KmG3PDHGLY + TRIO[t]/((1 + KeqTPI)*KmG3PDHDHAP))), v17 -> (KATPASE*(-SUMAXP + Prb[t] - 4*KeqAK*Prb[t] + (SUMAXP^2 - 2*SUMAXP*Prb[t] + 8*KeqAK*SUMAXP*Prb[t] + Prb[t]^2 - 4*KeqAK*Prb[t]^2)^0.5))/(compartment*(2 - 8*KeqAK)), v18 -> (KPHOS*(pT - PHOS[t]))/compartment, v2 -> (VmPGI*(-(F6P[t]/KeqPGI) + G6P[t]))/(compartment*KmPGIG6P*(1 + F6P[t]/KmPGIF6P + G6P[t]/KmPGIG6P)), v3 -> (KmGLYCOGEN*G6P[t])/compartment, v4 -> (KTREHALOSE*wt*G6P[t])/compartment, v5 -> (gR*VmPFK*F6P[t]*(-SUMAXP + Prb[t] - 4*KeqAK*Prb[t] + (SUMAXP^2 - 2*SUMAXP*Prb[t] + 8*KeqAK*SUMAXP*Prb[t] + Prb[t]^2 - 4*KeqAK*Prb[t]^2)^0.5)*(1 + F6P[t]/KmPFKF6P + (-SUMAXP + Prb[t] - 4*KeqAK*Prb[t] + (SUMAXP^2 - 2*SUMAXP*Prb[t] + 8*KeqAK*SUMAXP*Prb[t] + Prb[t]^2 - 4*KeqAK*Prb[t]^2)^0.5)/((2 - 8*KeqAK)*KmPFKATP) + (gR*F6P[t]*(-SUMAXP + Prb[t] - 4*KeqAK*Prb[t] + (SUMAXP^2 - 2*SUMAXP*Prb[t] + 8*KeqAK*SUMAXP*Prb[t] + Prb[t]^2 - 4*KeqAK*Prb[t]^2)^0.5))/((2 - 8*KeqAK)*KmPFKATP*KmPFKF6P)))/(compartment*(2 - 8*KeqAK)*KmPFKATP*KmPFKF6P*((L0*(1 + (CPFKF26BP*F26BP)/KPFKF26BP + (CPFKF16BP*F16P[t])/KPFKF16BP)^2*(1 + (2*CPFKAMP*KeqAK*(SUMAXP - (SUMAXP^2 - 2*SUMAXP*Prb[t] + 8*KeqAK*SUMAXP*Prb[t] + Prb[t]^2 - 4*KeqAK*Prb[t]^2)^0.5)^2)/((-1 + 4*KeqAK)*KPFKAMP*(SUMAXP - Prb[t] + 4*KeqAK*Prb[t] - (SUMAXP^2 - 2*SUMAXP*Prb[t] + 8*KeqAK*SUMAXP*Prb[t] + Prb[t]^2 - 4*KeqAK*Prb[t]^2)^0.5)))^2*(1 + (CiPFKATP*(-SUMAXP + Prb[t] - 4*KeqAK*Prb[t] + (SUMAXP^2 - 2*SUMAXP*Prb[t] + 8*KeqAK*SUMAXP*Prb[t] + Prb[t]^2 - 4*KeqAK*Prb[t]^2)^0.5))/((2 - 8*KeqAK)*KiPFKATP))^2*(1 + (CPFKATP*(-SUMAXP + Prb[t] - 4*KeqAK*Prb[t] + (SUMAXP^2 - 2*SUMAXP*Prb[t] + 8*KeqAK*SUMAXP*Prb[t] + Prb[t]^2 - 4*KeqAK*Prb[t]^2)^0.5))/((2 - 8*KeqAK)*KmPFKATP))^2)/((1 + F26BP/KPFKF26BP + F16P[t]/KPFKF16BP)^2*(1 + (2*KeqAK*(SUMAXP - (SUMAXP^2 - 2*SUMAXP*Prb[t] + 8*KeqAK*SUMAXP*Prb[t] + Prb[t]^2 - 4*KeqAK*Prb[t]^2)^0.5)^2)/((-1 + 4*KeqAK)*KPFKAMP*(SUMAXP - Prb[t] + 4*KeqAK*Prb[t] - (SUMAXP^2 - 2*SUMAXP*Prb[t] + 8*KeqAK*SUMAXP*Prb[t] + Prb[t]^2 - 4*KeqAK*Prb[t]^2)^0.5)))^2*(1 + (-SUMAXP + Prb[t] - 4*KeqAK*Prb[t] + (SUMAXP^2 - 2*SUMAXP*Prb[t] + 8*KeqAK*SUMAXP*Prb[t] + Prb[t]^2 - 4*KeqAK*Prb[t]^2)^0.5)/((2 - 8*KeqAK)*KiPFKATP))^2) + (1 + F6P[t]/KmPFKF6P + (-SUMAXP + Prb[t] - 4*KeqAK*Prb[t] + (SUMAXP^2 - 2*SUMAXP*Prb[t] + 8*KeqAK*SUMAXP*Prb[t] + Prb[t]^2 - 4*KeqAK*Prb[t]^2)^0.5)/((2 - 8*KeqAK)*KmPFKATP) + (gR*F6P[t]*(-SUMAXP + Prb[t] - 4*KeqAK*Prb[t] + (SUMAXP^2 - 2*SUMAXP*Prb[t] + 8*KeqAK*SUMAXP*Prb[t] + Prb[t]^2 - 4*KeqAK*Prb[t]^2)^0.5))/((2 - 8*KeqAK)*KmPFKATP*KmPFKF6P))^2)), v6 -> (VmALD*(F16P[t] - (KeqTPI*TRIO[t]^2)/(KeqALD*(1 + KeqTPI)^2)))/(compartment*KmALDF16P*(1 + F16P[t]/KmALDF16P + TRIO[t]/((1 + KeqTPI)*KmALDDHAP) + (KeqTPI*TRIO[t])/((1 + KeqTPI)*KmALDGAP) + (KeqTPI*F16P[t]*TRIO[t])/((1 + KeqTPI)*KmALDF16P*KmALDGAPi) + (KeqTPI*TRIO[t]^2)/((1 + KeqTPI)^2*KmALDDHAP*KmALDGAP))), v7 -> (-((VmGAPDHr*BPG[t]*NADH[t])/(KmGAPDHBPG*KmGAPDHNADH)) + (KeqTPI*VmGAPDHf*(NADt - NADH[t])*PHOS[t]*TRIO[t])/((1 + KeqTPI)*KmGAPDHGAP*KmGAPDHNAD))/(compartment*(1 + (NADt - NADH[t])/KmGAPDHNAD + NADH[t]/KmGAPDHNADH)*(1 + PHOS[t])*(1 + BPG[t]/KmGAPDHBPG + (KeqTPI*TRIO[t])/((1 + KeqTPI)*KmGAPDHGAP))), v8 -> (VmPGK*((KeqPGK*BPG[t]*(SUMAXP - (SUMAXP^2 - 2*SUMAXP*Prb[t] + 8*KeqAK*SUMAXP*Prb[t] + Prb[t]^2 - 4*KeqAK*Prb[t]^2)^0.5))/(1 - 4*KeqAK) - (P3G[t]*(-SUMAXP + Prb[t] - 4*KeqAK*Prb[t] + (SUMAXP^2 - 2*SUMAXP*Prb[t] + 8*KeqAK*SUMAXP*Prb[t] + Prb[t]^2 - 4*KeqAK*Prb[t]^2)^0.5))/(2 - 8*KeqAK)))/(compartment*KmPGKATP*KmPGKP3G*(1 + BPG[t]/KmPGKBPG + P3G[t]/KmPGKP3G)*(1 + (SUMAXP - (SUMAXP^2 - 2*SUMAXP*Prb[t] + 8*KeqAK*SUMAXP*Prb[t] + Prb[t]^2 - 4*KeqAK*Prb[t]^2)^0.5)/((1 - 4*KeqAK)*KmPGKADP) + (-SUMAXP + Prb[t] - 4*KeqAK*Prb[t] + (SUMAXP^2 - 2*SUMAXP*Prb[t] + 8*KeqAK*SUMAXP*Prb[t] + Prb[t]^2 - 4*KeqAK*Prb[t]^2)^0.5)/((2 - 8*KeqAK)*KmPGKATP))), v9 -> (VmPGM*(-(P2G[t]/KeqPGM) + P3G[t]))/(compartment*KmPGMP3G*(1 + P2G[t]/KmPGMP2G + P3G[t]/KmPGMP3G))
};

parameters = {
CPFKAMP -> 0.0845, CPFKATP -> 3.0, CPFKF16BP -> 0.397, CPFKF26BP -> 0.0174, CPFKF6P -> 0.0, CiPFKATP -> 100.0, F26BP -> 0.02, KATPASE -> 39.5, KEQAK -> 0.45, KEQHK -> 3800.0, KMHKADP -> 0.23, KMHKATP -> 0.15, KMHKG6P -> 30.0, KMHKGLCi -> 0.08, KPFKAMP -> 0.0995, KPFKF16BP -> 0.111, KPFKF26BP -> 0.000682, KPHOS -> 0.1, KSUCC -> 21.4, KTREHALOSE -> 2.32, KeqADH -> 6.9*^-05, KeqAK -> 0.45, KeqALD -> 0.069, KeqENO -> 6.7, KeqG3PDH -> 4300.0, KeqGLT -> 1.0, KeqPGI -> 0.314, KeqPGK -> 3200.0, KeqPGM -> 0.19, KeqPYK -> 6500.0, KeqTPI -> 0.045, KiADHACE -> 1.1, KiADHETOH -> 90.0, KiADHNAD -> 0.92, KiADHNADH -> 0.031, KiHKG6P -> 0.07, KiPFKATP -> 0.65, KmADHACE -> 1.11, KmADHETOH -> 17.0, KmADHNAD -> 0.17, KmADHNADH -> 0.11, KmALDDHAP -> 2.4, KmALDF16P -> 0.3, KmALDGAP -> 2.0, KmALDGAPi -> 10.0, KmENOP2G -> 0.04, KmENOPEP -> 0.5, KmG3PDHDHAP -> 0.4, KmG3PDHGLY -> 1.0, KmG3PDHNAD -> 0.93, KmG3PDHNADH -> 0.023, KmGAPDHBPG -> 0.0098, KmGAPDHGAP -> 0.21, KmGAPDHNAD -> 0.09, KmGAPDHNADH -> 0.06, KmGLTGLCi -> 1.1918, KmGLTGLCo -> 1.1918, KmGLYCOGEN -> 5.8, KmPDCPYR -> 4.33, KmPFKATP -> 0.71, KmPFKF6P -> 0.1, KmPGIF6P -> 0.3, KmPGIG6P -> 1.4, KmPGKADP -> 0.2, KmPGKATP -> 0.3, KmPGKBPG -> 0.003, KmPGKP3G -> 0.53, KmPGMP2G -> 0.08, KmPGMP3G -> 1.2, KmPYKADP -> 0.3, KmPYKATP -> 9.3, KmPYKF16P -> 0.2, KmPYKPEP -> 0.19, L0 -> 0.66, L0PYK -> 60000.0, NADt -> 1.59, SUMAXP -> 3.1, VMAXHK -> 226.452, VmADH -> 810.0, VmALD -> 322.258, VmENO -> 365.806, VmG3PDH -> 70.15, VmGAPDHf -> 1184.52, VmGAPDHr -> 6549.68, VmGLT -> 198.0, VmPDC -> 1062.58, VmPFK -> 182.903, VmPGI -> 339.677, VmPGK -> 1306.45, VmPGM -> 2525.81, VmPYK -> 1088.71, gR -> 5.12, nPDC -> 1.9, npyk -> 4.0, pT -> 10.0, wt -> 1.0, CO2 -> 1.0, ETOH -> 0.0, GLCo -> 110.0, GLY -> 0.15, Glyc -> 0.0, NAD -> 0.0, SUCC -> 0.0, Trh -> 0.0, X -> 0.0, compartment -> 1.0
};

assignments = {
ATP -> (-SUMAXP + (1 - 4*KEQAK)*Prb[t] + (4*(1 - 4*KEQAK)*KEQAK*Prb[t]^2 + (SUMAXP - (1 - 4*KEQAK)*Prb[t])^2)^0.5)/(2*(1 - 4*KEQAK))
};

events = {

};

speciesAnnotations = {

};

reactionAnnotations = {

};

units = {
{"time" -> "", "metabolite" -> "", "extent" -> ""}
};

(* Time evolution *)
odes = {
ACE'[t] == 1.0*v12 -2.0*v13 -1.0*v15, BPG'[t] == 1.0*v7 -1.0*v8, F16P'[t] == 1.0*v5 -1.0*v6, F6P'[t] == 1.0*v2 -1.0*v5, G6P'[t] == 1.0*v1 -1.0*v2 -1.0*v3 -2.0*v4, GLCi'[t] == 1.0*v14 -1.0*v1, NADH'[t] == 1.0*v7 +3.0*v13 -1.0*v15 -1.0*v16, P2G'[t] == 1.0*v9 -1.0*v10, P3G'[t] == 1.0*v8 -1.0*v9, PEP'[t] == 1.0*v10 -1.0*v11, PHOS'[t] == 2.0*v3 +3.0*v4 +1.0*v16 +1.0*v17 +1.0*v18 -1.0*v7, PYR'[t] == 1.0*v11 -1.0*v12, Prb'[t] == 1.0*v8 +1.0*v11 -1.0*v1 -1.0*v3 -1.0*v4 -1.0*v5 -1.0*v17, TRIO'[t] == 2.0*v6 -1.0*v7 -1.0*v16
};
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]}]