(* 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 = {
AHKph[t], ARF2[t], ARFIAA[t], ARRAm[t], ARRAp[t], ARRAph[t], ARRBph[t], ARm[t], ARp[t], Aux[t], AuxTIAA[t], AuxTIR1[t], CRm[t], CRp[t], Ck[t], IAAm[t], IAAp[t], IAAs[t], PINm[t], PINp[t]
};

initialValues = {
AHKph[0] == 1.0, ARF2[0] == 0.0, ARFIAA[0] == 0.0, ARRAm[0] == 0.0, ARRAp[0] == 0.0, ARRAph[0] == 0.0, ARRBph[0] == 0.0, ARm[0] == 0.0, ARp[0] == 0.0, Aux[0] == 1.0, AuxTIAA[0] == 0.0, AuxTIR1[0] == 0.0, CRm[0] == 0.0, CRp[0] == 0.0, Ck[0] == 1.0, IAAm[0] == 0.0, IAAp[0] == 0.0, IAAs[0] == 0.0, PINm[0] == 0.0, PINp[0] == 0.0
};

rates = {
r1, r10, r11, r12, r13, r14, r15, r16, r17, r18, r19, r2, r20, r3, r4, r5, r6, r7, r8, r9
};

rateEquations = {
r1 -> (F2 + F1*lambda1 + F3*lambda3)*phiIAAp - IAAm[t], r10 -> (deltaPINp*PINm[t] - PINp[t])/eps, r11 -> (F5b + F5a*lambda1)*phiARp - ARm[t], r12 -> (deltaARp*ARm[t] - ARp[t])/eps, r13 -> F4*phiCRp - CRm[t], r14 -> (deltaCRp*CRm[t] - CRp[t])/eps, r15 -> (CkAHKph*rd - ra*AHKph[t]*Ck[t])/eps, r16 -> muCk*(alphaCk - Ck[t]) - (etaCkPh*(-(CkAHKph*rd) + ra*AHKph[t]*Ck[t]))/eps, r17 -> (ARRBp*CkAHKph*ua - CkAHK*ud*ARRBph[t])/eps, r18 -> (CkAHKph*sa*ARRAp[t] - CkAHK*sd*ARRAph[t])/eps, r19 -> F6*phiARRAp - ARRAm[t], r2 -> etaARFIAA*(pd*ARFIAA[t] - ARF*pa*IAAp[t]) + (ld*AuxTIAA[t] + deltaIAAp*IAAm[t] - la*AuxTIR1[t]*IAAp[t])/eps, r20 -> (deltaARRAp*ARRAm[t] - ARRAp[t] + etaAHKph*(-(CkAHKph*sa*ARRAp[t]) + CkAHK*sd*ARRAph[t]))/eps, r3 -> (ka*TIR1*Aux[t] + (1 + ld)*AuxTIAA[t] - kd*AuxTIR1[t] - la*AuxTIR1[t]*IAAp[t])/eps, r4 -> (-((1 + ld)*AuxTIAA[t]) + la*AuxTIR1[t]*IAAp[t])/eps, r5 -> (AuxTIAA[t] - muIAAs*IAAs[t])/eps, r6 -> -(pd*ARFIAA[t]) + ARF*pa*IAAp[t], r7 -> ARF^2*qa - qd*ARF2[t], r8 -> muAux*(alphaAux - Aux[t]) - (etaAuxTIR1*(ka*TIR1*Aux[t] - kd*AuxTIR1[t]))/eps, r9 -> (F5b + F5a*lambda1)*phiPINp - PINm[t]
};

parameters = {
alphaAHK -> 1.0, alphaARF -> 1.0, alphaARRB -> 2.0, alphaAux -> 1.0, alphaCk -> 1.0, alphaPH -> 1.0, alphaTIR1 -> 1.0, deltaARRAp -> 1.0, deltaARp -> 1.0, deltaCRp -> 1.0, deltaIAAp -> 1.0, deltaPINp -> 1.0, eps -> 0.01, etaAHKph -> 1.0, etaARFIAA -> 1.0, etaAuxTIR1 -> 10.0, etaCkPh -> 1.0, ka -> 100.0, kd -> 1.0, la -> 0.5, lambda1 -> 0.1, lambda3 -> 0.02, ld -> 0.1, muAux -> 0.1, muCk -> 0.1, muIAAs -> 1.0, pa -> 10.0, pd -> 10.0, phiARRAp -> 100.0, phiARp -> 2.0, phiCRp -> 2.0, phiIAAp -> 100.0, phiPINp -> 100.0, psiARF -> 0.1, psiARFIAA -> 0.1, qa -> 1.0, qd -> 1.0, ra -> 1.0, rd -> 1.0, sa -> 1.0, sd -> 1.0, thARFIAA -> 0.1, thARRAph -> 0.1, thARRBph -> 0.1, thetaARF -> 0.1, thetaARF2 -> 0.01, thetaARp -> 0.1, thetaIAAp -> 0.1, ua -> 1.0, ud -> 1.0, cell -> 1.0
};

assignments = {
CkAHK -> alphaAHK - etaAHKph*(CkAHKph + AHKph[t]), CkAHKph -> alphaPH - AHKph[t] - ARRAph[t] - ARRBph[t], ARRBp -> alphaARRB - etaAHKph*ARRBph[t], ARF -> alphaARF - 2*ARF2[t] - ARFIAA[t], TIR1 -> alphaTIR1 - AuxTIAA[t] - AuxTIR1[t], F6 -> ARp[t]/(thetaARp*(1 + ARp[t]/thetaARp)), F5b -> (ARF^2/psiARF + ARF2[t]/thetaARF2)/(1 + ARF^2/psiARF + ARF/thetaARF + ARF2[t]/thetaARF2 + ARFIAA[t]/thARFIAA + (ARF*IAAp[t])/psiARFIAA), F5a -> ARF/(thetaARF*(1 + ARF^2/psiARF + ARF/thetaARF + ARF2[t]/thetaARF2 + ARFIAA[t]/thARFIAA + (ARF*IAAp[t])/psiARFIAA)), F4 -> ARRBph[t]/(thARRBph*(1 + ARRAph[t]/thARRAph + ARRBph[t]/thARRBph)), F3 -> ARRBph[t]/(thARRBph*(1 + ARF^2/psiARF + ARF/thetaARF + ARF2[t]/thetaARF2 + ARFIAA[t]/thARFIAA + ARRBph[t]/thARRBph + (ARF*IAAp[t])/psiARFIAA)), F2 -> (ARF^2/psiARF + ARF2[t]/thetaARF2)/(1 + ARF^2/psiARF + ARF/thetaARF + ARF2[t]/thetaARF2 + ARFIAA[t]/thARFIAA + ARRBph[t]/thARRBph + (ARF*IAAp[t])/psiARFIAA), F1 -> ARF/(thetaARF*(1 + ARF^2/psiARF + ARF/thetaARF + ARF2[t]/thetaARF2 + ARFIAA[t]/thARFIAA + ARRBph[t]/thARRBph + (ARF*IAAp[t])/psiARFIAA))
};

events = {

};

speciesAnnotations = {

};

reactionAnnotations = {

};

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

(* Time evolution *)
odes = {
AHKph'[t] == 1.0*r15 , ARF2'[t] == 1.0*r7 , ARFIAA'[t] == 1.0*r6 , ARRAm'[t] == 1.0*r19 , ARRAp'[t] == 1.0*r20 , ARRAph'[t] == 1.0*r18 , ARRBph'[t] == 1.0*r17 , ARm'[t] == 1.0*r11 , ARp'[t] == 1.0*r12 , Aux'[t] == 1.0*r8 , AuxTIAA'[t] == 1.0*r4 , AuxTIR1'[t] == 1.0*r3 , CRm'[t] == 1.0*r13 , CRp'[t] == 1.0*r14 , Ck'[t] == 1.0*r16 , IAAm'[t] == 1.0*r1 , IAAp'[t] == 1.0*r2 , IAAs'[t] == 1.0*r5 , PINm'[t] == 1.0*r9 , PINp'[t] == 1.0*r10 
};
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]}]