Symmetry model Section 5.5

BOP and OR 
3 October 2007.
AB - Last updated: Mar 26, 2010.
NS, March 31, 2011

SetDirectory[NotebookDirectory[]];
<<"GraphicsDirectives.m";

metabolites={"s","e","p","x","i","x1","x2","x3","x4"};
reactions={"v1","v2","v3","v4","v5","v6","v7","v8","v9","v10","v11"};
var={x1,x2,x3,x4,x5,x6,x7,x8,x9};
S=(({{-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, -1, 1, 0, 0, 0, 0, 0, 0, 0, 0}, {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 1, -1, -1, 1, 0, 0, 0, 0, 0, 0}, {0, 0, 0, -1, 1, -1, 1, -1, 1, -1, 1}, {0, 0, 0, 1, -1, -1, 1, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 1, -1, -1, 1, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 1, -1, -1, 1}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1}}));
MatrixForm[S,TableHeadings->{metabolites,reactions}]
MatrixForm[Table[{var[[i]],metabolites[[i]]},{i,9}],TableHeadings->{None,{"Notation", "Compounds"}}]

Left Null space

MatrixForm[(({{1, 1, 0}, {4, 3, 0}, {0, 0, 1}})).NullSpace[Transpose[S]].var]

Parameter Definition

ν=4;
k=1;
kl=100;
kli=100;
ki=2;
kii=2;

Setting up the rate laws and the dynamic mass balances

v1=k*x1[t]*x2[t];
v2=kl*x2[t];
v3=kli*x4[t];
v4=4*ki*x4[t]*x5[t];
v5=kii*x6[t];
v6=3*ki*x5[t]*x6[t];
v7=2*kii*x7[t];
v8=2*ki*x7[t]*x5[t];
v9=3*kii*x8[t];
v10=ki*x8[t]*x5[t];
v11=4*kii*x9[t];
v={v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11};
dxdt=S.v;

Simulation and Plotting

tfinal=50;
solution=NDSolve[{x1'[t]==dxdt[[1]],x2'[t]==dxdt[[2]],x3'[t]==dxdt[[3]],x4'[t]==dxdt[[4]],x5'[t]==dxdt[[5]],x6'[t]==dxdt[[6]],x7'[t]==dxdt[[7]],x8'[t]==dxdt[[8]],x9'[t]==dxdt[[9]],x1[0]==1,x2[0]==1,x3[0]==0,x4[0]==0,x5[0]==1,x6[0]==0,x7[0]==0,x8[0]==0,x9[0]==0},{x1,x2,x3,x4,x5,x6,x7,x8,x9},{t,0,tfinal}][[1]];

Plots

graph1=ParametricPlot[Evaluate[{x1[t],x2[t]}/.solution],{t,0,tfinal},FrameLabel->{"s","e"},Epilog->ppAnnotate[{x1,x2},solution,tfinal]]

graph2=ParametricPlot[Evaluate[{x2[t],x4[t]+x6[t]+x7[t]+x8[t]+x9[t]}/.solution],{t,0,10},FrameLabel->{"e","x+x1+x2+x3+x4"},Epilog->ppAnnotate[{x2,(x4[#]+x6[#]+x7[#]+x8[#]+x9[#])&},solution,tfinal]]

graph3=ParametricPlot[Evaluate[{x5[t],x6[t]+2 x7[t]+3 x8[t]+4 x9[t]}/.solution],{t,0,10},AxesLabel->{"i","x1+2x2+3x3+4x4"},Epilog->ppAnnotate[{x5,(x6[#]+2 x7[#]+3 x8[#]+4 x9[#])&},solution,tfinal]]

graph4=Plot[Evaluate[{x1[t],x2[t],x3[t],x4[t]+x6[t]+x7[t]+x8[t]+x9[t],x6[t]+2 x7[t]+3 x8[t]+4 x9[t]}/.solution],{t,0,0.1`},FrameLabel->{"Time",""},Epilog->createLegend[{"S","E","P","Xtot","Bound I"}]]

graph5=Plot[Evaluate[{x1[t],x2[t],x3[t],x4[t]+x6[t]+x7[t]+x8[t]+x9[t],x6[t]+2 x7[t]+3 x8[t]+4 x9[t]}/.solution],{t,0,10.},FrameLabel->{"Time",""},Epilog->createLegend[{"S","E","P","Xtot","Bound I"}]]

Step 7. Export Graphs.

(* Export["Section5.5.A.eps",graph1];
Export["Section5.5.B.eps",graph2];
Export["Section5.5.C.eps",graph3];
Export["Section5.5.D.eps",graph4];
Export["Section5.5.E.eps",graph5]; *)