(* 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 = {
cLc[t], cLm[t], cLn[t], cPn[t], cTc[t], cTm[t], cTn[t], cXc[t], cXm[t], cXn[t], cYc[t], cYm[t], cYn[t]
};

initialValues = {
cLc[0] == 0.015, cLm[0] == 0.539, cLn[0] == 0.0855, cPn[0] == 0.835, cTc[0] == 8.7, cTm[0] == 0.456, cTn[0] == 0.4, cXc[0] == 1.34, cXm[0] == 0.08, cXn[0] == 0.4, cYc[0] == 0.004, cYm[0] == 0.037, cYn[0] == 0.013
};

rates = {
four12, four15, four18, four21, one11, one12, one13, one14, one15, one16, one17, one18, one19, one20, one21, one22, one23, three11, three12, three13, three15, three16, three18, three19, three21, three22, three23, two11, two14, two17, two20, two23
};

rateEquations = {
four12 -> (compartment*m2*cLc[t])/(k2 + cLc[t]), four15 -> (compartment*cTc[t]*(m6 + m5*(1 - LD[t, dayLength])))/(k5 + cTc[t]), four18 -> (compartment*m10*cXc[t])/(k8 + cXc[t]), four21 -> (compartment*m13*cYc[t])/(k11 + cYc[t]), one11 -> compartment*q1*cPn[t]*LD[t, dayLength], one12 -> compartment*p1*cLm[t], one13 -> compartment*r1*cLc[t], one14 -> (compartment*g3^c*n2*cYn[t]^b)/((g3^c + cLn[t]^c)*(g2^b + cYn[t]^b)), one15 -> compartment*p2*cTm[t], one16 -> compartment*r3*cTc[t], one17 -> (compartment*n3*cTn[t]^d)/(g4^d + cTn[t]^d), one18 -> compartment*p3*cXm[t], one19 -> compartment*r5*cXc[t], one20 -> (compartment*g6^f*(q2*cPn[t]*LD[t, dayLength] + (g5^e*(n5 + n4*LD[t, dayLength]))/(g5^e + cTn[t]^e)))/(g6^f + cLn[t]^f), one21 -> compartment*p4*cYm[t], one22 -> compartment*r7*cYc[t], one23 -> compartment*p5*(1 - LD[t, dayLength]), three11 -> (compartment*m1*cLm[t])/(k1 + cLm[t]), three12 -> compartment*r2*cLn[t], three13 -> (compartment*m3*cLn[t])/(k3 + cLn[t]), three15 -> compartment*r4*cTn[t], three16 -> (compartment*cTn[t]*(m8 + m7*(1 - LD[t, dayLength])))/(k6 + cTn[t]), three18 -> compartment*r6*cXn[t], three19 -> (compartment*m11*cXn[t])/(k9 + cXn[t]), three21 -> compartment*r8*cYn[t], three22 -> (compartment*m14*cYn[t])/(k12 + cYn[t]), three23 -> compartment*q3*cPn[t]*LD[t, dayLength], two11 -> (compartment*n1*cXn[t]^a)/(g1^a + cXn[t]^a), two14 -> (compartment*m4*cTm[t])/(k4 + cTm[t]), two17 -> (compartment*m9*cXm[t])/(k7 + cXm[t]), two20 -> (compartment*m12*cYm[t])/(k10 + cYm[t]), two23 -> (compartment*m15*cPn[t])/(k13 + cPn[t])
};

parameters = {
CP -> 24.0, Lmax -> 1.0, Lmin -> 0.0, a -> 3.3064, b -> 1.0258, c -> 1.0258, d -> 1.4422, dayLength -> 12.0, e -> 3.6064, f -> 1.0237, g1 -> 0.876738488, g2 -> 0.036805783, g3 -> 0.26593318, g4 -> 0.538811228, g5 -> 1.17803247, g6 -> 0.064455137, k1 -> 1.817, k10 -> 1.7303, k11 -> 1.8258, k12 -> 1.8066, k13 -> 1.2, k2 -> 1.5644, k3 -> 1.2765, k4 -> 2.5734, k5 -> 2.7454, k6 -> 0.4033, k7 -> 6.5585, k8 -> 0.6632, k9 -> 17.1111, m1 -> 1.5283, m10 -> 0.2179, m11 -> 3.3442, m12 -> 4.297, m13 -> 0.1347, m14 -> 0.6114, m15 -> 1.2, m2 -> 20.44, m3 -> 3.6888, m4 -> 3.8231, m5 -> 0.0013, m6 -> 3.1741, m7 -> 0.0492, m8 -> 4.0424, m9 -> 10.1132, n1 -> 5.1694, n2 -> 3.0087, n3 -> 0.2431, n4 -> 0.0857, n5 -> 0.1649, p1 -> 0.8295, p2 -> 4.324, p3 -> 2.147, p4 -> 0.2485, p5 -> 0.5, q1 -> 2.4514, q2 -> 2.40178, q3 -> 1.0, r1 -> 16.8363, r2 -> 0.1687, r3 -> 0.3166, r4 -> 2.1509, r5 -> 1.0352, r6 -> 3.3017, r7 -> 2.2123, r8 -> 0.2002, two11\[LetterSpace]Fch\[LetterSpace]0 -> 8.0, four12\[LetterSpace]Fch\[LetterSpace]0 -> 8.0, compartment -> 1.0
};

assignments = {
LD[tod_,length_] -> Ceiling[Sin[0.001 + (Pi*tod)/length]/2]
};

events = {

};

speciesAnnotations = {
cLc[t]->"http://identifiers.org/uniprot/O81713", cLn[t]->"http://identifiers.org/uniprot/O81713", cTc[t]->"http://identifiers.org/uniprot/Q9LKL2", cTn[t]->"http://identifiers.org/uniprot/Q9LKL2"
};

reactionAnnotations = {

};

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

(* Time evolution *)
odes = {
cLc'[t] == 1.0*one12 +1.0*three12 -1.0*four12 -1.0*one13, cLm'[t] == 1.0*two11 +1.0*one11 -1.0*three11, cLn'[t] == 1.0*one13 -1.0*three12 -1.0*three13, cPn'[t] == 1.0*one23 -1.0*two23 -1.0*three23, cTc'[t] == 1.0*one15 +1.0*three15 -1.0*one16 -1.0*four15, cTm'[t] == 1.0*one14 -1.0*two14, cTn'[t] == 1.0*one16 -1.0*three15 -1.0*three16, cXc'[t] == 1.0*one18 +1.0*three18 -1.0*four18 -1.0*one19, cXm'[t] == 1.0*one17 -1.0*two17, cXn'[t] == 1.0*one19 -1.0*three18 -1.0*three19, cYc'[t] == 1.0*one21 +1.0*three21 -1.0*four21 -1.0*one22, cYm'[t] == 1.0*one20 -1.0*two20, cYn'[t] == 1.0*one22 -1.0*three21 -1.0*three22
};
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]}]