This Workbook sets up, simulates and graphs the  equations that describe a coupled set of simple reversible reactions.  It  corresponds to the Ch4.4 in SB2.

Notation:

Kinetic  parameters: k1, km1, k2, k3, km3.
Equilibrium constants: K1, K3.
Fluxes:  v1, v1m, v2,  v3 and v3m.
Exhanges fluxes: b1 and b2.

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

Step 1: define the equations

stoich=(({{-1, 1, 0, 0, 0}, {1, -1, -1, 0, 0}, {0, 0, 1, -1, 1}, {0, 0, 0, 1, -1}}));
v ={v1,vm1,v2,v3,vm3};
dxdt=stoich.v

step 2 : define the parmeter values and rate laws

k1=1;km1=1;k2=.2;k3=1;km3=1;
K1=k1/km1;K3=k3/km3;
v1=k1*x1[t];
vm1=km1*x2[t];
v2=k2*x2[t];
v3=k3*x3[t];
vm3=km3*x4[t];

step 3:  set up the conditions and simulate

tfinal=25;
solution=NDSolve[{x1'[t]==dxdt[[1]],x2'[t]==dxdt[[2]],x3'[t]==dxdt[[3]],x4'[t]==dxdt[[4]],x1[0]==1.0,x2[0]==0,x3[0]==0,x4[0]==0},{x1,x2,x3,x4},{t,0,tfinal}];

step 4:  graph the results

<<PlotLegends`
graph1=Plot[{x1[t]/.solution,x2[t]/.solution,x3[t]/.solution,x4[t]/.solution},
{t,0,tfinal},PlotLabel->"Concentrations vs time",PlotLegend->{"x1","x2","x3","x4"},
PlotStyle->{Hue[0.5],Hue[0.6],Hue[0.8],Hue[1.0]},PlotRange->{0,1.},AxesLabel->{"Time"," "},
LegendPosition->{1.1,-.4},LegendSize->{.35,.35}]
graph2=Plot[{x1[t]-x2[t]/K1/.solution,x3[t]-x4[t]/K3/.solution,x1[t]+x2[t]/.solution,x3[t]+x4[t]/.solution},{t,0,tfinal},PlotLabel->"Pools versus time",PlotLegend->{"P1","P2","P3","P4"},PlotStyle->{Hue[0.5],Hue[0.6],Hue[0.8],Hue[1.0]},AxesLabel->{"Time"," "},LegendPosition->{1.1,-.4},LegendSize->{.35,.35}];
graph3=ParametricPlot[Evaluate[{x1[t]-x2[t]/K1,x3[t]-x4[t]/K3}/.solution],{t,0,tfinal},PlotLabel->"The two dis-equilibrium quantities",PlotStyle->{Hue[.7]}];
graph4=ParametricPlot[Evaluate[{x1[t],x2[t]}/.solution],{t,0,tfinal},PlotLabel->"x1 vs.x2",AxesLabel->{"x1","x2"}, PlotStyle->{Hue[.7]}];
graph5=ParametricPlot[Evaluate[{x3[t],x4[t]}/.solution],{t,0,tfinal},PlotLabel->"x3 vs.x4",AxesLabel->{"x3","x4"}, PlotStyle->{Hue[.7]}];
x1[tfinal]/.solution;
x2[tfinal]/.solution;
x3[tfinal]/.solution;

graph6=Plot[{x1[t]/.solution,x2[t]/.solution,x3[t]/.solution,x4[t]/.solution},
{t,0,10},
PlotStyle->{Hue[0.5],Hue[0.6],Hue[0.8],Hue[1.0]},PlotRange->{0,1.},PlotLegend->{"x1","x2","x3","x4"},LegendPosition->{1.1,-.4},LegendSize->{.35,.35}]
graph7=Plot[{x1[t]-x2[t]/K1/.solution,x3[t]-x4[t]/K3/.solution,x1[t]+x2[t]/.solution,x3[t]+x4[t]/.solution},{t,0,tfinal},PlotStyle->{Hue[0.5],Hue[0.6],Hue[0.8],Hue[1.0]},PlotLegend->{"P1","P2","P3","P4"},LegendPosition->{1.1,-.4},LegendSize->{.35,.35}]
graph8=Plot[{x1[t]-x2[t]/K1/.solution,x3[t]-x4[t]/K3/.solution,x2[t]/.solution,x1[t]+x2[t]+x3[t]+x4[t]/.solution},{t,0,tfinal},PlotStyle->{Hue[0.5],Hue[0.6],Hue[0.8],Hue[1.0]},PlotLegend->{"P1","P2","x2","xtot"},LegendPosition->{1.1,-.4},LegendSize->{.35,.35}]

step 5: Export the graphs

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

Appendix: Pool formation

Clear[K1,k1,km1,k2,K3,k3,km3]
K1=k1/km1; K3=k3/km3;
fluxes={k1*x1,km1*x2,k2*x2,k3*x3,km3*x4}

MatrixForm[stoich=(({{-1, 1, 0, 0, 0}, {1, -1, -1, 0, 0}, {0, 0, 1, -1, 1}, {0, 0, 0, 1, -1}}))]
MatrixRank[stoich]
leftnull=MatrixForm[NullSpace[Transpose[stoich]]]
MatrixForm[pools=(({{1, -1/K1, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, -1/K3}, {1, 1, 1, 1}}))]
MatrixForm[p=pools.stoich]
MatrixForm[p.fluxes]
TeXForm[MatrixForm[p.fluxes]];
Null

Appendix: Gradient/Jacobian Work

MatrixForm[grad=(({{k1, 0, 0, 0}, {0, km1, 0, 0}, {0, k2, 0, 0}, {0, 0, k3, 0}, {0, 0, 0, km3}}))];
MatrixForm[jac=stoich.grad]
MatrixRank[jac];
NullSpace[Transpose[jac]];
NullSpace[jac];