(* 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 = {
A[t], D\[LetterSpace]IA[t], D\[LetterSpace]II[t], I[t], R[t]
};

initialValues = {
A[0] == 10.0, D\[LetterSpace]IA[0] == 10.0, D\[LetterSpace]II[0] == 10.0, I[0] == 10.0, R[0] == 1.0
};

rates = {
reac\[LetterSpace]DIA, reac\[LetterSpace]DIAxA, reac\[LetterSpace]DII, reac\[LetterSpace]DIIxI, reac\[LetterSpace]degrA, reac\[LetterSpace]degrDIA, reac\[LetterSpace]degrDII, reac\[LetterSpace]degrI, reac\[LetterSpace]degrR, reac\[LetterSpace]prodA, reac\[LetterSpace]prodI, reac\[LetterSpace]prodR
};

rateEquations = {
reac\[LetterSpace]DIA -> -(kprime\[LetterSpace]IA*D\[LetterSpace]IA[t]) + k\[LetterSpace]IA*A[t]*I[t], reac\[LetterSpace]DIAxA -> (aux*d\[LetterSpace]I*gamma\[LetterSpace]I*K\[LetterSpace]aux*D\[LetterSpace]IA[t])/(1 + aux*K\[LetterSpace]aux), reac\[LetterSpace]DII -> -(kprime\[LetterSpace]II*D\[LetterSpace]II[t]) + k\[LetterSpace]II*I[t]^2, reac\[LetterSpace]DIIxI -> (aux*d\[LetterSpace]I*gamma\[LetterSpace]I*K\[LetterSpace]aux*D\[LetterSpace]II[t])/(1 + aux*K\[LetterSpace]aux), reac\[LetterSpace]degrA -> d\[LetterSpace]A*A[t], reac\[LetterSpace]degrDIA -> d\[LetterSpace]IA*D\[LetterSpace]IA[t], reac\[LetterSpace]degrDII -> d\[LetterSpace]II*D\[LetterSpace]II[t], reac\[LetterSpace]degrI -> (aux*d\[LetterSpace]I*gamma\[LetterSpace]I*K\[LetterSpace]aux*I[t])/(1 + aux*K\[LetterSpace]aux), reac\[LetterSpace]degrR -> d\[LetterSpace]r*R[t], reac\[LetterSpace]prodA -> pi\[LetterSpace]A, reac\[LetterSpace]prodI -> pi\[LetterSpace]I*R[t], reac\[LetterSpace]prodR -> (1 + (f\[LetterSpace]c*A[t]*(1 + (f\[LetterSpace]A*w\[LetterSpace]A*A[t])/B\[LetterSpace]d))/B\[LetterSpace]d)/(1 + k\[LetterSpace]Am + (A[t]*(1 + (w\[LetterSpace]A*A[t])/B\[LetterSpace]d))/B\[LetterSpace]d + (w\[LetterSpace]D*D\[LetterSpace]IA[t])/B\[LetterSpace]d + (w\[LetterSpace]I*A[t]*I[t])/(B\[LetterSpace]d*K\[LetterSpace]IA))
};

parameters = {
B\[LetterSpace]d -> 100.0, K\[LetterSpace]IA -> 10.0, K\[LetterSpace]II -> 10.0, K\[LetterSpace]aux -> 1.0, d\[LetterSpace]A -> 0.003, d\[LetterSpace]I -> 0.05, d\[LetterSpace]IA -> 0.003, d\[LetterSpace]II -> 0.003, d\[LetterSpace]r -> 0.007, f\[LetterSpace]A -> 10.0, f\[LetterSpace]c -> 10.0, gamma\[LetterSpace]I -> 10.0, k\[LetterSpace]Am -> 10.0, k\[LetterSpace]IA -> 1.0, k\[LetterSpace]II -> 1.0, pi\[LetterSpace]A -> 1.0, pi\[LetterSpace]I -> 1.0, w\[LetterSpace]A -> 10.0, w\[LetterSpace]D -> 10.0, w\[LetterSpace]I -> 10.0, compartment\[LetterSpace]0000001 -> 1.0
};

assignments = {
aux\[LetterSpace]basal -> 1/((-1 + gamma\[LetterSpace]I)*K\[LetterSpace]aux), kprime\[LetterSpace]II -> k\[LetterSpace]II*K\[LetterSpace]II, kprime\[LetterSpace]IA -> k\[LetterSpace]IA*K\[LetterSpace]IA, aux -> 1.11 + 1.11*Sin[0.007854*(-200 + t)]
};

events = {

};

speciesAnnotations = {

};

reactionAnnotations = {

};

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

(* Time evolution *)
odes = {
A'[t] == 1.0*reac\[LetterSpace]DIAxA +1.0*reac\[LetterSpace]prodA -1.0*reac\[LetterSpace]DIA -1.0*reac\[LetterSpace]degrA, D\[LetterSpace]IA'[t] == 1.0*reac\[LetterSpace]DIA -1.0*reac\[LetterSpace]degrDIA -1.0*reac\[LetterSpace]DIAxA, D\[LetterSpace]II'[t] == 1.0*reac\[LetterSpace]DII -1.0*reac\[LetterSpace]degrDII -1.0*reac\[LetterSpace]DIIxI, I'[t] == 1.0*reac\[LetterSpace]DIIxI +1.0*reac\[LetterSpace]prodI -1.0*reac\[LetterSpace]DIA -1.0*reac\[LetterSpace]DII -1.0*reac\[LetterSpace]DII -1.0*reac\[LetterSpace]degrI, R'[t] == 1.0*reac\[LetterSpace]prodR -1.0*reac\[LetterSpace]degrR
};
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]}]