This Workbook sets up, simulates and graphs the  equations that describe a simple feedback loop. In this case the enzyme synthetic rate responds on a slow time scale (ie transcription/translation) to the buildup of the end product in ther reaction chain. It corresponds to Section 9.6.1 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
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=80;
k8=.1;
K7=10/0.05615528128088304;
k7=11*0.6403882032022076*k8;

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+k7/(1+K7*y5)-k8*y6==0,
 -k6*y5*y6+km6*y7==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];
v7=k7/(1.+K7*x5[t]);
v8=k8*x6[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+v7-v8,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,x5[t]/.solution,x6[t]+x7[t]/.solution},
{t,0,tfinal},PlotLabel->"Concentrations vs time",PlotLegend->{"x1","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->{.3,.4},ImageSize->600]

graph2=Plot[{k1*x1[t]*x6[t]/.solution,k7/(1+K7*x5[t])/.solution,k8*x6[t]/.solution},
{t,0,tfinal},PlotLabel->"Fluxes vs time",AxesLabel->{"Time"," "},PlotLegend->{"v1","v7","v8"},
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.5},
LegendPosition->{1.,-.4},LegendSize->{.3,.4},ImageSize->600]

graph3 =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.45}},AspectRatio->.7,ImageSize->480]

graph4=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,.8},{0,.9}},AspectRatio->.7,ImageSize->500]

(x6[tfinal]+x7[tfinal])/.solution;
x5[tfinal]/.solution;

This cell will export the graphs

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