(* 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 = {
NADPH[t], PGA[t]
};

initialValues = {
NADPH[0] == 0.21, PGA[0] == 2.4
};

rates = {
NADPH\[LetterSpace]prod, PGA\[LetterSpace]cons, PGA\[LetterSpace]prod\[LetterSpace]Vc, PGA\[LetterSpace]prod\[LetterSpace]Vo
};

rateEquations = {
NADPH\[LetterSpace]prod -> chloroplast*NADPH\[LetterSpace]production[j, NADP, Nt], PGA\[LetterSpace]cons -> chloroplast*PGA\[LetterSpace]consumption[PGA[t], Rp, Nt, kc, Epar, NADPH[t]], PGA\[LetterSpace]prod\[LetterSpace]Vc -> chloroplast*carboxylation[kc, CO2, Kc, O2, Ko, Epar, RuBP], PGA\[LetterSpace]prod\[LetterSpace]Vo -> chloroplast*oxygenation[phi, kc, CO2, Kc, O2, Ko, Epar, RuBP]
};

parameters = {
Epar -> 1.33846153846, Ipar -> 1000.0, Kc -> 460.0, Ko -> 330.0, Nt -> 0.5, Rp -> 3.2, f -> 0.23, kc -> 2.5, CO2 -> 230.0, O2 -> 210.0, RuBP -> 2.0, chloroplast -> 1.0
};

assignments = {
oxygenation[phi_,kc_,S2_,K_,R_,KI_,Epar_,S1_] -> (phi*((Epar*kc*S2)/(K*(1 + R/KI) + S2) + (kc*S1*S2)/(K*(1 + R/KI) + S2) - Abs[(Epar*kc*S2)/(K*(1 + R/KI) + S2) - (kc*S1*S2)/(K*(1 + R/KI) + S2)]))/2, PGA\[LetterSpace]consumption[S1_,Rp_,Nt_,kc_,Epar_,R_] -> (Epar*kc*R*S1)/(Nt*Rp), carboxylation[kc_,S2_,K_,R_,KI_,Epar_,S1_] -> ((Epar*kc*S2)/(K*(1 + R/KI) + S2) + (kc*S1*S2)/(K*(1 + R/KI) + S2) - Abs[(Epar*kc*S2)/(K*(1 + R/KI) + S2) - (kc*S1*S2)/(K*(1 + R/KI) + S2)])/2, NADPH\[LetterSpace]production[j_,S1_,Nt_] -> (j*S1)/(2*Nt), NADP -> Nt - NADPH[t], jp -> ((1 - f)*Ipar)/2, phi -> (0.21*Kc*O2)/(CO2*Ko), j -> (-Abs[-467 + jp] + Abs[467 + jp])/130
};

events = {

};

speciesAnnotations = {

};

reactionAnnotations = {

};

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

(* Time evolution *)
odes = {
NADPH'[t] == 1.0*NADPH\[LetterSpace]prod -2.0*PGA\[LetterSpace]prod\[LetterSpace]Vc -2.0*PGA\[LetterSpace]prod\[LetterSpace]Vo, PGA'[t] == 2.0*PGA\[LetterSpace]prod\[LetterSpace]Vc +1.5*PGA\[LetterSpace]prod\[LetterSpace]Vo -1.0*PGA\[LetterSpace]cons
};
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]}]