Functions used in parameter setting

double setkvals1(double K_eq, double halftime)
{
   return log(2.0)/(halftime*(1.0 + 1.0/K_eq));
}

//////////

double setkvals2(double K_eq, double halftime)
{
   return log(2.0)/(halftime*(1.0 + K_eq));
}

/////////

double setkvals1acid(double K_eq, double halftime, double normdiss)
{
   // here normdiss is the normal H+ concentration
   return log(2.0)/(halftime*(1.0 + normdiss/K_eq));
}

/////////

double setkvals2acid(double K_eq, double halftime, double normdiss)
{
   return log(2.0)/(halftime*(normdiss + K_eq));
}

/***************************** SPECIFIC FUNCTIONS *******************/
double k_ADPATP_eq_fcn()
{
   return k_ADPATP/k_nADPATP;
}

///////////

double _AMPn_fcn()
{
// at equilibrium
   return k_ADPATP_eq*pow(_ADPn, 2)/_ATPn;
}

///////////

double k_AMPAd_eq_fcn()
{
   return _Adn/_AMPn;
}

///////////

double k_AMPAd_fcn()
{
   return log(2.0)/((1.0 + 1.0/k_AMPAd_eq)*k_AMPAd_ht);
}

///////////

double k_nAMPAd_fcn()
{
   return log(2.0)/((1.0 + k_AMPAd_eq)*k_AMPAd_ht);
}

///////////

double _mADPn_fcn()
{
   double rat;
   rat = _ATPn/_ADPn/cmphos_ratio;
   return _mATPADPpool/(1.0 + rat);
}

///////////

double _mATPn_fcn()
{
   double rat;
   rat = _ATPn/_ADPn/cmphos_ratio;
   return rat*_mATPADPpool/(1.0 + rat);
}

///////////

double Vmax_ANT_fcn()
{
   double temp;
   temp = (30.0*CMRglucn + 15.0*CMRlacn)*Vol_BrW/(R_ATPshuttle1);
     if(temp <= 0.0){
     fprintf(stderr, "ERROR in Vmax_ANT_fcn: maximal rate is negative.\n");
     return 0.0;
   }
   else
     return temp;
}

///////////

double Km_BiC_fcn()
{
   return 1.0*_eBiCn;
}

///////////

double Vmax1_BiC_fcn()
{
   return _eBiCn*(R_csfout)/(R_BEdiffBiC1)/(1.0 - Rbf_BiC);
}

///////////

double Vmax2_BiC_fcn()
{
   return Rbf_BiC*Vmax1_BiC*(R_BEdiffBiC1)/(R_BEdiffBiC2);
}

///////////

double k_CCn_fcn()
{
   double temp;
   temp = netCO2prodn*Vol_BrW/(_eCO2n - CO2n);
   if (temp <= 0){
     fprintf(stderr, "ERROR in k_CCn_fcn - extracellular CO2 less that blood CO2");
     return 0.0;
   }
   else
     return temp;
}

///////////

double Km_KK_fcn()
{
   return 1.0*_eKn;
}

///////////

double Vmax1_KK_fcn()
{
   return _eKn*(R_csfout)/(R_BEdiffK1)/(1.0 - Rbf_KK);
}

///////////

double Vmax2_KK_fcn()
{
   return Rbf_KK*Vmax1_KK*(R_BEdiffK1)/(R_BEdiffK2);
}

///////////

double k_Na_fcn()
{
   return log(2.0)/k_Na_ht/(1.0/Vol_c + 1.0/Vol_ecs);
}

///////////

double k_OOn_fcn()
{
   return 1.1*fabs(CMRO2n*Vol_BrW/(O2n - _eO2n));
}

///////////

double k_BiCshut_eq_fcn()
{
   return _mBiCn/_BiCn;
}

///////////

double k_BiCshut_fcn()
{
   return log(2.0)/k_BiCshut_ht/(1.0/Vol_exm + 1.0/k_BiCshut_eq/Vol_inm);
}

///////////

double k_nBiCshut_fcn()
{
   return log(2.0)/k_BiCshut_ht/(k_BiCshut_eq/Vol_exm + 1.0/Vol_inm);
}

///////////

double k_Pbuf_fcn()
{
   return setkvals1acid(k_Pbuf_eq, k_Pbuf_ht, _Hyn);
}

///////////

double k_nPbuf_fcn()
{
   return setkvals2acid(k_Pbuf_eq, k_Pbuf_ht, _Hyn);
}

///////////

double k_Pbuf_eq_fcn()
{
   return Pbufn*_Hyn/PbufHn;
}

//////////

double k_Cain1n_fcn()
{
   // see section secCadyn
   return Vmax_Caout1*Ca_i1n/(Km_Caout + Ca_i1n);
}

///////////

double k_Cain2n_fcn()
{
   // see section secCadyn
   return Vmax_Caout2*Ca_i2n/(Km_Caout + Ca_i2n);
}

///////////

double Vmax_Caout1_fcn()
{
// see section secCadyn
   return (log(2.0)/k_Ca_ht)*(Ca_i1n + Km_Caout);
}

///////////

double Vmax_Caout2_fcn()
{
   // see section secCadyn
   return (log(2.0)/k_Ca_ht)*(Ca_i2n + Km_Caout);
}

///////////

double Km_Caout_fcn()
{
   // change in tandem with k_Ca_ht
   return Km_Caout_ratio*Ca_i1n;
}

///////////

double k_CY_eq_fcn()
{
   // from the normal arterial proportions
   return CY_an/(CO2_an*Y_an);
}

///////////

double k_CY_fcn()
{
   double Ptemp;
   double Qtemp;
   Ptemp = 15.0;
   // very approximate value of CO2n + Yn
   Qtemp = 1.0;
   return log(2.0)/((Ptemp + Qtemp/k_CY_eq)*k_CY_ht);
}

///////////

double k_nCY_fcn()
{
   double Ptemp;
   double Qtemp;
   Ptemp = 15.0;
   // very approximate value of CO2n + Yn
   Qtemp = 1.0;
   return log(2.0)/((k_CY_eq*Ptemp + Qtemp)*k_CY_ht);
}

///////////

