This Workbook sets up, simulates and graphs the  equations that describe a simple reversible reaction in an open system. It  corresponds to Section 6.2 of 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: k1, km1, k2.
Equilibrium constant: K1.
Fluxes: v1, v1m.
 Exhanges fluxes: b1 and b2.

BOP April 7, 2007
Addiel de Alba. April 2010

k1=1.`;km1=2.`;k2=0.1`;K1=(k1/km1);
v1=k1 x1[t];
vm1=km1 x2[t];
b1=0.02`;
b2=k2 x2[t];
tfinal=150;
x10=0.21`;
x20=0.1`;
solution=NDSolve[{x1'[t]==b1-v1+vm1,x2'[t]==v1-vm1-b2,x1[0]==x10,x2[0]==x20},{x1,x2},{t,0,tfinal}];
<<"PlotLegends`"

plot solution with time

graph1=Plot[{x1[t]/.solution,x2[t]/.solution},{t,0,tfinal},PlotLabel->"Concentrations vs time",PlotLegend->{"x1","x2"},AxesLabel->{"Time"," "},PlotStyle->{Hue[0.5`],Hue[0.6`]},PlotRange->{0,1.2`},LegendPosition->{1.1`,-0.4`},LegendSize->{0.35`,0.35`},ImageSize->600];
graph2=Plot[{x1[t]-((x2[t])/K1)/.solution,{x1[t]+x2[t]}/.solution},{t,0,tfinal},PlotLabel->"Pools versus time",PlotLegend->{"p1=x1-x2/K1","p2=x1+x2"},PlotStyle->{Hue[0.5`],Hue[0.6`]},AxesLabel->{"Time"," "},LegendPosition->{1.1`,-0.4`},LegendSize->{.5`,0.35`},ImageSize->600]

plot dynamic phase portraits

graph3=ParametricPlot[Evaluate[{x1[t],x2[t]}/.solution],{t,0,tfinal},PlotLabel->"Dynamic Phase portrait, x1 vs x2",PlotStyle->{Hue[0.5`]},AxesLabel->{"x1","x2"},PlotRange->{{0,1},{0,0.5`}}];
graph4=ParametricPlot[Evaluate[{x1[t]+x2[t],x1[t]-((x2[t])/K1)}/.solution],{t,0,tfinal},PlotLabel->"Dynamic Phase portrait, p2 vs p1",AxesLabel->{"p2","p1"},PlotStyle->{Hue[0.5`]},PlotRange->{{0.3,0.65`},{0,.02`}},AspectRatio->1/GoldenRatio]
x1f=x1[tfinal]/.solution;
x2f=x2[tfinal]/.solution;
Dev=1-(x2f/(x1f K1));
Options[ParametricPlot];

graph1=Plot[{x1[t]/.solution,x2[t]/.solution},{t,0,tfinal},PlotStyle->{Hue[0.5`],Hue[0.6`]},PlotRange->{0,0.5`}];
graph2=Plot[{x1[t]-((x2[t])/K1)/.solution,{x1[t]+x2[t]}/.solution},{t,0,tfinal},PlotStyle->{Hue[0.5`],Hue[0.6`]}];
graph3=ParametricPlot[Evaluate[{x1[t],x2[t]}/.solution],{t,0,tfinal},PlotRange->{{0,1},{0,0.5`}}];
graph4=ParametricPlot[Evaluate[{x1[t]+x2[t],x1[t]-((x2[t])/K1)}/.solution],{t,0,tfinal},PlotStyle->{Hue[0.5`]},PlotRange->{{0.3,0.65`},{0,0.02`}},AspectRatio->.7,ImageSize->400];

Here we export the graphs

Export["Section6.2.A.eps",graph2];
Export["Section6.2.B.eps",graph4];

Clear[b1,k1,km1,k2];
Solve[{b1-k1*x1+km1*x2==0, k1*x1-(km1+k2)x2==0},{x1,x2},{x1}];
TeXForm[Solve[{b1-k1*x1+km1*x2==0, k1*x1-(km1+k2)x2==0},{x1,x2},{x1}]];