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
NS, March 23, 2011

SetDirectory[NotebookDirectory[]];
<<"GraphicsDirectives.m"
SetOptions[Plot,AspectRatio->1/2];

Step 1: define the parameters

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];
x10=0.75;
x20=0.75;
x30=1.5;
x40=0.75;
x50=0.75;

step 2:  set up the conditions and simulate

tfinal=5;

A: k1 = 10 k2

BalanceX1=10(-v1+vm1);
BalanceX3=-v2+vm2;
solutionA=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}];

B: k1 =  k2

BalanceX1=-v1+vm1;
BalanceX3=-v2+vm2;
solutionB=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}];

C: k1 = 0.1 k2

BalanceX1=0.10(-v1+vm1);
BalanceX3=-v2+vm2;
solutionC=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}];

step 4:  graph the results

Pools

pool1=x1[t] x2[t]-((x3[t])/K1);
pool2=x3[t]-((x4[t] x5[t])/K2);
poolLabels={"p1 = x1[t] x2[t]-((x3[t])/K1)","p2 = x3[t]-((x4[t] x5[t])/K2)"};

SetOptions[Plot,AxesLabel->{"Time",""}];
legend=createLegend[poolLabels];
panelA=Plot[Evaluate[{pool1,pool2}/.solutionA],{t,0,1},PlotLabel->"k1 = 10 k2",Epilog->legend];
panelB=Plot[Evaluate[{pool1,pool2}/.solutionB],{t,0,1},PlotLabel->"k1 = k2",Epilog->legend];
panelC=Plot[Evaluate[{pool1,pool2}/.solutionC],{t,0,tfinal},PlotLabel->"k1 = .1 k2",Epilog->legend];
SlideView[{panelA,panelB,panelC}]

poolFunctions=Evaluate[{pool1,pool2}/.#&/@{solutionA,solutionB,solutionC}];
tooltipped=Thread[Tooltip[poolFunctions,{panelA,panelB,panelC},LabelStyle->{Background->White}]];
annotation=ppAnnotate[{x1[#] x2[#]-((x3[#])/K1)&,x3[#]-((x4[#] x5[#])/K2)&},#,tfinal]&/@{solutionA,solutionB,solutionC};
graph1=ParametricPlot[tooltipped,{t,0,tfinal},PlotStyle->{Directive[Black],Directive[Dotted,Black],Directive[Dashed,Black]},Epilog->annotation,FrameLabel->poolLabels,AspectRatio->1/2]

Linearized pools

{linPool1,linPool2}={x1[t]-1+x2[t]-1-x3[t]+1,x3[t]-1-x4[t]+1-x5[t]+1};

legend=createLegend[{"full solution","linearized solution"},{Directive[Black],Directive[Dashed,Black]},-10,-60];
insets=MapThread[Plot[Evaluate[{pool1,linPool1}/.#1],{t,0,tfinal},PlotStyle->{Directive[Black],Directive[Dashed,Black]},AspectRatio->1/2,ImageSize->500,Epilog->legend,AxesLabel->{"t","p1"},PlotLabel->#2]&,{{solutionA,solutionB,solutionC},{"k1 = 10 k2","k1 = k2","k1 = 0.1 k2"}}];
SlideView[insets]

tooltipped=Flatten[{Tooltip[{pool1,pool2}/.#[[2]],#[[1]],LabelStyle->{Background->White}],Tooltip[{linPool1,linPool2}/.#[[2]],#[[1]],LabelStyle->{Background->White}]}&/@Thread[List[insets,{solutionA,solutionB,solutionC}]],
{1,2}];

legend=createLegend[{"full solution","linearized solution"},{Directive[Black],Directive[Dashed,Black]},0,-10];
graph6=ParametricPlot[tooltipped,{t,0,tfinal},PlotStyle->{Directive[Black],Directive[Dashed,Black]}, PlotRange->All,AspectRatio->1/2,ImageSize->600,FrameLabel->{"p1","p2"},Epilog->legend]

step 5: Export the graphs

(* 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];*)

Appendix: Pool formation

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