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
NS, checked on March 15, 2011 OK

SetDirectory[NotebookDirectory[]];
<<"GraphicsDirectives.m"

tfinal=80;
parameters={k0->.5,k1->1,k2->1,k3->1,k4->1,k5->1,k6->10,km6->1,b1->0.1,k8->.1,K7->10/0.05615528128088304,k7->11*0.6403882032022076*k8,b1->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+k7/(1+K7*y5)-k8*y6==0,
 -k6*y5*y6+km6*y7==0
}//.parameters,{y1,y2,y3,y4,y5, y6, y7}];
{z1,z2,z3,z4,z5,z6,z7}={y1,y2,y3,y4,y5, y6, y7}/.ststsolution[[2]]

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];

parameters=parameters/.(b1->.1)->(b1->1.);
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}//.parameters,{x1,x2,x3,x4,x5, x6, x7},{t,0,tfinal}];

graph1=Plot[Evaluate[{x1[t],x5[t],x6[t]+x7[t]}/.solution],
{t,0,tfinal},PlotLabel->"Concentrations vs time",AxesLabel->{"Time"," "},Epilog->createLegend[{"x1","x2","x6+x7"}]]

graph2=Plot[Evaluate[{k1*x1[t]*x6[t],k7/(1+K7*x5[t]),k8*x6[t]}/.solution//.parameters],
{t,0,tfinal},PlotLabel->"Fluxes vs time",AxesLabel->{"Time"," "},Epilog->createLegend[{"v1","v7","v8"}],
PlotRange->{0,0.5}]

graph3 =Show[ParametricPlot[Evaluate[{k1*x1[t]*x6[t],k5*x5[t]}/.solution//.parameters],{t,0,tfinal},PlotLabel->"Dynamic Phase Portrait for v1 vs. v5",AxesLabel->{"v1","v5"},PlotRange->{{0,.6},{0,0.45}},Epilog->ppAnnotate[{k1*x1[#]*x6[#]&,k5*x5[#]&}//.parameters,solution,tfinal],FrameLabel->{"v1","v5"}],Graphics[{Black,Dashed,Line[{{0,0},{1.,1.}}]}]]

graph4=Show[ParametricPlot[Evaluate[{x6[t],x7[t]}/.solution],{t,0,tfinal},AxesLabel->{"x6","x7"},PlotLabel->"Dynamic Phase Portrait for x6 and x7",PlotRange->{{0,.8},{0,.9}},Epilog->ppAnnotate[{x6,x7},solution,tfinal],FrameLabel->{"x6","x7"}],Graphics[{Black,Dashed,Line[{{0,1},{1,0}}],Line[{{0,x6[tfinal]+x7[tfinal]},{x6[tfinal]+x7[tfinal],0}}]/.solution}]]

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

Manipulate[
parameters=parameters/.(b1->_)->(b1->b1param);
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}//.parameters,{x1,x2,x3,x4,x5, x6, x7},{t,0,tfinal}];

Column[{
Show[ParametricPlot[Evaluate[{k1*x1[t]*x6[t],k5*x5[t]}/.solution//.parameters],{t,0,tfinal},PlotLabel->"Dynamic Phase Portrait for v1 vs. v5",AxesLabel->{"v1","v5"},PlotRange->{{0,1.},{0,1.}},Epilog->ppAnnotate[{k1*x1[#]*x6[#]&,k5*x5[#]&}//.parameters,solution,tfinal],FrameLabel->{"v1","v5"}],Graphics[{Black,Dashed,Line[{{0,0},{1.,1.}}]}]],
Show[ParametricPlot[Evaluate[{x6[t],x7[t]}/.solution],{t,0,tfinal},AxesLabel->{"x6","x7"},PlotLabel->"Dynamic Phase Portrait for x6 and x7",PlotRange->{{0,.8},{0,.9}},Epilog->ppAnnotate[{x6,x7},solution,tfinal],FrameLabel->{"x6","x7"}],Graphics[{Black,Dashed,Line[{{0,1},{1,0}}],Line[{{0,x6[tfinal]+x7[tfinal]},{x6[tfinal]+x7[tfinal],0}}]/.solution}],Graphics[Text["x6+x7="<>ToString[x6[tfinal]+x7[tfinal]/.solution[[1]]],{x6[tfinal],x7[tfinal]}+{.1,0}/.solution[[1]]]]]
}
]
,{b1param,.1,2.,.1},{{tfinal,400},1,400,10}]

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];