double Y_tot_fcn()
{
   // the four terminal amino groups of haemoglobin
   return 4.0*X_tot;
}

///////////

double CY_an_fcn()
{
   // assuming a certain proportion of the dissolved CO2. Perhaps a little more than twice as much.
   return 2.5*CO2_an;
}

///////////

double CY_a_fcn()
{
// from the equilibrium
   return Y_tot*k_CY_eq*CO2_a/(1.0 + k_CY_eq*CO2_a);
}

///////////

double Y_a_fcn()
{
// from conservation
   return Y_tot - CY_a;
}

// To generate the initial conditions

///////////

double Y_an_fcn()
{
   // not combined with CO2
   double temp;
   temp = Y_tot - CY_an;
   if (temp <= 0){
     fprintf(stderr, "ERROR in Y_an_fcn - Y_an <= 0\n");
     return 0.00001;
     // so that we don't get an error in k_CY_eq
   }
   else
     return temp;
}

///////////

double Yn_fcn()
{
   // From conservation
   return Y_tot - CYn;
}

///////////

double CYn_fcn()
{
   // By conservation and setting dCY/dt = 0
   return (k_CY*CO2n*Y_tot + 2.0*qn/Vol_c*CY_an)/(k_CY*CO2n + k_nCY + 2.0*qn/Vol_c);
}

///////////

double k_CO2shut_fcn()
{
   double tempval;
   tempval = _mCO2n - _CO2n;
   if (tempval <= 0){
     fprintf(stderr, "ERROR in k_CO2shut_fcn, mitochondrial CO2 less than cytoplasmic.\n");
     return 0.0;
   }
   else
     return netCO2prodn*Vol_BrW/tempval;
}

///////////

double k_CH_fcn()
{
   return setkvals1acid(k_CH_eq, k_CH_ht, Hyn);
}

///////////

double k_nCH_fcn()
{
   return setkvals2acid(k_CH_eq, k_CH_ht, Hyn);
}

///////////

double k_CHe_fcn()
{
// the third quantity is the approximate bicarbonate concentration
   return setkvals1acid(k_CH_eq, k_CHe_ht, _eHyn);
}

///////////

double k_nCHe_fcn()
{
// the third quantity is the approximate bicarbonate concentration
   return setkvals2acid(k_CH_eq, k_CHe_ht, _eHyn);
}

///////////

double k_CHi_fcn()
{
// the third quantity is the approximate bicarbonate concentration
   return setkvals1acid(k_CH_eq, k_CHi_ht, _Hyn);
}

///////////

double k_nCHi_fcn()
{
// the third quantity is the approximate bicarbonate concentration
   return setkvals2acid(k_CH_eq, k_CHi_ht, _Hyn);
}

///////////

double k_CHm_fcn()
{
// the third quantity is the approximate bicarbonate concentration
   return setkvals1acid(k_CH_eq, k_CHm_ht, _mHyn);
}

///////////

double k_nCHm_fcn()
{
// the third quantity is the approximate bicarbonate concentration
   return setkvals2acid(k_CH_eq, k_CHm_ht, _mHyn);
}

///////////

double VmaxC_oxp_fcn()
{
   double temp;
   temp = Inh_oxp*Vol_BrW*(V_oxpn)/(R_Cortoxphos);
// temp = Inh_oxp*Vol_BrW*((1.0-Hleak_frac)*(2.5*N1_H + 0.5*N2_H - 1.0)*(2.0*CMRglucn + CMRlacn)/4.0)/(R_Cortoxphos);
   if(temp >=0.0){
     fprintf(stderr, "ERROR in VmaxC_oxp_fcn: normal rate has the wrong sign.\n");
     return 0.0;
   }
   else
     return temp;
}

///////////

double k_ecAd_fcn()
{
   return log(2.0)/k_ecAd_ht/(1.0/Vol_ecs + 1.0/Vol_exm);
}

///////////

double _eAdn_fcn()
{
   // because extracellular adenosine currently does not escape.
   return _Adn;
}

///////////

double k_ecCC_fcn()
{
   // based on the assumption that transport of BiC between ecs and ics is minimal.
   double temp;
   temp = netCO2prodn*Vol_BrW/(_CO2n - _eCO2n);
   if (temp <= 0){
     fprintf(stderr, "ERROR in k_ecCC_fcn: intracellular CO2 less than extracellular");
     return 0;
   }
   else
     return temp;
}

///////////

double k_ecOO_fcn()
{
   // conservation
   return CMRO2n*Vol_BrW/(_eO2n - _O2n);
   //log(2.0)/k_ecOO_ht/(1.0/Vol_ecs + 1.0/Vol_exm);
}

///////////

double k_ecPy0_fcn()
{
   return log(2.0)/k_ecPy0_ht/(1.0/Vol_ecs + 1.0/Vol_exm);
}

///////////

double gCa1_tot_fcn()
{
   // the total possible Calcium conductivity - entirely determined by what we want the normal conductivity to be and the normal probability to be
   double temp;
   temp = 1.0 + exp((V_Ca - PD_mem1)/kV_Ca);
   return gCa1n*temp;
}

///////////

double gCa1n_fcn()
{
   // normal Calcium conductivity determined by normal membrane potential and the normal requirement for inflow. Assumed to be the same for both segments in the resting state...
   return k_Cain1n/(PD_Ca1n - PD_mem1n);
}

///////////

double gCa2n_fcn()
{
   // normal Calcium conductivity determined by normal membrane potential and the normal requirement for inflow. Assumed to be the same for both segments in the resting state...
   return k_Cain2n/(PD_Ca2n - PD_mem2n);
}

///////////

double gCa2_tot_fcn()
{
   // the total possible Calcium conductivity - entirely determined by what we want the normal conductivity to be and the normal probability to be
   double temp;
   temp = 1.0 + exp((V_Ca - PD_mem2)/kV_Ca);
   return gCa2n*temp;
}

///////////

double gKATP1n_fcn()
{
   // just a guess - in faraci and heistad it is stated that these channels do not really affect basal tone
   double tempval;
   tempval = gKATP1n_frac*gK1n;
   if (tempval > 1.0 - gKpres1n){
     fprintf(stderr, "ERROR in gKATP1n_fcn()\n");
     return 0.0;
     }
   else
     return tempval;
}

