A workbook to analytically solve for the steady states of the PFK subnetwork

1. Define the structure of the subnetwork

variables = {R0,R0A,R0AF,R1,R1A,R1AF,R2,R2A,R2AF,R3,R3A,R3AF,R4,R4A,R4AF,T0,T1,T2,T3,T4};
varnames={"R","RA","RAF","R1","R1A","R1AF","R2","R2A","R2AF","R3","R3A","R3AF","R4","R4A","R4AF","T0","T1","T2","T3","T4"};
fluxnames={"vR01","vR02","vR03","vR10","vR11","vR12","vR13","vR20","vR21","vR22","vR23","vR30","vR31","vR32","vR33","vR40","vR41","vR42","vR43","vL","vT1","vT2","vT3","vT4"};

MatrixForm[stoich=(({{-1, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0}, {1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 1, -1, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 1, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1}}))];


MatrixForm[stoich,TableHeadings->{varnames,fluxnames}]
Dimensions[stoich,2]
MatrixRank[stoich]

Note that kap has the ATP conc in it as it is a constant and kfp has F6p in it -- thus the equations are linear in the concentrations of the enzyme states.  all of these can be solved simultaneously

kap = (k^+_ATP) ATP and kfp = (k^+_F6P) F6P

2. Solve for the enzyme forms analytically

fluxes={                              kap*R0-kam*R0A,kfp*R0A-kfm*R0AF,kpfk*R0AF,
4*kp*R0-1*km*R1,kap*R1-kam*R1A,kfp*R1A-kfm*R1AF,kpfk*R1AF,
3*kp*R1-2*km*R2,kap*R2-kam*R2A,kfp*R2A-kfm*R2AF,kpfk*R2AF,
2*kp*R2-3*km*R3,kap*R3-kam*R3A,kfp*R3A-kfm*R3AF,kpfk*R3AF,
1*kp*R3-4*km*R4,kap*R4-kam*R4A,kfp*R4A-kfm*R4AF,kpfk*R4AF,
T0-L*R0,
4*kip*T0-1*kim*T1,
3*kip*T1-2*kim*T2,
2*kip*T2-3*kim*T3,
1*kip*T3-4*kim*T4};
eqn1=R0+R0A+R0AF+R1+R1A+R1AF+R2+R2A+R2AF+R3+R3A+R3AF+R4+R4A+R4AF+T0+T1+T2+T3+T4-PFKtot;
eqn2=kpfk*(R0AF+R1AF+R2AF+R3AF+R4AF)-vpfk;
MatrixForm[stoich.fluxes,TableHeadings->{varnames}]

Solve for the concentrations using only flux as the constraint

solution=Solve[{stoich.fluxes==0, eqn2==0},{R0,R0A,R0AF,R1,R1A,R1AF,R2,R2A,R2AF,R3,R3A,R3AF,R4,R4A,R4AF,T0,T1,T2,T3,T4}];

