Supplemental Mathematica Notebook for: Joseph P. Dexter, Tathagata Dasgupta, and Jeremy Gunawardena, “Invariants reveal multiple forms of robustness in bifunctional enzyme systems.”

This Mathematica notebook is intended to provide details on the algebraic calculations that are omitted from the main paper. It should be read in conjunction with the paper. The notation is the same as in the main paper except that R is used in place of E to denote the bifunctional enzyme, and in certain reaction networks the complexes have been renumbered (as defined later). 

The papers analyzes robustness in 12 reaction networks, which are numbered consecutively as 0 to xi between Fig. 2, 4, and 5 in the main paper. The order of the notebook follows the ordering of these networks.  The biochemical assumptions of and type of robustness observed in each network are summarized in the following table (TC refers to the ternary complex). 

({{Network, Assumptions, Robustness}, {0, full network, robust upper bound}, {i, modification only TC, robust upper bound}, {ii, demodification only TC, robust upper bound}, {iii, ordered binding TC (Sp then S), hybrid robustness}, {iv, ordered binding (Sp then S), modification only TC, robust upper bound}, {v, ordered binding(Sp then S), demodification only TC, robust two-sided bounds on concentration ratio}, {vi, ordered binding TC(S then Sp), hybrid robustness}, {vii, ordered binding(S then Sp), modification only TC, robust two-sided bounds on concentration ratio}, {viii, ordered binding(S then Sp), demodification only TC, robust upper bound}, {ix, no ES, restricted TC formation, ACR}, {x, no ESp, restricted TC formation, ACR}, {xi, no TC formation, robust concentration ratio}})

Preliminaries

MakeCRN

MakeCRN is a custom Mathematica function to display a chemical reaction network, compute the deficiency, and generate the system of polynomial equations (Subscript[f, 1]),...,(Subscript[f, n])  that arise from steady-state mass-action analysis. MakeCRN is contained in the package “reaction-networks.m,” which can be accessed freely and publicly from http://vcp.med.harvard.edu/software.html. 

Ensure that the path to the directory containing “reaction-networks.m” is entered between the quotes in the next command. If the package is in “/home/user/foobar/”, then SetDirectory[“/home/user/foobar”].

SetDirectory[""]

<<reaction-networks.m

Replacement rules

Replacement rule to set all parameters to a random number between 0 and 1 (so that they are not listed by the Mathematica function Variables).

ignoreParameters ={(Subscript[κ, 1->2])->RandomReal[], (Subscript[κ, 2->1])->RandomReal[], (Subscript[κ, 2->3])->RandomReal[], (Subscript[κ, 3->4])->RandomReal[], (Subscript[κ, 4->3])->RandomReal[], (Subscript[κ, 4->1])->RandomReal[], (Subscript[κ, 5->6])->RandomReal[],(Subscript[κ, 6->5])->RandomReal[], (Subscript[κ, 6->9])->RandomReal[], (Subscript[κ, 9->6])->RandomReal[], (Subscript[κ, 6->7])->RandomReal[], (Subscript[κ, 6->8])->RandomReal[],  (Subscript[κ, 8->6])->RandomReal[], (Subscript[κ, 7->6])->RandomReal[], (Subscript[κ, 4->5])->RandomReal[], (Subscript[κ, 5->4])->RandomReal[], (Subscript[κ, 1->4])->RandomReal[], (Subscript[κ, 3->2])->RandomReal[],(Subscript[κ, 7->8])->RandomReal[],(Subscript[κ, 8->7])->RandomReal[]} ;

Replacement rule to set all parameters to 1 (for certain plots).

parametersTo1 ={(Subscript[κ, 1->2])->1, (Subscript[κ, 2->1])->1, (Subscript[κ, 2->3])->1, (Subscript[κ, 3->4])->1, (Subscript[κ, 4->3])->1, (Subscript[κ, 4->1])->1, (Subscript[κ, 5->6])->1,(Subscript[κ, 6->5])->1, (Subscript[κ, 6->9])->1, (Subscript[κ, 9->6])->1, (Subscript[κ, 6->7])->1, (Subscript[κ, 6->8])->1,  (Subscript[κ, 8->6])->1, (Subscript[κ, 7->6])->1, (Subscript[κ, 4->5])->1, (Subscript[κ, 5->4])->1, (Subscript[κ, 1->4])->1, (Subscript[κ, 3->2])->1};

Replacement rule to set S and Sp to 0.

substrateZero={S->0, Sp->0};

Function for inspecting the structure of invariants

substrateVars  = { (S^2)Sp, (Sp^2) S, S Sp, (S^2), (Sp^2), S, Sp};

coeffs[x] returns a 7 x 2 matrix in which the first entry on row i is the i-th monomial in the substrateVars list (above) and the second entry is the corresponding coefficient in the polynomial x.

cq[x_,i_]:={Part[substrateVars,i],FullSimplify[Coefficient[x,Part[substrateVars,i]]/.substrateZero]};
coeffs[x_]:=MatrixForm[Table[cq[x,i],{i,1,7}]];

Network 0

Summary of network 0 (full network)

MakeCRN[{{S + R, RS, r},{RS,R + Sp},{R + Sp,RSp, r},{RSp, S + R}, {RSp + S, RSpS, r}, {RSpS, RSp + Sp}, {RSpS, RS + S}, {RSpS, RS + Sp, r}}]

Total enzyme and total substrate are conserved

Enzyme conservation.

FullSimplify[(Subscript[f, 1])+(Subscript[f, 2])+(Subscript[f, 3])+(Subscript[f, 4])]

Substrate conservation.

FullSimplify[(Subscript[f, 2])+(Subscript[f, 3])+2*(Subscript[f, 4])+(Subscript[f, 5])+(Subscript[f, 6])]

Groebner calculations

We compute a lexicographically ordered Groebner basis on the polynomial system for network 0.

vars = {RSp,RSpS,RS,R,Sp,S};
basis = Table[(Subscript[f, i]), {i, 1, 6}];
grbasis=GroebnerBasis[basis,vars,MonomialOrder -> Lexicographic, CoefficientDomain -> RationalFunctions];
gr[i_]:= Simplify[Part[grbasis,i]];

The Groebner basis has five terms.

Length[grbasis]

There is no invariant in terms of just S or just Sp and the rate constants. Only the first Groebner basis term has no enzyme-substrate complex terms.

Variables[gr[1]/.ignoreParameters]
Variables[gr[2]/.ignoreParameters]
Variables[gr[3]/.ignoreParameters]
Variables[gr[4]/.ignoreParameters]
Variables[gr[5]/.ignoreParameters]

We obtain the cubic invariant by dividing the first term of the first Groebner basis by R.

cubic0=FullSimplify[gr[1]/R]

The structure of the cubic invariant can be readily seen using the coeffs custom function, which generates a table of monomial coefficients.

FullSimplify[coeffs[cubic0]]

Computation of robust upper bounds

Following the notation in Eq. (10) in the main paper, we refer to the seven coefficients as (Subscript[a, 1]), ..., (Subscript[a, 7]).

(Subscript[a, 1])=cq[cubic0,1][[2]]
(Subscript[a, 2])=cq[cubic0,2][[2]]
(Subscript[a, 3])=cq[cubic0,4][[2]]
(Subscript[a, 4])=cq[cubic0,5][[2]]
(Subscript[a, 5])=cq[cubic0,3][[2]]
(Subscript[a, 6])=cq[cubic0,6][[2]]
(Subscript[a, 7])=cq[cubic0,7][[2]]