///////////

double gKATP2n_fcn()
{
   // just a guess - in faraci and heistad it is stated that these channels do not really affect basal tone
   double tempval;
   tempval = gKATP2n_frac*gK2n;
   if (tempval > 1.0 - gKpres2n){
     fprintf(stderr, "error in gKATP2n_fcn()\n");
     return 0.0;
     }
   else
     return tempval;
}

///////////

double gKATP1tot_fcn()
{
   // saturating function
   return gKATP1n/gKATPfracn;
}

///////////

double gKATP2tot_fcn()
{
   return gKATP2n/gKATPfracn;
}

///////////

double m_pH_fcn()
{
   double tempval;
   tempval = gKATPfracn/(1.0 - gKATPfracn);
   return m_pH_frac*tempval;
}

///////////

double m_Ad_fcn()
{
   double tempval;
   tempval = gKATPfracn/(1.0 - gKATPfracn);
   return m_Ad_frac*tempval;
}

///////////

double c_pHAd_fcn()
{
   double tempval;
   tempval = (1.0 - m_pH_frac - m_Ad_frac)*gKATPfracn/(1.0 - gKATPfracn);
     if(tempval < 0){
       fprintf(stderr, "ERROR in c_pHAd_fcn\n");
       return 0;
     }
     else
       return tempval;
}

///////////

double gK1n_fcn()
{
   // solve for normal physiological gK1 as a function of normal membrane potential. Note that this is reduced because of the normal pressure
   return (gNa*(PD_mem1n - PD_Na) + gCl*(PD_mem1n - PD_Cl))/(PD_Kn - PD_mem1n);

}

///////////

double gK2n_fcn()
{
   // solve for normal physiological gK2 as a function of normal membrane potential. Note that this is reduced because of the normal pressure
   return (gNa*(PD_mem2n - PD_Na) + gCl*(PD_mem2n - PD_Cl))/(PD_Kn - PD_mem2n);

}

///////////

double gKpot1n_fcn()
{
   double tempval;
   tempval = gK1n - gKpres1n - gKATP1n - gKV1n;
   if (tempval < 0){
     fprintf(stderr, "error in gKpot1n_fcn()\n");
     return 0.0;
   }
   else
     return tempval;
}

///////////

double gKpot2n_fcn()
{
   double tempval;
   tempval = gK2n - gKpres2n - gKATP2n - gKV2n;
   if (tempval < 0){
     fprintf(stderr, "error in gKpot2n_fcn()\n");
     return 0.0;
   }
   else
     return tempval;
}

///////////

double gKpot1tot_fcn()
{
   // saturating function
   return gKpot1n/gKpotfracn;
}

///////////

double gKpot2tot_fcn()
{
   return gKpot2n/gKpotfracn;
}

///////////

double m_K_fcn()
{
   return gKpotfracn/(1.0 - gKpotfracn) - c_K;
}

///////////

double c_K_fcn()
{
   return c_K_frac*gKpotfracn/(1.0 - gKpotfracn);
}

///////////

double gKpres1n_fcn()
{
   // about half of normal conductivity is of pressure sensitive channels?
   return gKpres1n_frac*gK1n;
}

///////////

double gKpres2n_fcn()
{
   // a little less than half of normal conductivity is of pressure sensitive channels? Just a guess...
   return gKpres2n_frac*gK2n;
}

///////////

double gKpres1tot_fcn()
{
   // the total possible pressure sensitive potassium channel conductivity
   return gKpres1n/gKpresfracn;
}

///////////

double gKpres2tot_fcn()
{
   return gKpres2n/gKpresfracn;
}

///////////

double c_pres_fcn()
{
   return 1.0/gKpresfracn - 1.0;
}

///////////

double c_V_fcn()
{
   double tempval;
   tempval = 1.0 - m_NO - m_Ca;
   if (tempval<0.0){
     fprintf(stderr, "ERROR in c_V_fcn.\n");
     return 0.0;
   }
   else
     return tempval;
}

///////////

double gKV1n_fcn()
{
   // just a guess - contribution to basal tone
   double tempval;
   tempval = gKV1n_frac*gK1n;
   if (tempval > 1.0 - gKpres1n - gKATP1n){
     fprintf(stderr, "error in gKV1n_fcn()\n");
     return 0.0;
     }
   else
     return tempval;
}

///////////

double gKV2n_fcn()
{
   // just a guess
   double tempval;
   tempval = gKV2n_frac*gK2n;
   if (tempval > 1.0 - gKpres2n - gKATP2n){
     fprintf(stderr, "error in gKV2n_fcn()\n");
     return 0.0;
   }
   else
     return tempval;
}

///////////

double gKV1tot_fcn()
{
   return gKV1n/(1.0/(1.0 + exp((V_act - PD_mem1)/kV_act)))/(1.0/(1.0 + exp(-(V_inact - PD_mem1)/kV_inact)));
}

///////////

double gKV2tot_fcn()
{
   return gKV2n/(1.0/(1.0 + exp((V_act - PD_mem2)/kV_act)))/(1.0/(1.0 + exp(-(V_inact - PD_mem2)/kV_inact)));
}

///////////

double Vmax_gluc_fcn()
{
   // total transported into tissue = total used by tissue
   return CMRglucn*Vol_BrW*((_eglucn + K_gluc*(1.0 + p_gluc))*(glucn + K_gluc*(1.0 + p_gluc)) - pow(p_gluc, 2)*pow(K_gluc, 2))/(glucn - _eglucn);
}

///////////

double K_gluc_fcn()
{
   return Km_gluc*(1.0 + p_gluc)/(1.0 + 2*p_gluc);
}

///////////

double k_ecgluc_fcn()
{
   return CMRglucn*Vol_BrW/(_eglucn - _glucn);
}

///////////

double Km_glucNN_fcn()
{
// from Brenda for the enzyme Glyceraldehyde-3-P dehydrogenase
   return 1.0;
}

////////////

double Km_glucNA_fcn()
{
   return 0.2*_ADPn;
}

///////////

double Km_glucNP_fcn()
{
   return 0.2*_Phosn;
}

