Supplemental Mathematica Notebook for Joseph P. Dexter and Jeremy Gunawardena, “Dimerization and bifunctionality confer robustness to the isocitrate dehydrogenase regulatory system in Escherichia coli”

This Mathematica notebook is intended to provide details about the derivation of the key invariant, the numerical calculations with the invariant, and the analysis of alternative models that are omitted from the main paper. It should be read in parallel with the relevant sections of the paper. The notation is the same as in the main paper except that R is used in place of E to denote IDHKP and (Subscript[x, 1])...(Subscript[x, 10]) are used to refer to the species in the network, as defined below and in Table 1 in the main paper.

The notebook requires custom Mathematica functions which are provided in the accompanying file CRN-functions.m, which is obtainable from our website at http://vcp.med.harvard.edu/papers.html.

1) Reaction network for the IDH regulatory system

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

SetDirectory[""]

<<CRN-functions.m

Reaction network for the model

MakeCRN is a custom  Mathematica function to display a chemical reaction network, which is provided in the accompanying package CRN-functions.m.

MakeCRN[{{ Ioo + R, RIoo,  r}, {RIoo, R + Iop}, {R + Iop, RIop, r}, {RIop, Ioo + R}, {RIop, Ipp + R}, {R + Ipp, RIpp, r}, {RIpp, R + Iop}, {RIop + R, RIopR, r},{RIopR, RIpp + R}, {RIopR, RIoo + R}, {RIpp + R, RIppR, r}, {RIoo + R, RIooR, r}, {RIppR, RIop + R}, {RIooR, RIop +  R}}, True, True]

2) Derivation of the key invariant (Eq. (2) in the main paper)

Replacement rules

For later use, we introduce the following replacement rule to assign numerical values to the rate constants. Parameter values are from Table 2 in the main text.

parameterClasses = {(Subscript[k, 1])->(Subscript[k, on]),(Subscript[k, 2])->(Subscript[k, off]),(Subscript[k, 3])->(Subscript[k, cat]), (Subscript[k, 4])->(Subscript[k, on]), (Subscript[k, 5])->(Subscript[k, off]), (Subscript[k, 6])->(Subscript[κ, cat]),(Subscript[k, 7])->(Subscript[k, cat]),(Subscript[k, 8])->(Subscript[κ, on]), (Subscript[k, 9])->(Subscript[κ, off]), (Subscript[k, 10])->(Subscript[κ, cat]),(Subscript[k, 11])->(Subscript[k, on]),(Subscript[k, 12])->(Subscript[k, off]),(Subscript[k, 13])->(Subscript[k, cat]),(Subscript[k, 14])->(Subscript[κ, cat]),(Subscript[k, 15])->(Subscript[κ, on]), (Subscript[k, 16])->(Subscript[κ, off]), (Subscript[k, 17])->(Subscript[k, on]),(Subscript[k, 18]) ->(Subscript[k, off]),(Subscript[k, 19])->(Subscript[κ, cat]),(Subscript[k, 20])->(Subscript[k, cat])};

numericalValues = {(Subscript[k, on])->1*10^8, (Subscript[k, off])->23, (Subscript[k, cat])->9.8*10^-2,(Subscript[κ, on])->1*10^8 , (Subscript[κ, off])->40, (Subscript[κ, cat])->6.5*10^-2};

Replacement rule for (Subscript[x, 3]) and (Subscript[x, 10]) at the midpoint in the grid of values considered.

xValues = {(Subscript[x, 3])->0.3, (Subscript[x, 10])->17.5};

Groebner basis calculation on the polynomial system describing the reaction network

