This Workbook sets up, simulates and graphs the equations that describe adenosine system. It corresponds to Section 8.3.2 of Chapter 8 in SB2.
The reactions considered are: 

v,use :  ATP--> ADP
v,form: ADP-ATP
v, distr :  2ADP<=>ATP+AMP
v, buff ATP+Buff <=> ADP +BuffP

BOP April 24, 2007
Addiel de Alba. April 2010.

FInd the steady state

kform=10/0.4;
kuse=10/1.6;
kdistr=1000;kmdistr=1000; Atot=2.1; Btot= 10;
kbuff =1000;
kmbuff=1000;
Kbuff=kbuff/kmbuff;
b0=0.4*Btot/(Kbuff*1.6+.4);
ststsolution=NSolve[{-kuse*a3+kform*a2==0,kdistr*a2*a2-kmdistr*a3*a1==0,a3+a2+a1-Atot==0},{a3,a2,a1}];
{z1,z2,z3}={a1,a2,a3}/.ststsolution[[2]];
 kbuff*z3*b0-kmbuff*z2*(Btot-b0);

Solve dynamic behavior from changing the IC or the use rate

kuse=((1.5` 10)/(1.6`));
vuse=kuse atp[t];
vform=kform adp[t];
vdistr=kdistr adp[t] adp[t];
vmdistr=kmdistr atp[t] amp[t];
vbuff=kbuff buff[t] atp[t];
vmbuff=kmbuff buffp[t] adp[t];
BalanceATP=-vuse+vform+(vdistr-vmdistr)-(vbuff-vmbuff);
BalanceADP=vuse-vform-2 (vdistr-vmdistr)+(vbuff-vmbuff);
BalanceAMP=vdistr-vmdistr;
BalanceBuff=vbuff-vmbuff;
tfinal=1;
solution=NDSolve[{atp'[t]==BalanceATP,adp'[t]==BalanceADP,amp'[t]==BalanceAMP,buff'[t]==-BalanceBuff,buffp'[t]==BalanceBuff,atp[0]==1.6`,adp[0]==0.4`,amp[0]==0.1`,buff[0]==b0,buffp[0]==Btot-b0},{atp,adp,amp,buff,buffp},{t,0,tfinal}];
<<"PlotLegends`"

graph1=Plot[{((atp[t]+((adp[t])/2))/(atp[t]+adp[t]+amp[t]))/.solution,((buffp[t])/10)/.solution},{t,0,1},PlotLabel->"Energy and buffer charge vs time",PlotLegend->{"E.C.","B.C."},PlotStyle->{Hue[0.5`],Hue[0.6`],Hue[0.8`]},PlotRange->{0.7`,0.9`},LegendPosition->{1.1`,-0.4`},LegendSize->{0.5`,0.35`},ImageSize->800]

graph2=Plot[{atp[t]/.solution,adp[t]/.solution,amp[t]/.solution,buff[t]/.solution},{t,0,1},PlotLabel->"Concentrations vs time",PlotLegend->{"atp","adp","amp","buff"},PlotStyle->{Hue[0.5`],Hue[0.6`],Hue[0.8`],Hue[1.`]},PlotRange->{0,2.1`},LegendPosition->{1.1`,-0.4`},LegendSize->{0.45`,0.35`}];

graph3=Plot[{atp[t]/.solution,adp[t]/.solution,amp[t]/.solution,buff[t]/.solution},{t,0,0.01`},PlotLabel->"Concentrations vs time",PlotLegend->{"atp","adp","amp","buff"},PlotStyle->{Hue[0.5`],Hue[0.6`],Hue[0.8`],Hue[1.`]},PlotRange->{0,2.1`},LegendPosition->{1.1`,-0.4`},LegendSize->{0.45`,0.35`}];

graph4=ParametricPlot[Evaluate[{kuse atp[t],kform adp[t]}/.solution],{t,0,tfinal},PlotLabel->"ATP load vs ATP synthesis",PlotStyle->{Hue[0.7`]},PlotRange->{{9,15},{9,15}}];

(((((atp[0]/.solution) atp[tfinal]/.solution) adp[0]/.solution) adp[tfinal]/.solution) amp[0]/.solution) amp[tfinal]/.solution;

graph5=Plot[{((atp[t]+((adp[t])/2))/(atp[t]+adp[t]+amp[t]))/.solution,((buffp[t])/10)/.solution},{t,0,1},PlotStyle->{Hue[0.5`],Hue[0.6`]},PlotRange->{0.7`,0.9`}];
((atp[0]+((adp[0])/2))/(atp[0]+adp[0]+amp[0]))/.solution;
((atp[tfinal]+((adp[tfinal])/2))/(atp[tfinal]+adp[tfinal]+amp[tfinal]))/.solution;

((buffp[0])/10)/.solution;
((buffp[tfinal])/10)/.solution;
((2 atp[0]+adp[0]+buffp[0])/(2 (atp[0]+adp[0]+amp[0])+buff[0]+buffp[0]))/.solution;
((2 atp[tfinal]+adp[tfinal]+buffp[tfinal])/(2 (atp[tfinal]+adp[tfinal]+amp[tfinal])+buff[tfinal]+buffp[tfinal]))/.solution;

here we export the graphs

Export["Section8.3.2.eps",graph1];