solution ={{T4->(((kip^4) (km^4) (kfp kpfk+kam (kfm+kpfk)) L vpfk)/(kap kfp (kim^4) (((km+kp))^4) kpfk)),T3->((4 (kip^3) (km^4) (kfp kpfk+kam (kfm+kpfk)) L vpfk)/(kap kfp (kim^3) (((km+kp))^4) kpfk)),R0A->(((km^4) (kfm+kpfk) vpfk)/(kfp (((km+kp))^4) kpfk)),R1A->((4 (km^3) kp (kfm+kpfk) vpfk)/(kfp (((km+kp))^4) kpfk)),R0AF->(((km^4) vpfk)/((((km+kp))^4) kpfk)),R1AF->((4 (km^3) kp vpfk)/((((km+kp))^4) kpfk)),R2A->((6 (km^2) (kp^2) (kfm+kpfk) vpfk)/(kfp (((km+kp))^4) kpfk)),R2AF->((6 (km^2) (kp^2) vpfk)/((((km+kp))^4) kpfk)),R1->((4 (km^3) kp (kfp kpfk+kam (kfm+kpfk)) vpfk)/(kap kfp (((km+kp))^4) kpfk)),R0->(((km^4) (kfp kpfk+kam (kfm+kpfk)) vpfk)/(kap kfp (((km+kp))^4) kpfk)),T0->(((km^4) (kfp kpfk+kam (kfm+kpfk)) L vpfk)/(kap kfp (((km+kp))^4) kpfk)),T1->((4 kip (km^4) (kfp kpfk+kam (kfm+kpfk)) L vpfk)/(kap kfp kim (((km+kp))^4) kpfk)),T2->((6 (kip^2) (km^4) (kfp kpfk+kam (kfm+kpfk)) L vpfk)/(kap kfp (kim^2) (((km+kp))^4) kpfk)),R2->((6 (km^2) (kp^2) (kfp kpfk+kam (kfm+kpfk)) vpfk)/(kap kfp (((km+kp))^4) kpfk)),R3A->((4 km (kp^3) (kfm+kpfk) vpfk)/(kfp (((km+kp))^4) kpfk)),R3AF->((4 km (kp^3) vpfk)/((((km+kp))^4) kpfk)),R3->((4 km (kp^3) (kfp kpfk+kam (kfm+kpfk)) vpfk)/(kap kfp (((km+kp))^4) kpfk)),R4->(((kp^4) (kfp kpfk+kam (kfm+kpfk)) vpfk)/(kap kfp (((km+kp))^4) kpfk)),R4A->(((kp^4) (kfm+kpfk) vpfk)/(kfp (((km+kp))^4) kpfk)),R4AF->(((kp^4) vpfk)/((((km+kp))^4) kpfk))}};
enzymes=Simplify[solution];

enzymestates=Simplify[{R0,R0A,R0AF,R1,R1A,R1AF,R2,R2A,R2AF,R3,R3A,R3AF,R4,R4A,R4AF,T0,T1,T2,T3,T4}/.enzymes];
MatrixForm[Flatten[Simplify[enzymestates]],TableHeadings->{{R0,R0A,R0AF,R1,R1A,R1AF,R2,R2A,R2AF,R3,R3A,R3AF,R4,R4A,R4AF,T0,T1,T2,T3,T4},{}}]

Plugging in equilibrium constants

MatrixForm[Flatten[Simplify[enzymestates/.{kim->kip*Ki,km->kp*Ka,kam->kap*Katp,kfm->kfp*Kf6p}]],TableHeadings->{{R0,R0A,R0AF,R1,R1A,R1AF,R2,R2A,R2AF,R3,R3A,R3AF,R4,R4A,R4AF,T0,T1,T2,T3,T4},{}}]

Determine ratios between different enzyme complex forms, determine total enzyme

sol={{enzymestates[[1,1]],enzymestates[[1,1]]/enzymestates[[1,1]],
 enzymestates[[1,2]]/enzymestates[[1,1]],
 enzymestates[[1,3]]/enzymestates[[1,1]],1,1}, 

{enzymestates[[1,4]],enzymestates[[1,4]]/enzymestates[[1,4]],
 enzymestates[[1,5]]/enzymestates[[1,4]],
 enzymestates[[1,6]]/enzymestates[[1,4]],enzymestates[[1,4]]/enzymestates[[1,1]],enzymestates[[1,17]]/enzymestates[[1,16]]},  

{enzymestates[[1,7]],enzymestates[[1,7]]/enzymestates[[1,7]],
 enzymestates[[1,8]]/enzymestates[[1,7]],
 enzymestates[[1,9]]/enzymestates[[1,7]],enzymestates[[1,7]]/enzymestates[[1,1]],enzymestates[[1,18]]/enzymestates[[1,16]]},  

{enzymestates[[1,10]],enzymestates[[1,10]]/enzymestates[[1,10]],
 enzymestates[[1,11]]/enzymestates[[1,10]],
 enzymestates[[1,12]]/enzymestates[[1,10]],enzymestates[[1,10]]/enzymestates[[1,1]],enzymestates[[1,19]]/enzymestates[[1,16]]},

{enzymestates[[1,13]],enzymestates[[1,13]]/enzymestates[[1,13]],
 enzymestates[[1,14]]/enzymestates[[1,13]],
 enzymestates[[1,15]]/enzymestates[[1,13]],enzymestates[[1,13]]/enzymestates[[1,1]],enzymestates[[1,20]]/enzymestates[[1,16]]}};

