(* 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 = {
ATP[t], AcCoA[t], Cit[t], KG[t], NAD[t], OAA[t], Pyr[t], DeltaPsi[t]
};

initialValues = {
DeltaPsi[0] == 150, ATP[0] == 3.536, AcCoA[0] == 0.063, Cit[0] == 0.44, KG[0] == 0.225, NAD[0] == 0.856, OAA[0] == 0.005, Pyr[0] == 0.154
};

rates = {
v1, v2, v3, v4, v5, v6, v7, v8, vANT, vATP, vleak, vresp
};

rateEquations = {
v1 -> mitochondrion*v1\[LetterSpace]k1, v2 -> mitochondrion*v2\[LetterSpace]k2*NAD[t]*Pyr[t], v3 -> mitochondrion*v3\[LetterSpace]k3*AcCoA[t]*OAA[t], v4 -> mitochondrion*v4\[LetterSpace]k4*Cit[t]*NAD[t], v5 -> mitochondrion*v5\[LetterSpace]k5*(At - ATP[t])*KG[t]*NAD[t], v6 -> mitochondrion*v6\[LetterSpace]k6*(-(KG[t]/v6\[LetterSpace]Keq) + OAA[t]), v7 -> mitochondrion*v7\[LetterSpace]k7*ATP[t]*Pyr[t], v8 -> mitochondrion*v8\[LetterSpace]k8*OAA[t], vANT -> JANT*mitochondrion, vATP -> JATP*mitochondrion, vleak -> Jleak*mitochondrion, vresp -> Jresp*mitochondrion
};

parameters = {
At -> 4.16, C -> 6.75*^-06, DeltaPsim -> 150.0, F -> 96485.0, K -> 2.0, Kapp -> 4.4*^-09, Nt -> 1.07, R -> 8314.0, T -> 298.0, a -> 0.1, b -> 0.004, kANT -> 0.05387, kATP -> 131.9, kleak -> 0.000426, kresp -> 2.5, H -> 1.0, H2O -> 1.0, He -> 1.0, O2 -> 1.0, iP -> 2.44, v1\[LetterSpace]k1 -> 0.038, v2\[LetterSpace]k2 -> 0.152, v3\[LetterSpace]k3 -> 57.142, v4\[LetterSpace]k4 -> 0.053, v5\[LetterSpace]k5 -> 0.082361, v6\[LetterSpace]k6 -> 0.0032, v6\[LetterSpace]Keq -> 0.3975, v7\[LetterSpace]k7 -> 0.04, v8\[LetterSpace]k8 -> 3.6, cytoplasm -> 1.0, mitochondrion -> 1.0
};

assignments = {
NADH -> Nt - NAD[t], ADP -> At - ATP[t], DeltaGtransport -> 1.2*F*DeltaPsi[t], ATPcrit -> At/(1 + 1/(E^((3*DeltaGtransport)/(R*T))*iP*Kapp)), JATP -> (-1 + 2/(1 + E^(b*(-ATPcrit + ATP[t]))))*kATP, Jresp -> (kresp*(Nt - NAD[t]))/((1 + E^(a*(-DeltaPsim + DeltaPsi[t])))*(K + Nt - NAD[t])), Jleak -> kleak*DeltaPsi[t], JANT -> kANT*ATP[t]
};

events = {

};

speciesAnnotations = {

};

reactionAnnotations = {

};

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

(* Time evolution *)
odes = {
ATP'[t] == 1.0*v5 +1.0*vATP -1.0*v7 -1.0*vANT, AcCoA'[t] == 1.0*v2 -1.0*v3, Cit'[t] == 1.0*v3 -1.0*v4, KG'[t] == 1.0*v4 +1.0*v6 -1.0*v5, NAD'[t] == 1.0*vresp -1.0*v2 -1.0*v4 -2.0*v5, OAA'[t] == 1.0*v5 +1.0*v7 -1.0*v3 -1.0*v6 -1.0*v8, Pyr'[t] == 1.0*v1 -1.0*v2 -1.0*v7, DeltaPsi'[t] == (-JANT - 3*JATP - Jleak + 10*Jresp)/C
};
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]}]