(* 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 = {
FeII[t], FeIII[t], FeIII\[LetterSpace]NO[t], FeIII\[LetterSpace]star[t], FeII\[LetterSpace]NO[t], FeII\[LetterSpace]O2[t], FeII\[LetterSpace]star[t], FeII\[LetterSpace]star\[LetterSpace]O2[t]
};

initialValues = {
FeII[0] == 0.0, FeIII[0] == 1.0, FeIII\[LetterSpace]NO[0] == 0.0, FeIII\[LetterSpace]star[0] == 0.0, FeII\[LetterSpace]NO[0] == 0.0, FeII\[LetterSpace]O2[0] == 0.0, FeII\[LetterSpace]star[0] == 0.0, FeII\[LetterSpace]star\[LetterSpace]O2[0] == 0.0
};

rates = {
r1, r10, r2, r3, r4, r5, r6, r9, rF, rG
};

rateEquations = {
r1 -> cytosol*k1*FeIII[t], r10 -> cytosol*k10*O2*FeII\[LetterSpace]NO[t], r2 -> cytosol*k2*O2*FeII[t], r3 -> cytosol*k3*FeII\[LetterSpace]O2[t], r4 -> cytosol*k4*FeIII\[LetterSpace]star[t], r5 -> cytosol*k5*O2*FeII\[LetterSpace]star[t], r6 -> cytosol*k6*FeII\[LetterSpace]star\[LetterSpace]O2[t], r9 -> cytosol*k9*FeII\[LetterSpace]NO[t], rF -> cytosol*k7*FeIII\[LetterSpace]NO[t], rG -> cytosol*k8*FeIII\[LetterSpace]NO[t]
};

parameters = {
k1 -> 2.6, k10 -> 0.0013, k2 -> 0.9, k3 -> 26.0, k4 -> 2.6, k5 -> 0.9, k6 -> 26.0, k7 -> 5.0, k8 -> 2.6, k9 -> 0.0001, NADPH -> 40.0, NADPplus -> 0.0, NO -> 0.0, NO3 -> 0.0, O2 -> 140.0, citrulline -> 0.0, cytosol -> 1.0
};

assignments = {
FeIII\[LetterSpace]t -> FeIII[t] + FeIII\[LetterSpace]star[t]
};

events = {

};

speciesAnnotations = {
NO[t]->"http://identifiers.org/chebi/CHEBI:16480", NO[t]->"http://identifiers.org/kegg.compound/C00533", NO3[t]->"http://identifiers.org/chebi/CHEBI:17632", NO3[t]->"http://identifiers.org/kegg.compound/C00244"
};

reactionAnnotations = {

};

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

(* Time evolution *)
odes = {
FeII'[t] == 1.0*r1 +1.0*r9 -1.0*r2, FeIII'[t] == 1.0*rF +1.0*r10 -1.0*r1, FeIII\[LetterSpace]NO'[t] == 1.0*r6 -1.0*rF -1.0*rG, FeIII\[LetterSpace]star'[t] == 1.0*r3 -1.0*r4, FeII\[LetterSpace]NO'[t] == 1.0*rG -1.0*r9 -1.0*r10, FeII\[LetterSpace]O2'[t] == 1.0*r2 -1.0*r3, FeII\[LetterSpace]star'[t] == 1.0*r4 -1.0*r5, FeII\[LetterSpace]star\[LetterSpace]O2'[t] == 1.0*r5 -1.0*r6
};
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]}]