(* 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 = {
E1F[t], E1Fp[t], E1W[t], E2F[t], E2Fp[t], E2W[t], MF[t], MW[t], PF[t], PFp[t], PW[t], PWL[t]
};

initialValues = {
E1F[0] == 0.43076, E1Fp[0] == 0.45583, E1W[0] == 5.84748, E2F[0] == 0.10647, E2Fp[0] == 0.09872, E2W[0] == 5.70265, MF[0] == 0.6935, MW[0] == 1.2689, PF[0] == 0.06565, PFp[0] == 0.07719, PW[0] == 26.4393, PWL[0] == 0.0
};

rates = {
E1Fdeg, E1Fpdeg, E1Fptrl, E1Ftrl, E1Wdeg, E1Wtrl, E2Fdeg, E2Fpdeg, E2Fptrl, E2Ftrl, E2Wdeg, E2Wtrl, MFdeg, MFtrn, MWdeg, MWtrn, PFdeg, PFpdeg, PFptrl, PFtrl, PWLdeg, PWdeg, PWtrl, PWtrs
};

rateEquations = {
E1Fdeg -> gam1*E1F[t], E1Fpdeg -> gam1p*E1Fp[t], E1Fptrl -> a3p*MF[t], E1Ftrl -> a3*MF[t], E1Wdeg -> gam2*E1W[t], E1Wtrl -> MW[t]*(a6 + a7*(PF[t] + PFp[t])), E2Fdeg -> gam1*E2F[t], E2Fpdeg -> gam1p*E2Fp[t], E2Fptrl -> f1p*E1Fp[t], E2Ftrl -> f1*E1F[t], E2Wdeg -> gam2*E2W[t], E2Wtrl -> f2*E1W[t], MFdeg -> (d1*MF[t])/(b5 + MF[t]), MFtrn -> (a2*PW[t]^m)/((1 + (PF[t] + PFp[t])/b3)*(b4^m + PW[t]^m)) + (a1*PWL[t]^n)/((1 + (PF[t] + PFp[t])/b1)*(b2^n + PWL[t]^n)), MWdeg -> (d3*MW[t])/(b8 + MW[t]), MWtrn -> a4 + (a5*PWL[t]^k)/(b7^k + PWL[t]^k), PFdeg -> d2*PF[t], PFpdeg -> d2p*PFp[t], PFptrl -> f1p*E2Fp[t], PFtrl -> f1*E2F[t], PWLdeg -> (d5*PWL[t])/(b10 + PWL[t]), PWdeg -> (d4*PW[t])/(b9 + PW[t]), PWtrl -> f2*E2W[t], PWtrs -> -(r2*PWL[t]) + (amp*r1*PW[t]*(1 + Tanh[24*(-dawn + t - 24*Floor[t/24])])*(1 - Tanh[24*(-dusk + t - 24*Floor[t/24])]))/4
};

parameters = {
a1 -> 24.9736, a2 -> 3.59797, a3 -> 0.2834, a3p -> 0.34578, a4 -> 0.46227, a5 -> 0.02917, a6 -> 0.20695, a7 -> 3.02856, amp -> 0.0, b1 -> 0.00209, b10 -> 93.03664, b2 -> 2.13476, b3 -> 0.08039, b4 -> 0.45859, b5 -> 0.13056, b6 -> 0.0, b7 -> 2.96739, b8 -> 0.11167, b9 -> 81.10381, d1 -> 1.43549, d2 -> 0.21251, d2p -> 0.18191, d3 -> 0.50309, d4 -> 3.36641, d5 -> 0.41085, dawn -> 12.0, dusk -> 24.0, f1 -> 0.09292, f1p -> 0.09588, f2 -> 0.14979, gam1 -> 0.34603, gam1p -> 0.40119, gam2 -> 0.00028, k -> 2.18234, m -> 1.34979, n -> 1.02419, r1 -> 2.71574, r2 -> 35.40005, cytosol -> 1.0, nucleus -> 1.0
};

assignments = {
WC1\[LetterSpace]tot -> E1W[t] + E2W[t] + PW[t] + PWL[t], Frq\[LetterSpace]tot -> lFrq\[LetterSpace]tot + sFrq\[LetterSpace]tot, lFrq\[LetterSpace]tot -> E1Fp[t] + E2Fp[t] + PFp[t], sFrq\[LetterSpace]tot -> E1F[t] + E2F[t] + PF[t]
};

events = {

};

speciesAnnotations = {

};

reactionAnnotations = {

};

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

(* Time evolution *)
odes = {
E1F'[t] == 1.0*E1Ftrl -1.0*E1Fdeg -1.0*E2Ftrl, E1Fp'[t] == 1.0*E1Fptrl -1.0*E1Fpdeg -1.0*E2Fptrl, E1W'[t] == 1.0*E1Wtrl -1.0*E1Wdeg -1.0*E2Wtrl, E2F'[t] == 1.0*E2Ftrl -1.0*E2Fdeg -1.0*PFtrl, E2Fp'[t] == 1.0*E2Fptrl -1.0*E2Fpdeg -1.0*PFptrl, E2W'[t] == 1.0*E2Wtrl -1.0*E2Wdeg -1.0*PWtrl, MF'[t] == 1.0*MFtrn -1.0*MFdeg, MW'[t] == 1.0*MWtrn -1.0*MWdeg, PF'[t] == 1.0*PFtrl -1.0*PFdeg, PFp'[t] == 1.0*PFptrl -1.0*PFpdeg, PW'[t] == 1.0*PWtrl -1.0*PWdeg -1.0*PWtrs, PWL'[t] == 1.0*PWtrs -1.0*PWLdeg
};
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]}]