This Workbook sets up, simulates and graphs the  equations that describe a simple reversible bi-linear reaction.  It  corresponds to the Ch 4.3 in SB2.

Notation: 

Kinetic parameters: k1,  km1.
Equilibrium constant: K1
Fluxes:  v1, v1m.


BOP April 7, 2007
AB - Last updated: Mar 26, 2010.

Step 1: define the equations

dxdt={-v1+vm1,-v1+vm1,v1-vm1};

step 2 : define the parmeter values and rate laws

k1=1;km1=1.;K1=k1/km1;
v1=k1*x1[t]*x2[t];
vm1=km1*x3[t];

step 3:  set up the conditions and simulate

tfinal=2;
solution=NDSolve[{x1'[t]==dxdt[[1]],x2'[t]==dxdt[[2]],x3'[t]==dxdt[[3]],x1[0]==3.0,x2[0]==2.0,x3[0]==0},{x1,x2,x3},{t,0,tfinal}];

step 4:  graph the results

<<"PlotLegends`"

graph1=Plot[{x1[t]/.solution,x2[t]/.solution,x3[t]/.solution},{t,0,tfinal},PlotLabel->"Concentrations vs time",PlotLegend->{"x1","x2","x3"},PlotStyle->{Hue[0.5`],Hue[0.6`],Hue[0.7`]},PlotRange->{0,3.2`},AxesLabel->{"Time"," "},LegendPosition->{.6`,0.2`},LegendSize->{0.35`,0.30`},ImageSize->{600}]

graph2=Plot[{x1[t] x2[t]-((x3[t])/K1)/.solution,{0.73` (x1[t]-1.73`)+1.73` (x2[t]-0.73`)-(x3[t]-1.27`)}/.solution,{x1[t]+x3[t]}/.solution,{x2[t]+x3[t]}/.solution},{t,0,tfinal},PlotLabel->"Pools versus time",PlotLegend->{"p1","p1-linear","p2","p3"},AxesLabel->{"Time"," "},PlotStyle->{Hue[0.5`],Hue[0.6`],Hue[0.7`],Hue[0.8`]},LegendPosition->{.6`,0.2`},LegendSize->{.4`,0.2`},ImageSize->{600}]

graph3=ParametricPlot[Evaluate[{x1[t],x2[t]}/.solution],{t,0,tfinal},PlotLabel->"Phase portrait, x1 vs x2",AxesLabel->{"x1","x2"},PlotStyle->{Hue[0.5`]},PlotRange->{{0,3},{0,3}},ImageSize->{500},AspectRatio->.7]

graph4=ParametricPlot[Evaluate[{x1[t],x3[t]}/.solution],{t,0,tfinal},PlotLabel->"Phase portrait, x1 vs x3",PlotStyle->{Hue[0.5`]},AxesLabel->{"x1","x3"},PlotRange->{{0,3},{0,3}},ImageSize->{500},AspectRatio->.7]

graph5=ParametricPlot[Evaluate[{x2[t],x3[t]}/.solution],{t,0,tfinal},PlotLabel->"Phase portrait, x2 vs x3",PlotStyle->{Hue[0.5`]},AxesLabel->{"x2","x3"},PlotRange->{{0,3},{0,3}},ImageSize->{500},AspectRatio->.7]
graph6=Plot[{x1[t] x2[t]-((x3[t])/K1)/.solution,{0.73` (x1[t]-1.73`)+1.73` (x2[t]-0.73`)-(x3[t]-1.27`)}/.solution},{t,0,tfinal},PlotLabel->"Pool versus time",PlotLegend->{"p1","p1-linear"},AxesLabel->{"Time"," "},PlotStyle->{Hue[0.5`],Hue[0.6`]},LegendPosition->{.6`,0.2`},LegendSize->{0.35`,0.30`},ImageSize->{500},AspectRatio->.7]
x1[tfinal]/.solution;
x2[tfinal]/.solution;

Step 5:  Export the graphs

(* Export["Section4.3A.eps",graph1];
Export["Section4.3B.eps",graph2]; *)

Appendix: Pool formation

Clear[K1,k1,km1]
K1=k1/km1

MatrixForm[stoich=(({{-1, 1}, {-1, 1}, {1, -1}}))]
MatrixRank[stoich]
leftnull=MatrixForm[(({{1, 0}, {1, 1}})).NullSpace[Transpose[stoich]]]
MatrixForm[pools=(({{.73, 1.73, -1}, {1, 0, 1}, {0, 1, 1}}))]
TeXForm[MatrixForm[pools=(({{.73, 1.73, -1}, {1, 0, 1}, {0, 1, 1}}))]];
TeXForm[MatrixForm[pools.stoich]];
Null

Solve[{x1*x2-x3/K==0, x1+x3==a, x2+x3==b},{x1,x2,x3},{x2,x3}]