As derived in the paper, there is the following robust upper bound on [S] when (Subscript[κ, 6->8]) < (Subscript[κ, 6->7]).

bound0a=FullSimplify[-(Subscript[a, 4])/(Subscript[a, 2])]

There is the following robust upper bound on [Sp] when (Subscript[κ, 6->7]) < (Subscript[κ, 6->8]).

bound0b=FullSimplify[-(Subscript[a, 3])/(Subscript[a, 1])]

Plots of the invariant

Plot of the invariant when (Subscript[κ, 6->8]) < (Subscript[κ, 6->7]), indicating a robust upper bound on [S]. The red line is the robust upper bound obtained from the cubic invariant. This is Fig. 3 in the main paper.

numericalValuesNetwork0a ={(Subscript[κ, 1->2])->1, (Subscript[κ, 2->1])->1, (Subscript[κ, 2->3])->1, (Subscript[κ, 3->4])->1, (Subscript[κ, 4->3])->1, (Subscript[κ, 4->1])->1, (Subscript[κ, 5->6])->1,(Subscript[κ, 6->5])->1, (Subscript[κ, 6->9])->1, (Subscript[κ, 9->6])->1, (Subscript[κ, 6->7])->2, (Subscript[κ, 6->8])->1,  (Subscript[κ, 8->6])->1, (Subscript[κ, 7->6])->1, (Subscript[κ, 4->5])->1, (Subscript[κ, 5->4])->1, (Subscript[κ, 1->4])->1, (Subscript[κ, 3->2])->1};

Show[ContourPlot[{(((cubic0/.numericalValuesNetwork0a)))== 0,Sp==0,S==0}, {S,-15,15},{Sp,-15,15},WorkingPrecision -> 100,PerformanceGoal-> "Quality",FrameTicks->{{{-20,-10,-5,0,5,10,20},None},{{-20,-10, -5,0,5,10,20},None}}
,ContourStyle->{{Blue},{Black,Thick},{Black,Thick}}, Axes->True, FrameLabel->{"[S]","[(Subscript[S, P])]"}, RotateLabel->False, LabelStyle->Directive[Bold]],
ContourPlot[{S==(bound0a/.numericalValuesNetwork0a)},{S,-15,15},{Sp,-15,15},ContourStyle->Directive[Red]]]

Plot of the invariant when (Subscript[κ, 6->7]) < (Subscript[κ, 6->8]), indicating a robust upper bound on [Sp]. The red line is the robust upper bound obtained from the cubic invariant.

numericalValuesNetwork0b ={(Subscript[κ, 1->2])->1, (Subscript[κ, 2->1])->1, (Subscript[κ, 2->3])->1, (Subscript[κ, 3->4])->1, (Subscript[κ, 4->3])->1, (Subscript[κ, 4->1])->1, (Subscript[κ, 5->6])->1,(Subscript[κ, 6->5])->1, (Subscript[κ, 6->9])->1, (Subscript[κ, 9->6])->1, (Subscript[κ, 6->7])->1, (Subscript[κ, 6->8])->2,  (Subscript[κ, 8->6])->1, (Subscript[κ, 7->6])->1, (Subscript[κ, 4->5])->1, (Subscript[κ, 5->4])->1, (Subscript[κ, 1->4])->1, (Subscript[κ, 3->2])->1};

Show[ContourPlot[{(((cubic0/.numericalValuesNetwork0b)))== 0,Sp==0,S==0}, {S,-15,15},{Sp,-15,15},WorkingPrecision -> 100,PerformanceGoal-> "Quality",FrameTicks->{{{-20,-10,-5,0,5,10,20},None},{{-20,-10, -5,0,5,10,20},None}}
,ContourStyle->{{Blue},{Black,Thick},{Black,Thick}}, Axes->True,FrameLabel->{"[S]","[(Subscript[S, P])]"}, RotateLabel->False, LabelStyle->Directive[Bold]],
ContourPlot[{Sp==(bound0b/.numericalValuesNetwork0b)},{S,-15,15},{Sp,-15,15},ContourStyle->Directive[Red]]]

Network i

The calculations for the subsequent networks are generally similar to those for network 0 and follow the format in the section above. As such, detailed explanations of each step are not repeated.

For convenience with certain Mathematica functions, the nodes in some networks are numbered differently in the notebook than in the main paper. Tables indicating notational correspondences are included before each renumbered network.

Summary of network i (random binding, modification only TC)

({{node, paper, notebook}, {1, E+S, E+S}, {2, ES, ES}, {3, E+Sp, E+Sp}, {4, ESp, ESp}, {5, ESp+S, ESp+S}, {6, ESpS, ESpS}, {7, ESp+Sp, ESp+Sp}, {8, --, ES+Sp}, {9, ES+Sp, --}})    ({{paper, notebook}, {(Subscript[κ, 1->2]), (Subscript[κ, 1->2])}, {(Subscript[κ, 2->1]), (Subscript[κ, 2->1])}, {(Subscript[κ, 2->3]), (Subscript[κ, 2->3])}, {(Subscript[κ, 3->4]), (Subscript[κ, 3->4])}, {(Subscript[κ, 4->3]), (Subscript[κ, 4->3])}, {(Subscript[κ, 4->1]), (Subscript[κ, 4->1])}, {(Subscript[κ, 5->6]), (Subscript[κ, 5->6])}, {(Subscript[κ, 6->5]), (Subscript[κ, 6->5])}, {(Subscript[κ, 6->7]), (Subscript[κ, 6->7])}, {(Subscript[κ, 6->9]), (Subscript[κ, 6->8])}, {(Subscript[κ, 9->6]), (Subscript[κ, 8->6])}})

MakeCRN[{{S + R, RS, r},{RS,R + Sp},{R + Sp,RSp, r},{RSp, S + R}, {RSp + S, RSpS, r}, {RSpS, RSp + Sp}, {RS+Sp, RSpS, r}}]

Total enzyme and total substrate are conserved

FullSimplify[(Subscript[f, 1])+(Subscript[f, 2])+(Subscript[f, 3])+(Subscript[f, 4])]

FullSimplify[(Subscript[f, 2])+(Subscript[f, 3])+2*(Subscript[f, 4])+(Subscript[f, 5])+(Subscript[f, 6])]

Groebner calculations

vars = {RSp,RSpS,RS,R,Sp,S};
basis = Table[(Subscript[f, i]), {i, 1, 6}];
grbasis=GroebnerBasis[basis,vars,MonomialOrder -> Lexicographic, CoefficientDomain -> RationalFunctions];
gr[i_]:= Simplify[Part[grbasis,i]];

Length[grbasis]

Variables[gr[1]/.ignoreParameters]
Variables[gr[2]/.ignoreParameters]
Variables[gr[3]/.ignoreParameters]
Variables[gr[4]/.ignoreParameters]
Variables[gr[5]/.ignoreParameters]

cubic1=FullSimplify[gr[1]/R]

FullSimplify[coeffs[cubic1]]

Computation of robust upper bound

Following the notation in Eq. (10) in the main paper, we refer to the seven coefficients as (Subscript[a, 1]), ..., (Subscript[a, 7]).