///////////

double Vmax_glucNn_fcn()
{
   return Vol_BrW/Vol_exm*CMRglucn/(R_glycol);
}

//////////

double k_Bohr_fcn()
{
   // see sechaemO2
   return k_XO4H_eq/k_XH_eq;
}

///////////

double k_XO4_eq_fcn()
{
   // see sechaemO2
   return (1.0 + k_XH_eq*3.981e-5)/(pow((26.0*henry_O2), 2.5)*(1.0 + k_Bohr*k_XH_eq*3.981e-5));
}

///////////

double k_XO4_fcn()
{
   double Ptemp;
   double Qtemp;
   Ptemp = pow(O2n, 1.5)*(O2n + 10.0*Xn);
   Qtemp = 1.0;
   return log(2.0)/((Ptemp + Qtemp/k_XO4_eq)*k_XO4_ht);
}

///////////

double k_nXO4_fcn()
{
   double Ptemp;
   double Qtemp;
   Ptemp = pow(O2n, 1.5)*(O2n + 10.0*Xn);
   Qtemp = 1.0;
   return log(2.0)/((k_XO4_eq*Ptemp + Qtemp)*k_XO4_ht);
}

///////////

double k_XH_eq_fcn()
{
   // although on p7 of Bhagavan, a pK value of about 6.5 is given, at http://www.hmc.psu.edu/scaduto/springteaching/htmlslides/page7.htm a pK value of 7.9 is given which we use.
   return 1.0/0.000013;
}

///////////

double k_XH_fcn()
{
   double Ptemp;
   double Qtemp;
   Ptemp = n_H*Xn;
   Qtemp = 1.0;
   return log(2.0)/((Ptemp + Qtemp/k_XH_eq)*k_XH_ht);
}

///////////

double k_nXH_fcn()
{
   double Ptemp;
   double Qtemp;
   Ptemp = n_H*Xn;
   Qtemp = 1.0;
   return log(2.0)/((k_XH_eq*Ptemp + Qtemp)*k_XH_ht);
}

///////////

double k_XHO_eq_fcn()
{
   // see sechaemO2
   return k_Bohr*k_XO4_eq;
}

///////////

double k_XHO_fcn()
{
   double Ptemp;
   double Qtemp;
   Ptemp = pow(O2n, 1.5)*(O2n + 10.0*X_Hn);
   Qtemp = 1.0;
   return log(2.0)/((Ptemp + Qtemp/k_XHO_eq)*k_XHO_ht);
}

///////////

double k_nXHO_fcn()
{
   double Ptemp;
   double Qtemp;
   Ptemp = pow(O2n, 1.5)*(O2n + 10.0*X_Hn);
   Qtemp = 1.0;
   return log(2.0)/((k_XHO_eq*Ptemp + Qtemp)*k_XHO_ht);
}

///////////

double k_XO4H_eq_fcn()
{
// at http://www.hmc.psu.edu/scaduto/springteaching/htmlslides/page7.htm a pK value of 6.7 is given which we use
   return 1.0/0.0002;
}

///////////

double k_XO4H_fcn()
{
   double Ptemp;
   double Qtemp;
   Ptemp = n_H*X_O4n;
   Qtemp = 1.0;
   return log(2.0)/((Ptemp + Qtemp/k_XO4H_eq)*k_XO4H_ht);
}

///////////

double k_nXO4H_fcn()
{
   double Ptemp;
   double Qtemp;
   Ptemp = n_H*X_O4n;
   Qtemp = 1.0;
   return log(2.0)/((k_XO4H_eq*Ptemp + Qtemp)*k_XO4H_ht);
}

///////////

double X_a_fcn()
{
// from conservation of haemoglobin and the equilibria
   return X_tot/(1.0 + k_XH_eq*Hy_a + pow(O2_a, 2.5)*k_XO4_eq + Hy_a*pow(O2_a, 2.5)*k_XH_eq*k_XHO_eq);
}

///////////

double X_an_fcn()
{
   // from conservation of haemoglobin and the equilibria
   return X_totn/(1.0 + k_XH_eq*Hy_an + pow(O2_an, 2.5)*k_XO4_eq + Hy_an*pow(O2_an, 2.5)*k_XH_eq*k_XHO_eq);
}

///////////

double Xn_fcn()
{
// conservation of haemoglobin and the equilibria
   return X_totn/(1.0 + k_XH_eq*Hyn + pow(O2n, 2.5)*k_XO4_eq + Hyn*pow(O2n, 2.5)*k_XH_eq*k_XHO_eq);
}

///////////

double X_O4_a_fcn()
{
   return k_XO4_eq*X_a*pow(O2_a, 2.5);
}

///////////

double X_O4_an_fcn()
{
   return k_XO4_eq*X_an*pow(O2_an, 2.5);
}

///////////

double X_O4n_fcn()
{
   return k_XO4_eq*Xn*pow(O2n, 2.5);
}

///////////

double X_O4_H_a_fcn()
{
   return k_XHO_eq*X_H_a*pow(O2_a, 2.5);
}

///////////

double X_O4_H_an_fcn()
{
   return k_XHO_eq*X_H_an*pow(O2_an, 2.5);
}

///////////

double X_O4_Hn_fcn()
{
   return k_XHO_eq*X_Hn*pow(O2n, 2.5);
}

///////////

double X_H_a_fcn()
{
   return k_XH_eq*X_a*Hy_a;
}

///////////

double X_H_an_fcn()
{
   return k_XH_eq*X_an*Hy_an;
}

///////////

double X_Hn_fcn()
{
   return k_XH_eq*Xn*Hyn;
}

///////////

double SaO2_fcn()
{
   return (X_O4_a + X_O4_H_a)/(X_O4_a + X_O4_H_a + X_a + X_H_a);
}

///////////

double SaO2n_fcn()
{
   return (X_O4_an + X_O4_H_an)/(X_O4_an + X_O4_H_an + X_an + X_H_an);
}

///////////

double O2n_fcn()
{
   // This is needed to calculate k_OOn. It is set heuristically at the moment, but is very important and can make a big difference.
   return 0.068;
   //return 0.4*(O2_an + O2_vn);
}

///////////

