(* 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 = {
Fvar[t], S[t], a0[t], a1[t], a10[t], a11[t], a12[t], a13[t], a14[t], a15[t], a16[t], a17[t], a18[t], a19[t], a2[t], a20[t], a3[t], a4[t], a5[t], a6[t], a7[t], a8[t], a9[t], c1[t], c10[t], c11[t], c12[t], c13[t], c14[t], c15[t], c16[t], c17[t], c18[t], c19[t], c2[t], c20[t], c3[t], c4[t], c5[t], c6[t], c7[t], c8[t], c9[t]
};

initialValues = {
Fvar[0] == 0.0, S[0] == 1.0, a0[0] == 0.0, a1[0] == 0.0, a10[0] == 0.0, a11[0] == 0.0, a12[0] == 0.0, a13[0] == 0.0, a14[0] == 0.0, a15[0] == 0.0, a16[0] == 0.0, a17[0] == 0.0, a18[0] == 0.0, a19[0] == 0.0, a2[0] == 0.0, a20[0] == 0.0, a3[0] == 0.0, a4[0] == 0.0, a5[0] == 0.0, a6[0] == 0.0, a7[0] == 0.0, a8[0] == 0.0, a9[0] == 0.0, c1[0] == 0.0, c10[0] == 0.0, c11[0] == 0.0, c12[0] == 0.0, c13[0] == 0.0, c14[0] == 0.0, c15[0] == 0.0, c16[0] == 0.0, c17[0] == 0.0, c18[0] == 0.0, c19[0] == 0.0, c2[0] == 0.0, c20[0] == 0.0, c3[0] == 0.0, c4[0] == 0.0, c5[0] == 0.0, c6[0] == 0.0, c7[0] == 0.0, c8[0] == 0.0, c9[0] == 0.0
};

rates = {
v\[LetterSpace]1, v\[LetterSpace]10, v\[LetterSpace]11, v\[LetterSpace]12, v\[LetterSpace]13, v\[LetterSpace]14, v\[LetterSpace]15, v\[LetterSpace]16, v\[LetterSpace]17, v\[LetterSpace]18, v\[LetterSpace]19, v\[LetterSpace]2, v\[LetterSpace]20, v\[LetterSpace]21, v\[LetterSpace]22, v\[LetterSpace]23, v\[LetterSpace]24, v\[LetterSpace]25, v\[LetterSpace]26, v\[LetterSpace]27, v\[LetterSpace]28, v\[LetterSpace]29, v\[LetterSpace]3, v\[LetterSpace]30, v\[LetterSpace]31, v\[LetterSpace]32, v\[LetterSpace]33, v\[LetterSpace]34, v\[LetterSpace]35, v\[LetterSpace]36, v\[LetterSpace]37, v\[LetterSpace]38, v\[LetterSpace]39, v\[LetterSpace]4, v\[LetterSpace]40, v\[LetterSpace]41, v\[LetterSpace]42, v\[LetterSpace]5, v\[LetterSpace]6, v\[LetterSpace]7, v\[LetterSpace]8, v\[LetterSpace]9
};

rateEquations = {
v\[LetterSpace]1 -> (2*Dvar*kVol*(-a0[t] + S[t]))/(lambda + Ls), v\[LetterSpace]10 -> kVol*Pdiff*(Aprota*a4[t] - Aprotc*c5[t]), v\[LetterSpace]11 -> kVol*(Pdiff*(-(Aprota*a5[t]) + Aprotc*c5[t]) + PPin*(-(Ba*a5[t]) + Bc*c5[t])), v\[LetterSpace]12 -> kVol*Pdiff*(Aprota*a5[t] - Aprotc*c6[t]), v\[LetterSpace]13 -> kVol*(Pdiff*(-(Aprota*a6[t]) + Aprotc*c6[t]) + PPin*(-(Ba*a6[t]) + Bc*c6[t])), v\[LetterSpace]14 -> kVol*Pdiff*(Aprota*a6[t] - Aprotc*c7[t]), v\[LetterSpace]15 -> kVol*(Pdiff*(-(Aprota*a7[t]) + Aprotc*c7[t]) + PPin*(-(Ba*a7[t]) + Bc*c7[t])), v\[LetterSpace]16 -> kVol*Pdiff*(Aprota*a7[t] - Aprotc*c8[t]), v\[LetterSpace]17 -> kVol*(Pdiff*(-(Aprota*a8[t]) + Aprotc*c8[t]) + PPin*(-(Ba*a8[t]) + Bc*c8[t])), v\[LetterSpace]18 -> kVol*Pdiff*(Aprota*a8[t] - Aprotc*c9[t]), v\[LetterSpace]19 -> kVol*(Pdiff*(-(Aprota*a9[t]) + Aprotc*c10[t]) + PPin*(-(Ba*a9[t]) + Bc*c9[t])), v\[LetterSpace]2 -> kVol*Pdiff*(Aprota*a0[t] - Aprotc*c1[t]), v\[LetterSpace]20 -> kVol*Pdiff*(Aprota*a9[t] - Aprotc*c10[t]), v\[LetterSpace]21 -> kVol*(Pdiff*(-(Aprota*a10[t]) + Aprotc*c10[t]) + PPin*(-(Ba*a10[t]) + Bc*c10[t])), v\[LetterSpace]22 -> kVol*Pdiff*(Aprota*a10[t] - Aprotc*c11[t]), v\[LetterSpace]23 -> kVol*(Pdiff*(-(Aprota*a11[t]) + Aprotc*c11[t]) + PPin*(-(Ba*a11[t]) + Bc*c11[t])), v\[LetterSpace]24 -> kVol*Pdiff*(Aprota*a11[t] - Aprotc*c12[t]), v\[LetterSpace]25 -> kVol*(Pdiff*(-(Aprota*a12[t]) + Aprotc*c12[t]) + PPin*(-(Ba*a12[t]) + Bc*c12[t])), v\[LetterSpace]26 -> kVol*Pdiff*(Aprota*a12[t] - Aprotc*c13[t]), v\[LetterSpace]27 -> kVol*(Pdiff*(-(Aprota*a13[t]) + Aprotc*c13[t]) + PPin*(-(Ba*a13[t]) + Bc*c13[t])), v\[LetterSpace]28 -> kVol*Pdiff*(Aprota*a13[t] - Aprotc*c14[t]), v\[LetterSpace]29 -> kVol*(Pdiff*(-(Aprota*a14[t]) + Aprotc*c14[t]) + PPin*(-(Ba*a14[t]) + Bc*c14[t])), v\[LetterSpace]3 -> kVol*(Pdiff*(-(Aprota*a1[t]) + Aprotc*c1[t]) + PPin*(-(Ba*a1[t]) + Bc*c1[t])), v\[LetterSpace]30 -> kVol*Pdiff*(Aprota*a14[t] - Aprotc*c15[t]), v\[LetterSpace]31 -> kVol*(Pdiff*(-(Aprota*a15[t]) + Aprotc*c15[t]) + PPin*(-(Ba*a15[t]) + Bc*c15[t])), v\[LetterSpace]32 -> kVol*Pdiff*(Aprota*a15[t] - Aprotc*c16[t]), v\[LetterSpace]33 -> kVol*(Pdiff*(-(Aprota*a16[t]) + Aprotc*c16[t]) + PPin*(-(Ba*a16[t]) + Bc*c16[t])), v\[LetterSpace]34 -> kVol*Pdiff*(Aprota*a16[t] - Aprotc*c17[t]), v\[LetterSpace]35 -> kVol*(Pdiff*(-(Aprota*a17[t]) + Aprotc*c17[t]) + PPin*(-(Ba*a17[t]) + Bc*c17[t])), v\[LetterSpace]36 -> kVol*Pdiff*(Aprota*a17[t] - Aprotc*c18[t]), v\[LetterSpace]37 -> kVol*(Pdiff*(-(Aprota*a18[t]) + Aprotc*c18[t]) + PPin*(-(Ba*a18[t]) + Bc*c18[t])), v\[LetterSpace]38 -> kVol*Pdiff*(Aprota*a18[t] - Aprotc*c19[t]), v\[LetterSpace]39 -> kVol*(Pdiff*(-(Aprota*a19[t]) + Aprotc*c19[t]) + PPin*(-(Ba*a19[t]) + Bc*c19[t])), v\[LetterSpace]4 -> kVol*Pdiff*(Aprota*a1[t] - Aprotc*c2[t]), v\[LetterSpace]40 -> kVol*Pdiff*(Aprota*a19[t] - Aprotc*c20[t]), v\[LetterSpace]41 -> kVol*(Pdiff*(-(Aprota*a20[t]) + Aprotc*c20[t]) + PPin*(-(Ba*a20[t]) + Bc*c20[t])), v\[LetterSpace]42 -> (2*Dvar*kVol*(a20[t] - Fvar[t]))/(lambda + Ls), v\[LetterSpace]5 -> kVol*(Pdiff*(-(Aprota*a2[t]) + Aprotc*c2[t]) + PPin*(-(Ba*a2[t]) + Bc*c2[t])), v\[LetterSpace]6 -> kVol*Pdiff*(Aprota*a2[t] - Aprotc*c3[t]), v\[LetterSpace]7 -> kVol*(Pdiff*(-(Aprota*a3[t]) + Aprotc*c3[t]) + PPin*(-(Ba*a3[t]) + Bc*c3[t])), v\[LetterSpace]8 -> kVol*Pdiff*(Aprota*a3[t] - Aprotc*c4[t]), v\[LetterSpace]9 -> kVol*(Pdiff*(-(Aprota*a4[t]) + Aprotc*c4[t]) + PPin*(-(Ba*a4[t]) + Bc*c4[t]))
};

parameters = {
Dvar -> 6.7*^-10, F -> 96485.3399, L -> 0.002, Ls -> 0.002, PPin -> 3.3*^-06, Pdiff -> 5.6*^-07, R -> 8.314472, T -> 295.15, V -> -0.12, ext -> 0.0, kVol -> 500.0, lambda -> 5*^-07, lcyt -> 0.0001, pHa -> 5.3, pHc -> 7.2, pK -> 4.8, w -> 1*^-05, default\[LetterSpace]compartment -> 1.0
};

assignments = {
Aprotc -> (1 + 10^(pHc - pK))^(-1), Aaniona -> 1 - Aprota, Bc -> -((Aanionc*subval)/(-1 + 2.71828^(-subval))), Aprota -> (1 + 10^(pHa - pK))^(-1), Ba -> (Aaniona*subval)/(-1 + 2.71828^subval), subval -> -((F*V)/(R*T)), Aanionc -> 1 - Aprotc
};

events = {

};

speciesAnnotations = {

};

reactionAnnotations = {

};

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

(* Time evolution *)
odes = {
Fvar'[t] == 1.0*v\[LetterSpace]42 , S'[t] ==  -1.0*v\[LetterSpace]1, a0'[t] == 40000.0*v\[LetterSpace]1 -40000.0*v\[LetterSpace]2, a1'[t] == 40000.0*v\[LetterSpace]3 -40000.0*v\[LetterSpace]4, a10'[t] == 40000.0*v\[LetterSpace]21 -40000.0*v\[LetterSpace]22, a11'[t] == 40000.0*v\[LetterSpace]23 -40000.0*v\[LetterSpace]24, a12'[t] == 40000.0*v\[LetterSpace]25 -40000.0*v\[LetterSpace]26, a13'[t] == 40000.0*v\[LetterSpace]27 -40000.0*v\[LetterSpace]28, a14'[t] == 40000.0*v\[LetterSpace]29 -40000.0*v\[LetterSpace]30, a15'[t] == 40000.0*v\[LetterSpace]31 -40000.0*v\[LetterSpace]32, a16'[t] == 40000.0*v\[LetterSpace]33 -40000.0*v\[LetterSpace]34, a17'[t] == 40000.0*v\[LetterSpace]35 -40000.0*v\[LetterSpace]36, a18'[t] == 40000.0*v\[LetterSpace]37 -40000.0*v\[LetterSpace]38, a19'[t] == 40000.0*v\[LetterSpace]39 -40000.0*v\[LetterSpace]40, a2'[t] == 40000.0*v\[LetterSpace]5 -40000.0*v\[LetterSpace]6, a20'[t] == 40000.0*v\[LetterSpace]41 -40000.0*v\[LetterSpace]42, a3'[t] == 40000.0*v\[LetterSpace]7 -40000.0*v\[LetterSpace]8, a4'[t] == 40000.0*v\[LetterSpace]9 -40000.0*v\[LetterSpace]10, a5'[t] == 40000.0*v\[LetterSpace]11 -40000.0*v\[LetterSpace]12, a6'[t] == 40000.0*v\[LetterSpace]13 -40000.0*v\[LetterSpace]14, a7'[t] == 40000.0*v\[LetterSpace]15 -40000.0*v\[LetterSpace]16, a8'[t] == 40000.0*v\[LetterSpace]17 -40000.0*v\[LetterSpace]18, a9'[t] == 40000.0*v\[LetterSpace]19 -40000.0*v\[LetterSpace]20, c1'[t] == 20.0*v\[LetterSpace]2 -20.0*v\[LetterSpace]3, c10'[t] == 20.0*v\[LetterSpace]20 -20.0*v\[LetterSpace]21, c11'[t] == 20.0*v\[LetterSpace]22 -20.0*v\[LetterSpace]23, c12'[t] == 20.0*v\[LetterSpace]24 -20.0*v\[LetterSpace]25, c13'[t] == 20.0*v\[LetterSpace]26 -20.0*v\[LetterSpace]27, c14'[t] == 20.0*v\[LetterSpace]28 -20.0*v\[LetterSpace]29, c15'[t] == 20.0*v\[LetterSpace]30 -20.0*v\[LetterSpace]31, c16'[t] == 20.0*v\[LetterSpace]32 -20.0*v\[LetterSpace]33, c17'[t] == 20.0*v\[LetterSpace]34 -20.0*v\[LetterSpace]35, c18'[t] == 20.0*v\[LetterSpace]36 -20.0*v\[LetterSpace]37, c19'[t] == 20.0*v\[LetterSpace]38 -20.0*v\[LetterSpace]39, c2'[t] == 20.0*v\[LetterSpace]4 -20.0*v\[LetterSpace]5, c20'[t] == 20.0*v\[LetterSpace]40 -20.0*v\[LetterSpace]41, c3'[t] == 20.0*v\[LetterSpace]6 -20.0*v\[LetterSpace]7, c4'[t] == 20.0*v\[LetterSpace]8 -20.0*v\[LetterSpace]9, c5'[t] == 20.0*v\[LetterSpace]10 -20.0*v\[LetterSpace]11, c6'[t] == 20.0*v\[LetterSpace]12 -20.0*v\[LetterSpace]13, c7'[t] == 20.0*v\[LetterSpace]14 -20.0*v\[LetterSpace]15, c8'[t] == 20.0*v\[LetterSpace]16 -20.0*v\[LetterSpace]17, c9'[t] == 20.0*v\[LetterSpace]18 -20.0*v\[LetterSpace]19
};
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]}]