(Subscript[a, 1])=cq[cubic1,1][[2]]
(Subscript[a, 2])=cq[cubic1,2][[2]]
(Subscript[a, 3])=cq[cubic1,4][[2]]
(Subscript[a, 4])=cq[cubic1,5][[2]]
(Subscript[a, 5])=cq[cubic1,3][[2]]
(Subscript[a, 6])=cq[cubic1,6][[2]]
(Subscript[a, 7])=cq[cubic1,7][[2]]

As derived in the paper, there is the following robust upper bound on [S].

bound1=FullSimplify[-(Subscript[a, 4])/(Subscript[a, 2])]

Plot of the invariant

Show[ContourPlot[{(((cubic1/.parametersTo1)))== 0,Sp==0,S==0}, {S,-15,15},{Sp,-15,15},WorkingPrecision -> 100,PerformanceGoal-> "Quality",FrameTicks->{{{-20,-10,-5,0,5,10,20},None},{{-20,-10, -5,0,5,10,20},None}}
,ContourStyle->{{Blue},{Black,Thick},{Black,Thick}}, Axes->True, FrameLabel->{"[S]","[(Subscript[S, P])]"}, RotateLabel->False, LabelStyle->Directive[Bold]],
ContourPlot[{S==(bound1/.parametersTo1)},{S,-15,15},{Sp,-15,15},ContourStyle->Directive[Red]]]

Network ii

Summary of network ii (random binding, demodification only TC).

({{node, paper, notebook}, {1, E+S, E+S}, {2, ES, ES}, {3, E+Sp, E+Sp}, {4, ESp, ESp}, {5, ESp+S, ESp+S}, {6, ESpS, ESpS}, {7, --, ES+Sp}, {8, ES+S, ES+S}, {9, ES+Sp, --}})    ({{paper, notebook}, {(Subscript[κ, 1->2]), (Subscript[κ, 1->2])}, {(Subscript[κ, 2->1]), (Subscript[κ, 2->1])}, {(Subscript[κ, 2->3]), (Subscript[κ, 2->3])}, {(Subscript[κ, 3->4]), (Subscript[κ, 3->4])}, {(Subscript[κ, 4->3]), (Subscript[κ, 4->3])}, {(Subscript[κ, 4->1]), (Subscript[κ, 4->1])}, {(Subscript[κ, 5->6]), (Subscript[κ, 5->6])}, {(Subscript[κ, 6->5]), (Subscript[κ, 6->5])}, {(Subscript[κ, 6->9]), (Subscript[κ, 6->7])}, {(Subscript[κ, 9->6]), (Subscript[κ, 7->6])}, {(Subscript[κ, 6->8]), (Subscript[κ, 6->8])}})

MakeCRN[{{S + R, RS, r},{RS,R + Sp},{R + Sp,RSp, r},{RSp, S + R}, {RSp + S, RSpS, r}, {RS+Sp, RSpS, r}, {RSpS, RS + S}}]

Total enzyme and total substrate are conserved

FullSimplify[(Subscript[f, 1])+(Subscript[f, 2])+(Subscript[f, 3])+(Subscript[f, 4])]

FullSimplify[(Subscript[f, 2])+(Subscript[f, 3])+2*(Subscript[f, 4])+(Subscript[f, 5])+(Subscript[f, 6])]

Groebner calculations

vars = {RSp,RSpS,RS,R,Sp,S};
basis = Table[(Subscript[f, i]), {i, 1, 6}];
grbasis=GroebnerBasis[basis,vars,MonomialOrder -> Lexicographic, CoefficientDomain -> RationalFunctions];
gr[i_]:= Simplify[Part[grbasis,i]];

Length[grbasis]

Variables[gr[1]/.ignoreParameters]
Variables[gr[2]/.ignoreParameters]
Variables[gr[3]/.ignoreParameters]
Variables[gr[4]/.ignoreParameters]
Variables[gr[5]/.ignoreParameters]

cubic2=FullSimplify[gr[1]/R]

FullSimplify[coeffs[cubic2]]

Computation of robust upper bound

Following the notation in Eq. (10) in the main paper, we refer to the seven coefficients as (Subscript[a, 1]), ..., (Subscript[a, 7]).

(Subscript[a, 1])=cq[cubic2,1][[2]]
(Subscript[a, 2])=cq[cubic2,2][[2]]
(Subscript[a, 3])=cq[cubic2,4][[2]]
(Subscript[a, 4])=cq[cubic2,5][[2]]
(Subscript[a, 5])=cq[cubic2,3][[2]]
(Subscript[a, 6])=cq[cubic2,6][[2]]
(Subscript[a, 7])=cq[cubic2,7][[2]]

As derived in the paper, there is the following robust upper bound on [Sp].

bound2=FullSimplify[-(Subscript[a, 3])/(Subscript[a, 1])]

Plot of the invariant

Show[ContourPlot[{(((cubic2/.parametersTo1)))== 0,Sp==0,S==0}, {S,-15,15},{Sp,-15,15},WorkingPrecision -> 100,PerformanceGoal-> "Quality",FrameTicks->{{{-20,-10,-5,0,5,10,20},None},{{-20,-10, -5,0,5,10,20},None}}
,ContourStyle->{{Blue},{Black,Thick},{Black,Thick}}, Axes->True, FrameLabel->{"[S]","[(Subscript[S, P])]"}, RotateLabel->False, LabelStyle->Directive[Bold]],
ContourPlot[{Sp==(bound2/.parametersTo1)},{S,-15,15},{Sp,-15,15},ContourStyle->Directive[Red]]]

Network iii

Summary of network iii (ordered binding, modification and demodification TC).

MakeCRN[{{S + R, RS, r},{RS,R + Sp},{R + Sp,RSp, r},{RSp, S + R}, {RSp + S, RSpS, r}, {RSpS, RSp+Sp}, {RSpS, RS + S}}]

Total enzyme and total substrate are conserved

FullSimplify[(Subscript[f, 1])+(Subscript[f, 2])+(Subscript[f, 3])+(Subscript[f, 4])]

FullSimplify[(Subscript[f, 2])+(Subscript[f, 3])+2*(Subscript[f, 4])+(Subscript[f, 5])+(Subscript[f, 6])]

Groebner calculations

vars = {RSp,RSpS,RS,R,Sp,S};
basis = Table[(Subscript[f, i]), {i, 1, 6}];
grbasis=GroebnerBasis[basis,vars,MonomialOrder -> Lexicographic, CoefficientDomain -> RationalFunctions];
gr[i_]:= Simplify[Part[grbasis,i]];

Length[grbasis]

Variables[gr[1]/.ignoreParameters]
Variables[gr[2]/.ignoreParameters]
Variables[gr[3]/.ignoreParameters]
Variables[gr[4]/.ignoreParameters]

quad3=FullSimplify[gr[1]/R]

FullSimplify[coeffs[quad3]]

Computation of bounds

Following the notation in Eq. (20) in the main paper, we refer to the four non-zero coefficients as (Subscript[a, 11]), ..., (Subscript[a, 14]).

(Subscript[a, 11])=-cq[quad3,6][[2]]
(Subscript[a, 12])=-cq[quad3,4][[2]]
(Subscript[a, 13])=cq[quad3,7][[2]]
(Subscript[a, 14])=cq[quad3,3][[2]]

