This Workbook sets up, simulates and graphs the  equations that describe feedback inhibition of pathways. It corresponds to Section 9.5 in SB2.

We 1) define the numerical values of the kinetic  parameters, 2) then the reaction rate laws, 3) then we integrate and, 4)  finally graph the solution.

Notation:

Kinetic parameters: k0, k1, k2,  k3, k4, k5, k6, km6.
External flux: b1.
Chemical Fluxes: v0, v1, v2, v3,  v4, v5, v6

BOP April 7, 2007
BOP to version 6.0, April 27 2009
Addiel de Alba. April 2010

k0=.5;k1=1;k2=1;k3=1;k4=1;k5=1;k6=10;km6=1;b1=0.1;tfinal=15;etot=1.0;

ststsolution=NSolve[{b1-k1*y1*y6-k0*y1==0,k1*y1*y6-k2*y2==0,k2*y2-k3*y3==0,k3*y3-k4*y4==0,k4*y4-k5*y5-(k6*y5*y6-km6*y7)==0,k6*y5*y6-km6*y7==0, y6+y7-etot==0},{y1,y2,y3,y4,y5, y6, y7}];
{z1,z2,z3,z4,z5,z6,z7}={y1,y2,y3,y4,y5, y6, y7}/.ststsolution[[2]];

b1=1.0;
v0=k0*x1[t];
v1=k1*x1[t]*x6[t];
v2=k2*x2[t];
v3=k3*x3[t];
v4=k4*x4[t];
v5=k5*x5[t];
v6=k6*x5[t]*x6[t]-km6*x7[t];

solution=NDSolve[{x1'[t]==b1-v1-v0,x2'[t]==v1-v2,x3'[t]==v2-v3,x4'[t]==v3-v4,x5'[t]==v4-v5-v6,x6'[t]==-v6,x7'[t]==v6,x1[0]==z1,x2[0]==z2,x3[0]==z3,x4[0]==z4,x5[0]==z5,x6[0]==z6,x7[0]==z7},{x1,x2,x3,x4,x5, x6, x7},{t,0,tfinal}];

<<"PlotLegends`"

graph1=Plot[{x1[t]/.solution,x2[t]/.solution,x3[t]/.solution,x4[t]/.solution,x5[t]/.solution,x6[t]/.solution,x7[t]/.solution},
{t,0,tfinal},PlotStyle->{Hue[0.5],Hue[0.6],Hue[0.7],Hue[0.8],Hue[0.9],Hue[1.0],Hue[1.1]},PlotRange->{0,1.5}]

graph2=Plot[{k0*x1[t]/.solution,k1*x1[t]*x6[t]/.solution,k2*x2[t]/.solution,k3*x3[t]/.solution,k4*x4[t]/.solution,k5*x5[t]/.solution,k6*x5[t]x6[t]-km6*x7[t]/.solution},
{t,0,tfinal},PlotStyle->{Hue[0.4],Hue[0.5],Hue[0.6],Hue[0.7],Hue[0.8],Hue[0.9],Hue[1.0]},PlotRange->{0,0.8}]

graph3 =ParametricPlot[{Evaluate[{k1*x1[t]*x6[t],k5*x5[t]}/.solution]},{t,0,tfinal},PlotStyle->{Hue[.7]},PlotRange->{{0,0.6},{0,0.6}},AspectRatio->.7,ImageSize->400]

graph4=ParametricPlot[Evaluate[{x6[t],x7[t]}/.solution],{t,0,tfinal},PlotStyle->{Hue[.7]},PlotRange->{{0,1},{0,1}},AspectRatio->.7,ImageSize->400]

graph5=Plot[{x1[t]/.solution,x2[t]/.solution,x3[t]/.solution,x4[t]/.solution,x5[t]/.solution,x6[t]/.solution,x7[t]/.solution},
{t,0,tfinal},PlotLabel->"Concentrations vs time",PlotLegend->{"x1","x2","x3","x4","x5","x6","x7"},
PlotStyle->{Hue[0.5],Hue[0.6],Hue[0.7],Hue[0.8],Hue[0.9],Hue[1.0],Hue[1.1]},PlotRange->{0,1.5},AxesLabel->{"Time"," "},
LegendPosition->{1.,-.4},LegendSize->{.4,.7}];

graph6=Plot[{k0*x1[t]/.solution,k1*x1[t]*x6[t]/.solution,k2*x2[t]/.solution,k3*x3[t]/.solution,k4*x4[t]/.solution,k5*x5[t]/.solution,k6*x5[t]x6[t]-km6*x7[t]/.solution},
{t,0,tfinal},PlotLabel->"Fluxes vs time",AxesLabel->{"Time"," "},PlotLegend->{"v0","v1","v2","v3","v4","v5","v6"},
PlotStyle->{Hue[0.4],Hue[0.5],Hue[0.6],Hue[0.7],Hue[0.8],Hue[0.9],Hue[1.0]},PlotRange->{0,0.8},
LegendPosition->{1.,-.4},LegendSize->{.4,.7}];

graph7 =ParametricPlot[{Evaluate[{k1*x1[t]*x6[t],k5*x5[t]}/.solution]},{t,0,tfinal},PlotLabel->"Dynamic Phase Portrait for v1 vs. v5",PlotStyle->{Hue[.7]},AxesLabel->{"v1","v5"},PlotRange->{{0,0.6},{0,0.6}}];

graph8=Plot[{x7[t]/(x6[t]+x7[t])/.solution},
{t,0,tfinal},PlotLabel->"inhibition fraction vs time",AxesLabel->{"Time"," fraction "},
PlotStyle->{Hue[0.4]},PlotRange->{0.0,1.0}];

graph9=ParametricPlot[Evaluate[{x6[t],x7[t]}/.solution],{t,0,tfinal},AxesLabel->{"x6","x7"},PlotLabel->"Dynamic Phase Portrait for x6 and x7", PlotStyle->{Hue[.7]},PlotRange->{{0,1},{0,1}}];
x6[tfinal]/.solution;
Null;

graph8=Plot[{k6*x5[t]x6[t]-km6*x7[t]/.solution},
{t,0,tfinal},PlotStyle->{Hue[0.4],Hue[0.5],Hue[0.6],Hue[0.7],Hue[0.8],Hue[0.9],Hue[1.0]},PlotRange->{-0.2,0.2}];

This cell will export the graphs

Export["Section9.5.A.eps",graph1];
Export["Section9.5.B.eps",graph2];
Export["Section9.5.A2.eps",graph3];
Export["Section9.5.B2.eps",graph4];