(Subscript[f, 1])=-(Subscript[x, 2]) (Subscript[x, 1]) (Subscript[k, 1])+(Subscript[x, 4]) (Subscript[k, 2])+(Subscript[x, 4]) (Subscript[k, 3])-(Subscript[x, 3]) (Subscript[x, 1]) (Subscript[k, 4])+(Subscript[x, 6]) (Subscript[k, 6])+(Subscript[x, 6]) (Subscript[k, 5])+(Subscript[x, 6]) (Subscript[k, 7])-(Subscript[x, 10]) (Subscript[x, 1]) (Subscript[k, 8])+(Subscript[x, 5]) (Subscript[k, 10])+(Subscript[x, 5]) (Subscript[k, 9])-(Subscript[x, 1]) (Subscript[x, 6]) (Subscript[k, 11])+(Subscript[x, 8]) (Subscript[k, 12])+(Subscript[x, 8]) (Subscript[k, 13])+(Subscript[x, 8]) (Subscript[k, 14])-(Subscript[x, 1]) (Subscript[x, 5]) (Subscript[k, 15])-(Subscript[x, 1]) (Subscript[x, 4]) (Subscript[k, 17])+(Subscript[x, 9]) (Subscript[k, 19])+(Subscript[x, 9]) (Subscript[k, 16])+(Subscript[x, 7]) (Subscript[k, 20])+(Subscript[x, 7]) (Subscript[k, 18]);
(Subscript[f, 2])=-(Subscript[x, 2]) (Subscript[x, 1]) (Subscript[k, 1])+(Subscript[x, 4]) (Subscript[k, 2])+(Subscript[x, 6]) (Subscript[k, 6]);
(Subscript[f, 3])=(Subscript[x, 4]) (Subscript[k, 3])-(Subscript[x, 3]) (Subscript[x, 1]) (Subscript[k, 4])+(Subscript[x, 6]) (Subscript[k, 5])+(Subscript[x, 5]) (Subscript[k, 10]);
(Subscript[f, 4])=(Subscript[x, 2]) (Subscript[x, 1]) (Subscript[k, 1])-(Subscript[x, 4]) ((Subscript[k, 2])+(Subscript[k, 3]))+(Subscript[x, 8]) (Subscript[k, 14])-(Subscript[x, 1]) (Subscript[x, 4]) (Subscript[k, 17])+(Subscript[x, 7]) (Subscript[k, 18]);
(Subscript[f, 5])=(Subscript[x, 10]) (Subscript[x, 1]) (Subscript[k, 8])-(Subscript[x, 5]) ((Subscript[k, 10])+(Subscript[k, 9]))+(Subscript[x, 8]) (Subscript[k, 13])-(Subscript[x, 1]) (Subscript[x, 5]) (Subscript[k, 15])+(Subscript[x, 9]) (Subscript[k, 16]);
(Subscript[f, 6])=(Subscript[x, 3]) (Subscript[x, 1]) (Subscript[k, 4])-(Subscript[x, 6]) ((Subscript[k, 6])+(Subscript[k, 5])+(Subscript[k, 7]))-(Subscript[x, 1]) (Subscript[x, 6]) (Subscript[k, 11])+(Subscript[x, 8]) (Subscript[k, 12])+(Subscript[x, 9]) (Subscript[k, 19])+(Subscript[x, 7]) (Subscript[k, 20]);
(Subscript[f, 7])=(Subscript[x, 1]) (Subscript[x, 4]) (Subscript[k, 17])-(Subscript[x, 7]) ((Subscript[k, 20])+(Subscript[k, 18]));
(Subscript[f, 8])=(Subscript[x, 1]) (Subscript[x, 6]) (Subscript[k, 11])-(Subscript[x, 8]) ((Subscript[k, 12])+(Subscript[k, 13])+(Subscript[k, 14]));
(Subscript[f, 9])=(Subscript[x, 1]) (Subscript[x, 5]) (Subscript[k, 15])-(Subscript[x, 9]) ((Subscript[k, 19])+(Subscript[k, 16]));
(Subscript[f, 10])=(Subscript[x, 6]) (Subscript[k, 7])-(Subscript[x, 10]) (Subscript[x, 1]) (Subscript[k, 8])+(Subscript[x, 5]) (Subscript[k, 9]);