The behavior of the network depends on the sign of (Subscript[a, 14]).  The ratio [Sp]/[S] is trapped between the following two bounds when((Subscript[κ, 2->1])+(Subscript[κ, 2->3])) (Subscript[κ, 6->7]) < (Subscript[κ, 2->1]) (Subscript[κ, 6->8]).

bound3a=FullSimplify[(Subscript[a, 11])/(Subscript[a, 13])]

bound3b=FullSimplify[(Subscript[a, 12])/(Subscript[a, 14])]

When this condition is met, the ratio [Sp]/[S] is monotonically decreasing.

FullSimplify[(((Subscript[a, 11])/(Subscript[a, 13]))/((Subscript[a, 12])/(Subscript[a, 14])))>1]

When this condition is met, the ratio [Sp]/[S] is monotonically increasing.

FullSimplify[(((Subscript[a, 11])/(Subscript[a, 13]))/((Subscript[a, 12])/(Subscript[a, 14])))<1]

When (Subscript[κ, 2->1]) (Subscript[κ, 6->8]) < ((Subscript[κ, 2->1])+(Subscript[κ, 2->3])) (Subscript[κ, 6->7]), we have the following robust upper bound on [S].

bound3c=FullSimplify[-(Subscript[a, 13])/(Subscript[a, 14])]

Plots of the invariant

sol3=FullSimplify[Solve[quad3==0,Sp]];

numericalValuesNetworkIIIa ={(Subscript[κ, 1->2])->1, (Subscript[κ, 2->1])->1, (Subscript[κ, 2->3])->1, (Subscript[κ, 3->4])->1, (Subscript[κ, 4->3])->1, (Subscript[κ, 4->1])->0.1, (Subscript[κ, 5->6])->1,(Subscript[κ, 6->5])->1, (Subscript[κ, 6->9])->1, (Subscript[κ, 9->6])->1, (Subscript[κ, 6->7])->1, (Subscript[κ, 6->8])->5,  (Subscript[κ, 8->6])->1, (Subscript[κ, 7->6])->1, (Subscript[κ, 4->5])->1, (Subscript[κ, 5->4])->1, (Subscript[κ, 1->4])->1, (Subscript[κ, 3->2])->1};

numericalValuesNetworkIIIb ={(Subscript[κ, 1->2])->1, (Subscript[κ, 2->1])->1, (Subscript[κ, 2->3])->1, (Subscript[κ, 3->4])->1, (Subscript[κ, 4->3])->1, (Subscript[κ, 4->1])->1, (Subscript[κ, 5->6])->1,(Subscript[κ, 6->5])->1, (Subscript[κ, 6->9])->1, (Subscript[κ, 9->6])->1, (Subscript[κ, 6->7])->1, (Subscript[κ, 6->8])->3,  (Subscript[κ, 8->6])->1, (Subscript[κ, 7->6])->1, (Subscript[κ, 4->5])->1, (Subscript[κ, 5->4])->1, (Subscript[κ, 1->4])->1, (Subscript[κ, 3->2])->1};

Monotonically decreasing ratio.

Plot[{(Sp/.sol3)/S/.numericalValuesNetworkIIIa,(bound3a/.numericalValuesNetworkIIIa),(bound3b/.numericalValuesNetworkIIIa)},{S,-10,10},AxesLabel->{"[S]","[(Subscript[S, P])]/[S]"},Ticks->{{-10,-5,5,10},{2,4,6,8,10}},LabelStyle->Directive[Bold]]

Monotonically increasing ratio.

Plot[{(Sp/.sol3)/S/.numericalValuesNetworkIIIb,(bound3a/.numericalValuesNetworkIIIb),(bound3b/.numericalValuesNetworkIIIb)},{S,-40,40},AxesLabel->{"[S]","[(Subscript[S, P])]/[S]"},Ticks->{{-30,-20,-10,10,20,30},{1,2,3,4,5,6}},LabelStyle->Directive[Bold]]

Robust upper bound on [S].

Show[ContourPlot[{(((quad3/.parametersTo1)))== 0,Sp==0,S==0}, {S,-40,40},{Sp,-40,40},WorkingPrecision -> 100,PerformanceGoal-> "Quality",FrameTicks->{{{-40,-30,-20,-10,-5,0,5,10,20,30,40},None},{{-40,-30,-20,-10, -5,0,5,10,20,30,40},None}}
,ContourStyle->{{Blue},{Black,Thick},{Black,Thick}}, Axes->True, FrameLabel->{"[S]","[(Subscript[S, P])]"}, RotateLabel->False, LabelStyle->Directive[Bold]],
ContourPlot[{S==(bound3c/.parametersTo1)},{S,-40,40},{Sp,-40,40},ContourStyle->Directive[Red]]]

Network iv

Summary of network iv (ordered binding, modification only TC).

MakeCRN[{{S + R, RS, r},{RS,R + Sp},{R + Sp,RSp, r},{RSp, S + R}, {RSp + S, RSpS, r}, {RSpS, RSp+Sp}}]

Total enzyme and total substrate are conserved

FullSimplify[(Subscript[f, 1])+(Subscript[f, 2])+(Subscript[f, 3])+(Subscript[f, 4])]

FullSimplify[(Subscript[f, 2])+(Subscript[f, 3])+2*(Subscript[f, 4])+(Subscript[f, 5])+(Subscript[f, 6])]

Groebner calculations

vars = {RSp,RSpS,RS,R,Sp,S};
basis = Table[(Subscript[f, i]), {i, 1, 6}];
grbasis=GroebnerBasis[basis,vars,MonomialOrder -> Lexicographic, CoefficientDomain -> RationalFunctions];
gr[i_]:= Simplify[Part[grbasis,i]];

Length[grbasis]

Variables[gr[1]/.ignoreParameters]
Variables[gr[2]/.ignoreParameters]
Variables[gr[3]/.ignoreParameters]
Variables[gr[4]/.ignoreParameters]

quad4=FullSimplify[gr[1]/R]

FullSimplify[coeffs[quad4]]

Computation of robust upper bound

Following the notation in Eq. (18) in the main paper, we refer to the three coefficients as (Subscript[a, 8]), ..., (Subscript[a, 10]).

(Subscript[a, 8])=cq[quad4,6][[2]]
(Subscript[a, 9])=cq[quad4,7][[2]]
(Subscript[a, 10])=cq[quad4,3][[2]]

As derived in the paper, there is the following robust upper bound on [S].

bound4=FullSimplify[-(Subscript[a, 9])/(Subscript[a, 10])]

Plot of the invariant

Show[ContourPlot[{(((quad4/.parametersTo1)))== 0,Sp==0,S==0}, {S,-15,15},{Sp,-15,15},WorkingPrecision -> 100,PerformanceGoal-> "Quality",FrameTicks->{{{-20,-10,-5,0,5,10,20},None},{{-20,-10, -5,0,5,10,20},None}}
,ContourStyle->{{Blue},{Black,Thick},{Black,Thick}}, Axes->True, FrameLabel->{"[S]","[(Subscript[S, P])]"}, RotateLabel->False, LabelStyle->Directive[Bold]],
ContourPlot[{S==(bound4/.parametersTo1)},{S,-15,15},{Sp,-15,15},ContourStyle->Directive[Red]]]

Network v

