Hb Notebook

Off[General::spell]; Off[General::spell1];
SetDirectory["c:/users/aabordba.alanine0/desktop/website"];  
<<functions.m

Section I: 	System Definition and Setting up the Equations

Define network structure/topology and the nameas of the concentrations and variables;  
	there are m=8 compounds and n=8 reactions.  
	the rank of S is r=7, 
	there is one time invariant pools, Dim[Left Null]]  =  m-r = 8-7 = 1
	there is one steady-state pathways, Dim[Null]] =  n-r = 8-7= 1
	there are two compounds that can exchange with system boundary; they are o2 and 23dpg
	hemoglobin stays confined with in the system

variables = {dpg23,o2,hb,hbo2,hbo22,hbo23,hbo24,dhb};
varnames={"dpg23","o2","hb","hbo2","hbo22","hbo23","hbo24","dhb"};
elemnames={"C","H","O","P","Hb"};
MatrixForm[elem=(({{3, 0, 0, 0, 0, 0, 0, 3}, {3, 0, 0, 0, 0, 0, 0, 3}, {10, 2, 0, 2, 4, 6, 8, 10}, {2, 0, 0, 0, 0, 0, 0, 2}, {0, 0, 1, 1, 1, 1, 1, 1}}))];
MatrixForm[elem,TableHeadings->{elemnames,varnames}]
TeXForm[MatrixForm[elem,TableHeadings->{elemnames,varnames}]];

fluxes={vdpgm,vdpgase,vhbo1,vhbo2,vhbo3,vhbo4,vhbdpg,vo2};
fluxnames={"vdpgm","vdpgase","vhbo1","vhbo2","vhbo3","vhbo4","vhbdpg","vo2"};
knames={"kdpgm","kdpgase","khbo1","khbo2","khbo3","khbo4","khbdpg","ko2"};
Knames={"Kdpgm","Kdpgase","Khbo1","Khbo2","Khbo3","Khbo4","Khbdpg","Ko2"};
MatrixForm[stoich=(({{1, -1, 0, 0, 0, 0, -1, 0}, {0, 0, -1, -1, -1, -1, 0, -1}, {0, 0, -1, 0, 0, 0, -1, 0}, {0, 0, 1, -1, 0, 0, 0, 0}, {0, 0, 0, 1, -1, 0, 0, 0}, {0, 0, 0, 0, 1, -1, 0, 0}, {0, 0, 0, 0, 0, 1, 0, 0}, {0, 0, 0, 0, 0, 0, 1, 0}}))];
Print["The Stoichiometric Matrix is: "]
MatrixForm[stoich,TableHeadings->{varnames,fluxnames}]
TeXForm[MatrixForm[stoich,TableHeadings->{varnames,fluxnames}]];
dim=Dimensions[stoich,2]
rank=MatrixRank[stoich]

The null spaces:  determining the pathways and the pools

null=NullSpace[stoich];
leftnull=NullSpace[Transpose[stoich]];

MatrixForm[null,TableHeadings->{{"RL-shunt"},fluxnames}]
TeXForm[MatrixForm[null,TableHeadings->{{"p1"},fluxnames}]];

Print["The left null space is: "];
MatrixForm[leftnull,TableHeadings->{{"hb"},varnames}]
TeXForm[
MatrixForm[leftnull,TableHeadings->{{"hb"},varnames}]];

MatrixForm[elem.stoich,TableHeadings->{elemnames,fluxnames}]
TeXForm[MatrixForm[elem.stoich,TableHeadings->{elemnames,fluxnames}]];

Define equilibrium constants
--use of thermodynamic information
--the irreversible ones are set very high;  arbitrarily at 1,000,000

infty=1000000;
pO2conv=2.8684*10^-4;
Kdpgm=infty;Kdpgase=infty;
Khbo1=0.012/pO2conv;Khbo2=0.021/pO2conv;Khbo3=0.051/pO2conv;Khbo4=0.37/pO2conv;Khbdpg=1/4;Ko2=1;

