(* 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 = {
C1[t], C2[t], C3[t], IC2[t], IC3[t], IF[t], IM1[t], IM2[t], O[t]
};

initialValues = {
C1[0] == 0.0, C2[0] == 0.0, C3[0] == 1.0, IC2[0] == 0.0, IC3[0] == 0.0, IF[0] == 0.0, IM1[0] == 0.0, IM2[0] == 0.0, O[0] == 0.0
};

rates = {
reaction\[LetterSpace]0000001, reaction\[LetterSpace]0000002, reaction\[LetterSpace]0000003, reaction\[LetterSpace]0000004, reaction\[LetterSpace]0000005, reaction\[LetterSpace]0000006, reaction\[LetterSpace]0000007, reaction\[LetterSpace]0000008, reaction\[LetterSpace]0000009, reaction\[LetterSpace]0000010, reaction\[LetterSpace]0000011
};

rateEquations = {
reaction\[LetterSpace]0000001 -> cell*(-(b11*IC2[t]) + a11*IC3[t]), reaction\[LetterSpace]0000002 -> cell*(a12*IC2[t] - b12*IF[t]), reaction\[LetterSpace]0000003 -> cell*(a4*IF[t] - b4*IM1[t]), reaction\[LetterSpace]0000004 -> cell*(a5*IM1[t] - b5*IM2[t]), reaction\[LetterSpace]0000005 -> cell*(b3*C3[t] - a3*IC3[t]), reaction\[LetterSpace]0000006 -> cell*(b11*C2[t] - a11*C3[t]), reaction\[LetterSpace]0000007 -> cell*(b3*C2[t] - a3*IC2[t]), reaction\[LetterSpace]0000008 -> cell*(b12*C1[t] - a12*C2[t]), reaction\[LetterSpace]0000009 -> cell*(b3*C1[t] - a3*IF[t]), reaction\[LetterSpace]0000010 -> cell*(b2*IF[t] - a2*O[t]), reaction\[LetterSpace]0000011 -> cell*(-(a13*C1[t]) + b13*O[t])
};

parameters = {
Fara -> 96485.0, Gna -> 23.5, Rk -> 8314.0, Temp -> 310.0, nai -> 15.0, nao -> 140.0, scale -> 1000.0, vhold -> -80.0, vtest -> -30.0, cell -> 1.0
};

assignments = {
a11 -> 3.802/(0.1027/E^(v/17) + 0.2/E^(v/150)), a12 -> 3.802/(0.1027/E^(v/15) + 0.23/E^(v/150)), a2 -> 9.178*E^(0.03369272237196765*v), b13 -> 0.22/E^(0.04926108374384236*(-10 + v)), b3 -> 0.0084 + v/50000, a4 -> 0.09178000000000001*E^(0.03369272237196765*v), b4 -> 3.7933*^-7/E^(0.12987012987012986*v), b5 -> 7.5866*^-9/E^(0.12987012987012986*v), a5 -> 0.00009661052631578948*E^(0.03369272237196765*v), Ena -> (Rk*Temp*Log[nao/nai])/Fara, b11 -> 0.1917/E^(0.04926108374384236*v), b2 -> (0.000060166489970363645*E^(0.04926108374384236*(-10 + v) - 0.0961774074981622*v))/((0.1027/E^(v/12) + 0.25/E^(v/150))*(0.0084 + v/50000)), b12 -> 0.2/E^(0.04926108374384236*(-5 + v)), Ina -> (Gna*(v - (Rk*Temp*Log[nao/nai])/Fara)*O[t])/(scale*(C1[t] + C2[t] + C3[t] + IC2[t] + IC3[t] + IF[t] + IM1[t] + IM2[t] + O[t])), a3 -> 3.7933*^-7/E^(0.12987012987012986*v), a13 -> 3.802/(0.1027/E^(v/12) + 0.25/E^(v/150)), v -> Piecewise[{{vtest, 5 < t && t <= 300}}, vhold]
};

events = {

};

speciesAnnotations = {

};

reactionAnnotations = {

};

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

(* Time evolution *)
odes = {
C1'[t] == 1.0*reaction\[LetterSpace]0000011 -1.0*reaction\[LetterSpace]0000008 -1.0*reaction\[LetterSpace]0000009, C2'[t] == 1.0*reaction\[LetterSpace]0000008 -1.0*reaction\[LetterSpace]0000006 -1.0*reaction\[LetterSpace]0000007, C3'[t] == 1.0*reaction\[LetterSpace]0000006 -1.0*reaction\[LetterSpace]0000005, IC2'[t] == 1.0*reaction\[LetterSpace]0000001 +1.0*reaction\[LetterSpace]0000007 -1.0*reaction\[LetterSpace]0000002, IC3'[t] == 1.0*reaction\[LetterSpace]0000005 -1.0*reaction\[LetterSpace]0000001, IF'[t] == 1.0*reaction\[LetterSpace]0000002 +1.0*reaction\[LetterSpace]0000009 -1.0*reaction\[LetterSpace]0000003 -1.0*reaction\[LetterSpace]0000010, IM1'[t] == 1.0*reaction\[LetterSpace]0000003 -1.0*reaction\[LetterSpace]0000004, IM2'[t] == 1.0*reaction\[LetterSpace]0000004 , O'[t] == 1.0*reaction\[LetterSpace]0000010 -1.0*reaction\[LetterSpace]0000011
};
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]}]