Summary of network v (ordered binding, demodification only TC).

({{node, paper, notebook}, {1, E+S, E+S}, {2, ES, ES}, {3, E+Sp, E+Sp}, {4, ESp, ESp}, {5, ESp+S, ESp+S}, {6, ESpS, ESpS}, {7, --, ES+Sp}, {8, ES+S, ES+S}, {9, --, --}})    ({{paper, notebook}, {(Subscript[κ, 1->2]), (Subscript[κ, 1->2])}, {(Subscript[κ, 2->1]), (Subscript[κ, 2->1])}, {(Subscript[κ, 2->3]), (Subscript[κ, 2->3])}, {(Subscript[κ, 3->4]), (Subscript[κ, 3->4])}, {(Subscript[κ, 4->3]), (Subscript[κ, 4->3])}, {(Subscript[κ, 4->1]), (Subscript[κ, 4->1])}, {(Subscript[κ, 5->6]), (Subscript[κ, 5->6])}, {(Subscript[κ, 6->5]), (Subscript[κ, 6->5])}, {(Subscript[κ, 6->8]), (Subscript[κ, 6->7])}})

MakeCRN[{{S + R, RS, r},{RS,R + Sp},{R + Sp,RSp, r},{RSp, S + R}, {RSp + S, RSpS, r}, {RSpS, RS + S}}]

Total enzyme and total substrate are conserved

FullSimplify[(Subscript[f, 1])+(Subscript[f, 2])+(Subscript[f, 3])+(Subscript[f, 4])]

FullSimplify[(Subscript[f, 2])+(Subscript[f, 3])+2*(Subscript[f, 4])+(Subscript[f, 5])+(Subscript[f, 6])]

Groebner calculations

vars = {RSp,RSpS,RS,R,Sp,S};
basis = Table[(Subscript[f, i]), {i, 1, 6}];
grbasis=GroebnerBasis[basis,vars,MonomialOrder -> Lexicographic, CoefficientDomain -> RationalFunctions];
gr[i_]:= Simplify[Part[grbasis,i]];

Length[grbasis]

Variables[gr[1]/.ignoreParameters]
Variables[gr[2]/.ignoreParameters]
Variables[gr[3]/.ignoreParameters]
Variables[gr[4]/.ignoreParameters]

quad5=FullSimplify[gr[1]/R]

FullSimplify[coeffs[quad5]]

Computation of bounds

Following the notation in Eq. (20) the main paper, we refer to the four non-zero coefficients as (Subscript[a, 11]), ..., (Subscript[a, 14]).

(Subscript[a, 11])=-cq[quad5,6][[2]]
(Subscript[a, 12])=-cq[quad5,4][[2]]
(Subscript[a, 13])=cq[quad5,7][[2]]
(Subscript[a, 14])=cq[quad5,3][[2]]

The ratio [Sp]/[S] is trapped between the following two bounds.

bound5a=FullSimplify[(Subscript[a, 11])/(Subscript[a, 13])]

bound5b=FullSimplify[(Subscript[a, 12])/(Subscript[a, 14])]

For this region of parameter space, the ratio [Sp]/[S] is monotonically decreasing.

FullSimplify[(((Subscript[a, 11])/(Subscript[a, 13]))/((Subscript[a, 12])/(Subscript[a, 14])))>1]

For this region of parameter space, the ratio [Sp]/[S] is monotonically increasing.

FullSimplify[(((Subscript[a, 11])/(Subscript[a, 13]))/((Subscript[a, 12])/(Subscript[a, 14])))<1]

Plots of the invariant

sol5=FullSimplify[Solve[quad5==0,Sp]];

numericalValuesNetworkVa ={(Subscript[κ, 1->2])->1, (Subscript[κ, 2->1])->2, (Subscript[κ, 2->3])->1, (Subscript[κ, 3->4])->1, (Subscript[κ, 4->3])->1, (Subscript[κ, 4->1])->1, (Subscript[κ, 5->6])->1,(Subscript[κ, 6->5])->1, (Subscript[κ, 6->9])->1, (Subscript[κ, 9->6])->1, (Subscript[κ, 6->7])->1, (Subscript[κ, 6->8])->1,  (Subscript[κ, 8->6])->1, (Subscript[κ, 7->6])->1, (Subscript[κ, 4->5])->1, (Subscript[κ, 5->4])->1, (Subscript[κ, 1->4])->1, (Subscript[κ, 3->2])->1};

numericalValuesNetworkVb ={(Subscript[κ, 1->2])->1, (Subscript[κ, 2->1])->1, (Subscript[κ, 2->3])->2, (Subscript[κ, 3->4])->1, (Subscript[κ, 4->3])->1, (Subscript[κ, 4->1])->1, (Subscript[κ, 5->6])->1,(Subscript[κ, 6->5])->1, (Subscript[κ, 6->9])->1, (Subscript[κ, 9->6])->1, (Subscript[κ, 6->7])->1, (Subscript[κ, 6->8])->1,  (Subscript[κ, 8->6])->1, (Subscript[κ, 7->6])->1, (Subscript[κ, 4->5])->1, (Subscript[κ, 5->4])->1, (Subscript[κ, 1->4])->1, (Subscript[κ, 3->2])->1};

Monotonically decreasing ratio.

Plot[{(Sp/.sol5)/S/.numericalValuesNetworkVa,(bound5a/.numericalValuesNetworkVa),(bound5b/.numericalValuesNetworkVa)},{S,-20,20}, Ticks->{{-20,-10,0,10,20},{0.3,0.4,0.5,0.6,0.7,0.8,0.9}},AxesLabel->{"[S]","[(Subscript[S, P])]/[S]"},LabelStyle->Directive[Bold]]

Monotonically increasing ratio.

Plot[{(Sp/.sol5)/S/.numericalValuesNetworkVb,(bound5a/.numericalValuesNetworkVb),(bound5b/.numericalValuesNetworkVb)},{S,-20,20},Ticks->{{-20,-10,0,10,20},{0.5,1,1.5,2,2.5,3}},AxesLabel->{"[S]","[(Subscript[S, P])]/[S]"},LabelStyle->Directive[Bold]]

Network vi

Summary of network vi (ordered binding, modification and demodification TC).

({{node, paper, notebook}, {1, E+S, E+S}, {2, ES, ES}, {3, E+Sp, E+Sp}, {4, ESp, ESp}, {5, --, ES+Sp}, {6, ESpS, ESpS}, {7, ESp+Sp, ESp+Sp}, {8, ES+S, ES+S}, {9, ES+Sp, --}})    ({{paper, notebook}, {(Subscript[κ, 1->2]), (Subscript[κ, 1->2])}, {(Subscript[κ, 2->1]), (Subscript[κ, 2->1])}, {(Subscript[κ, 2->3]), (Subscript[κ, 2->3])}, {(Subscript[κ, 3->4]), (Subscript[κ, 3->4])}, {(Subscript[κ, 4->3]), (Subscript[κ, 4->3])}, {(Subscript[κ, 4->1]), (Subscript[κ, 4->1])}, {(Subscript[κ, 6->7]), (Subscript[κ, 6->7])}, {(Subscript[κ, 6->8]), (Subscript[κ, 6->8])}, {(Subscript[κ, 6->9]), (Subscript[κ, 6->5])}, {(Subscript[κ, 9->6]), (Subscript[κ, 5->6])}})

