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

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

BOP April 24, 2007
Addiel de Alba. April 2010
NS, checked on March 15, 2011

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

FInd the steady state

parameters={kform->10/0.4,kuse->10/1.6,kdist->1000,kmdist->1000, Atot->2.1};

vuse=kuse*atp[t];
vform=kform*adp[t];
vdist=kdist*adp[t]*adp[t];
vmdist=kmdist*atp[t]*amp[t];

BalanceATP=-vuse+vform+(vdist-vmdist);
BalanceADP=vuse-vform-2*(vdist-vmdist);
BalanceAMP=(vdist-vmdist);


ststsolution=NSolve[{-kuse*atp+kform*adp==0,kdist*adp*adp-kmdist*atp*amp==0, amp+adp+atp-Atot==0}/.parameters,{amp,adp,atp}]

{fixpoint1,fixpoint2}=ststsolution;

tfinal=10;
solution=NDSolve[{atp'[t]==BalanceATP,adp'[t]==BalanceADP,amp'[t]==BalanceAMP,atp[0]==(atp/.fixpoint2),adp[0]==(adp/.fixpoint2),amp[0]==(amp/.fixpoint2)}/.parameters,{atp,adp,amp},{t,0,tfinal}];

graph1=Plot[Evaluate[{ atp[t],adp[t],amp[t]}/.solution],
{t,0,1},PlotLabel->"Concentrations vs time",PlotRange->{0,2.1},Epilog->createLegend[{"atp","adp","amp"}]]

test={kuse*atp[t],kform*adp[t]}/.parameters/.solution;
ststPlot=ParametricPlot[test[[1]],{t,0,tfinal},PlotLabel->"ATP load vs ATP synthesis",PlotRange->{{9,16},{9,16}},Epilog->ppAnnotate[{kuse*atp[#]&,kform*adp[#]&}/.parameters,solution,tfinal,endColor->Black],FrameLabel->{"rate of use","rate of formation"}]

Solve dynamic behavior from changing the IC or the use rate

perturbation=1.5;
parameters={kform->10/0.4,kuse->perturbation*10/1.6,kdist->1000,kmdist->1000, Atot->2.1};

tfinal=10;
solution=NDSolve[{atp'[t]==BalanceATP,adp'[t]==BalanceADP,amp'[t]==BalanceAMP,atp[0]==1.60,adp[0]==0.40,amp[0]==0.1}/.parameters,{atp,adp,amp},{t,0,tfinal}];

graph1=Plot[Evaluate[{atp[t],adp[t],amp[t]}/.parameters/.solution],
{t,0,1},PlotLabel->"Concentrations vs time",AxesLabel->{"t","Concentration"},PlotRange->{0,2.1},Epilog->createLegend[{"atp","adp","amp"}]]

graph2=ParametricPlot[Evaluate[{kuse*atp[t],kform*adp[t]}/.parameters/.solution],{t,0,tfinal},PlotLabel->"ATP load vs ATP synthesis",PlotRange->All,FrameLabel->{"rate of use","rate of formation"} ,Epilog->ppAnnotate[{kuse*atp[#]&,kform*adp[#]&}/.parameters,solution,tfinal]]

graph3=Plot[{(atp[t]+adp[t]/2)/(atp[t]+adp[t]+amp[t])/.solution},
{t,0,1},PlotLabel->"Energy charge vs time",Epilog->createLegend[{"E.C"}],PlotRange->{0.7,0.9},AxesLabel->{"t",None}]


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

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

Manipulate[

parameters={kform->10/0.4,kuse->perturbation*10/1.6,kdist->1000,kmdist->1000, Atot->2.1};

solution=NDSolve[{atp'[t]==BalanceATP,adp'[t]==BalanceADP,amp'[t]==BalanceAMP,atp[0]==atpInit,adp[0]==adpInit,amp[0]==ampInit}/.Dispatch[parameters],{atp,adp,amp},{t,0,tfinal+1}];

graph1=Quiet@LogLinearPlot[Evaluate[{atp[t],adp[t],amp[t]}/.parameters/.solution],
{t,1*^-4,tfinal},PlotLabel->"Concentrations vs time",PlotRange->{All,{0,All}},Epilog->createLegend[{"atp","adp","amp"}],AxesLabel->{"t","Concentration"}];

graph2=ParametricPlot[Evaluate[{kuse*atp[t],kform*adp[t]}/.parameters/.solution],{t,0,tfinal},PlotLabel->"ATP load vs ATP synthesis",PlotRange->{{0,30},{0,30}} ,Epilog->ppAnnotate[{kuse*atp[#]&,kform*adp[#]&}/.parameters,solution,tfinal],FrameLabel->{"rate of use","rate of formation"}];

GraphicsRow[{graph1,Show[graph2,Graphics[{Dashed,Black,Line[{{0,0},{50,50}}]}]]},ImageSize->600*2]
,{{perturbation,1.,"kuse perturbation"},.1,2.9,.1},{{atpInit,1.6,"ATP init. conc."},0.,2.,.1},{{adpInit,0.4,"ADP init. conc."},0.,1.,.1},{{ampInit,0.1,"AMP init. conc."},0.,1.,.1},{{tfinal,10},.1,50,1}]

here we export the graphs

Export["Section8.3.1.A.eps",graph2];
Export["Section8.3.1.B.eps",graph1];
Export["Section8.3.1.C.eps",graph3];