vars = { (Subscript[x, 4]), (Subscript[x, 5]),(Subscript[x, 6]), (Subscript[x, 7]), (Subscript[x, 8]), (Subscript[x, 9]), (Subscript[x, 1]), (Subscript[x, 3]), (Subscript[x, 10]),(Subscript[x, 2])};
basis = Table[(Subscript[f, i]), {i, 1, 10}];
grbasis=GroebnerBasis[basis,vars,MonomialOrder -> Lexicographic, CoefficientDomain -> RationalFunctions];
gr[i_]:= Simplify[Part[grbasis,i]];

The Groebner basis has nine terms.

Length[grbasis]

There is no invariant in terms of just (Subscript[x, 2]) and the rate constants.

Variables[gr[1]/.parameterClasses/.numericalValues]
Variables[gr[2]/.parameterClasses/.numericalValues]
Variables[gr[3]/.parameterClasses/.numericalValues]
Variables[gr[4]/.parameterClasses/.numericalValues]
Variables[gr[5]/.parameterClasses/.numericalValues]
Variables[gr[6]/.parameterClasses/.numericalValues]
Variables[gr[7]/.parameterClasses/.numericalValues]
Variables[gr[8]/.parameterClasses/.numericalValues]
Variables[gr[9]/.parameterClasses/.numericalValues]

The first three terms contain (Subscript[x, 1]), (Subscript[x, 2]), (Subscript[x, 3]),  and (Subscript[x, 10]), but (Subscript[x, 1]) can only be factored out of the first term. As such, we focus on this invariant for the rest of the notebook.

Collect[gr[1],(Subscript[x, 1]),Simplify]
Collect[gr[2],(Subscript[x, 1]),Simplify]
Collect[gr[3],(Subscript[x, 1]),Simplify]

Eq. (2) is obtained by dividing the first term of the Groebner basis by (Subscript[x, 1])

FullSimplify[gr[1]/(Subscript[x, 1])]

Coefficients of the terms in the invariant

The invariant has the form (Subscript[x, 2]^2) + (Subscript[a, 2])(Subscript[x, 3]^2) + (Subscript[a, 3])(Subscript[x, 10]^2) + (Subscript[a, 4])(Subscript[x, 2])*(Subscript[x, 3]) + (Subscript[a, 5])(Subscript[x, 2])*(Subscript[x, 10]) + (Subscript[a, 6])(Subscript[x, 3])*(Subscript[x, 10]) = 0 when normalized by (Subscript[a, 1]) (the coefficient of (Subscript[x, 2]^2)).

invar1=Expand[gr[1]/(Subscript[x, 1])];

norm=Simplify[Coefficient[invar1,(Subscript[x, 2]),2]];

num1=Simplify[Coefficient[invar1, (Subscript[x, 2]),2]];
num2=Simplify[Coefficient[invar1, (Subscript[x, 3]),2]];
num3=Simplify[Coefficient[invar1, (Subscript[x, 10]),2]];
num4=Simplify[Coefficient[invar1, (Subscript[x, 2])(Subscript[x, 3]),1]];
num5=Simplify[Coefficient[invar1, (Subscript[x, 2])(Subscript[x, 10]),1]];
num6=Simplify[Coefficient[invar1, (Subscript[x, 3])(Subscript[x, 10]),1]];

coefficients={a1->Simplify[num1/norm],
a2->Simplify[num2/norm],
a3->Simplify[num3/norm],
a4->Simplify[num4/norm],
a5->Simplify[num5/norm],
a6->Simplify[num6/norm]}

Coefficients in terms of the six parameter classes

For clarity, the invariant will be written in terms of the six classes of parameters (rather than the underlying 20 rate constants) throughout the remainder of the notebook.

Simplify[a1/.coefficients/.parameterClasses]
Simplify[a2/.coefficients/.parameterClasses]
Simplify[a3/.coefficients/.parameterClasses]
Simplify[a4/.coefficients/.parameterClasses]
Simplify[a5/.coefficients/.parameterClasses]
Simplify[a6/.coefficients/.parameterClasses]

Numerical values for the coefficients

At the reference parameter values, the coefficients have the following numerical values (given in Table 3 in the main text).