MakeCRN[{{S + R, RS, r},{RS,R + Sp},{R + Sp,RSp, r},{RSp, S + R}, {RS + Sp, RSpS, r}, {RSpS, RSp+Sp},{RSpS, RS + S}}]

Total enzyme and total substrate are conserved

FullSimplify[(Subscript[f, 1])+(Subscript[f, 2])+(Subscript[f, 3])+(Subscript[f, 4])]

FullSimplify[(Subscript[f, 2])+(Subscript[f, 3])+2*(Subscript[f, 4])+(Subscript[f, 5])+(Subscript[f, 6])]

Groebner calculations

vars = {RSp,RSpS,RS,R,Sp,S};
basis = Table[(Subscript[f, i]), {i, 1, 6}];
grbasis=GroebnerBasis[basis,vars,MonomialOrder -> Lexicographic, CoefficientDomain -> RationalFunctions];
gr[i_]:= Simplify[Part[grbasis,i]];

Length[grbasis]

Variables[gr[1]/.ignoreParameters]
Variables[gr[2]/.ignoreParameters]
Variables[gr[3]/.ignoreParameters]
Variables[gr[4]/.ignoreParameters]

quad6=FullSimplify[gr[1]/R]

FullSimplify[coeffs[quad6]]

Computation of bounds

We refer to the four non-zero coefficients as (Subscript[a, 11]), ..., (Subscript[a, 14]).

(Subscript[a, 11])=-cq[quad6,6][[2]]
(Subscript[a, 12])=cq[quad6,5][[2]]
(Subscript[a, 13])=cq[quad6,7][[2]]
(Subscript[a, 14])=-cq[quad6,3][[2]]

The behavior of the network depends on the sign of (Subscript[a, 14]).  The ratio [Sp]/[S] is trapped between the following two bounds when ((Subscript[κ, 4->1])+(Subscript[κ, 4->3])) (Subscript[κ, 6->8])< (Subscript[κ, 4->3]) (Subscript[κ, 6->7]) .

bound6a=FullSimplify[(Subscript[a, 11])/(Subscript[a, 13])]

bound6b=FullSimplify[(Subscript[a, 14])/(Subscript[a, 12])]

When this condition is met, the ratio [Sp]/[S] is monotonically decreasing.

FullSimplify[(((Subscript[a, 11])/(Subscript[a, 13]))/((Subscript[a, 14])/(Subscript[a, 12])))>1]

When this condition is met, the ratio [Sp]/[S] is monotonically increasing.

FullSimplify[(((Subscript[a, 11])/(Subscript[a, 13]))/((Subscript[a, 14])/(Subscript[a, 12])))<1]

When (Subscript[κ, 4->3]) (Subscript[κ, 6->7])<((Subscript[κ, 4->1])+(Subscript[κ, 4->3])) (Subscript[κ, 6->8]), we have the following robust upper bound on [Sp].

bound6c=FullSimplify[-(Subscript[a, 11])/(Subscript[a, 14])]

Plots of the invariant

sol6=FullSimplify[Solve[quad6==0,S]];

numericalValuesNetworkVIa ={(Subscript[κ, 1->2])->1, (Subscript[κ, 2->1])->0.5, (Subscript[κ, 2->3])->1, (Subscript[κ, 3->4])->1, (Subscript[κ, 4->3])->1, (Subscript[κ, 4->1])->1, (Subscript[κ, 5->6])->1,(Subscript[κ, 6->5])->1, (Subscript[κ, 6->9])->1, (Subscript[κ, 9->6])->1, (Subscript[κ, 6->7])->3, (Subscript[κ, 6->8])->1,  (Subscript[κ, 8->6])->1, (Subscript[κ, 7->6])->1, (Subscript[κ, 4->5])->1, (Subscript[κ, 5->4])->1, (Subscript[κ, 1->4])->1, (Subscript[κ, 3->2])->1};

numericalValuesNetworkVIb ={(Subscript[κ, 1->2])->1, (Subscript[κ, 2->1])->0.5, (Subscript[κ, 2->3])->1, (Subscript[κ, 3->4])->1, (Subscript[κ, 4->3])->1, (Subscript[κ, 4->1])->1, (Subscript[κ, 5->6])->1,(Subscript[κ, 6->5])->1, (Subscript[κ, 6->9])->1, (Subscript[κ, 9->6])->1, (Subscript[κ, 6->7])->3, (Subscript[κ, 6->8])->1,  (Subscript[κ, 8->6])->1, (Subscript[κ, 7->6])->1, (Subscript[κ, 4->5])->1, (Subscript[κ, 5->4])->1, (Subscript[κ, 1->4])->1, (Subscript[κ, 3->2])->1};

Monotonically decreasing ratio.

Plot[{1/((S/.sol6)/Sp)/.numericalValuesNetworkVIa,(bound6a/.numericalValuesNetworkVIa),(bound6b/.numericalValuesNetworkVIa)},{Sp,-30,30},AxesLabel->{"[S]","[(Subscript[S, P])]/[S]"},Ticks->{{-30,-20,-10,10,20,30},{0.5,1,1.5,2}},LabelStyle->Directive[Bold]]

Monotonically increasing ratio.

Plot[{(((S/.sol6)/Sp))/.numericalValuesNetworkVIb,(1/bound6a/.numericalValuesNetworkVIb),(1/bound6b/.numericalValuesNetworkVIb)},{Sp,-30,30},AxesLabel->{"[S]","[(Subscript[S, P])]/[S]"},Ticks->{{-30,-20,-10,10,20,30},{1,2,3,4,5,6}},LabelStyle->Directive[Bold]]

Robust upper bound on [Sp].

Show[ContourPlot[{(((quad6/.parametersTo1)))== 0,Sp==0,S==0}, {S,-40,40},{Sp,-40,40},WorkingPrecision -> 100,PerformanceGoal-> "Quality",FrameTicks->{{{-40,-30,-20,-10,-5,0,5,10,20,30,40},None},{{-40,-30,-20,-10, -5,0,5,10,20,30,40},None}}
,ContourStyle->{{Blue},{Black,Thick},{Black,Thick}}, Axes->True, FrameLabel->{"[S]","[(Subscript[S, P])]"}, RotateLabel->False, LabelStyle->Directive[Bold]],
ContourPlot[{Sp==(bound6c/.parametersTo1)},{S,-40,40},{Sp,-40,40
},ContourStyle->Directive[Red]]]

Network vii

Summary of network vii (ordered binding, modification only TC).

({{node, paper, notebook}, {1, E+S, E+S}, {2, ES, ES}, {3, E+Sp, E+Sp}, {4, ESp, ESp}, {5, --, ES+Sp}, {6, ESpS, ESpS}, {7, ESp+Sp, ESp+Sp}, {8, --, --}, {9, ES+Sp, --}})    ({{paper, notebook}, {(Subscript[κ, 1->2]), (Subscript[κ, 1->2])}, {(Subscript[κ, 2->1]), (Subscript[κ, 2->1])}, {(Subscript[κ, 2->3]), (Subscript[κ, 2->3])}, {(Subscript[κ, 3->4]), (Subscript[κ, 3->4])}, {(Subscript[κ, 4->3]), (Subscript[κ, 4->3])}, {(Subscript[κ, 4->1]), (Subscript[κ, 4->1])}, {(Subscript[κ, 6->7]), (Subscript[κ, 6->7])}, {(Subscript[κ, 6->9]), (Subscript[κ, 6->5])}, {(Subscript[κ, 9->6]), (Subscript[κ, 5->6])}})

