// Simulation of oxyhemoglobin (HbO2) and carbomino hemoglobin (HbCO2)
// dissociation curves and computation of total O2 and CO2 contents in 
// whole blood: A new approach

JSim v1.1
math SHbCO2 {		// CO2 saturation of hemoglobin
   realDomain pCO2;	// CO2 partial pressure, mmHg
   pCO2.min=0; pCO2.max=100; pCO2.delta=0.5;

// State variables those are computed as functions of the independent
// variable pCO2 and physiological parameters pO2, pH, DPG and Temp
   real 
	CO2(pCO2),	// CO2 concentration in RBCs, M
	K4p(pCO2),	// equilibrium constant for Hb-O2 binding, 1/M
			// = K4dp*([O2]/[O20])^n0*([H+]0/[H+])^n1*
			// ([CO2]0/[CO2])^n2*([DPG]0/[DPG])^n3*(T0/T)^n4
	KHbO2(pCO2),	// Hill coefficient for SHbO2, 1/M
	KHbCO2(pCO2),	// Hill coefficient for SHbCO2, 1/M
	SHbO2(pCO2),	// O2 saturation of hemoglobin, unitless
	SHbCO2(pCO2),	// CO2 saturation of hemoglobin, unitless
	O2free,		// free O2 concentration in blood, M
	O2bound(pCO2),	// Hb-bound O2 concentration in blood, M
	O2total(pCO2),	// total O2 concentration in blood, M
	O2content(pCO2),// total O2 content in blood (ml O2/100 ml blood)
	CO2free(pCO2),	// free CO2 concentration in blood, M
	CO2bicarb(pCO2),// CO2 in bicarbonate (HCO3-) form, M
	CO2bound(pCO2),	// Hb-bound CO2 concentration in blood, M
	CO2total(pCO2),	// total CO2 concentration in blood, M
	CO2content(pCO2),// total CO2 content in blood (ml O2/100 ml blood)
	HbNH2(pCO2),	// HbNH2 concentration in RBCs, M
	HbNH3p(pCO2),	// HbNH3+ concentration in RBCs, M
	O2HbNH2(pCO2),	// O2HbNH2 concentration in RBCs, M
	O2HbNH3p(pCO2),	// O2HbNH3+ concentration in RBCs, M
	HbNHCOOH(pCO2),	// HbNHCOOH concentration in RBCs, M
	HbNHCOOm(pCO2),	// HbNHCOO- concentration in RBCs, M
	O2HbNHCOOH(pCO2),// O2HbNHCOOH concentration in RBCs, M
	O2HbNHCOOm(pCO2),// O2HbNHCOO- concentration in RBCs, M
	SHbO2kin(pCO2),	// O2 saturation of hemoglobin from kinetics, unitless
	SHbCO2kin(pCO2);// CO2 saturation of hemoglobin from kinetics, unitless

// State variables those represent the O2 and CO2 contents in region-wise
   real O2freepl1, O2freepl2, O2freerbc1, O2freerbc2, O2boundrbc1(pCO2), 
	O2boundrbc2(pCO2), CO2freepl1(pCO2), CO2freepl2(pCO2), 
	CO2freerbc1(pCO2), CO2freerbc2(pCO2), CO2bicarbpl1(pCO2), 
	CO2bicarbpl2(pCO2), CO2bicarbrbc1(pCO2), CO2bicarbrbc2(pCO2), 
	CO2boundrbc1(pCO2), CO2boundrbc2(pCO2);

// Parameters those are privately fixed in the model (i.e., water fractions,
// hemoglobin concentration in RBCs, equilibrium constants, and Hill index)
   real
	Wpl = 0.94,	// fractional water space in plasma, unitless
	Wrbc = 0.65,	// fractional water space in RBCs, unitless
	Rrbc = 0.69,	// Gibbs-Donnan ratio across RBC membrane, unitless
	Hbrbc = 5.18e-3,// hemoglobin concentration in RBCs, M
	K2 = 2.95e-5,	// CO2 + HbNH2 equilibrium constant, unitless
	K2dp = 1.0e-6,	// HbNHCOOH dissociation constant, M
	K2p = K2/K2dp,	// kf2p/kb2p, 1/M
	K3 = 2.51e-5,	// CO2 + O2HbNH2 equilibrium constant, unitless
	K3dp = 1.0e-6,	// O2HbNHCOOH dissociation constant, M
	K3p = K3/K3dp,	// kf3p/kb3p, 1/M
	K5dp = 2.63e-8,	// HbNH3+ dissociation constant, M
	K6dp = 1.91e-8,	// O2HbNH3+ dissociation constant, M
	nhill = 2.7,	// Hill index, unitless
	n0 = nhill-1.0;	// deviation of Hill index or cooperativity from
			// the stochiometry of O2 for each heme site

// Variables those are privately fixed in the model with the standard 
// physiological values (i.e., pH0, pO2, pCO20, DPG0, Temp0)
   private real 
	pO20 = 100.0,		// standard O2 partial pressure in blood, mmHg
	pCO20 = 40.0,		// standard CO2 partial pressure in blood, mmHg
	pH0 = 7.24,		// standard pH in RBCs, unitless
	DPG0 = 4.65e-3,		// standard 2,3-DPG concentration in RBCs, M
	Temp0 = 37.0,		// standard temperature in blood; degC
	fact = 1.0e-6/Wpl,	// a multiplicative factor, M/mmHg
	alphaO20 = fact*1.37,	// solubility of O2 in water at 37 C, M/mmHg
	alphaCO20 = fact*30.7,	// solubility of CO2 in water at 37 C, M/mmHg
	O20 = alphaO20*pO20,	// standard O2 concentration in RBCs, M
	CO20 = alphaCO20*pCO20,	// standard CO2 concentration in RBCs, M
	Hp0 = 10^(-pH0),	// standard H+ concentration in RBCs, M
	pHpl0 = pH0-log(Rrbc),	// standard pH in plasma, unitless
	P500 = 26.8,		// standard pO2 at 50% SHbO2, mmHg
	C500 = alphaO20*P500;	// standard O2 concentration at 50% SHbO2, M

// Variables those are varied under various physiological conditions and
// those are computed as functions of pH, pCO2, DPG and Temp
   real
	Hct = 0.45,	// hematocrit, unitless
	pO2 = 100.0,	// O2 partial pressure in blood, mmHg
	pH = 7.24, 	// pH in RBCs, unitless
	DPG = 4.65e-3,	// 2,3-DPG concentration in blood, M
	Temp = 37.0,	// temperature in blood; degC
	Wbl,		// fractional water space in blood, unitless
	alphaO2,	// solubility of O2 in water, M/mmHg
	alphaCO2,	// solubility of CO2 in water, M/mmHg
	pK1,		// pK of CO2 hydration reaction, unitless
	K1,		// equilibrium constant of CO2 hydration reaction, M
	pHpl,		// pH in RBCs:pHpl = pHrbc-log(Rrbc), unitless
	Hppl,		// H+ concentration in RBCs:Hppl = Rrbc*Hprbc, M
	Hp,		// H+ concentration in RBCs, M
	O2,		// O2 concentration in RBCs, M
	P501,		// pO2 at 50% SHbO2 as a function of pH, mmHg
	P502(pCO2),	// pO2 at 50% SHbO2 as a function of pCO2, mmHg
	P503,		// pO2 at 50% SHbO2 as a function of DPG, mmHg
	P504,		// pO2 at 50% SHbO2 as a function of Temp, mmHg
	C501,		// O2 concentration at 50% SHbO2 as function of pH, M
	C502(pCO2),	// O2 concentration at 50% SHbO2 as function of pCO2, M
	C503,		// O2 concentration at 50% SHbO2 as function of DPG, M
	C504,		// O2 concentration at 50% SHbO2 as function of Temp, M
	K4dp,		// proportionality constant for K4p, 1/M
	n1,		// index on Hp0/Hp as a function of pH, unitless
	n2(pCO2),	// index on CO20/CO2 as a function of pCO2, unitless
	n3,		// index on DPG0/DPG as a function of DPG, unitless
	n4;		// index on Temp0/Temp as a function of Temp, unitless

// Other intermediate variables used in the computations of state variables
   real pCO2diff(pCO2), pHdiff, DPGdiff, Tempdiff, pHpldiff, Term1, Term2, 
	Term3, Term4, Term5(pCO2), Term10, Term20, Term30, Term40, K4tp, 
	Kratio10, Kratio11, Kratio12(pCO2), Kratio20, Kratio21, Kratio22(pCO2);

// Calculation of all the intermediate variables and the indices n1, n2, n3 
// and n4 in the computations of SHbO2 & SHbCO2
	Wbl = (1-Hct)*Wpl + Hct*Wrbc;
	pHpl = pH-log(Rrbc);
	pHpldiff = pHpl-pHpl0;
	pHdiff = pH-pH0;
	pCO2diff = pCO2-pCO20;
	DPGdiff = DPG-DPG0;
	Tempdiff = Temp-Temp0;
	alphaO2 = fact*(1.37 - 0.0137*Tempdiff + 0.00058*Tempdiff^2);
	alphaCO2 = fact*(30.7 - 0.57*Tempdiff + 0.02*Tempdiff^2);
	pK1 = 6.091 - 0.0434*pHpldiff + 0.0014*Tempdiff*pHpldiff;
	K1 = 10^(-pK1);
	O2 = alphaO2*pO2;
	CO2 = alphaCO2*pCO2;
	Hp = 10^(-pH);
	Hppl = 10^(-pHpl);

	Term1 = K2p*(1+K2dp/Hp);
	Term2 = K3p*(1+K3dp/Hp);
	Term3 = (1+Hp/K5dp);
	Term4 = (1+Hp/K6dp);
	Term10 = K2p*(1+K2dp/Hp0);
	Term20 = K3p*(1+K3dp/Hp0);
	Term30 = (1+Hp0/K5dp);
	Term40 = (1+Hp0/K6dp);
	Term5 = (Hp0/Hp)^n1*(CO20/CO2)^n2*(DPG0/DPG)^n3*(Temp0/Temp)^n4;
	Kratio10 = (Term10*CO20+Term30)/(Term20*CO20+Term40);
	Kratio11 = (Term1*CO20+Term3)/(Term2*CO20+Term4);
	Kratio12 = (Term10*alphaCO20*pCO2+Term30)/
		   (Term20*alphaCO20*pCO2+Term40);
	K4dp = Kratio10*O20^n0/C500^nhill;
	K4tp = K4dp/O20^n0;
	Kratio20 = Kratio10/K4tp;
	Kratio21 = Kratio11/K4tp;
	Kratio22 = Kratio12/K4tp;

	P501 = 26.765 - 21.279*pHdiff + 8.872*pHdiff^2;
	P502 = 26.80 + 0.0428*pCO2diff + 3.64e-5*pCO2diff^2;
	P503 = 26.78 + 795.633533*DPGdiff - 19660.8947*DPGdiff^2;
	P504 = 26.75 + 1.4945*Tempdiff + 0.04335*Tempdiff^2 
	     + 0.0007*Tempdiff^3;
	C501 = alphaO20*P501;
	C502 = alphaO20*P502;
	C503 = alphaO20*P503;
	C504 = alphaO2*P504;
	n1 = if (abs(pH-pH0) < 1.0e-6) 1.0  	// can be any arbitrary value
	     else (log(Kratio21)-nhill*log(C501))/(pH-pH0);
	n2 = if (abs(pCO2-pCO20) < 1.0e-6) 1.0 	// can be any arbitrary value
	     else (log(Kratio22)-nhill*log(C502))/(log(CO20)-log(CO2));
	n3 = if (abs(DPG-DPG0) < 1.0e-6) 1.0 	// can be any arbitrary value
	     else (log(Kratio20)-nhill*log(C503))/(log(DPG0)-log(DPG));
	n4 = if (abs(Temp-Temp0) < 1.0e-6) 1.0	// can be any arbitrary value
	     else (log(Kratio20)-nhill*log(C504))/(log(Temp0)-log(Temp));

// Calculation of Hill coefficients (KHbO2 and KHbCO2), O2 and CO2 saturations
// of hemoglobin (SHbO2 and SHbCO2), and O2 and CO2 contents in blood. These 
// are computed as functions of pO2 and pCO2. Also compute the concentrations
// of all the components of HbO2 and HbCO2.
	K4p = K4dp*(O2/O20)^n0*Term5;
	KHbO2 = K4p*(Term2*CO2+Term4)/(Term1*CO2+Term3);
	KHbCO2 = (Term1+Term2*K4p*O2)/(Term3+Term4*K4p*O2);
	SHbO2 = KHbO2*O2/(1+KHbO2*O2);
	SHbCO2 = KHbCO2*CO2/(1+KHbCO2*CO2);
	O2free = Wbl*alphaO2*pO2;
	O2bound = 4*Hct*Hbrbc*SHbO2;
	O2total = O2free+O2bound;
	O2content = 2225.6*O2total;
	CO2free = Wbl*alphaCO2*pCO2;
	CO2bicarb = ((1-Hct)*Wpl+Hct*Wrbc*Rrbc)*(K1*CO2/Hppl);
	CO2bound = 4*Hct*Hbrbc*SHbCO2;
	CO2total = CO2free+CO2bicarb+CO2bound;
	CO2content = 2225.6*CO2total;
	HbNH2 = Hbrbc/((Term1*CO2+Term3)+ K4p*O2*(Term2*CO2+Term4));
	HbNH3p = HbNH2*Hp/K5dp;
	O2HbNH2 = K4p*O2*HbNH2;
	O2HbNH3p = O2HbNH2*Hp/K6dp;
	HbNHCOOH = K2p*CO2*HbNH2;
	HbNHCOOm = K2dp*HbNHCOOH/Hp;
	O2HbNHCOOH = K3p*CO2*O2HbNH2;
	O2HbNHCOOm = K3dp*O2HbNHCOOH/Hp;
	SHbO2kin = (O2HbNH2+O2HbNH3p+O2HbNHCOOH+O2HbNHCOOm)/Hbrbc;
	SHbCO2kin = (HbNHCOOH+HbNHCOOm+O2HbNHCOOH+O2HbNHCOOm)/Hbrbc;

	O2freepl1 = Wpl*(1-Hct)*alphaO2*pO2;
	O2freepl2 = 2225.6*O2freepl1;
	O2freerbc1 = Wrbc*Hct*alphaO2*pO2;
	O2freerbc2 = 2225.6*O2freerbc1;
	O2boundrbc1 = 4*Hct*Hbrbc*SHbO2;
	O2boundrbc2 = 2225.6*O2boundrbc1;
	CO2freepl1 = Wpl*(1-Hct)*alphaCO2*pCO2;
	CO2freepl2 = 2225.6*CO2freepl1;
	CO2freerbc1 = Wrbc*Hct*alphaCO2*pCO2;
	CO2freerbc2 = 2225.6*CO2freerbc1;
	CO2bicarbpl1 = Wpl*(1-Hct)*K1*alphaCO2*pCO2/Hppl;
	CO2bicarbpl2 = 2225.6*CO2bicarbpl1;
	CO2bicarbrbc1 = Wrbc*Hct*Rrbc*K1*alphaCO2*pCO2/Hppl;
	CO2bicarbrbc2 = 2225.6*CO2bicarbrbc1;
	CO2boundrbc1 = 4*Hct*Hbrbc*SHbCO2;
	CO2boundrbc2 = 2225.6*CO2boundrbc1;
}
//-----------------------------------------------------------------------------
// The equations for O2 and CO2 saturations of hemoglobin (SHbO2 and SHbCO2)  
// are derived by considering the various kinetic reactions involving the
// binding of O2 and CO2 with hemoglobin in RBCs:
//
//            kf1p       K1dp
// 1. CO2+H2O <--> H2CO3 <--> HCO3- + H+;  K1=(kf1p/kb1p)*K1dp
//            kb1p		K1 = 7.43e-7 M, K1dp = 5.5e-4 M
//
//              kf2p          K2dp
// 2. CO2+HbNH2 <--> HbNHCOOH <--> HbNHCOO- + H+;  K2=(kf2p/kb2p)*K2dp
//              kb2p		K2 = 2.95e-5, K2dp = 1.0e-6 M
//
//                kf3p            K3dp
// 3. CO2+O2HbNH2 <--> O2HbNHCOOH <--> O2HbNHCOO- + H+; K3=(kf3p/kb3p)*K3dp
//                kb3p		K3 = 2.51e-5, K3dp = 1.0e-6 M
//
//              kf4p          
// 4. O2+HbNH2 <--> O2HbNH2;  K4p=K4dp*func([O2],[H+],[CO2],[DPG],T)
//              kb4p		K4dp and K4p are to be determined
//
//    func = ([O2]/[O2]s)^n0*([H+]s/[H+])^n1*([CO2]s/[CO2])^n2*
//           ([DPG]s/[DPG])^n3*(Temps/Temp)^n4
//
//           K5dp
// 5. HbNH3+ <--> HbNH2 + H+;  K5 = 2.63e-8 M
//
//             K6dp
// 6. O2HbNH3+ <--> O2HbNH2 + H+;  K6 = 1.91e-9 M
//
// The association and dissociation rate constants of O2 with hemoglobin is
// assumed to be dependent on [O2], [H+], [CO2], [DPG] and temperature (Temp)
// such that the equilibrium constant K4p is proportional to ([O2]/[O2]s)^n0,
// ([H+]s/[H+])^n1, ([CO2]s/[CO2])^n2, ([DPG]s/[DPG])^n3, and (Temps/Temp)^n4.
// The problem is to estimate the values of the proportionality constant K4dp 
// and the indices n0, n1, n2, n3 and n4 such that SHbO2 is 50% at pO2 = 26.8 
// mmHg, pH = 7.24, pCO2 = 40 mmHg, [DPG] = 4.65 mM and Temp = 37 C in RBCs 
// and the HbO2 dissociation curve shifts appropriately w.r.t. pH and pCO2.
//-----------------------------------------------------------------------------