double Hyn_fcn()
{
   // very tentative! simple exponential function
   return sqrt(Hy_an*Hy_vn);
}

///////////

double Vmax_PyN_fcn()
{
   return Vol_BrW/Vol_inm*netPyusen/(R_Krebs1);
}

///////////

double Km_PyNP_fcn()
{
   return 0.1*_mPy0n;
}

///////////

double Km_PyNN_fcn()
{
   return 0.1*_mNADn;
}

///////////

double Vmax_AKGN_fcn()
{
   // note that the volume terms are important.
   return Vol_BrW/Vol_inm*netPyusen/(R_Krebs2);
}

///////////

double Km_AKGNA_fcn()
{
   return 3.0*_AKGn;
}

///////////

double Km_AKGNN_fcn()
{
   // note that the volume terms are important.
   return 3.0*_mNADn;
}

///////////

double Km_AKGNAD_fcn()
{
   return 0.1*_mADPn;
}

///////////

double Km_AKGNP_fcn()
{
   return 0.2*_mPhosn;
}

///////////

double Vmax_sucN_fcn()
{
   // note that the volume terms are important.
   return Vol_BrW/Vol_inm*netPyusen*(Km_sucNs + _sucn)/_sucn*(Km_sucNN + _mNADn)/_mNADn*(Km_sucNF + _mFADn)/_mFADn;
}

///////////

double Km_sucNs_fcn()
{
   // note that the volume terms are important.
   return 3.0*_sucn;
}

///////////

double Km_sucNN_fcn()
{
   // note that the volume terms are important.
   return 0.05*_mNADn;
}

///////////

double Km_sucNF_fcn()
{
   return 1.0*_mFADn;
}

///////////

double k_LH_fcn()
{
   return setkvals1acid(k_LH_eq, k_LH_ht, _Hyn);
}

///////////

double k_nLH_fcn()
{
   return setkvals2acid(k_LH_eq, k_LH_ht, _Hyn);
}

///////////

double K_L0H_fcn()
{
   return Km_L0*_eHyn*(1.0 + p_L0H)/(1.0 + 2.0*p_L0H);
}

///////////

double Vmax_L0H_fcn()
{
   // net produced by tissue = total transported out
   double tempval1;
   tempval1 = -(1.0 - L_trans_frac)*CMRlacn*Vol_BrW - _eL0n*(R_csfout);
   if(tempval1 < 0.0){
     fprintf(stderr, "ERROR in Vmax_L0H_fcn: net lactate production not sufficient to match outflow in CSF.\n");
     return 0.0;
   }
   else
     return tempval1*((_eL0n*_eHyn + K_L0H*(1.0 + p_L0H))*(L0n*Hyn + K_L0H*(1.0 + p_L0H)) - pow(p_L0H, 2)*pow(K_L0H, 2))/(_eL0n*_eHyn - L0n*Hyn);
}

///////////

double K_ecL0H_fcn()
{
   return Km_L0*_Hyn*(1.0 + p_ecL0H)/(1.0 + 2.0*p_ecL0H);
}

///////////

double Vmax_ecL0H_fcn()
{
   // net produced by tissue = total transported out
   return (1.0 - L_trans_frac)*CMRlacn*Vol_BrW/(R_ECdiffL0H);
}

///////////

double k_ecL_fcn()
{
   return L_trans_frac*CMRlacn*Vol_BrW/(_eLn - _lacn);
}

///////////

double k_LL_fcn()
{
   return L_trans_frac*CMRlacn*Vol_BrW/(Ln - _eLn);
}

///////////

double k_nMAshut_fcn()
{
   return k_MAshut/Keq_MAshut;
}

///////////

double CMRglucn_fcn()
{
   return Mass_Br*CMRglucnw/Vol_BrW;
}

///////////

double CMRlacn_fcn()
{
   return Mass_Br*CMRlacnw/Vol_BrW;
}

///////////

double CMRO2n_fcn()
{
   return 6.0*CMRglucn + 3.0*CMRlacn;
}

///////////

double netATPuse_fcn()
{
// assuming that only a certain proportion of protons re-entering the matrix are involved in the production of ATP
   return (Z_ATPH*N2_H/2.0 + 5.0*Z_ATPH*N1_H/2.0 + 1.0)*CMRlacn + (Z_ATPH*N2_H + 5.0*Z_ATPH*N1_H + 4.0)*CMRglucn;
// return 4.0*CMRglucn + 1.0*CMRlacn + (1.0 - Hleak_frac)*(2.5*N1_H + 0.5*N2_H - 1.0)*(2.0*CMRglucn + CMRlacn)/4.0;
}

///////////

double netPyusen_fcn()
{
   return 2.0*CMRglucn + CMRlacn;
}

///////////

double netCO2prodn_fcn()
{
   return CMRO2n;
}

///////////

double Z_ATPH_fcn()
{
   return ((1.0 - Hleak_frac)*(N2_H + 5.0*N1_H) - 2.0)/(4.0*N2_H + 20.0*N1_H);
}

///////////

double V_oxpn_fcn()
{
   return (Z_ATPH*N2_H/2.0 + 5.0*Z_ATPH*N1_H/2.0)*(CMRlacn + 2.0*CMRglucn);
}

///////////

double V_leakn_fcn()
{
   return Hleak_frac/(1.0 - Hleak_frac)*(2.0*Z_ATPH*N2_H + 10.0*Z_ATPH*N1_H + 1.0)*(CMRlacn + 2.0*CMRglucn);
}

///////////

double k_MLC1inac0_fcn()
{
   return k_MLC1inac0_frac*k_MLC1inacn;
}

///////////

double k_MLC2inac0_fcn()
{
   return k_MLC1inac0_frac*k_MLC1inacn;
}

///////////

double K1_NO1_fcn()
{
   return (1.0 - k_MLC1inac0_frac)*k_MLC1inacn*(1.0 + pow(K2_NO, n_NO));
}

///////////

double K1_NO2_fcn()
{
   return (1.0 - k_MLC2inac0_frac)*k_MLC2inacn*(1.0 + pow(K2_NO, n_NO));
}

///////////