MakeCRN[{{S + R, RS, r},{RS,R + Sp},{R + Sp,RSp, r},{RSp, S + R}, {RS + Sp, RSpS, r}, {RSpS, RSp+Sp}}]

Total enzyme and total substrate are conserved

FullSimplify[(Subscript[f, 1])+(Subscript[f, 2])+(Subscript[f, 3])+(Subscript[f, 4])]

FullSimplify[(Subscript[f, 2])+(Subscript[f, 3])+2*(Subscript[f, 4])+(Subscript[f, 5])+(Subscript[f, 6])]

Groebner calculations

vars = {RSp,RSpS,RS,R,Sp,S};
basis = Table[(Subscript[f, i]), {i, 1, 6}];
grbasis=GroebnerBasis[basis,vars,MonomialOrder -> Lexicographic, CoefficientDomain -> RationalFunctions];
gr[i_]:= Simplify[Part[grbasis,i]];

Length[grbasis]

Variables[gr[1]/.ignoreParameters]
Variables[gr[2]/.ignoreParameters]
Variables[gr[3]/.ignoreParameters]
Variables[gr[4]/.ignoreParameters]

quad7=FullSimplify[gr[1]/R]

FullSimplify[coeffs[quad7]]

We refer to the four non-zero coefficients as (Subscript[a, 11]), ..., (Subscript[a, 14]).

(Subscript[a, 11])=-cq[quad7,6][[2]]
(Subscript[a, 12])=cq[quad7,5][[2]]
(Subscript[a, 13])=cq[quad7,7][[2]]
(Subscript[a, 14])=-cq[quad7,3][[2]]

The ratio [Sp]/[S] is trapped between the following two bounds.

bound7a=FullSimplify[(Subscript[a, 11])/(Subscript[a, 13])]

bound7b=FullSimplify[(Subscript[a, 14])/(Subscript[a, 12])]

When this condition is met, the ratio [Sp]/[S] is monotonically decreasing.

FullSimplify[(((Subscript[a, 11])/(Subscript[a, 14]))/((Subscript[a, 13])/(Subscript[a, 12])))>1]

When this condition is met, the ratio [Sp]/[S] is monotonically increasing.

FullSimplify[(((Subscript[a, 11])/(Subscript[a, 14]))/((Subscript[a, 13])/(Subscript[a, 12])))<1]

Plots of the invariant

numericalValuesNetworkVIIa ={(Subscript[κ, 1->2])->1, (Subscript[κ, 2->1])->1, (Subscript[κ, 2->3])->1, (Subscript[κ, 3->4])->1, (Subscript[κ, 4->3])->1, (Subscript[κ, 4->1])->2, (Subscript[κ, 5->6])->1,(Subscript[κ, 6->5])->1, (Subscript[κ, 6->9])->1, (Subscript[κ, 9->6])->1, (Subscript[κ, 6->7])->1, (Subscript[κ, 6->8])->1,  (Subscript[κ, 8->6])->1, (Subscript[κ, 7->6])->1, (Subscript[κ, 4->5])->1, (Subscript[κ, 5->4])->1, (Subscript[κ, 1->4])->1, (Subscript[κ, 3->2])->1};

numericalValuesNetworkVIIb ={(Subscript[κ, 1->2])->1, (Subscript[κ, 2->1])->2, (Subscript[κ, 2->3])->1, (Subscript[κ, 3->4])->1, (Subscript[κ, 4->3])->1, (Subscript[κ, 4->1])->1, (Subscript[κ, 5->6])->1,(Subscript[κ, 6->5])->1, (Subscript[κ, 6->9])->1, (Subscript[κ, 9->6])->1, (Subscript[κ, 6->7])->1, (Subscript[κ, 6->8])->1,  (Subscript[κ, 8->6])->1, (Subscript[κ, 7->6])->1, (Subscript[κ, 4->5])->1, (Subscript[κ, 5->4])->1, (Subscript[κ, 1->4])->1, (Subscript[κ, 3->2])->1};

sol7=FullSimplify[Solve[quad7==0,S]];

Monotonically decreasing ratio.

Plot[{1/((S/.sol7)/Sp)/.numericalValuesNetworkVIIa,(bound7a/.numericalValuesNetworkVIIa),(bound7b/.numericalValuesNetworkVIIa)},{Sp,-20,20},AxesLabel->{"[(Subscript[S, P])]","[(Subscript[S, P])]/[S]"},LabelStyle->Directive[Bold]]

Monotonically increasing ratio.

Plot[{1/((S/.sol7)/Sp)/.numericalValuesNetworkVIIb,(bound7a/.numericalValuesNetworkVIIb),(bound7b/.numericalValuesNetworkVIIb)},{Sp,-20,20},AxesLabel->{"[(Subscript[S, P])]","[(Subscript[S, P])]/[S]"},LabelStyle->Directive[Bold]]

Network viii

Summary of network viii (ordered binding, demodification only TC)

({{node, paper, notebook}, {1, E+S, E+S}, {2, ES, ES}, {3, E+Sp, E+Sp}, {4, ESp, ESp}, {5, --, ES+Sp}, {6, ESpS, ESpS}, {7, --, ES+S}, {8, ES+S, --}, {9, ES+Sp, --}})    ({{paper, notebook}, {(Subscript[κ, 1->2]), (Subscript[κ, 1->2])}, {(Subscript[κ, 2->1]), (Subscript[κ, 2->1])}, {(Subscript[κ, 2->3]), (Subscript[κ, 2->3])}, {(Subscript[κ, 3->4]), (Subscript[κ, 3->4])}, {(Subscript[κ, 4->3]), (Subscript[κ, 4->3])}, {(Subscript[κ, 4->1]), (Subscript[κ, 4->1])}, {(Subscript[κ, 6->8]), (Subscript[κ, 6->7])}, {(Subscript[κ, 6->9]), (Subscript[κ, 6->5])}, {(Subscript[κ, 9->6]), (Subscript[κ, 5->6])}})

MakeCRN[{{S + R, RS, r},{RS,R + Sp},{R + Sp,RSp, r},{RSp, S + R}, {RS + Sp, RSpS, r}, {RSpS, RS+S}}]

Total enzyme and total substrate are conserved

FullSimplify[(Subscript[f, 1])+(Subscript[f, 2])+(Subscript[f, 3])+(Subscript[f, 4])]

FullSimplify[(Subscript[f, 2])+(Subscript[f, 3])+2*(Subscript[f, 4])+(Subscript[f, 5])+(Subscript[f, 6])]

Groebner calculations

vars = {RSp,RSpS,RS,R,Sp,S};
basis = Table[(Subscript[f, i]), {i, 1, 6}];
grbasis=GroebnerBasis[basis,vars,MonomialOrder -> Lexicographic, CoefficientDomain -> RationalFunctions];
gr[i_]:= Simplify[Part[grbasis,i]];

Length[grbasis]

Variables[gr[1]/.ignoreParameters]
Variables[gr[2]/.ignoreParameters]
Variables[gr[3]/.ignoreParameters]
Variables[gr[4]/.ignoreParameters]

