(* 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], h1[t], h2[t]
};

initialValues = {
c1[0] == 0.3, c2[0] == 0.1, h1[0] == 0.8, h2[0] == 0.1
};

rates = {
Inactivated\[LetterSpace]to\[LetterSpace]S\[LetterSpace]Cell1, Inactivated\[LetterSpace]to\[LetterSpace]S\[LetterSpace]Cell2, Jleak\[LetterSpace]Cell1, Jleak\[LetterSpace]Cell2, Jpump\[LetterSpace]Cell1, Jpump\[LetterSpace]Cell2, Jreceptor\[LetterSpace]Cell1, Jreceptor\[LetterSpace]Cell2, Open\[LetterSpace]to\[LetterSpace]Inactivated\[LetterSpace]Cell1, Open\[LetterSpace]to\[LetterSpace]Inactivated\[LetterSpace]Cell2, diffusion
};

rateEquations = {
Inactivated\[LetterSpace]to\[LetterSpace]S\[LetterSpace]Cell1 -> compartment*Phi3\[LetterSpace]c1*(1 - h1[t]), Inactivated\[LetterSpace]to\[LetterSpace]S\[LetterSpace]Cell2 -> compartment*Phi3\[LetterSpace]c2*(1 - h2[t]), Jleak\[LetterSpace]Cell1 -> compartment*Jleak\[LetterSpace]Cell1\[LetterSpace]Jleak, Jleak\[LetterSpace]Cell2 -> compartment*Jleak\[LetterSpace]Cell2\[LetterSpace]Jleak, Jpump\[LetterSpace]Cell1 -> (compartment*Jpump\[LetterSpace]Cell1\[LetterSpace]Vp*c1[t]^2)/(Jpump\[LetterSpace]Cell1\[LetterSpace]Kp^2 + c1[t]^2), Jpump\[LetterSpace]Cell2 -> (compartment*Jpump\[LetterSpace]Cell2\[LetterSpace]Vp*c2[t]^2)/(Jpump\[LetterSpace]Cell2\[LetterSpace]Kp^2 + c2[t]^2), Jreceptor\[LetterSpace]Cell1 -> (compartment*Jreceptor\[LetterSpace]Cell1\[LetterSpace]kf*Jreceptor\[LetterSpace]Cell1\[LetterSpace]p^4*Phi1\[LetterSpace]c1^4*h1[t]^4)/(Jreceptor\[LetterSpace]Cell1\[LetterSpace]p*Phi1\[LetterSpace]c1 + Phi\[LetterSpace]minus1\[LetterSpace]c1)^4, Jreceptor\[LetterSpace]Cell2 -> (compartment*Jreceptor\[LetterSpace]Cell2\[LetterSpace]kf*Jreceptor\[LetterSpace]Cell2\[LetterSpace]p^4*Phi1\[LetterSpace]c2^4*h2[t]^4)/(Jreceptor\[LetterSpace]Cell2\[LetterSpace]p*Phi1\[LetterSpace]c2 + Phi\[LetterSpace]minus1\[LetterSpace]c2)^4, Open\[LetterSpace]to\[LetterSpace]Inactivated\[LetterSpace]Cell1 -> (compartment*Open\[LetterSpace]to\[LetterSpace]Inactivated\[LetterSpace]Cell1\[LetterSpace]p*Phi1\[LetterSpace]c1*Phi2\[LetterSpace]c1*h1[t])/(Open\[LetterSpace]to\[LetterSpace]Inactivated\[LetterSpace]Cell1\[LetterSpace]p*Phi1\[LetterSpace]c1 + Phi\[LetterSpace]minus1\[LetterSpace]c1), Open\[LetterSpace]to\[LetterSpace]Inactivated\[LetterSpace]Cell2 -> (compartment*Open\[LetterSpace]to\[LetterSpace]Inactivated\[LetterSpace]Cell2\[LetterSpace]p*Phi1\[LetterSpace]c2*Phi2\[LetterSpace]c2*h2[t])/(Open\[LetterSpace]to\[LetterSpace]Inactivated\[LetterSpace]Cell2\[LetterSpace]p*Phi1\[LetterSpace]c2 + Phi\[LetterSpace]minus1\[LetterSpace]c2), diffusion -> compartment*diffusion\[LetterSpace]D*(-c1[t] + c2[t])
};

