(* 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 = {
D\[LetterSpace]1[t], E\[LetterSpace]1[t], RS\[LetterSpace]1[t], R\[LetterSpace]1[t], X\[LetterSpace]1[t]
};

initialValues = {
D\[LetterSpace]1[0] == 0.1, E\[LetterSpace]1[0] == 0.6, RS\[LetterSpace]1[0] == 1.0, R\[LetterSpace]1[0] == 0.5, X\[LetterSpace]1[0] == 0.7
};

rates = {
\[LetterSpace]1, \[LetterSpace]7, cycECDK2degradation\[LetterSpace]1, cyclebreak\[LetterSpace]1, cycleprogression\[LetterSpace]1, cyclinCDK4degradation\[LetterSpace]1, cyclin\[LetterSpace]1, pRBE2Fcomplexassociation\[LetterSpace]1, pRBE2FcomplexdeassociationviacycDCDK4\[LetterSpace]1, pRBpdephosphorylation\[LetterSpace]1
};

rateEquations = {
\[LetterSpace]1 -> (aD\[LetterSpace]1*GF\[LetterSpace]1*k\[LetterSpace]1)/(1 + GF\[LetterSpace]1*k\[LetterSpace]1), \[LetterSpace]7 -> (pE\[LetterSpace]1*E\[LetterSpace]1[t]*RS\[LetterSpace]1[t])/(qE\[LetterSpace]1 + E\[LetterSpace]1[t] + RS\[LetterSpace]1[t]), cycECDK2degradation\[LetterSpace]1 -> dE\[LetterSpace]1*E\[LetterSpace]1[t]*X\[LetterSpace]1[t], cyclebreak\[LetterSpace]1 -> dX\[LetterSpace]1*X\[LetterSpace]1[t], cycleprogression\[LetterSpace]1 -> E2F\[LetterSpace]1*f\[LetterSpace]1 + aX\[LetterSpace]1*E\[LetterSpace]1[t] + g\[LetterSpace]1*E\[LetterSpace]1[t]*X\[LetterSpace]1[t]^2, cyclinCDK4degradation\[LetterSpace]1 -> dD\[LetterSpace]1*D\[LetterSpace]1[t]*E\[LetterSpace]1[t], cyclin\[LetterSpace]1 -> aE\[LetterSpace]1*(1 + af\[LetterSpace]1*E2F\[LetterSpace]1), pRBE2Fcomplexassociation\[LetterSpace]1 -> E2F\[LetterSpace]1*pS\[LetterSpace]1*R\[LetterSpace]1[t], pRBE2FcomplexdeassociationviacycDCDK4\[LetterSpace]1 -> (pD\[LetterSpace]1*D\[LetterSpace]1[t]*RS\[LetterSpace]1[t])/(qD\[LetterSpace]1 + D\[LetterSpace]1[t] + RS\[LetterSpace]1[t]), pRBpdephosphorylation\[LetterSpace]1 -> (pX\[LetterSpace]1*RP\[LetterSpace]1*X\[LetterSpace]1[t])/(qX\[LetterSpace]1 + RP\[LetterSpace]1 + X\[LetterSpace]1[t])
};

parameters = {
GF\[LetterSpace]1 -> 6.3, RT\[LetterSpace]1 -> 2.5, aD\[LetterSpace]1 -> 0.4, aE\[LetterSpace]1 -> 0.16, aX\[LetterSpace]1 -> 0.08, af\[LetterSpace]1 -> 0.9, dD\[LetterSpace]1 -> 0.4, dE\[LetterSpace]1 -> 0.2, dX\[LetterSpace]1 -> 1.04, f\[LetterSpace]1 -> 0.2, g\[LetterSpace]1 -> 0.528, k\[LetterSpace]1 -> 0.05, pD\[LetterSpace]1 -> 0.48, pE\[LetterSpace]1 -> 0.096, pS\[LetterSpace]1 -> 0.6, pX\[LetterSpace]1 -> 0.48, qD\[LetterSpace]1 -> 0.6, qE\[LetterSpace]1 -> 0.6, qX\[LetterSpace]1 -> 0.8, theta\[LetterSpace]1 -> 1.5, cell\[LetterSpace]1 -> 1.0
};

assignments = {
RP\[LetterSpace]1 -> RT\[LetterSpace]1 - RS\[LetterSpace]1[t] - R\[LetterSpace]1[t], E2F\[LetterSpace]1 -> theta\[LetterSpace]1 - RS\[LetterSpace]1[t], unpho\[LetterSpace]RB -> RS\[LetterSpace]1[t] + R\[LetterSpace]1[t]
};

events = {

};

speciesAnnotations = {

};

reactionAnnotations = {

};

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

(* Time evolution *)
odes = {
D\[LetterSpace]1'[t] == 1.0*\[LetterSpace]1 -1.0*cyclinCDK4degradation\[LetterSpace]1, E\[LetterSpace]1'[t] == 1.0*cyclin\[LetterSpace]1 -1.0*cycECDK2degradation\[LetterSpace]1, RS\[LetterSpace]1'[t] == 1.0*pRBE2Fcomplexassociation\[LetterSpace]1 -1.0*pRBE2FcomplexdeassociationviacycDCDK4\[LetterSpace]1 -1.0*\[LetterSpace]7, R\[LetterSpace]1'[t] == 1.0*pRBpdephosphorylation\[LetterSpace]1 -1.0*pRBE2Fcomplexassociation\[LetterSpace]1, X\[LetterSpace]1'[t] == 1.0*cycleprogression\[LetterSpace]1 -1.0*cyclebreak\[LetterSpace]1
};
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]}]