(* 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 = {
A3c[t], A3m[t], N2[t], S1[t], aco[t], aps[t], cys[t], eth[t], hyd[t], oah[t], oxy[t], pap[t], sul[t]
};

initialValues = {
A3c[0] == 1.5, A3m[0] == 1.5, N2[0] == 2.0, S1[0] == 1.5, aco[0] == 0.3, aps[0] == 0.5, cys[0] == 0.3, eth[0] == 4.0, hyd[0] == 0.5, oah[0] == 1.5, oxy[0] == 7.0, pap[0] == 0.4, sul[0] == 0.4
};

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

rateEquations = {
v1 -> (c0*k\[LetterSpace]v0)/(1 + (cys[t]/Kc)^n), v10 -> c0*k\[LetterSpace]v10, v11a -> (c2*k11*N2[t]*oxy[t])/((1 + (hyd[t]/Kh)^m)*(a*N2[t] + oxy[t])), v11a2 -> (c2*k11*N2[t]*oxy[t])/((1 + (hyd[t]/Kh)^m)*(a*N2[t] + oxy[t])), v11b -> (3*A2m*c2*k11*N2[t]*oxy[t])/((A2m + Ka)*(1 + (hyd[t]/Kh)^m)*(a*N2[t] + oxy[t])), v12 -> c1*k12*A3c[t], v13 -> c0*k\[LetterSpace]v13, v14 -> c2*k14*oxy[t], v15 -> c1*k15*aco[t], v16 -> A2c*c2*k16*A3m[t], v17 -> c1*k17*hyd[t], v18 -> c1*k18*oah[t], v2 -> c1*k2*A3c[t]*sul[t], v3 -> c1*k3*A3c[t]*aps[t], v4 -> c1*k4*N2[t]*pap[t], v5 -> c1*k5*hyd[t]*oah[t], v6 -> c1*k6*cys[t], v7 -> c1*k7*N1*eth[t], v8 -> c2*k8*S2*aco[t], v9 -> c2*k9*N1*S1[t], vLEAK -> 0
};

parameters = {
Ac -> 2.0, Am -> 2.0, Ka -> 1.0, Kc -> 0.1, Kh -> 0.5, N -> 2.0, S -> 2.0, a -> 0.1, k11 -> 10.0, k12 -> 5.0, k14 -> 10.0, k15 -> 5.0, k16 -> 10.0, k17 -> 0.02, k18 -> 1.0, k2 -> 0.2, k3 -> 0.2, k4 -> 0.2, k5 -> 0.1, k6 -> 0.12, k7 -> 10.0, k8 -> 10.0, k9 -> 10.0, k\[LetterSpace]v0 -> 1.6, k\[LetterSpace]v10 -> 80.0, k\[LetterSpace]v13 -> 4.0, m -> 4.0, n -> 4.0, C1 -> 0.0, C2 -> 0.0, H2O -> 0.0, Hm -> 0.0, Ho -> 0.0, PPi -> 0.0, eth\[LetterSpace]ex -> 0.0, oxy\[LetterSpace]ex -> 0.0, sul\[LetterSpace]ex -> 0.0, c0 -> 1.0, c1 -> 1.0, c2 -> 1.0
};

assignments = {
S2 -> S - S1[t], N1 -> N - N2[t], A2m -> Am - A3m[t], A2c -> Ac - A3c[t]
};

events = {

};

speciesAnnotations = {
A2c[t]->"http://identifiers.org/obo.chebi/CHEBI:16761", A2m[t]->"http://identifiers.org/obo.chebi/CHEBI:16761", A3c[t]->"http://identifiers.org/obo.chebi/CHEBI:15422", A3m[t]->"http://identifiers.org/obo.chebi/CHEBI:15422", H2O[t]->"http://identifiers.org/obo.chebi/CHEBI:15377", Hm[t]->"http://identifiers.org/obo.chebi/CHEBI:24636", Ho[t]->"http://identifiers.org/obo.chebi/CHEBI:24636", N1[t]->"http://identifiers.org/obo.chebi/CHEBI:15846", N2[t]->"http://identifiers.org/obo.chebi/CHEBI:16908", PPi[t]->"http://identifiers.org/obo.chebi/CHEBI:18361", aco[t]->"http://identifiers.org/obo.chebi/CHEBI:15351", aps[t]->"http://identifiers.org/obo.chebi/CHEBI:17709", cys[t]->"http://identifiers.org/obo.chebi/CHEBI:17561", eth[t]->"http://identifiers.org/obo.chebi/CHEBI:16236", eth\[LetterSpace]ex[t]->"http://identifiers.org/obo.chebi/CHEBI:16236", hyd[t]->"http://identifiers.org/obo.chebi/CHEBI:16136", oah[t]->"http://identifiers.org/obo.chebi/CHEBI:16288", oxy[t]->"http://identifiers.org/obo.chebi/CHEBI:15379", oxy\[LetterSpace]ex[t]->"http://identifiers.org/obo.chebi/CHEBI:15379", pap[t]->"http://identifiers.org/obo.chebi/CHEBI:17980", sul[t]->"http://identifiers.org/obo.chebi/CHEBI:16189", sul\[LetterSpace]ex[t]->"http://identifiers.org/obo.chebi/CHEBI:16189"
};

reactionAnnotations = {
v11a->"http://identifiers.org/obo.go/GO:0015990", v11a->"http://identifiers.org/obo.go/GO:0042775", v11a->"http://identifiers.org/obo.go/GO:0002082", v15->"http://identifiers.org/ec-code/2.3.1.31", v16->"http://identifiers.org/obo.go/GO:0005471", v2->"http://identifiers.org/ec-code/2.7.7.4", v3->"http://identifiers.org/ec-code/2.7.1.25", v3->"http://identifiers.org/obo.go/GO:0004020"
};

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

(* Time evolution *)
odes = {
A3c'[t] == 1.0*v16 -1.0*v2 -1.0*v3 -1.0*v12, A3m'[t] == 1.0*v11b -1.0*v16, N2'[t] == 2.0*v7 +4.0*v9 -3.0*v4 -1.0*v11a, S1'[t] == 1.0*v8 -1.0*v9, aco'[t] == 1.0*v7 -1.0*v15 -1.0*v8, aps'[t] == 1.0*v2 -1.0*v3, cys'[t] == 1.0*v5 -1.0*v6, eth'[t] == 1.0*v13 -1.0*v7, hyd'[t] == 1.0*v4 -1.0*v5 -1.0*v17, oah'[t] == 1.0*v15 -1.0*v5 -1.0*v18, oxy'[t] == 1.0*v10 -1.0*v14 -1.0*v11a2, pap'[t] == 1.0*v3 -1.0*v4, sul'[t] == 1.0*v1 -1.0*v2
};
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]}]