(* 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 = {
c[t], cer[t], V[t], a[t], n[t], rIR[t], nIR[t]
};

initialValues = {
V[0] == -60.0, a[0] == 0.46, n[0] == 0.01, rIR[0] == 0.282, nIR[0] == 0.008, c[0] == 0.1, cer[0] == 100.0
};

rates = {
reaction\[LetterSpace]0000005, reaction\[LetterSpace]0000006
};

rateEquations = {
reaction\[LetterSpace]0000005 -> cell*fcyt*(jer + jmem), reaction\[LetterSpace]0000006 -> -(cell*fer*jer*vcytver)
};

parameters = {
alpha -> 4.5*^-06, cm -> 5300.0, dact -> 0.35, dinact -> 0.4, dip3 -> 0.5, fcyt -> 0.01, fer -> 0.01, gca -> 1000.0, gir -> 5.0, gk -> 1400.0, gkatp -> 500.0, gkca -> 900.0, ip3 -> 0.0, kd -> 0.3, kpmca -> 0.2, kserca -> 0.4, pleak -> 0.0005, r -> 0.14, sa -> 0.1, sm -> 12.0, sn -> 5.0, taua -> 300000.0, taun -> 16.0, vca -> 50.0, vcytver -> 5.0, vir -> -75.0, vk -> -75.0, vm -> -20.0, vn -> -16.0, cell -> 1.0
};

assignments = {
nIRinf -> (1 + betaIRn/alphaIRn)^(-1), betaIRn -> 0.00035*E^(0.07*(25 + V[t])), alphaIRn -> 0.09/(1 + E^(0.11*(100 + V[t]))), oinf -> (dinact*ip3*c[t])/((dip3 + ip3)*(dact + c[t])*(dinact + c[t])), jleak -> pleak*(-c[t] + cer[t]), jserca -> kserca*c[t], jmem -> -(alpha*ica) - kpmca*c[t], ikca -> gkca*w*(-vk + V[t]), ainf -> (1 + E^((r - c[t])/sa))^(-1), ninf -> (1 + E^((vn - V[t])/sn))^(-1), ica -> gca*minf*(-vca + V[t]), ikatp -> gkatp*a[t]*(-vk + V[t]), w -> c[t]^5/(kd^5 + c[t]^5), ik -> gk*n[t]*(-vk + V[t]), minf -> (1 + E^((vm - V[t])/sm))^(-1), tauIRr -> (alphaIRr + betaIRr)^(-1), rIRinf -> (1 + betaIRr/alphaIRr)^(-1), betaIRr -> 0.15/(1 + E^(-0.05*(120 + V[t]))), alphaIRr -> 30/(1 + E^(0.04*(230 + V[t]))), tauIRn -> (alphaIRn + betaIRn)^(-1), iir -> gir*nIR[t]*rIR[t]*(-vir + V[t]), jer -> jip3 + jleak - jserca, jip3 -> oinf*(-c[t] + cer[t])
};

events = {

};

speciesAnnotations = {
c[t]->"http://identifiers.org/chebi/CHEBI:29108", c[t]->"http://identifiers.org/kegg.compound/C00076", cer[t]->"http://identifiers.org/chebi/CHEBI:29108", cer[t]->"http://identifiers.org/kegg.compound/C00076"
};

reactionAnnotations = {

};

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

(* Time evolution *)
odes = {
c'[t] == 1.0*reaction\[LetterSpace]0000005 , cer'[t] == 1.0*reaction\[LetterSpace]0000006 , V'[t] == -((ica + iir + ik + ikatp + ikca)/cm), a'[t] == (ainf - a[t])/taua, n'[t] == (ninf - n[t])/taun, rIR'[t] == (rIRinf - rIR[t])/tauIRr, nIR'[t] == (nIRinf - nIR[t])/tauIRn
};
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]}]