aValues=List[Simplify[a1/.coefficients/.parameterClasses/.numericalValues],
Simplify[a2/.coefficients/.parameterClasses/.numericalValues],
Simplify[a3/.coefficients/.parameterClasses/.numericalValues],
Simplify[a4/.coefficients/.parameterClasses/.numericalValues],
Simplify[a5/.coefficients/.parameterClasses/.numericalValues],
Simplify[a6/.coefficients/.parameterClasses/.numericalValues]];

Grid[{{"(Subscript[a, 1])","(Subscript[a, 2])","(Subscript[a, 3])","(Subscript[a, 4])","(Subscript[a, 5])","(Subscript[a, 6])"},aValues},Alignment->Left,Spacings->{2,1},Frame->All,ItemStyle->"Text"]

3) Prediction of robustness from exact calculations with the invariant

Solving the invariant for (Subscript[x, 2])

The invariant is quadratic in (Subscript[x, 2]) and thus yields two solutions. Only the second solution is positive at the reference parameter values.

invar1a = (Subscript[x, 2]^2) + a2*(Subscript[x, 3]^2) + a3*(Subscript[x, 10]^2) + a4*(Subscript[x, 2])*(Subscript[x, 3]) + a5*(Subscript[x, 2])*(Subscript[x, 10]) + a6*(Subscript[x, 3])*(Subscript[x, 10]);

solall=Solve[invar1a ==0, (Subscript[x, 2])];
sol1 = ((Subscript[x, 2])/.solall[[1]]);
sol2 =((Subscript[x, 2])/.solall[[2]]);
(sol1/.coefficients/.parameterClasses/.numericalValues/.xValues)
(sol2/.coefficients/.parameterClasses/.numericalValues/.xValues)

The second solution in terms of the coefficients.

FullSimplify[sol2]

Plot of (Subscript[x, 2]) in terms of (Subscript[x, 3]) and (Subscript[x, 10])

Plot3D[(sol2/.coefficients/.parameterClasses/.numericalValues/.xValues),  {(Subscript[x, 10]), 0, 35}, {(Subscript[x, 3]), 0, 1}, AxesLabel -> Automatic, BaseStyle -> {FontWeight -> "Bold", FontSize -> 16}]

Partial derivative of (Subscript[x, 2]) with respect to (Subscript[x, 10])

The following expression is computed by differentiating the second solution to the invariant.

d2 =FullSimplify[ D[sol2,(Subscript[x, 10])] ]

Plot of the partial derivative in terms of (Subscript[x, 3]) and (Subscript[x, 10])

The flatness of this plot indicates robust behavior.

Plot3D[(d2/.coefficients/.parameterClasses/.numericalValues/.xValues),  {(Subscript[x, 10]), 0, 35}, {(Subscript[x, 3]), 0, 1}, AxesLabel -> Automatic, PlotRange -> {-0.5, 0.1},  BaseStyle -> {FontWeight -> "Bold", FontSize -> 16}]

4) Explanation for robustness using the approximate conservation law

Elimination of (Subscript[x, 10 ])or (Subscript[x, 2]) using the approximate conservation law for total IDH

Replacement rules from the conservation law and for fixing (Subscript[x, 3]) at ϵ.

conservation1 = {(Subscript[x, 10])->T -(Subscript[x, 2])-(Subscript[x, 3])};
conservation2 ={(Subscript[x, 2])->T -(Subscript[x, 3])-(Subscript[x, 10])};
small={(Subscript[x, 3])->ϵ};

Expression for (Subscript[x, 2]) in terms of ϵ and T.

solx2=Simplify[Solve[(sol2/.conservation1/.small)==(Subscript[x, 2]), (Subscript[x, 2])]];
(Subscript[x, 2])/.solx2[[2]]

Plot of (Subscript[x, 2]) in terms of total IDH (T) at a particular ϵ value

epsilonValue={ϵ->0.5};

Plot[((Subscript[x, 2])/.solx2[[2]]/.coefficients/.parameterClasses/.numericalValues/.epsilonValue), {T, 0, 10},PlotRange->{0,3}, AxesLabel->Automatic]