double K1_CaM1_fcn()
{
   return k_MLC1acn*(1.0 + pow(K2_CaM, n_Ca));
}

///////////

double K1_CaM2_fcn()
{
   return k_MLC2acn*(1.0 + pow(K2_CaM, n_Ca));
}

///////////

double MLC1n_fcn()
{
   // see secMyos
   return MLC_tot*MLC1_inac_fracn;
}

///////////

double MLC2n_fcn()
{
   // see secMyos
   return MLC_tot*MLC2_inac_fracn;
}

///////////

double MLCp1n_fcn()
{
   // see secMyos
   return MLC_tot - MLC1n;
}

///////////

double MLCp2n_fcn()
{
   // see secMyos
   return MLC_tot - MLC2n;
}

///////////

double k_MLC1_eq_fcn()
{
   // the normal equilibrium
   return MLC1n/MLCp1n;
}

///////////

double k_MLC2_eq_fcn()
{
   // the normal equilibrium
   return MLC2n/MLCp2n;
}

///////////

double k_MLC1inacn_fcn()
{
   // see secMyos
   return setkvals1(k_MLC1_eq, k_MLC_ht);
}

///////////

double k_MLC1acn_fcn()
{
   // see secMyos
   return setkvals2(k_MLC1_eq, k_MLC_ht);
}

///////////

double k_MLC2inacn_fcn()
{
   // see secMyos
   return setkvals1(k_MLC2_eq, k_MLC_ht);
}

///////////

double k_MLC2acn_fcn()
{
   // see secMyos
   return setkvals2(k_MLC2_eq, k_MLC_ht);
}

///////////

double T_max1n_fcn()
{
   return T_max01*(k_aut1*M_1n + (1.0 - k_aut1));
}

///////////

double T_max2n_fcn()
{
   return T_max02*(k_aut2*M_2n + (1.0 - k_aut2));
}

///////////

double K1_MLC_fcn()
{
// maximum possible force
   return pow(K2_MLC, n_MLC) + 1.0;
}

///////////

double _NADn_fcn()
{
   return _NADpool/(1.0 + 1.0/nadnadh1_ratio);
}

///////////

double _NADHn_fcn()
{
   return _NADn/nadnadh1_ratio;
}

///////////

double _mNADn_fcn()
{
   return _mNADpool/(1.0 + 1.0/nadnadh2_ratio);
}

///////////

double _mNADHn_fcn()
{
   return _mNADn/nadnadh2_ratio;
}

///////////

double k_ecNa_fcn()
{
   return k_ecKK_mod*3.0*Vol_BrW*netATPuse*(1.0 - ATP_frac)/(_eNan - _Nan);
}

///////////

double k_ecKK_fcn()
{
   return k_ecKK_mod*2.0*Vol_BrW*netATPuse*(1.0 - ATP_frac)/(_Kn - _eKn);
}

///////////

double Vmax_KATP_fcn()
{
   return Vol_BrW*netATPuse*(1.0 - ATP_frac)/(R_potpump);
}

///////////

double Km_KATPK_fcn()
{
   return 0.7*(_eKn/_Kn);
}

///////////

double Km_KATPNa_fcn()
{
   return 2.0*(_Nan/_eNan);
}

///////////

double Km_ATPuse_fcn()
{
   return Km_ATPuse_frac*_ATPn;
}

///////////

double Vmax_ATPuse_fcn()
{
   return ATP_frac*netATPuse*Vol_BrW/Vol_exm/(R_ATPuse);
}

///////////

double Km_NQN_fcn()
{
   return 2.0*pow((_mNADHn/_mNADn), 1.0);
}

///////////

double Km_NQH_fcn()
{
// The way mitochondrial NADH levels change as O2 drops depends to a significant level on this parameter. If it is large, then the change in pH outweighs the drop in O2, and the rate of oxidative phosphorylation increases despite the drop in O2. If it is very small, then there is significant shift in NADH levels (rise) but also a significant drop in pH.
   return 1.0*_mHyn/_Hyn;
}

///////////

double Vmax_NQ_fcn()
{
   // not valid for isolated mitochondria
   // to ensure that oxygen is used at the rate it is supplied.
   return Inh_etc*Vol_BrW*5.0/6.0*CMRO2n/(R_Noxy);
}

///////////

double Km_FQF_fcn()
{
   return 2.0*(_mFADHn/_mFADn);
}

///////////

double Km_FQO_fcn()
{
   return 0.05*_mO2n;
}

///////////

double Km_FQH_fcn()
{
   return 1.0*_mHyn/_Hyn;
}

///////////

double Vmax_FQ_fcn()
{
   return Inh_etc*Vol_BrW*1.0/6.0*CMRO2n/(R_Foxy);
}

///////////

double k_NOout1_fcn()
{
   // see secNOdyn
   return log(2.0)/k_NO_ht;
}

///////////

double k_NOout2_fcn()
{
   // see secNOdyn
   return log(2.0)/k_NO_ht;
}

///////////

double k_NOprod1n_fcn()
{
   // see secNOdyn
   return k_NO1_eq*k_NOout1;
}

///////////

double k_NOprod2n_fcn()
{
   // see secNOdyn
   return k_NO2_eq*k_NOout2;
}

///////////

double k_NO1_eq_fcn()
{
   return NO1n;
}

///////////

double k_NO2_eq_fcn()
{
   return NO2n;
}

///////////

double k_O2shut_fcn()
{
// Oxygen shuttling in must equal total O2 usage
   double tempval;
   tempval = _O2n - _mO2n;
   if (tempval < 0){
     fprintf(stderr, "ERROR in k_O2shut_fcn, mitochondrial O2 greater than cytoplasmic.\n");
     return 0.00001;
   }
   else
     return CMRO2n*Vol_BrW/(tempval);
}

///////////

double pH_mn_fcn()
{
   return pH_tn + pH_mc_diff;
}

///////////

double _ePy0n_fcn()
{
   return _Py0n;
}

///////////

double _mPy0n_fcn()
{
// just a guess. The pH gradient contributes to the flow of Py0 across the mitochondrial membrane.
   return 0.8*_Py0n*_Hyn/_mHyn;
}

///////////

double L_a_fcn()
{
   return L0_a*Hy_a/k_LH_eq;
}