MatrixForm[sol,TableHeadings->{{0,1,2,3,4},{"Ri0","Ri/Ri0","RiA/Ri0","RiAF/Ri0","Ri0/R10","Ti/T0"}}]
TeXForm[MatrixForm[sol,TableHeadings->{{0,1,2,3,4},{"Ri0","Ri/Ri0","RiA/Ri0","RiAF/Ri0","Ri0/R10","Ti/T0"}}]];

rtot=Simplify[Sum[enzymestates[[1,i]],{i,1,15}]]
ttot=Simplify[Sum[enzymestates[[1,i]],{i,16,20}]]

pfktot=Simplify[rtot+ttot]

sumRiAF=Simplify[enzymestates[[1,3]]+enzymestates[[1,6]]+enzymestates[[1,9]]+enzymestates[[1,12]]+enzymestates[[1,15]]]

3. Estimated Numerical Solution

The enzyme is inhibited between 10-15% under physiological conditions --> Ponce et al. Biochimica et Biophysica Acta 1971 250(1):63-74

Keq={kim->kip*Ki,km->kp*Ka,kam->kap*Katp,kfm->kfp*Kf6p};
estimates={Ki->0.1/1.6,Ka->0.033/0.08672812499999999`,Katp->0.068/1.6,Kf6p->0.1/0.0198,vpfk->1.12,L->0.0011};

x=MatrixForm[Flatten[Simplify[enzymestates/.Keq/.estimates]],TableHeadings->{{R0,R0A,R0AF,R1,R1A,R1AF,R2,R2A,R2AF,R3,R3A,R3AF,R4,R4A,R4AF,T0,T1,T2,T3,T4},{}}]

Sum[x[[1,i]],{i,1,20}]

Plugging in all the concentrations and flux constraints reduces the system to only three dimensions, this space can be sampled/FVA if need be
Here, I just systematically solve each PERC by making sure that the denominator stays greater than 0 for the solution to be positive
This provides us with bounds for which the PERCs must lie for positive values
I set the PERCs within these bounds and tried to keep them all near each other
Another interesting constraint would be to set a ratio of Ttot/PFKtot if it exists in literature, reducing the dimension to 2
	Also, the L value is only a small portion of calculating Ttot/PFKtot, as can be seen above, Ttot is also related to kap, kfp, and kpfk
		This can only be seen in the combined (R and T) version of the module
If kfp is large and kap is small, Ttot/PFKtot increases, tried to manually find a value that would give a large ratio (33%)

kpfk1=Solve[{((1.7138509118570582`)/kap)+((1.1928386637539248`)/kfp)+((7.144437695726893`)/kpfk)==33*10^-6},{kpfk,kap,kfp}]