Kvalues={Kdpgm,Kdpgase,Khbo1,Khbo2,Khbo3,Khbo4,Khbdpg,Ko2};
Print["The values for the equilibrium constants are:"]
MatrixForm[Kvalues,TableHeadings->{Knames}]
Length[Kvalues];
TeXForm[MatrixForm[Kvalues,TableHeadings->{Knames}]];

Define exchange fluxes and loads on cofactors--one flux in each null space pathway vector--then compute the stst

rlsshunt=0.441;
givenststfluxes=rlsshunt;
MatrixForm[null];
solution=givenststfluxes*Flatten[null];
MatrixForm[solution,TableHeadings->{fluxnames}]

Define  external concentrations--this is an arbitrary number that must be scaled relative to the binding constant for hemoglobin

o2plasma=((100-40)/2+40)*pO2conv;

Define steady state concentration values; ie use aarash's K values from the literature, Kdhb unknown to me

dpg23stst=3.1;
o2stst=o2plasma;
hbtot=7.3;

eqn1 =Khbo1*hb0*o2stst-hb1;
eqn2 =Khbo2*hb1*o2stst-hb2;
eqn3 =Khbo3*hb2*o2stst-hb3;
eqn4 =Khbo4*hb3*o2stst-hb4;
eqn5 =Khbdpg*hb0*dpg23stst-hbd; (* 23dpg binding to hemoglobin independent of the o2 binding state *)
eqn6=hb0+hb1+hb2+hb3+hb4+hbd-hbtot; (* total mass balance on hemoglobin *)


ststsolution=NSolve[{eqn1==0,eqn2==0,eqn3==0,eqn4==0,eqn5==0,eqn6==0},{hb0,hb1,hb2,hb3,hb4,hbd}];
{hbstst,hbo2stst,hbo22stst,hbo23stst,hbo24stst,dhbstst}={hb0,hb1,hb2,hb3,hb4,hbd}/.ststsolution[[1]];

varstst={dpg23stst,o2stst,hbstst,hbo2stst,hbo22stst,hbo23stst,hbo24stst,dhbstst};

MatrixForm[varstst,TableHeadings->{varnames}]
TeXForm[MatrixForm[varstst,TableHeadings->{varnames}]];

Compute Pools: hard and soft pools
	use the left null space for 'hard pools'
	'soft pools' defined based on intuition
	compute pool sizes based on the steady state concentrations

totalpoolnames={"O2 on Hb","tot Hb cap","dhb amount","tot 23dpg"};
MatrixForm[pools=(({{0, 0, 0, 1, 2, 3, 4, 0}, {0, 0, 4, 4, 4, 4, 4, 4}, {0, 0, 0, 0, 0, 0, 0, 1}, {1, 0, 0, 0, 0, 0, 0, 1}}))];
Dimensions[pools,2];
MatrixForm[pools,TableHeadings->{totalpoolnames,varnames}]
TeXForm[MatrixForm[pools,TableHeadings->{totalpoolnames,varnames}]];
MatrixForm[poolsize=pools.varstst,TableHeadings->{totalpoolnames}]
BarChart[poolsize]

Calculate mass action ratios--uninformative here as all reactions except two are at equilibrium

Gvalues={0,0,Khbo1,Khbo2,Khbo3,Khbo4,Khbdpg,Ko2};

MatrixForm[Transpose[{0.001*Round[1000*Gvalues],0.001*Round[1000*Kvalues],0.001*Round[1000*Gvalues/Kvalues]}],TableHeadings->{fluxnames,{"G","K","G/K"}}];

TeXForm[MatrixForm[Transpose[{0.001*Round[1000*Gvalues],0.001*Round[1000*Kvalues],0.001*Round[1000*Gvalues/Kvalues]}],TableHeadings->{fluxnames,{"G","K","G/K"}}]];

Determine kinetic constants from the given fluxes and concentration values

kdpgase=solution[[2]]/(dpg23stst);
kdpgm=1;
khbo1=506935.270702259;
khbo2=511077.050923776;
khbo3=509243.459567699;
khbo4=501595.340624411;
khbdpg=519612.560391792;
ko2=509725.707725914;