Plot of (Subscript[x, 10]) in terms of T  at a fixed (Subscript[x, 3]) value

solx10=Simplify[Solve[(sol2/.conservation2/.small)==((Subscript[x, 2])/.conservation2), (Subscript[x, 10])]];

Plot[((Subscript[x, 10])/.solx10[[1]]/.coefficients/.small/.parameterClasses/.numericalValues/.epsilonValue), {T, 0, 10}, PlotRange -> {0, 10}, AxesLabel->Automatic]

5) Dominant terms in the invariant

Minimum and maximum terms in the coefficients

Each monomial in the numerator of (Subscript[a, 2])...(Subscript[a, 6]) is written below.

mon2=MonomialList[Coefficient[invar1,(Subscript[x, 3]),2]/.parameterClasses]
mon3=MonomialList[Coefficient[invar1,(Subscript[x, 10]),2]/.parameterClasses]
mon4=MonomialList[Coefficient[invar1,(Subscript[x, 2])(Subscript[x, 3]),1]/.parameterClasses]
mon5=MonomialList[Coefficient[invar1,(Subscript[x, 2])(Subscript[x, 10]),1]/.parameterClasses]
mon6=MonomialList[Coefficient[invar1,(Subscript[x, 3])(Subscript[x, 10]),1]/.parameterClasses]

At the reference parameter values, the largest and smallest monomials in each coefficient have the following (normalized) values.

max2=(Max[mon2]/norm)/.parameterClasses/.numericalValues;
min2=(Min[mon2]/norm)/.parameterClasses/.numericalValues;
max3=(Max[mon3]/norm)/.parameterClasses/.numericalValues;
min3=(Min[mon3]/norm)/.parameterClasses/.numericalValues;
max4=(Max[mon4]/norm)/.parameterClasses/.numericalValues;
min4=(Min[mon4]/norm)/.parameterClasses/.numericalValues;
max5=(Max[mon5]/norm)/.parameterClasses/.numericalValues;
min5=(Min[mon5]/norm)/.parameterClasses/.numericalValues;
max6=(Max[mon6]/norm)/.parameterClasses/.numericalValues;
min6=(Min[mon6]/norm)/.parameterClasses/.numericalValues;
max=List[max2,max3,max4,max5,max6];
min=List[min2,min3,min4,min5,min6];

Grid[{{"","(Subscript[a, 2])","(Subscript[a, 3])","(Subscript[a, 4])","(Subscript[a, 5])","(Subscript[a, 6])"},Prepend[max,"max"],Prepend[min,"min"]},Alignment->Left,Spacings->{2,1},Frame->All,ItemStyle->"Text"]

Dominant terms for (Subscript[a, 5]) and (Subscript[a, 6])

(Subscript[a, 5]) and (Subscript[a, 6]) are approximated to 99% accuarcy by their two largest terms .

approxa5=Simplify[max5+min5]
approxa6=Simplify[max6+min6]
approxa5/a5/.coefficients/.parameterClasses/.numericalValues
approxa6/a6/.coefficients/.parameterClasses/.numericalValues

The symbolic expressions for the dominant terms .

Position[mon5/.parameterClasses/.numericalValues, Max[mon5/.parameterClasses/.numericalValues]]
Position[mon5/.parameterClasses/.numericalValues, Min[mon5/.parameterClasses/.numericalValues]]
Position[mon6/.parameterClasses/.numericalValues, Max[mon6/.parameterClasses/.numericalValues]]
Position[mon6/.parameterClasses/.numericalValues, Min[mon6/.parameterClasses/.numericalValues]]

mon5[[37]]
mon5[[35]]
mon6[[32]]
mon6[[34]]

These are the dominant monomials, as shown in Eq. (6)-(8) in the main paper.

Excluding the min and max terms, there are no terms in (Subscript[a, 5]) and (Subscript[a, 6]) outside of [-200, 200].

