This Workbook sets up, simulates and graphs the equations that describe adenosine system in an open system. It corresponds to the adenosine phosphate system with changes in the total AXP inventory Section 8.3.3 of Ch 8 in SB2.

We 1) define the numerical values of the kinetic parameters, 2) then the reaction rate laws, 3) then we integrate and, 4) finally graph the solution.

BOP May 16, 2007
Addiel de Alba. April 2010
NS, checked on March 15, 2011 OK

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

tfinal=500;
parameters={ka1->6.25`,ka2->25.`,ka5->1000,kam5->1000,b1->0.03`,ka4->(b1/(0.1`))};

FInd the steady state

va1=ka1 atp[t];
va2=ka2 adp[t];
va3=b1;
va4=ka4 amp[t];
va5=ka5 adp[t] adp[t];
vam5=kam5 atp[t] amp[t];

ststsolution=NSolve[{-ka1 atp+ka2 adp+(ka5 adp adp-kam5 atp amp)==0,ka1 atp-ka2 adp-2 (ka5 adp adp-kam5 atp amp)==0,b1-ka4 amp+(ka5 adp adp-kam5 atp amp)==0}//.parameters,{atp,adp,amp}]
{fixpoint1,fixpoint2}=ststsolution;

Perturb the system

perturbParam={ka1->1.5*6.25`,ka2->25.`,ka5->1000,kam5->1000,b1->0.03`,ka4->(b1/(0.1`)),ka1->1.5*6.25};

va1=ka1 atp[t];
BalanceATP=-va1+va2+(va5-vam5);
BalanceADP=va1-va2-2 (va5-vam5);
BalanceAMP=va3-va4+(va5-vam5);

solution=NDSolve[{atp'[t]==BalanceATP,adp'[t]==BalanceADP,amp'[t]==BalanceAMP,atp[0]==1.6`,adp[0]==0.4`,amp[0]==0.1`}//.perturbParam,{atp,adp,amp},{t,0,tfinal}];

graph1=ParametricPlot[Evaluate[{ka1 atp[t],ka2 adp[t]}/.solution//.perturbParam],{t,0,tfinal},PlotLabel->"ATP load vs ATP synthesis",PlotPoints->1000,Epilog->ppAnnotate[{ka1 atp[#]&,ka2 adp[#]&}/.perturbParam,solution,tfinal],FrameLabel->{"ka1 * atp","ka2 * adp"}]

graph2=ParametricPlot[Evaluate[{ka1 atp[t],ka4 amp[t]}/.solution//.perturbParam],{t,0,tfinal},PlotLabel->"ATP load vs AMP drain",Epilog->ppAnnotate[{ka1 atp[#]&,ka4 amp[#]&}//.perturbParam,solution,tfinal],FrameLabel->{"ka1 * atp","ka4 * amp"}]

graph3=ParametricPlot[Evaluate[{2 atp[t]+adp[t],atp[t]+adp[t]+amp[t]}/.solution//.perturbParam],{t,0,tfinal},PlotLabel->"occupancy vs capacity",PlotRange->{{1,4},{1.`,2.4`}},Epilog->ppAnnotate[{2 atp[#]+adp[#]&,atp[#]+adp[#]+amp[#]&},solution,tfinal],FrameLabel->{"2 atp + adp","atp + adp + amp"}]

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

Quiet@LogLinearPlot[{((atp[t]+((adp[t])/2))/(atp[t]+adp[t]+amp[t]))/.solution},{t,1*^-7,tfinal},PlotLabel->"Energy charge vs time",Epilog->createLegend[{"E.C"}],AxesLabel->{"t",None}]

graph5=Plot[Evaluate[{atp[t],adp[t],amp[t]}/.solution],{t,0,1},PlotLabel->"Concentrations vs time",Epilog->createLegend[{"atp","adp","amp"}],PlotRange->{0,1.6`}]
graph6=Plot[Evaluate[{atp[t],adp[t],amp[t]}/.solution],{t,0,tfinal},PlotLabel->"Concentrations vs time",Epilog->createLegend[{"atp","adp","amp"}],PlotRange->{0,1.6`}]
graph6=Quiet@LogLinearPlot[Evaluate[{atp[t],adp[t],amp[t]}/.solution],{t,1*^-7,tfinal},PlotLabel->"Concentrations vs time",Epilog->createLegend[{"atp","adp","amp"}],PlotRange->{0,1.6`}]

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;

here we export the graphs

Export["Section8.3.3.A.eps",graph1];
Export["Section8.3.3.B.eps",graph2];
Export["Section8.3.3.A2.eps",graph3];
Export["Section8.3.3.B2.eps",graph4];
Export["Section8.3.3.C2.eps",graph5];