(* 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], cAMP[t], cer[t], h[t], inh[t], V[t], n[t]
};

initialValues = {
h[0] == 0.0, inh[0] == 1.0, V[0] == -60.0, n[0] == 0.0, c[0] == 0.3, cAMP[0] == 1.0, cer[0] == 260.0
};

rates = {
reaction\[LetterSpace]0000001, reaction\[LetterSpace]0000002, reaction\[LetterSpace]000003
};

rateEquations = {
reaction\[LetterSpace]0000001 -> cell*f*(jertot + jmemtot), reaction\[LetterSpace]0000002 -> -(cell*fer*jertot*sigmav), reaction\[LetterSpace]000003 -> (cell*ETswitch*(cAMPlow - cAMP[t]))/taudir
};

parameters = {
ETswitch -> 0.0, IP3 -> 0.0, alpha -> 4.5*^-06, cAMPlow -> 0.2, cm -> 5300.0, dact -> 0.35, dinh -> 0.4, dip3 -> 0.5, f -> 0.01, fer -> 0.01, gca -> 2000.0, girk -> 1000.0, gk -> 3500.0, kc -> 0.15, ki -> 0.5, kserca -> 0.4, lambda -> 1.25, perl -> 0.0005, sh -> 70.0, sigmav -> 10.0, sm -> 12.0, sn -> 5.0, taudir -> 20000.0, tauh -> 20.0, taun -> 20.0, vca -> 25.0, vh -> -20.0, vk -> -75.0, vm -> -20.0, vn -> -16.0, cell -> 1.0
};

assignments = {
igirk -> girk*h[t]*(-vk + V[t]), ica -> gca*minf*(-vca + V[t]), jerp -> kserca*c[t], jerip3 -> o*(-c[t] + cer[t]), hinfer -> (1 + c[t]/dinh)^(-1), ainf -> (1 + dact/c[t])^(-1), ik -> gk*n[t]*(-vk + V[t]), perl\[LetterSpace]inf -> (c[t]^4*cAMP[t]*inh[t])/(ki^4 + c[t]^4), jertot -> jerip3 + jerleak - kserca*c[t], jerleak -> perl*(-c[t] + cer[t]), jmemtot -> -(kc*c[t]) - (alpha*gca*(-vca + V[t]))/(1 + E^((vm - V[t])/sm)), o -> ainf^3*binf^3*hinfer^3, binf -> IP3/(dip3 + IP3), hinf -> (1 + E^((vh - V[t])/sh))^(-1), ninf -> (1 + E^((vn - V[t])/sn))^(-1), minf -> (1 + E^((vm - V[t])/sm))^(-1)
};

events = {

};

speciesAnnotations = {
c[t]->"http://identifiers.org/obo.chebi/CHEBI:29108", c[t]->"http://identifiers.org/kegg.compound/C00076", cAMP[t]->"http://identifiers.org/obo.chebi/CHEBI:17489", cAMP[t]->"http://identifiers.org/kegg.compound/C00575", cer[t]->"http://identifiers.org/obo.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]0000001 , cAMP'[t] == 1.0*reaction\[LetterSpace]000003 , cer'[t] == 1.0*reaction\[LetterSpace]0000002 , h'[t] == (hinf - h[t])/tauh, inh'[t] == (ETswitch*(0.2 - inh[t]))/taudir, V'[t] == (-ica - igirk - ik)/cm, n'[t] == (lambda*(ninf - n[t]))/taun
};
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]}]