kap1=Solve[-1.728623785100803`*^29 kap-2.483658134478425`*^29 kfp+4.78225485488576`*^24 kap kfp==0,{kap,kfp}]

kfp1=Solve[-1.728623785100803`*^29+4.78225485488576`*^24 kfp==0,{kfp}]

kap1/.kfp->1000000
kpfk1/.kfp->1000000/.kap->200000

Final concentrations calculated from these values as well as the final concentration

Keq={kim->kip*Ki,km->kp*Ka,kam->kap*Katp,kfm->kfp*Kf6p};
estimates={Ki->0.1/1.6,Ka->0.033/0.08672812499999999`,Katp->0.068/1.6,Kf6p->0.1/0.0198,vpfk->1.12,L->0.0011,kpfk->307447.5581772384`,kfp->1000000,kap->200000};

x=MatrixForm[Flatten[Simplify[enzymestates/.Keq/.estimates]],TableHeadings->{{R0,R0A,R0AF,R1,R1A,R1AF,R2,R2A,R2AF,R3,R3A,R3AF,R4,R4A,R4AF,T0,T1,T2,T3,T4},{}}]

Sum[x[[1,i]],{i,1,20}]

MatrixForm[{Sum[x[[1,i]],{i,1,3}],Sum[x[[1,i]],{i,4,6}],Sum[x[[1,i]],{i,7,9}],Sum[x[[1,i]],{i,10,12}],Sum[x[[1,i]],{i,13,15}],Sum[x[[1,i]],{i,16,20}]},TableHeadings->{{R0i,R1i,R2i,R3i,R4i,T},{}}]
Rfrac=Sum[x[[1,i]],{i,1,15}]/0.000033
Rcatfrac=(x[[1,3]]+x[[1,6]]+x[[1,9]]+x[[1,12]]+x[[1,15]])/Rfrac
MatrixForm[{Sum[x[[1,i]],{i,1,3}],Sum[x[[1,i]],{i,4,6}],Sum[x[[1,i]],{i,7,9}],Sum[x[[1,i]],{i,10,12}],Sum[x[[1,i]],{i,13,15}]}/Rfrac,TableHeadings->{{v0,v1,v2,v3,v4},{}}]

Check Flux through catalytic reaction

eqn2/.estimates/.R0AF->2.10241288178717`*^-8/.R1AF->2.210161541978762`*^-7/.R2AF->8.712871203769399`*^-7/.R3AF->1.5265676421604298`*^-6/.R4AF->1.0030026461382197`*^-6

T total/ PFK ratio

ttot/pfktot/.Keq/.estimates

fluxes/.{R0->enzymestates[[1,1]],R0A->enzymestates[[1,2]],R0AF->enzymestates[[1,3]],R1->enzymestates[[1,4]],R1A->enzymestates[[1,5]],R1AF->enzymestates[[1,6]],R2->enzymestates[[1,7]],R2A->enzymestates[[1,8]],R2AF->enzymestates[[1,9]],R3->enzymestates[[1,10]],R3A->enzymestates[[1,11]],R3AF->enzymestates[[1,12]],R4->enzymestates[[1,13]],R4A->enzymestates[[1,14]],R4AF->enzymestates[[1,15]],T0->enzymestates[[1,16]],T1->enzymestates[[1,17]],T2->enzymestates[[1,18]],T3->enzymestates[[1,19]],T4->enzymestates[[1,20]]}/.Keq/.estimates

4. examining T fraction

the last figure in the new notebook has the hue set to the ratio of ttot/pfktot.
the color is scaled so the greener it is, the larger the fraction and the more inhibited the system is. the green region has about a 33% inhibiton and the redder regions approach 0%
it seems that the lower the kap value is the more inhibited the system is.
aarash

ttot1=(((((kim+kip))^4) (km^4) (kfp z+kam (kfm+z)) L vpfk)/(kap kfp (kim^4) (((km+kp))^4) z));
pfktot1=(1/(kap kfp z))(kfp kpfk+kam (kfm+z)+kap (kfm+kfp+z)+(((((kim+kip))^4) (km^4) (kfp z+kam (kfm+z)) L)/((kim^4) (((km+kp))^4)))) vpfk;
Keq1={kim->kip*Ki,km->kp*Ka,kam->kap*Katp,kfm->kfp*Kf6p};
estimates1={Ki->0.1/1.6,Ka->0.033/0.08672812499999999`,Katp->0.068/1.6,Kf6p->0.1/0.0198,vpfk->1.12,L->0.0011};

ttot1/pfktot1/.Keq1/.estimates1

surface=Plot3D[Log[10,((7.144437695726893`)/((33/(10^6))-((1.7138509118570582`)/kap)-((1.1928386637539248`)/kfp)))],{kap,10^4,10^7},{kfp,10^4,10^7},ColorFunction->Function[{kap,kfp,z},Hue[(0.5302240284438019` (kfp z+0.0425` kap (5.05050505050505` kfp+z)))/(kfp z+0.0425` kap (5.05050505050505` kfp+z)+kap (6.05050505050505` kfp+z)+0.5302240284438019` (kfp z+0.0425` kap (5.05050505050505` kfp+z)))]],PlotRange->{{10^4,10^7},{10^4,10^7},{5.32,5.5}},PlotPoints->40,PlotStyle->Directive[Opacity[0.7]]];

point=ListPointPlot3D[{{2*^5,1*^6,Log[10,307447.5581772384`]}},PlotStyle->{PointSize[0.02],Hue[1]}];
Show[surface,point]

Density plot shows 2d version of the above 3d plot, as values increase in the z direction (kpfk), the lighter the plot becomes

DensityPlot[Log[10,((7.144437695726893`)/((33/(10^6))-((1.7138509118570582`)/kap)-((1.1928386637539248`)/kfp)))],{kap,10^4,10^7},{kfp,10^4,10^7},ColorFunction->"BlueGreenYellow",PlotPoints->40]