(Max[Delete[mon5, {{37},{35}}]]/norm)/.parameterClasses/.numericalValues
(Min[Delete[mon5, {{37},{35}}]]/norm)/.parameterClasses/.numericalValues
(Max[Delete[mon6, {{32},{34}}]]/norm)/.parameterClasses/.numericalValues
(Min[Delete[mon6, {{32},{34}}]]/norm)/.parameterClasses/.numericalValues

6) Sensitivity analysis

Generation of random parameter values according to constraints

We first randomly generate 200,000 combinations of parameters. Each of the six parameters is picked independently from the interval [v/2, 2v], where v is the reference value.

randomParameters={(Subscript[k, on])->RandomReal[{0.5*1*10^8,2*1*10^8}, {200000}],
(Subscript[k, off])->RandomReal[{0.5*23, 2*23}, {200000}],
(Subscript[k, cat])->RandomReal[{0.0098*0.5, 2*0.0098}, {200000}],
(Subscript[κ, on])->RandomReal[{0.5*1*10^8,2*1*10^8}, {200000}],
(Subscript[κ, off])->RandomReal[{0.5*40, 2*40}, {200000}],
(Subscript[κ, cat])->RandomReal[{0.0065*0.5, 2*0.0065}, {200000}]};

The kinase and phosphatase Michaelis-Menten constants.

(Subscript[K, M])=((Subscript[k, off])+(Subscript[k, cat]))/(Subscript[k, on]);
(Subscript[κ, M])=((Subscript[κ, off])+(Subscript[κ, cat]))/(Subscript[κ, on]);

Parameter combinations are excluded if they correspond to one of the Michaelis-Menten constants being outside of [0.8v, 1.2v], where v is the reference value. The following code generates the list of the parameter combinations that lead to unacceptable Michaelis-Menten values.  The reference Michaelis-Menten values are experimentally determined and given in Table 2 in the main paper.

ex1=((Subscript[K, M])/.parameterClasses/.randomParameters)/.x_/;x>1.2*2.3*10^-7->t;
ex2=((Subscript[K, M])/.parameterClasses/.randomParameters)/.x_/;x<0.8*2.3*10^-7->t;
ex3 =((Subscript[κ, M])/.parameterClasses/.randomParameters)/.x_/;x>1.2*4*10^-7->t;
ex4=((Subscript[κ, M])/.parameterClasses/.randomParameters)/.x_/;x<0.8*4*10^-7->t;
pos1=Position[ex1, t];
pos2=Position[ex2, t];
pos3 = Position[ex3, t];
pos4 = Position[ex4, t];
combined = Union[pos1, pos2, pos3, pos4];

Removing those leaves around 18,000 acceptable combinations.

Simplify[200000-Length[combined]]

We check that the Michaelis-Menten constants for the remaining combinations are indeed restricted to [0.8v, 1.2v].

Max[Delete[((Subscript[K, M])/.parameterClasses/.randomParameters),combined]]
Min[Delete[((Subscript[K, M])/.parameterClasses/.randomParameters),combined]]
Max[Delete[((Subscript[κ, M])/.parameterClasses/.randomParameters),combined]]
Min[Delete[((Subscript[κ, M])/.parameterClasses/.randomParameters),combined]]

One positive solution to the invariant is obtained for all parameter combinations

sol1mod and sol2mod are matrices with five columns; each column contains the solutions for all parameter combinations at one of the (Subscript[x, 3]) and (Subscript[x, 10]) values. Column 1: ((Subscript[x, 3]), (Subscript[x, 10])) = (0.3, 17.5), Column 2: ((Subscript[x, 3]), (Subscript[x, 10])) = (0.1, 5), Column 3: ((Subscript[x, 3]), (Subscript[x, 10])) = (0.1, 30), Column 4: ((Subscript[x, 3]), (Subscript[x, 10])) = (0.5, 5), Column 5: ((Subscript[x, 3]), (Subscript[x, 10])) = (0.5, 30).

xGrid = {(Subscript[x, 3])->{0.3, 0.1,0.1, 0.5,0.5},(Subscript[x, 10])->{17.5, 5,30,5, 30}};