///////////

double _Kn_fcn()
{
// rather arbitrary
   return 180.0;
}

///////////

double _eL0n_fcn()
{
// average
   return 0.9*_L0n + 0.1*L0n;
}

///////////

double _eglucn_fcn()
{
// rapidly crosses the membrane?
   return 0.8*_glucn + 0.2*glucn;
}

///////////

double BiC_an_fcn()
{
   return k_CH_eq*CO2_an/Hy_an;
}

///////////

double BiC_a_fcn()
{
   return BiC_an;
}

///////////

double BiCn_fcn()
{
   // from dCO2/dt + dBiC/dt = 0 and the assumption of equilibrium in the blood.
   return 0.95*k_CH_eq*CO2n/Hyn;
}

///////////

double _eLn_fcn()
{
   return _eL0n*_eHyn/k_LH_eq;
}

///////////

double _lacn_fcn()
{
   // the normal concentration of undissociated lactate (at normal pH)
   return _Hyn*_L0n/k_LH_eq;
}

///////////

double Ln_fcn()
{
   // necessarily at equilibrium - L only participates in this one reaction
   return L0n*Hyn/k_LH_eq;
}

///////////

double Hy_a_fcn()
{
// from the equilibrium
   return k_CH_eq*CO2_a/BiC_a;
}

///////////

double _eHyn_fcn()
{
   return 1000.0*pow(10.0, -pH_en);
}

///////////

double _Hyn_fcn()
{
   return 1000.0*pow(10.0, -pH_tn);
}

///////////

double _mHyn_fcn()
{
   return 1000.0*pow(10.0, -pH_mn);
}

///////////

double CO2n_fcn()
{
   // From some basic analysis we can conclude that that BiC levels don't change very much in the blood, so change in CO2 is roughly equal to that added by metabolism. The 0.5 is currently rather randomly added!
   return CO2_an + 0.2*netCO2prodn*Vol_BrW/2.0/qn;
}

///////////

double _eCO2n_fcn()
{
   // this must be less than the tissue values since the primary means of transfer out of the tissue is as CO2 rather than BiC
   return 0.5*_CO2n + 0.5*CO2n;
}

///////////

double _CO2n_fcn()
{
   return PtCO2n*henry_CO2t;
}

///////////

double _mCO2n_fcn()
{
// a little higher than outside the mitochondria?
   return 1.2*_CO2n;
}

///////////

double _eBiCn_fcn()
{
// equilibrium
   return k_CH_eq*_eCO2n/_eHyn;
}

///////////

double _BiCn_fcn()
{
// equilibrium
   return k_CH_eq*_CO2n/_Hyn;
}

///////////

double _mBiCn_fcn()
{
   // equilibrium
   return k_CH_eq*_mCO2n/_mHyn;
}

///////////

double glucn_fcn()
{
// zeroing the dynamics...
   return gluc_an - CMRglucn*Vol_BrW/2.0/qn;
}

///////////

double O2_a_fcn()
{
   return PaO2*henry_O2;
}

///////////

double CO2_an_fcn()
{
   return Pa_CO2n*henry_CO2;
}

///////////

double Hy_an_fcn()
{
   return 1000.0*pow(10.0, -pH_an);
}

///////////

double O2_an_fcn()
{
   return PaO2n*henry_O2;
}

///////////

double Hy_vn_fcn()
{
   return 1000.0*pow(10.0, -pH_vn);
}

///////////

double O2_vn_fcn()
{
   return PvO2n*henry_O2;
}

///////////

double CO2_vn_fcn()
{
   // appears to give a very suspect value at normal production rates - however me must realise that a lot of the CO2 produced is carried out as bicarbonate
   return PvCO2n*henry_CO2;
}

///////////

double CO2_a_fcn()
{
   return Pa_CO2*henry_CO2;
}

///////////

double _eO2n_fcn()
{
   return 0.8*_O2n + 0.2*O2n;
}

///////////

double _O2n_fcn()
{
   return Mass_Br*(_O2nw/Vol_BrW);
}

///////////

double _mO2n_fcn()
{
   return 0.5*_O2n;
}

///////////

double X_vn_fcn()
{
   return X_totn - X_O4_vn;
}

///////////

double X_O4_vn_fcn()
{
   return -1.0/4.0/qn*k_OOn*(O2n - _O2n) + 1/4.0*((O2_an - _O2n) + 4.0*X_O4_an);
}

///////////

double pH_bn_fcn()
{
   return -log10(Hyn/1000.0);
}

///////////

double PMFn_fcn()
{
   return PD_mitn + C_elec*log(_Hyn/_mHyn);
}

///////////

double k1_lk_fcn()
{
// check this out
   return Vol_BrW*V_leakn/(exp(k2_lk*PMFn) - 1.0);
}

///////////

double k_PCrATP_eq_fcn()
{
   // asuming the reaction is normally at equilibrium... This is an acid dependent reaction
   return _Crn*_ATPn/(_PCrn*_ADPn*_Hyn);
}

///////////

double k_PCrATP_fcn()
{
   double Ptemp;
   double Qtemp;
   Ptemp = _PCrn*_ADPn;
   Qtemp = _Crn + _ATPn;
   return log(2.0)/((Ptemp + Qtemp/k_PCrATP_eq)*k_PCrATP_ht);

}

///////////

double k_nPCrATP_fcn()
{
   double Ptemp;
   double Qtemp;
   Ptemp = _PCrn*_ADPn;
   Qtemp = _Crn + _ATPn;
   return log(2.0)/((k_PCrATP_eq*Ptemp + Qtemp)*k_PCrATP_ht);
}

///////////

double _PCrn_fcn()
{
   return _PCrnw*Mass_Br/Vol_BrW;
}

///////////

double _Crn_fcn()
{
   return _Crnw*Mass_Br/Vol_BrW;
}

///////////

double PD_Ca1n_fcn()
{
   return C_elec/2.0*log(Ca_en/Ca_i1n);
}

///////////

double PD_Ca2n_fcn()
{
   return C_elec/2.0*log(Ca_en/Ca_i2n);
}

///////////

double PD_Kn_fcn()
{
   // constant field theory
   return C_elec*log(_eKn/_K_musn);
}