kvalues={kdpgm,kdpgase,khbo1,khbo2,khbo3,khbo4,khbdpg,ko2};

MatrixForm[Transpose[{0.001*Round[1000*Gvalues],0.001*Round[1000*Kvalues],0.001*Round[1000*Gvalues/Kvalues],kvalues}],TableHeadings->{fluxnames,{"G","K","G/K","k"}}]
TeXForm[MatrixForm[kvalues,TableHeadings->{knames}]];

MatrixForm[pools.stoich,TableHeadings->{totalpoolnames,kvalues}];

Define the rate laws

vo2 = ko2*(o2-o2plasma);

vdpgm=0.441;
vdpgase=kdpgase*dpg23;
vhbo1=khbo1*(hb*o2-hbo2/Khbo1);
vhbo2=khbo2*(hbo2*o2-hbo22/Khbo2);
vhbo3=khbo3*(hbo22*o2-hbo23/Khbo3);
vhbo4=khbo4*(hbo23*o2-hbo24/Khbo4);
vhbdpg=khbdpg*(hb*dpg23-dhb/Khbdpg);



fluxstst= Flatten[fluxes/.sublist[variables,varstst]];

Print["The steady state fluxes are:  "]
MatrixForm[0.001*Round[1000*fluxstst],TableHeadings->{fluxnames}]
BarChart[fluxstst]

TeXForm[
MatrixForm[MatrixForm[0.001*Round[1000*fluxstst],TableHeadings->{fluxnames}]]];

MatrixForm[0.01*Round[100*(solution-fluxstst)], TableHeadings->{fluxnames}];

Calculate the flux balances

Print["Checking the flux balances:  S.v = "]
Chop[MatrixForm[0.0001*Round[10000*stoich.fluxstst],TableHeadings-> {varnames}]]

Section II: 	Defining pools and ratios

Calculate the reaction that move the pools

MatrixForm[pstoich=pools.stoich,TableHeadings->{totalpoolnames,fluxnames}]

TeXForm[MatrixForm[pools.stoich,TableHeadings->{totalpoolnames,fluxnames}]];
MatrixForm[pools.stoich.fluxnames,TableHeadings->{totalpoolnames}];
MatrixForm[pools.stoich.fluxstst,TableHeadings->{totalpoolnames}];

MatrixForm[pools.stoich.fluxnames,TableHeadings->{totalpoolnames,{"reactions"}}]

Calculate the response times of the pools

ppstoich=pstoich;
Do[Do[ppstoich[[i,j]]=If[pstoich[[i,j]]>0,pstoich[[i,j]],0],{j,1,8}],{i,1,4}]

pnstoich=pstoich;
Do[Do[pnstoich[[i,j]]=If[pstoich[[i,j]]<0,pstoich[[i,j]],0],{j,1,8}],{i,1,4}]MatrixForm[tau = {0,Infinity,0}];


MatrixForm[tau = {0,Infinity,0,0}];
Do[tau[[i]]=poolsize[[i]]/(ppstoich[[i]].fluxstst),{i,1,1}]
Do[tau[[i]]=poolsize[[i]]/(ppstoich[[i]].fluxstst),{i,3,4}]
MatrixForm[fluxnet ={0,".",0,0}];
Do[fluxnet[[i]]=ppstoich[[i]].fluxstst,{i,1,1}]
Do[fluxnet[[i]]=ppstoich[[i]].fluxstst,{i,3,4}]
MatrixForm[Transpose[{poolsize,fluxnet,tau}],TableHeadings->{totalpoolnames,{"size","net flux","tau"}}]
TeXForm[MatrixForm[Transpose[{poolsize,fluxnet,tau}],TableHeadings->{totalpoolnames,{"size","net flux","tau"}}]];

Compute charge state of pools