sol1mod=Delete[sol1/.coefficients/.parameterClasses/.randomParameters/.xGrid, combined];
sol2mod=Delete[sol2/.coefficients/.parameterClasses/.randomParameters/.xGrid, combined];

For all parameter combinations, only one positive (Subscript[x, 2]) solution is obtained at the various (Subscript[x, 3]) and (Subscript[x, 10]) values.

count1=BinCounts[sol1mod[[All,1]], {{0, Infinity}} ];
count2=BinCounts[sol2mod[[All,1]], {{0, Infinity}}];
count3=BinCounts[sol1mod[[All,2]], {{0, Infinity}} ];
count4=BinCounts[sol2mod[[All,2]], {{0, Infinity}}];
count5=BinCounts[sol1mod[[All,3]], {{0, Infinity}} ];
count6=BinCounts[sol2mod[[All,3]], {{0, Infinity}}];
count7=BinCounts[sol1mod[[All,4]], {{0, Infinity}} ];
count8=BinCounts[sol2mod[[All,4]], {{0, Infinity}}];
count9=BinCounts[sol1mod[[All,5]], {{0, Infinity}} ];
count10=BinCounts[sol2mod[[All,5]], {{0, Infinity}}];
negsol=List[count1,count3,count5,count7,count9];
possol=List[count2,count4,count6,count8,count10];

Grid[{{(Subscript[x, 3]),0.3,0.1,0.1,0.5,0.5},{(Subscript[x, 10]),17.5,5,30,5,30},Prepend[negsol,"positive sol1"],Prepend[possol,"positive sol2"]},Alignment->Left,Spacings->{2,1},Frame->All,ItemStyle->"Text"]

Partial derivative values indicate parameter-independent robustness

Additionally, the mean partial derivative values are on the order of 10^-4, confirming that robustness is maintained.  (d2, the partial derivative of the positive solution to the invariant, is defined earlier in the notebook.)

m=Mean[Delete[d2/.coefficients/.parameterClasses/.randomParameters/.xGrid, combined]];
s=StandardDeviation[Delete[d2/.coefficients/.parameterClasses/.randomParameters/.xGrid, combined]];

Grid[{{(Subscript[x, 3]),0.3,0.1,0.1,0.5,0.5},{(Subscript[x, 10]),17.5,5,30,5,30},Prepend[m,"mean"],Prepend[s,"SD"]},Alignment->Left,Spacings->{2,1},Frame->All,ItemStyle->"Text"]

7) Analysis of hypothetical model of IDH regulation

Groebner basis calculation for network with monomeric IDH and bifunctional IDHKP

(Subscript[f, 1]) =-Io R (Subscript[k, 1])+RIo (Subscript[k, 2])+RIp (Subscript[k, 6]);
(Subscript[f, 2]) = RIo (Subscript[k, 3])-Ip R (Subscript[k, 4])+RIp (Subscript[k, 5]);
(Subscript[f, 3]) = -Io R (Subscript[k, 1])+RIo (Subscript[k, 2])+RIo (Subscript[k, 3])-Ip R (Subscript[k, 4])+RIp (Subscript[k, 6])+RIp (Subscript[k, 5]);
(Subscript[f, 4 ])= Io R (Subscript[k, 1])-RIo ((Subscript[k, 2])+(Subscript[k, 3]));
(Subscript[f, 5 ])= Ip R (Subscript[k, 4])-RIp ((Subscript[k, 6])+(Subscript[k, 5]));

vars = {R, Io, Ip, RIo, RIp};
basis = Table[(Subscript[f, i]), {i, 1, 5}];
grbasis=GroebnerBasis[basis,vars,MonomialOrder -> Lexicographic, CoefficientDomain -> RationalFunctions];
gr[i_]:= Simplify[Part[grbasis,i]]

The Groebner basis has four terms.

Length[grbasis]

The robust ratio can be derived from Groebner basis terms 3 and 4

t1=Solve[gr[3]==0, RIp]

t2=Solve[gr[4]==0, RIp]

FullSimplify[Solve[(RIp/.t1)==(RIp/.t2), Io]]