parameters = {
R1 -> 6.0, R3 -> 50.0, R5 -> 1.6, k1 -> 44.0, k2 -> 26.5, k3 -> 1.6, r2 -> 100.0, r4 -> 20.0, Jreceptor\[LetterSpace]Cell1\[LetterSpace]kf -> 28.0, Jreceptor\[LetterSpace]Cell1\[LetterSpace]p -> 0.2778, Jpump\[LetterSpace]Cell1\[LetterSpace]Vp -> 1.2, Jpump\[LetterSpace]Cell1\[LetterSpace]Kp -> 0.18, Jleak\[LetterSpace]Cell1\[LetterSpace]Jleak -> 0.2, Open\[LetterSpace]to\[LetterSpace]Inactivated\[LetterSpace]Cell1\[LetterSpace]p -> 0.2778, Jreceptor\[LetterSpace]Cell2\[LetterSpace]kf -> 28.0, Jreceptor\[LetterSpace]Cell2\[LetterSpace]p -> 0.2778, Jpump\[LetterSpace]Cell2\[LetterSpace]Vp -> 1.2, Jpump\[LetterSpace]Cell2\[LetterSpace]Kp -> 0.18, Jleak\[LetterSpace]Cell2\[LetterSpace]Jleak -> 0.2, Open\[LetterSpace]to\[LetterSpace]Inactivated\[LetterSpace]Cell2\[LetterSpace]p -> 0.2778, diffusion\[LetterSpace]D -> 0.01, compartment -> 1.0
};

assignments = {
Phi3\[LetterSpace]c2 -> k3/(R5 + c2[t]), Phi2\[LetterSpace]c2 -> (k2 + r4*c2[t])/(R3 + c2[t]), Phi\[LetterSpace]minus1\[LetterSpace]c2 -> k1/(R3 + c2[t]), Phi1\[LetterSpace]c2 -> (r2*c2[t])/(R1 + c2[t]), Phi3\[LetterSpace]c1 -> k3/(R5 + c1[t]), Phi2\[LetterSpace]c1 -> (k2 + r4*c1[t])/(R3 + c1[t]), Phi\[LetterSpace]minus1\[LetterSpace]c1 -> k1/(R3 + c1[t]), Phi1\[LetterSpace]c1 -> (r2*c1[t])/(R1 + c1[t])
};

events = {

};

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

reactionAnnotations = {
Jpump\[LetterSpace]Cell1->"http://identifiers.org/go/GO:0005388", Jpump\[LetterSpace]Cell2->"http://identifiers.org/go/GO:0005388", Jreceptor\[LetterSpace]Cell1->"http://identifiers.org/go/GO:0005220", Jreceptor\[LetterSpace]Cell1->"http://identifiers.org/go/GO:0007204", Jreceptor\[LetterSpace]Cell2->"http://identifiers.org/go/GO:0005220", Jreceptor\[LetterSpace]Cell2->"http://identifiers.org/go/GO:0007204"
};

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

(* Time evolution *)
odes = {
c1'[t] == 1.0*Jreceptor\[LetterSpace]Cell1 +1.0*Jleak\[LetterSpace]Cell1 +1.0*diffusion -1.0*Jpump\[LetterSpace]Cell1, c2'[t] == 1.0*Jreceptor\[LetterSpace]Cell2 +1.0*Jleak\[LetterSpace]Cell2 -1.0*Jpump\[LetterSpace]Cell2 -1.0*diffusion, h1'[t] == 1.0*Inactivated\[LetterSpace]to\[LetterSpace]S\[LetterSpace]Cell1 -1.0*Open\[LetterSpace]to\[LetterSpace]Inactivated\[LetterSpace]Cell1, h2'[t] == 1.0*Inactivated\[LetterSpace]to\[LetterSpace]S\[LetterSpace]Cell2 -1.0*Open\[LetterSpace]to\[LetterSpace]Inactivated\[LetterSpace]Cell2
};
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]}]