This Workbook sets up, simulates and graphs the  equations that describe the dynamic behavior of simple connected reversible bi-linear reactions in terms of the pools. It corresponds to Section 4.5 in Ch 4 of SB2.

We will procedure  as follows:

1) define the numerical values of the kinetic parameters.
2)  then define the reaction rate laws.
3) then we integrate the set of  differential equations.
 4) finally graph the solution.

Notation: 

 Kinetic parameters: k1, km1, k2, km2.
 Equilibrium constants: K1, K2.
 Fluxes:  v1, v1m, v2, v2m.

BOP May 21, 2008
Addiel de Alba. April 2010

k1=1;km1=1;k2=1;km2=1;K1=(k1/km1);K2=(k2/km2);
v1=k1 x1[t] x2[t];
vm1=km1 x3[t];
v2=k2 x3[t];
vm2=km2 x4[t] x5[t];
tfinal=5;
x10=0.75;
x20=0.75;
x30=1.5;
x40=0.75;
x50=0.75;

BalanceX1=-v1+vm1;
BalanceX3=-v2+vm2;
solution=NDSolve[{x1'[t]==BalanceX1,x2'[t]==BalanceX1,x3'[t]==-BalanceX1+BalanceX3,x4'[t]==-BalanceX3,x5'[t]==-BalanceX3,x1[0]==x10,x2[0]==x20,x3[0]==x30,x4[0]==x40,x5[0]==x50},{x1,x2,x3,x4,x5},{t,0,tfinal}];

BalanceX1=10(-v1+vm1);
BalanceX3=-v2+vm2;
solution1=NDSolve[{x1'[t]==BalanceX1,x2'[t]==BalanceX1,x3'[t]==-BalanceX1+BalanceX3,x4'[t]==-BalanceX3,x5'[t]==-BalanceX3,x1[0]==x10,x2[0]==x20,x3[0]==x30,x4[0]==x40,x5[0]==x50},{x1,x2,x3,x4,x5},{t,0,tfinal}];

BalanceX1=0.10(-v1+vm1);
BalanceX3=-v2+vm2;
solution2=NDSolve[{x1'[t]==BalanceX1,x2'[t]==BalanceX1,x3'[t]==-BalanceX1+BalanceX3,x4'[t]==-BalanceX3,x5'[t]==-BalanceX3,x1[0]==x10,x2[0]==x20,x3[0]==x30,x4[0]==x40,x5[0]==x50},{x1,x2,x3,x4,x5},{t,0,tfinal}];

graph1=ParametricPlot[Evaluate[{{x1[t] x2[t]-((x3[t])/K1),x3[t]-((x4[t] x5[t])/K2)}/.solution,{x1[t] x2[t]-((x3[t])/K1),x3[t]-((x4[t] x5[t])/K2)}/.solution1,{x1[t] x2[t]-((x3[t])/K1),x3[t]-((x4[t] x5[t])/K2)}/.solution2}],{t,0,tfinal},PlotStyle->{Directive[Dotted,Black],Directive[Black],Directive[Dashed,Black]}, PlotRange->All,AspectRatio->1/2,ImageSize->500]

graph2=Plot[{x1[t] x2[t]-((x3[t])/K1)/.solution2,x3[t]-((x4[t] x5[t])/K2)/.solution2},{t,0,tfinal},PlotStyle->{GrayLevel[0.0],GrayLevel[0.0]}, PlotRange->All,AspectRatio->1/2,ImageSize->500]

graph3=Plot[{x1[t] x2[t]-((x3[t])/K1)/.solution,x3[t]-((x4[t] x5[t])/K2)/.solution},{t,0,1},PlotStyle->{GrayLevel[0.0],GrayLevel[0.0]}, PlotRange->All,AspectRatio->1/2,ImageSize->500]

graph4=Plot[{x1[t] x2[t]-((x3[t])/K1)/.solution1,x3[t]-((x4[t] x5[t])/K2)/.solution1},{t,0,1},PlotStyle->{GrayLevel[0.0],GrayLevel[0.0]}, PlotRange->All,AspectRatio->1/2,ImageSize->500]

Options[Plot];
graph5=Plot[{x1[t] x2[t]-((x3[t])/K1)/.solution2,x1[t]-1+x2[t]-1-x3[t]+1/.solution2},{t,0,tfinal},PlotStyle->{Directive[Black],Directive[Dashed,Black]}, PlotRange->All,AspectRatio->1/2,ImageSize->500]

graph6=ParametricPlot[Evaluate[{{x1[t] x2[t]-((x3[t])/K1),x3[t]-((x4[t] x5[t])/K2)}/.solution,{x1[t]-1+x2[t]-1-x3[t]+1,x3[t]-1-x4[t]+1-x5[t]+1}/.solution,{x1[t] x2[t]-((x3[t])/K1),x3[t]-((x4[t] x5[t])/K2)}/.solution1,{x1[t]-1+x2[t]-1-x3[t]+1,x3[t]-1-x4[t]+1-x5[t]+1}/.solution1,{x1[t] x2[t]-((x3[t])/K1),x3[t]-((x4[t] x5[t])/K2)}/.solution2,{x1[t]-1+x2[t]-1-x3[t]+1,x3[t]-1-x4[t]+1-x5[t]+1}/.solution2}],{t,0,tfinal},PlotStyle->{Directive[Black],Directive[Dashed,Black]}, PlotRange->All,AspectRatio->1/2,ImageSize->500]

Export["Section4.5C.eps",graph1];
Export["Section4.5D.eps",graph2];
Export["Section4.5E.eps",graph3];
Export["Section4.5F.eps",graph4];
Export["Section4.5G.eps",graph5];
Export["Section4.5H.eps",graph6];

Clear[K1, K2,k1,km1,k2,km2];
MatrixForm[stoich=(({{-1, 0}, {-1, 0}, {1, -1}, {0, 1}, {0, 1}}))];
MatrixForm[v=(({{k1.p1}, {k2.p2}}))];
MatrixRank[stoich];
leftnull=MatrixForm[(({{1, 0, 0}, {0, 1, 1}, {1, 0, 1}})).NullSpace[Transpose[stoich]]];

MatrixForm[pools=(({{x2, x1, -1/K1, 0, 0}, {0, 0, 1, -x5/K2, -x4/K2}, {1, 0, 1, 0, 1}, {0, 1, 1, 1, 0}, {0, 1, 1, 0, 1}}))];
MatrixForm[pools.stoich];
MatrixForm[pools.stoich.v];
TeXForm[MatrixForm[pools.stoich.v]];