rationames={"hb_usage","dhb ratio"};
MatrixForm[rnumerator=({{1, 0, 0, 0}, {0, 0, 1, 0}})];
MatrixForm[rdenominator=({{0, 1, 0, 0}, {0, 0, 0, 1}})];
ratiovars={0,0};
MatrixForm[Do[ratiovars[[i]]=(rnumerator[[i]].poolsize/rdenominator[[i]].poolsize),{i,1,2}]];
MatrixForm[Transpose[{rnumerator.totalpoolnames,rnumerator.poolsize,rdenominator.totalpoolnames,rdenominator.poolsize,ratiovars}],
TableHeadings->{rationames,{"numerator","size(mM)","denominator","size(mM)","charge"}}]

Section III: 	Simulation of the Non-Linear Mass Balance Equations (oscillatory)

Set perturbations

tfinal=1;
o2plasma=(30*Cos[120*Pi*t]+70)*pO2conv;
vo2=ko2*(o2-o2plasma);
massbal=stoich.fluxes;

Compute the solution

ic=varstst;

solution=solveODE[massbal,ic,variables,0,tfinal,20000];

Section IV: 	Plots and Interpretation (oscillatory)

Fluxes

Examine the beginning, interim and the end flux map

tHigh=0.018;
tLow= 0.027;
tflux=fluxes/.tsub[variables];
tset1=Table[Flatten[{tflux[[i]]/.t->0,tflux[[i]]/.t->tHigh,tflux[[i]]/.t->tLow,tflux[[i]]/.t->tfinal}/.solution],{i,1,Length[fluxes]}];
MatrixForm[tset1,TableHeadings->{fluxnames,{0,tHigh,tLow,tfinal}}]
BarChart[Transpose[tset1][[1]],Transpose[tset1][[2]],Transpose[tset1][[3]],Transpose[tset1][[4]],BarLabels->fluxnames]

Selected flux phase planes

graph1vp=MASSPhasePlaneFlux[{vdpgase,vhbo1},{"vdpgase","vhbo1"},variables,solution,0,tfinal];
graph2vp=MASSPhasePlaneFlux[{vhbo4,vo2},{"vhbo4","vo2"},variables,solution,tLow,tHigh];

{graph1vp,graph2vp}

Select time responses of fluxes

graph1v=MASSTimeProfileFlux[{vdpgm,vdpgase},{"vdpgm","vdpgase"},variables,solution,0,0.1];
graph2v=MASSTimeProfileFlux[{vhbo2,vhbo3,vhbo4},{"vhbo2","vhbo3","vhbo4"},variables,solution,0,0.1];
graph3v=MASSTimeProfileFlux[{vhbo1},{"vhbo1","vhbdpgo1"},variables,solution,0,0.1];
graph4v=MASSTimeProfileFlux[{vo2},{"vo2"},variables,solution,0,0.1];
graph5v=MASSTimeProfileFlux[{vhbdpg},{"vhbdpg"},variables,solution,0,0.1];
{graph1v,graph2v,graph3v,graph4v,graph5v}

Concentrations

Get the start and end concentration states of the simulation

tvar=variables/.tsub[variables];
tset=Table[Flatten[{tvar[[i]]/.t->0,tvar[[i]]/.t->tHigh,tvar[[i]]/.t->tLow,tvar[[i]]/.t->tfinal}/.solution],{i,1,Length[variables]}];
MatrixForm[tset,TableHeadings->{variables,{0,tHigh,tLow,tfinal}}]
BarChart[Transpose[tset][[1]],Transpose[tset][[2]],Transpose[tset][[3]],Transpose[tset][[4]],BarLabels->varnames]

Selected time responses of concentrations

graph1c=MASSTimeProfileMet[{dpg23},{"o2"},solution,0,0.25];
graph2c=MASSTimeProfileMet[{dpg23,dhb,hbo24},{"dpg23","dhb","hbo24"},solution,0,0.1];
{graph1c,graph2c}

Pool fluxes

calculate the fluxes in and out of the pools

tpoolfluxin = ppstoich . tflux;
tset3 = Table[Flatten[{tpoolfluxin[[i]] /. t -> 0, tpoolfluxin[[i]] /. t -> tHigh, tpoolfluxin[[i]] /. t -> tLow,tpoolfluxin[[i]] /. t -> tfinal} /. solution], {i, 1, Length[tpoolfluxin]}]; 
matin=MatrixForm[tset3, TableHeadings -> {totalpoolnames, {0, tHigh,tLow, tfinal}}];barin=BarChart[Transpose[tset3][[1]], Transpose[tset3][[2]], Transpose[tset3][[3]],Transpose[tset3][[4]], BarLabels -> totalpoolnames];