quad8=FullSimplify[gr[1]/R]

FullSimplify[coeffs[quad8]]

Computation of robust upper bound

We refer to the three non-zero coefficients as (Subscript[a, 8]), ..., (Subscript[a, 10]).

(Subscript[a, 8])=cq[quad8,6][[2]]
(Subscript[a, 9])=-cq[quad8,7][[2]]
(Subscript[a, 10])=-cq[quad8,3][[2]]

There is the following robust upper bound on [Sp].

bound8=FullSimplify[(Subscript[a, 8])/(Subscript[a, 10])]

Plot of the invariant

Show[ContourPlot[{(((quad8/.parametersTo1)))== 0,Sp==0,S==0}, {S,-15,15},{Sp,-15,15},WorkingPrecision -> 100,PerformanceGoal-> "Quality",FrameTicks->{{{-20,-10,-5,0,5,10,20},None},{{-20,-10, -5,0,5,10,20},None}}
,ContourStyle->{{Blue},{Black,Thick},{Black,Thick}}, Axes->True, FrameLabel->{"[S]","[(Subscript[S, P])]"}, RotateLabel->False, LabelStyle->Directive[Bold]],
ContourPlot[{Sp==(bound8/.parametersTo1)},{S,-15,15},{Sp,-15,15},ContourStyle->Directive[Red]]]

Network ix

Summary of network ix (no ES formation, modification only TC)

({{node, paper, notebook}, {1, E+S, E+Sp}, {2, --, ESp}, {3, E+Sp, E+S}, {4, ESp, ESp+S}, {5, ESp+S, ESpS}, {6, ESpS, ESp+Sp}, {7, ESp+Sp, --}, {8, --, --}, {9, --, --}})    ({{paper, notebook}, {(Subscript[κ, 3->4]), (Subscript[κ, 1->2])}, {(Subscript[κ, 4->3]), (Subscript[κ, 2->1])}, {(Subscript[κ, 4->1]), (Subscript[κ, 2->3])}, {(Subscript[κ, 5->6]), (Subscript[κ, 4->5])}, {(Subscript[κ, 6->5]), (Subscript[κ, 5->4])}, {(Subscript[κ, 6->7]), (Subscript[κ, 5->6])}})

MakeCRN[{{R + Sp,RSp, r},{RSp, S + R}, {RSp + S, RSpS, r}, {RSpS, RSp+Sp}}]

Total enzyme and total substrate are conserved

FullSimplify[(Subscript[f, 2])+2*(Subscript[f, 3])+(Subscript[f, 4])+(Subscript[f, 5])]

FullSimplify[(Subscript[f, 1])+(Subscript[f, 2])+(Subscript[f, 3])]

Groebner calculations

vars = {RSp,RSpS,R,Sp,S};
basis = Table[(Subscript[f, i]), {i, 1, 6}];
grbasis=GroebnerBasis[basis,vars,MonomialOrder -> Lexicographic, CoefficientDomain -> RationalFunctions];
gr[i_]:= Simplify[Part[grbasis,i]];

Length[grbasis]

Variables[gr[1]/.ignoreParameters]
Variables[gr[2]/.ignoreParameters]
Variables[gr[3]/.ignoreParameters]
Variables[gr[4]/.ignoreParameters]
Variables[gr[5]/.ignoreParameters]

ACR

ACR in S can be deduced from the first term of the Groebner basis.

FullSimplify[Solve[gr[1]==0,S]]

Network x

Summary of network x (no ESp formation, demodification only TC)

({{node, paper, notebook}, {1, E+S, E+Sp}, {2, ES, ESp}, {3, E+Sp, E+S}, {4, --, ESp+S}, {5, --, ESpS}, {6, ESpS, ESp+Sp}, {7, --, --}, {8, ES+S, --}, {9, ES+Sp, --}})    ({{paper, notebook}, {(Subscript[κ, 1->2]), (Subscript[κ, 1->2])}, {(Subscript[κ, 2->1]), (Subscript[κ, 2->1])}, {(Subscript[κ, 2->3]), (Subscript[κ, 2->3])}, {(Subscript[κ, 6->8]), (Subscript[κ, 5->6])}, {(Subscript[κ, 6->9]), (Subscript[κ, 5->4])}, {(Subscript[κ, 9->6]), (Subscript[κ, 4->5])}})

MakeCRN[{{R + S,RS, r},{RS, Sp + R}, {RS + Sp, RSpS, r}, {RSpS, RS+S}}]

Total enzyme and total substrate are conserved

FullSimplify[(Subscript[f, 2])+2*(Subscript[f, 3])+(Subscript[f, 4])+(Subscript[f, 5])]

FullSimplify[(Subscript[f, 1])+(Subscript[f, 2])+(Subscript[f, 3])]

Groebner calculations

vars = {RSpS,RS,R,Sp,S};
basis = Table[(Subscript[f, i]), {i, 1, 6}];
grbasis=GroebnerBasis[basis,vars,MonomialOrder -> Lexicographic, CoefficientDomain -> RationalFunctions];
gr[i_]:= Simplify[Part[grbasis,i]];

Length[grbasis]

Variables[gr[1]/.ignoreParameters]
Variables[gr[2]/.ignoreParameters]
Variables[gr[3]/.ignoreParameters]
Variables[gr[4]/.ignoreParameters]
Variables[gr[5]/.ignoreParameters]

ACR

ACR in Sp can be deduced from the first term of the Groebner basis.

FullSimplify[Solve[gr[1]==0,Sp]]

Network xi

Summary of network xi (no TC formation)

MakeCRN[{{R + S,RS, r},{RS, Sp + R}, {Sp+R, RSp, r}, {RSp, R+S}}]

Total enzyme and total substrate are conserved

FullSimplify[(Subscript[f, 2])+(Subscript[f, 3])+(Subscript[f, 4])+(Subscript[f, 5])]

FullSimplify[(Subscript[f, 1])+(Subscript[f, 2])+(Subscript[f, 3])]

Groebner calculations

vars = {RSp,RS,R,Sp,S};
basis = Table[(Subscript[f, i]), {i, 1, 6}];
grbasis=GroebnerBasis[basis,vars,MonomialOrder -> Lexicographic, CoefficientDomain -> RationalFunctions];
gr[i_]:= Simplify[Part[grbasis,i]];

Length[grbasis]

Variables[gr[1]/.ignoreParameters]
Variables[gr[2]/.ignoreParameters]
Variables[gr[3]/.ignoreParameters]
Variables[gr[4]/.ignoreParameters]

Robust ratio

The first term of the Groebner basis is linear and implies a robust ratio in S and Sp.

linear11=gr[1];

FullSimplify[Solve[linear11==0,S]]

Plot of the invariant

Show[ContourPlot[{(((linear11/.parametersTo1)))== 0,Sp==0,S==0}, {S,-25,25},{Sp,-15,15},WorkingPrecision -> 100,PerformanceGoal-> "Quality",FrameTicks->{{{-20,-10,-5,0,5,10,20},None},{{-20,-10, -5,0,5,10,20},None}}
,ContourStyle->{{Blue},{Black,Thick},{Black,Thick}}, Axes->True, FrameLabel->{"[S]","[(Subscript[S, P])]"}, RotateLabel->False, LabelStyle->Directive[Bold]]]