///////////

double _K_musn_fcn()
{
// rather arbitrary
   return _Kn;
}

///////////

double K_Phosshut_fcn()
{
   return 5.0*_Phosn*_Hyn*(1.0 + p_Phosshut)/(1.0 + 2*p_Phosshut);
}

///////////

double Vmax_Phosshut_fcn()
{
// assuming that total phosphate flowing in is equal to total ATP production (apart from that produced by glycolysis)
   return (1.0 + (1.0 - Hleak_frac)*(2.5*N1_H + 0.5*N2_H - 1.0)/4.0)*(2.0*CMRglucn + CMRlacn)*Vol_BrW/(R_Phosshuttle);
}

///////////

double K_Pyshut_fcn()
{
   return Km_Pyshut*_Hyn*(1.0 + p_Pyshut)/(1.0 + 2*p_Pyshut);
}

///////////

double Vmax_Pyshut_fcn()
{
   return netPyusen*Vol_BrW/(R_pyrshuttle);
}

///////////

double Km_PyLH_fcn()
{
   return Km_PyLH_frac*_Hyn;
}

///////////

double VR_PyL_fcn()
{
   return Rbf_PyL*(1.0 + Km_PyLL/_L0n + Km_PyLN/_NADn + Km_PyLLN/_L0n/_NADn)/(1.0 + Km_PyLPy/_Py0n + Km_PyLNH/_NADHn + Km_PyLPN/_Py0n/_NADHn + Km_PyLH/_Hyn);
}

///////////

double Vmax_PyL_fcn()
{
   double tempval;
   tempval = (R_PytoL);
   if(tempval < 0){
     fprintf(stderr, "ERROR in Vmax_PyL_fcn - normal production of lactate = %.4f - less than zero.\n", tempval);
     return 0.0;
}
   return -CMRlacn*Vol_BrW/Vol_exm/(tempval);
}

///////////

double P_icn_fcn()
{
   // from the normal outflow of CSF
   return q_csfn/G_o + P_vsn;
}

//////////

double P_cn_fcn()
{
   // from the normal production of CSF
   return q_csfn/G_f + P_icn;
}

//////////

double P_vn_fcn()
{
   // from conservation of flow
   return (qn - q_csfn)/G1_vs + P_icn;
}

//////////

double P_vsn_fcn()
{
   return qn/G_ve + P_cv;
}

//////////

double r_1n_fcn()
{
   return pow(G_1n/K_G1, 0.25);
}

//////////

double r_2n_fcn()
{
   return pow(G_2n/K_G2, 0.25);
}

//////////

double G_1n_fcn()
{
   return G_2n*G_ratn;
}

//////////

double G_2n_fcn()
{
   // from the ratio, and the fact that q = G_1G_2/(G_1 + G_2)*(P_a - P_c)
   return qn*(1.0 + G_ratn)/G_ratn/(P_an - P_cn);
}

//////////

double V_1n_fcn()
{
   return K_V1*r_1n*r_1n;

}

//////////

double V_2n_fcn()
{
   return K_V2*r_2n*r_2n;
}

//////////

double P_1n_fcn()
{
   // from q = 2*G_1*(P_a - P_1)
   return P_an - qn/2.0/G_1n;
}

///////////

double P_2n_fcn()
{
   // from q = 2*G_2*(P_2 - P_c)
   return P_cn + qn/2.0/G_2n;
}

///////////

double h_1n_fcn()
{
   return -r_1n + sqrt(r_1n*r_1n + 2.0*r_01*h_01 + h_01*h_01);
}

///////////

double h_2n_fcn()
{
   return -r_2n + sqrt(r_2n*r_2n + 2.0*r_02*h_02 + h_02*h_02);
}

///////////

double sigma_e1n_fcn()
{
   return sigma_e01*(exp(K_sigma1*(r_1n-r_01)/r_01) - 1.0) - sigma_coll1;
}

///////////

double sigma_e2n_fcn()
{
   return sigma_e02*(exp(K_sigma2*(r_2n-r_02)/r_02) - 1.0) - sigma_coll2;
}

///////////

double T_e1n_fcn()
{
   return sigma_e1n*h_1n;
}

///////////

double T_e2n_fcn()
{
   return sigma_e2n*h_2n;
}

///////////

double T_m1n_fcn()
{
   return T_max1n*exp(-pow(fabs((r_1n-r_m1)/(r_t1-r_m1)), n_m1));
}

///////////

double T_m2n_fcn()
{
   return T_max2n*exp(-pow(fabs((r_2n-r_m2)/(r_t2-r_m2)), n_m2));
}

///////////

double V_mcan_fcn()
{
   return qn/3.0*(1.0/(Pi*r_mcan*r_mcan));
}

//////////

double SvO2n_fcn()
{
   return SaO2n - CMRO2n*Vol_BrW/qn/4/X_totn;
}

///////////

double Vol_BrW_fcn()
{
   return Vol_BrW_frac*Vol_Br;
}

//////////

double Vol_ecs_fcn()
{
   return Vol_ecs_frac*Vol_BrW;
}

///////////

double Vol_inm_fcn()
{
   return Vol_inm_frac*Vol_BrW;
}

//////////

double Vol_exm_fcn()
{
   return Vol_exm_mod*(1.0-Vol_inm_frac-Vol_ecs_frac)*Vol_BrW;
}

//////////

double Vol_blood_fcn()
{
   return Vol_blood_frac*Vol_BrW;
}

///////////

double Vol_c_fcn()
{
   return Vol_c_frac*Vol_blood;
}

///////////

double Vol_artn_fcn()
{
   // must be compatible with Ursino model
   return V_1n + V_2n;
}

///////////

double Vol_venn_fcn()
{
   return Vol_blood - Vol_c - Vol_artn;
}

///////////

double Vol_ven00_fcn()
{
   double temp;
   temp = Vol_venn - (C_vin*(P_vn-P_icn) + C_ve*P_vsn);
   if(temp < 0.0){
     fprintf(stderr, "ERROR in Vol_ven00_fcn: Vol_ven00 is less than 0.\n");
     return 0.0;
   }
   return temp;
}

///////////