tpoolfluxout = -pnstoich . tflux; 
tset4 = Table[Flatten[{tpoolfluxout[[i]] /. t -> 0, tpoolfluxout[[i]] /. t -> tHigh, tpoolfluxout[[i]] /. t -> tLow, tpoolfluxout[[i]] /. t -> tfinal} /. solution], {i, 1, Length[tpoolfluxout]}]; 
matout=MatrixForm[tset4, TableHeadings -> {totalpoolnames, {0, tHigh,tLow, tfinal}}];
barout=BarChart[Transpose[tset4][[1]], Transpose[tset4][[2]], Transpose[tset4][[3]],Transpose[tset4][[4]], BarLabels -> totalpoolnames];

{matin,matout}
{barin,barout}

flux phase planes for pools

Table[MASSPhasePlanePoolFlux[i,totalpoolnames[[i]],fluxes,variables,ppstoich,pnstoich,solution,0,tfinal],{i,1,4}]

Time responses of total pool fluxes (can break up into in/out if of interest)

Table[MASSTimeProfilePoolFlux[i,totalpoolnames[[i]],fluxes,variables,ppstoich,pnstoich,solution,0,tfinal],{i,1,4}]

Pools

Calculate pools at start and end of simulations

tpool=pools.tvar;
tset2=Table[Flatten[{tpool[[i]]/.t->0,tpool[[i]]/.t->tHigh,tpool[[i]]/.t->tLow,tpool[[i]]/.t->tfinal}/.solution],{i,1,Length[tpool]}];
MatrixForm[tset2,TableHeadings->{totalpoolnames,{0,tHigh,tLow,tfinal}}]
BarChart[Transpose[tset2][[1]],Transpose[tset2][[2]],Transpose[tset2][[3]],Transpose[tset2][[4]],BarLabels->totalpoolnames]
tpool[[1]]/.t->1/.solution;

Selected pool phase planes

graph1pp=MASSPhasePlanePool[{pools[[1]],pools[[2]]},{totalpoolnames[[1]],totalpoolnames[[2]]},variables,solution,0,tfinal];
graph2pp=MASSPhasePlanePool[{pools[[2]],pools[[3]]},{totalpoolnames[[2]],totalpoolnames[[3]]},variables,solution,0,tfinal];
graph3pp=MASSPhasePlanePool[{pools[[1]],pools[[3]]},{totalpoolnames[[1]],totalpoolnames[[3]]},variables,solution,0,tfinal];
{graph1pp,graph2pp,graph3pp}

Time responses of pools

graph1p=MASSTimeProfilePool[{pools[[1]],pools[[2]]},{totalpoolnames[[1]],totalpoolnames[[2]]},variables,solution,0,tfinal];
graph2p=MASSTimeProfilePool[{pools[[3]]},{totalpoolnames[[3]]},variables,solution,0,tfinal];

{graph1p,graph2p}

Ratios

Compute the ratios at the beginning and the end

tratio={0,0};
Do[tratio[[i]]=(rnumerator[[i]].tpool/rdenominator[[i]].tpool),{i,1,2}]
tset5=Table[Flatten[{tratio[[i]]/.t->0,tratio[[i]]/.t->tHigh,tratio[[i]]/.t->tLow,tratio[[i]]/.t->tfinal}/.solution],{i,1,Length[tratio]}];
MatrixForm[tset5,TableHeadings->{rationames,{0,tHigh,tLow,tfinal}}]
BarChart[Transpose[tset5][[1]],Transpose[tset5][[2]],Transpose[tset5][[3]],Transpose[tset5][[4]],BarLabels->rationames]

Compute the time responses of the ratios

Table[MASSTimeProfileRatios[{i},rationames[[i]],variables,pools,rnumerator,rdenominator,solution,0,tfinal],{i,1,2}]