function dc = reactormodel_co2_5(t,C,rateFlag,growthKinetics,metabolismKinetics,Cinit)


Ceglc = C(1);
Ccglc = C(2);
Ccg6p = C(3);
Ccf6p = C(4);
Ccfbp = C(5);
Ccdhap = C(6);
Ccgap = C(7);
Cc13p2g = C(8);
Cc3pg = C(9);
Cc2pg = C(10);
Ccpep = C(11);
Ccpyr = C(12);
Cmpyr = C(13);
Cclac = C(14);
Celac = C(15);
Cc6pg = C(16);
Cnadp = C(17);
Ccgsh = C(18);
Ccru5p = C(19);
Ccxyl5p = C(20);
Ccr5p = C(21);
Ccsh7p = C(22);
Cce4p = C(23);
Cmaccoa = C(24);
Cmicit = C(25);
Cmakg = C(26);
Cnad = C(27);
Cmscoa = C(28);
Cmsuc = C(29);
Cmfum = C(30);
Cmmal = C(31);
Cmoaa = C(32);
Cmasp = C(33);
Cmglu = C(34);
Ccoaa = C(35);
Ccmal = C(36);
Ccasp = C(37);
Ccakg = C(38);
Ccglu = C(39);
Cmcit = C(40);
Cccit = C(41);
Ccf26p = C(42);
Ccell = C(43);
CeI = C(44);
Ceosm = C(45);
Cgco2 = C(46);
Ceco2 = C(47);
Ceh2co3 = C(48);
CeH = C(49);
Cehco3 = C(50);
Ceco3 = C(51);
pHI = C(52);


% ------------------------------------------------------------------------
%             Cell growth and extracellular metabolites
% ------------------------------------------------------------------------

% Adjustable growth parameters
mumax = growthKinetics(1);
kdmax = growthKinetics(2);
Kelac = growthKinetics(3);
CIexp = growthKinetics(4);
qI = growthKinetics(5);
Kos = growthKinetics(6);
Cos = growthKinetics(7);
cellosmo = growthKinetics(8);
kla = growthKinetics(9);

% Cell growth rate
mu = (mumax) * Kelac/(Kelac + Celac^2) * 1/(1 + CeI^CIexp) * Ceosm/(Ceosm + Kos * (Ceosm-Cos)^2);

% Death rate
kd = kdmax * (mumax - mu)/mumax;

% Cell concentration balance
rcell = (mu-kd)*Ccell;

% ------------------------------------------------------------------------
%             Interaction of environment with signaling (KBP)
% ------------------------------------------------------------------------

% Adjustable metabolic parameters (for signaling)
qcell = metabolismKinetics(1);
Kampk = metabolismKinetics(2);
Kstress = metabolismKinetics(3);

% Growth effect on Akt
Akt = 100*(mu/(mu+mumax));

% Stress effect of AMPK
inhibition = (Ceosm + Kstress * (Ceosm-Cos)^2)/Ceosm * (1 + CeI^CIexp);
AMPK = inhibition/(Kampk + inhibition);

% ------------------------------------------------------------------------
%             Enzyme variables
% ------------------------------------------------------------------------

% Adjustable metabolic parameters
E0perm = metabolismKinetics(4);
E0hk = metabolismKinetics(5);
E0pgi = metabolismKinetics(6);
E0pfkfb = metabolismKinetics(7);
E0pfk = metabolismKinetics(8);
E0ald = metabolismKinetics(9);
E0tpi = metabolismKinetics(10);
E0gapd = metabolismKinetics(11);
E0pgk = metabolismKinetics(12);
E0pgm = metabolismKinetics(13);
E0en = metabolismKinetics(14);
E0pk = metabolismKinetics(15);
E0ldh = metabolismKinetics(16);
E0mct = metabolismKinetics(17);
E0g6pd = metabolismKinetics(18);
E0pdhc = metabolismKinetics(19);
Ncd = metabolismKinetics(20);
KBP = metabolismKinetics(21);


% ------------------------------------------------------------------------
%             Model Parameters
% ------------------------------------------------------------------------

Kx11=1;
Kx12=0.15;
Kx14=0.17;
Kx15=0.197;
Kx16=1.1989;
Kx17=1.0106;
Kx18=1.1308;
Kx19=2.4439;
Kx110=0.1;
Kx111=0.10211;
Kx115=3;
Kx116=1.6;

fct2=9.9989*2.25;
fct3=0.048;
fct5=600.5;

Ccpi=2.5;
Cmpi=2.5;

Mperm=1;
Mpgm=1;
Mg6pd=1;
Mep=1;
Mrpi=1;
M6pgd=1;
Mgssgr=1;
Mpfk=1;
   
% Parameters
Vc=(0.75)*1.5*10^-12; %L
Vm=(0.25)*1.5*10^-12; %L
CcH = 1*10^(-7.3); %molar
Mg = 0.7; %units?
pi = 2.5;
Co2 = 1.2;

Cnadh = Ncd - Cnad;
Ndp = 0.065132;
Cnadph = Ndp - Cnadp;
Gsn = 3.15;
Gssg = Gsn - Ccgsh;
Ccamp = 0.03;
Ccatp= 0.311;
Ccadn = 0.847;
Ccadp = Ccadn - Ccatp;
KeqMgAtp = 0.081;
KeqMgAdp = 0.81;
MgAtp = Mg*Ccatp/KeqMgAtp;
MgAdp = Mg*Ccadp/KeqMgAdp;
Ccg16p=0.1;
Kisglum = 1;
Cccoash=0.02;
Ccaccoa=0.001;
Cmgtp=0.1;
Cmgdp=0.1;
Cc23p2g=3.1;
Cmatp=0.1;
Cmadp=0.1;
Cmamp=0.1;
Cmco2=21.4;
Cmcoash=0.04;
Ccala=1;

KAkt=50;


% Rate Constants and Vmax for all the enzymes considered in the modeling
% Glycolysis

% Glucose transporter
rmaxperm = 7.27*(0.2+0.8/(1+KAkt/Akt));
rmaxperm2 = 7.27*(0.2+0.8/(1+KAkt/Akt));
Kglc = 1.5;

%Hexokinase
Khkmgatp=1.0;
Khkimgatp=1.0;
Khkiglc=0.1;
Khkig6p=0.47;
Khkmgadp=1.0;
Khkimgadp=1.0;
Khki23p2g1=4.0;
Khkig16p1=0.03;
Khkig6p1=0.02;
KhkiGSH1=3.0;
Ehk=0.047;
khkf=3600*180*Ehk*(0.8+0.2/(1+0/Akt));
khkr=3600*1.16*Ehk*(0.8+0.2/(1+0/Akt));


% Glucose Phosphate Isomerase
Egpi = 0.95;
k2gpi = 1760*3600;
k3gpi = 1470*3600;
Kfmgpi = 7*10^(-1);
rfmgpi = Egpi*k2gpi;
rrmgpi = Egpi*k3gpi;
Kigpifbp=100000000000000000000007;
Krmgpi=123.89*10^(-3)*(1+(Ccfbp/Kigpifbp));

% Phosphofructokinase 2 (modeled as 2 enzymes PFK2 and F26BPase)
% Phosphofructokinase 2 parameters
Kmpfk2atp= 0.15;
Kmpfk2f6p= 0.032;
Kmpfk2f26p= 0.008*(0.2+0.8/(1+KAkt/Akt));
Kmpfk2adp= 0.062;
Kipfk2atp= 0.150;
Kipfk2f6p= 0.001;
Kipfk2f26p= 0.020;
Kipfk2adp= 0.230;
Keqpfk2= 16;
Kipfk2pep=0.013;
Vfpfk2= 300*(0.2+0.8/(1+KAkt/Akt));

% F26BPase parameters
Vff26bpase=13.86;
Kmf26bpasef26p=0.001;
Kif26bpasef6p=0.025;

% Phosphofructokinase
Epfk=0.0002;
kpfkf=822*3600*Epfk*(0.2+0.8/(1+0/Akt));
kpfkr=36*3600*Epfk*(0.2+0.8/(1+0/Akt));
Kpfkf6p =0.06;
Kpfkmgatp = 0.068;
Kpfkatp = 0.1;
Kpfkamp = 0.3;
Kpfkmg = 0.2;
Lpfk = 0.002;
Kpfkf26p=2*0.0027544;
Kpfkpi=30;
Kpfkg16p=0.1;
Kpfk23p2g=0.5;
Kpfkfbp=0.3;
Kpfkmgadp=0.54;
Kpfkilac=30;

% Aldolase
Eald=0.0053;
kaldf=244800*Eald;
kaldr=842400*Eald;
Kaldfbp=0.05;
Kaldifbp=.0198;
Kalddhap=.035;
Kaldgap=.189;
Kaldidhap=.011;
Kaldi23p2g=1.5;

% Triose-phosphate Isomerase
Etpi = 0.0097;
k2tpi = 5.26*10^7;
k3tpi = 4.75*10^6;
Kmftpi = 1.62*10^(-1);
Kmrtpi = 4.3*10^(-1);
rmftpi = k2tpi*Etpi;
rmrtpi = k3tpi*Etpi;

% Glyceraldehyde Phosphate Dehydrogenase
Kgapdnad=0.310;0.045;
Kgapdinad=0.045;
Kgapdipi=2.5;3.16;
Kgapdgap=0.095;
Kgapdigap1=0.031;
Kgapdigap=1.59*10^(-19)*10^3;
Kgapdnadh=0.0033;
Kgapdinadh=0.010;
Kgapd13p2g=0.000671;0.00671;
Kgapdi13p2g=1.52*10^(-21)*1000;
Kgapdi13p2g1=0.001;
Egap=.0045;
kgapdf=232*3600*Egap;
kgapdr=171*3600*Egap;

% Phosphoglycerate Kinase
Kpgkimgadp=.08;
Kpgk13p2g=.002;
Kpgki13p2g=1.6;
Kpgkimgatp=.186;
Kpgk3pg=1.1;
Kpgki3pg=.205;
Epgk=0.025;
kpgkf=2290*3600*Epgk;
kpgkr=917*3600*Epgk;

% Bisphosphoglycerate Mutase
k13dpg = 2.75*10^5;
k23dpg = 4*10^-2;

% Bisphosphoglycerate Phosphatase
rmdpgase = 0.52;
Kmdpgase = 0.2;

% Phosphoglycerate Mutase
Epgm=0.057;
kpgmf=795*3600*Epgm;
kpgmr=714*3600*Epgm;
Kpgm3pg=.168;
Kpgm2pg=.0256;

% Enolase
Een=0.0009;
kenf=190*3600*Een;
kenr=50*3600*Een;
Ken2pg=.14;
Keni2pg=.14;
Kenpep=.11;
Kenipep=.11;
Kenimg=.046;

% Pyruvate Kinase
Epk=0.027;
rmpkf = 4989600*Epk;
rmpkr = 11736*Epk;
Kpkpep = 0.225;
Kpkmgadp = 0.474;
Kpkpyr=4;
Kpkmgatp=3.0;
Kpkatp = 3.39;
Kpkg16p=0.1;
Kpkfdp =0.3;0.1;
Lpk = 0.398;
Kpkala=0.02;

% Lactate Dehydrogenase
Kldhnadh=0.00844;
Kldhinadh=0.00245;
Kldhpyr=0.2;0.137;
Kldhipyr=0.228;
Kldhi2pyr=0.101;
Kldhnad=0.107;
Kldhinad=0.503;
Kldhlac=10.0701;
Kldhilac=30;7.33;
Eldh=0.00343;
kldhf=1648800*Eldh;
kldhr=115*3600*Eldh;

% Lactate export
hyd = 7.0;
CeHexpt = 10^(-hyd); %%% Note this was removed so that pH can affect transport
KmcHmct= 0.0001;
Kmclacmct= 2.5;
KicHmct= 0.0002;
KieHmct= 0.0002;

% Pentose Phosphate Pathway

% Glutathione Oxidation
kgshox = 0.260;

% Glutathione Reductase
Egssgr = 1.25*10^-7;
k1gssgr = 3.6*8.5*10^7;
k2gssgr = 3600*5.1*10^2;
k3gssgr = 3.6*1*10^8;
k4gssgr = 3600*7.2*10^3;
k5gssgr = 3600*8.1*10^2;
k6gssgr = 3600*1*10^3;
k7gssgr = 3600*1*10^6;
k8gssgr = 3.6*5*10^7;
k9gssgr = 3600*1*10^6;
k10gssgr = 3.6*5*10^7;
k11gssgr = 3600*7*10^3;
k12gssgr = 3.6*1*10^8;
N1gssgr = k1gssgr*k3gssgr*k5gssgr*k7gssgr*k9gssgr*k11gssgr;
N2gssgr = k2gssgr*k4gssgr*k6gssgr*k8gssgr*k10gssgr*k12gssgr;
D1gssgr = k2gssgr*k9gssgr*k11gssgr*(k4gssgr*k6gssgr + k4gssgr*k7gssgr + k5gssgr*k7gssgr);
D2gssgr = k1gssgr*k9gssgr*k11gssgr*(k4gssgr*k6gssgr + k4gssgr*k7gssgr + k5gssgr*k7gssgr);
D3gssgr = k3gssgr*k5gssgr*k7gssgr*k9gssgr*k11gssgr;
D4gssgr = k2gssgr*k4gssgr*k6gssgr*k8gssgr*k11gssgr;
D5gssgr = k2gssgr*k9gssgr*k12gssgr*(k4gssgr*k6gssgr + k4gssgr*k7gssgr + k5gssgr*k7gssgr);
D6gssgr = k1gssgr*k3gssgr*(k5gssgr*k9gssgr*k11gssgr + k6gssgr*k9gssgr*k11gssgr + k7gssgr*k9gssgr*k11gssgr + k5gssgr*k7gssgr*k9gssgr + k5gssgr*k7gssgr*k11gssgr);
D7gssgr = k1gssgr*k4gssgr*k6gssgr*k8gssgr*k11gssgr;
D8gssgr = k3gssgr*k5gssgr*k7gssgr*k9gssgr*k12gssgr;
D9gssgr = k2gssgr*k4gssgr*k6gssgr*k8gssgr*k10gssgr;
D10gssgr = k2gssgr*k4gssgr*k6gssgr*k8gssgr*k12gssgr;
D11gssgr = k2gssgr*k10gssgr*k12gssgr*(k4gssgr*k6gssgr + k4gssgr*k7gssgr + k5gssgr*k7gssgr);
D12gssgr = k1gssgr*k3gssgr*k8gssgr*k11gssgr*(k5gssgr + k6gssgr);
D13gssgr = k1gssgr*k3gssgr*k5gssgr*k7gssgr*k10gssgr;
D14gssgr = k1gssgr*k4gssgr*k6gssgr*k8gssgr*k10gssgr;
D15gssgr = k3gssgr*k5gssgr*k7gssgr*k10gssgr*k12gssgr;
D16gssgr = k8gssgr*k10gssgr*k12gssgr*(k2gssgr*k4gssgr + k2gssgr*k5gssgr + k2gssgr*k6gssgr + k6gssgr*k4gssgr);
D17gssgr = k1gssgr*k3gssgr*k8gssgr*k10gssgr*(k5gssgr + k6gssgr);
D18gssgr = k3gssgr*k8gssgr*k10gssgr*k12gssgr*(k5gssgr + k6gssgr);
Detgssgr = D1gssgr + D2gssgr*Cnadph + D3gssgr*Gssg + D4gssgr*Ccgsh + D5gssgr*Cnadp + D6gssgr*Cnadph*Gssg + D7gssgr*Cnadph*Ccgsh +...
    D8gssgr*Gssg*Cnadp + D9gssgr*Ccgsh*Ccgsh + D10gssgr*Ccgsh*Cnadp + D11gssgr*Ccgsh*Cnadp + D12gssgr*Cnadph*Gssg*Ccgsh +...
    D13gssgr*Cnadph*Gssg*Ccgsh + D14gssgr*Cnadph*Ccgsh*Ccgsh + D15gssgr*Gssg*Ccgsh*Cnadp + D16gssgr*Ccgsh*Ccgsh*Cnadp + D17gssgr*Cnadph*Gssg*Ccgsh*Ccgsh + D18gssgr*Gssg*Ccgsh*Cnadp*Ccgsh;

% Glucose 6-Phosphate Dehydrogenase
Eg6pd = 9.3*10^-8;
k1g6pd = 3.6*1.1*10^8;
k2g6pd = 3600*8.7*10^2;
k3g6pd = 3.6*2.6*10^7;
k4g6pd = 3600*3*10^2;
k5g6pd = 3600*7.5*10^2;
k6g6pd = 3600*2*10^3;
k7g6pd = 3600*2.2*10^5;
k8g6pd = 3.6*1.1*10^9;
k9g6pd = 3600*1*10^4;
k10g6pd = 3.6*1.4*10^9;
N1g6pd = k1g6pd*k3g6pd*k5g6pd*k7g6pd*k9g6pd;
N2g6pd = k2g6pd*k4g6pd*k6g6pd*k8g6pd*k10g6pd;
D1g6pd = k2g6pd*k9g6pd*(k4g6pd*k6g6pd + k5g6pd*k6g6pd + k5g6pd*k7g6pd);
D2g6pd = k1g6pd*k9g6pd*(k4g6pd*k6g6pd + k5g6pd*k6g6pd + k5g6pd*k7g6pd);
D3g6pd = k3g6pd*k5g6pd*k7g6pd*k9g6pd;
D4g6pd = k2g6pd*k4g6pd*k6g6pd*k8g6pd;
D5g6pd = k2g6pd*k10g6pd*(k4g6pd*k6g6pd + k5g6pd*k6g6pd + k5g6pd*k7g6pd);
D6g6pd = k1g6pd*k3g6pd*(k5g6pd*k7g6pd + k5g6pd*k9g6pd + k6g6pd*k9g6pd + k7g6pd*k9g6pd);
D7g6pd = k1g6pd*k4g6pd*k6g6pd*k8g6pd;
D8g6pd = k3g6pd*k5g6pd*k7g6pd*k10g6pd;
D9g6pd = k8g6pd*k10g6pd*(k2g6pd*k4g6pd + k2g6pd*k5g6pd + k2g6pd*k6g6pd + k4g6pd*k6g6pd);
D10g6pd = k1g6pd*k3g6pd*k8g6pd*(k5g6pd + k6g6pd);
D11g6pd = k3g6pd*k8g6pd*k10g6pd*(k5g6pd + k6g6pd);
Detg6pd = D1g6pd + D2g6pd*Cnadp+ D3g6pd*Ccg6p+ D4g6pd*Cc6pg+ D5g6pd*Cnadph + D6g6pd*Cnadp*Ccg6p+ D7g6pd*Cnadp*Cc6pg + D8g6pd*Ccg6p*Cnadph +...
    D9g6pd*Cc6pg*Cnadph + D10g6pd*Cnadp*Ccg6p*Cc6pg+ D11g6pd*Ccg6p*Cc6pg*Cnadph;

% 6-Phosphogluconate Dehydrogenase
E6pgd = 2.1*10^-6;
k16pgd = 3.6*1.2*10^6;
k26pgd = 3600*4.1*10^2;
k36pgd = 3.6*1*10^9;
k46pgd = 3600*2.6*10^4;
k56pgd = 3600*4.8*10^1;
k66pgd = 3600*3*10^1;
k76pgd = 3600*6.3*10^2;
k86pgd = 3.6*3.6*10^4;
k96pgd = 3600*8*10^2;
k106pgd = 3.6*4.5*10^5;
k116pgd = 3600*2.8*3*10^2;
k126pgd = 3.6*9.9*10^6;
N16pgd = k16pgd*k36pgd*k56pgd*k76pgd*k96pgd*k116pgd;
N26pgd = k26pgd*k46pgd*k66pgd*k86pgd*k106pgd*k126pgd;
D16pgd = k26pgd*k96pgd*k116pgd*(k46pgd*k66pgd + k46pgd*k76pgd + k56pgd*k76pgd);
D26pgd = k16pgd*k96pgd*k116pgd*(k46pgd*k66pgd + k46pgd*k76pgd + k56pgd*k76pgd);
D36pgd = k36pgd*k56pgd*k76pgd*k96pgd*k116pgd;
D46pgd = k26pgd*k46pgd*k66pgd*k86pgd*k116pgd;
D56pgd = k26pgd*k96pgd*k126pgd*(k46pgd*k66pgd + k46pgd*k76pgd + k56pgd*k76pgd);
D66pgd = k16pgd*k36pgd*(k56pgd*k96pgd*k116pgd + k66pgd*k96pgd*k116pgd + k76pgd*k96pgd*k116pgd + k56pgd*k76pgd*k96pgd + k56pgd*k76pgd*k116pgd);
D76pgd = k16pgd*k46pgd*k66pgd*k86pgd*k116pgd;
D86pgd = k36pgd*k56pgd*k76pgd*k96pgd*k126pgd;
D96pgd = k26pgd*k46pgd*k66pgd*k86pgd*k106pgd;
D106pgd = k26pgd*k46pgd*k66pgd*k86pgd*k126pgd;
D116pgd = k26pgd*k106pgd*k126pgd*(k46pgd*k66pgd + k46pgd*k76pgd + k56pgd*k76pgd);
D126pgd = k16pgd*k36pgd*k86pgd*k116pgd*(k56pgd + k66pgd);
D136pgd = k16pgd*k36pgd*k56pgd*k76pgd*k106pgd;
D146pgd = k16pgd*k46pgd*k66pgd*k86pgd*k106pgd;
D156pgd = k36pgd*k56pgd*k76pgd*k106pgd*k126pgd;
D166pgd = k86pgd*k106pgd*k126pgd*(k26pgd*k46pgd + k26pgd*k56pgd + k26pgd*k66pgd + k66pgd*k46pgd);
D176pgd = k16pgd*k36pgd*k86pgd*k106pgd*(k56pgd + k66pgd);
D186pgd = k36pgd*k86pgd*k106pgd*k126pgd*(k56pgd + k66pgd);
Det6pgd = D16pgd + D26pgd*Cnadp+ D36pgd*Cc6pg+ D46pgd*Co2 + D56pgd*Cnadph + D66pgd*Cnadp*Cc6pg + D76pgd*Cnadp*Co2 + D86pgd*Cc6pg*Cnadph +...
D96pgd*Co2*Ccru5p + D106pgd*Co2*Cnadph + D116pgd*Ccru5p*Cnadph + D126pgd*Cnadp*Cc6pg*Co2 + D136pgd*Cnadp*Cc6pg*Ccru5p + D146pgd*Cnadp*Co2*Ccru5p +...
D156pgd*Cc6pg*Ccru5p*Cnadph + D166pgd*Co2*Ccru5p*Cnadph + D176pgd*Cnadp*Cc6pg*Co2*Ccru5p+ D186pgd*Cc6pg*Co2*Ccru5p*Cnadph;

% Ribose Phosphate Epimerase
Eep = 4.22*10^-3;
k2ep = 1.58*10^6;
k3ep = 1.1*10^6;
Kfmep = 1.9*10^-1;
Krmep = 5*10^-1;
rmfep = Eep*k3ep;
rmrep = Eep*k2ep;
          
% Ribose Phosphate Isomerase
Eki = 1.42*10^-2;
k2ki = 5.11*10^4;
k3ki = 1.2*10^5;
Kfmki = 7.8*10^-1;
Krmki = 2.2;
rmfki = Eki*k3ki;
rmrki = Eki*k2ki;

% Phosphorybosylpyrophosphate Synthase
rmprpps = 1.1;
Kprppsatp = 0.01;
Kprppsr5p = 0.57;

% Transketalose 1
Etk1 = 3.3*10^-4;
k1tk1 = 7.78*10^5;
k2tk1 = 1.37*10^5;
k3tk1 = 1.22*10^5;
k4tk1 = 5.62*10^5;
k5tk1 = 1.18*10^6;
k6tk1 = 6.30*10^5;
k7tk1 = 0.5*1.44*10^5;
k8tk1 = 1.61*10^5;

% Transaldolose
Eta = 6.9*10^-4;
k1ta = 2.09*10^6;
k2ta = 1.63*10^5;
k3ta = 5.87*10^4;
k4ta = 3.64*10^6;
k5ta = 1.76*10^6;
k6ta = 2.16*10^5;
k7ta = 6.12*10^4;
k8ta = 0.2*2.84*10^5;
N1ta = k1ta*k3ta*k5ta*k7ta;
N2ta = k2ta*k4ta*k6ta*k8ta;
D1ta = k1ta*k3ta*(k6ta + k7ta);
D2ta = k5ta*k7ta*(k2ta + k3ta);
D3ta = k2ta*k4ta*(k6ta + k7ta);
D4ta = k6ta*k8ta*(k2ta + k3ta);
D5ta = k1ta*k5ta*(k3ta + k7ta);
D6ta = k4ta*k8ta*(k2ta + k6ta);
D7ta = k5ta*k8ta*(k2ta + k3ta);
D8ta = k1ta*k4ta*(k6ta + k7ta);

% Transketalose 2
k5tk2 = 8.06*10^6;
k6tk2 = 6.3*10^5;
k7tk2 = 1.44*10^5;
k8tk2 = 7.67*10^4;


%TCA Cycle concentrations and parameters
CmH = 1*10^-8;
CmK = 1.2*10^-1;
CmMg = 1*10^-3;
Cmnh3= 0.45;
KHatp = 2.57*10^-7;
KMgatp = 1.51*10^-4;
KKatp = 1.35*10^-2;
Patp = 1 + CmH/KHatp + CmMg/KMgatp + CmK/KKatp;
KHadp = 3.80*10^-7;
KMgadp = 1.62*10^-3;
KKadp = 2.95*10^-2;
Padp = 1 + CmH/KHadp + CmMg/KMgadp + CmK/KKadp;
KHamp = 6.03*10^-7;
KMgamp = 1.38*10^-2;
KKamp = 8.91*10^-2;
Pamp = 1 + CmH/KHamp + CmMg/KMgamp + CmK/KKamp;
KHcoash = 7.41*10^-9;
Pcoash = 1 + CmH/KHcoash;
KMgoaa = 9.9*10^-1;
Poaa = 1 + CmMg/KMgoaa;
KHcit = 2.34*10^-6;
KMgcit = 4.27*10^-4;
KKcit = 4.58*10^-1;
Pcit = 1 + CmH/KHcit + CmMg/KMgcit + CmK/KKcit;
KHicit = 2.29*10*-6;
KMgicit = 3.47*10^-3;
Picit = 1 + CmH/KHicit + CmMg/KMgicit;
KHscoa = 1.10^10^-4;
Pscoa = 1 + CmH/KHscoa;
KHsuc = 7.41*10^-6;
KMgsuc = 6.76*10^-2;
KKsuc = 3.14*10^-1;
Psuc = 1 + CmH/KHsuc + CmMg/KMgsuc + CmK/KKsuc;
KHfum = 7.94*10^-5;
Pfum = 1 + CmH/KHfum;
KHmal = 1.78*10^-5;
KMgmal = 2.82*10^-2;
KKmal = 1.28;
Pmal = 1 + CmH/KHmal + CmMg/KMgmal + CmK/KKmal;
KHasp = 2.24*10-4;
KMgasp = 4.79*10^-3;
Pasp = 1 + CmH/KHasp + CmMg/KMgasp;
KHglu = 8.71*10^-5;
KMgglu = 1.51*10^-2;
Pglu = 1 + CmH/KHglu + CmMg/KMgglu;
KHco2 = 1.78*10^-10;
Pco2 = 1 + CmH*10/KHco2;
Paccoa = 1;
Pnadh = 1;
Pnad = 1;
Pakg = 1;
Pgtp = Patp;
Pgdp = Padp;
Pqh2 = 1;
Pcoq = 1;
KMgpyr = 0.382;
Ppyr = 1+CmMg/KMgpyr;
Pfatp= 1 + CmH*10/KHatp + CmK/KKatp;
Pfadp = 1 + CmH/KHadp + CmK/KKadp;
Pfamp = 1 + CmH/KHamp + CmK/KKamp;
Pfgtp = Pfatp;
Pfgdp = Pfadp;
Cmnadtot = 2.97;
Cmnadh=0.1;
Cmnad = Cmnadtot - Cmnadh;
Cmqtot = 1.35;
Cmqh2=0.2*Cmqtot;
Cmcoq = Cmqtot - Cmqh2;

% Pyruvate dehydrogenase
kpdhc = 1.22*10^-1*3.6*10^6;
Kpdhcnad = 60.7*0.001;
Kpdhccoa = 9.9*0.001;
Kpdhcpyr = 38.8*0.001;
Kpdhcinadh = 40.2*0.001;
Kpdhciaccoa = 40.0*0.001;
Kpdhceq0 = 5.02*10^-7;
Kpdhceq = Kpdhceq0*Pco2*Paccoa*Pnadh/(CmH*Ppyr*Pcoash*Pnad);
Kpdhcai1 = 1 + Cmaccoa/Kpdhciaccoa;
Kpdhcai2 = 1 + Cmnadh/Kpdhcinadh;

%citrate synthase
kcs = 1*11.6*3.6*10^6;
Kcseq0 = 7.34*10^-5;      
Kcseq = Kcseq0*Pcoash*Pcit/((CmH^2)*Poaa*Paccoa);
Kcsaccoa = 1.4*0.001;
Kcsoaa = 4*0.001;
Kcsia = 3.33*0.001;
Kcsicit = 1600*0.001;
Kcsiatp = 900*0.001;
Kcsiadp = 1800*0.001;
Kcsiamp = 6000*0.001;
Kcsicoa = 67*0.001;
Kcsiscoa = 140*0.001;
Kcsai1 = 1 + (Cmcit/Pcit)/Kcsicit;
Kcsai2 = 1 + (Cmatp/Patp)/Kcsiatp + (Cmadp/Padp)/Kcsiadp + (Cmamp/Pamp)/Kcsiamp + Cmcoash/Kcsicoa + Cmscoa/Kcsiscoa;

% Aconitase
kacon = 3.21*10^-2*3.6*10^6;
Kaconeq0 = 7.59*10^-2;
Kaconeq = Kaconeq0*Picit/Pcit;
Kaconicit = 434*0.001;
Kaconcit = 1161*0.001;

% Isocitrate dehydrogenase
kidh = 0.425*3.6*10^6;
Kidheq0 = 3.50*10^-16;
Kidheq = Kidheq0*Pakg*Pnadh*Pco2/((CmH^2)*50*Pnad*Picit);
Kidhicit = 183*0.001;
Kidhnad = 74*0.001;
Kidhib = 23.8*0.001;
Kidhiq = 29*0.001;
nh = 3;
Kidhaadp = 50*0.001;
Kidhiatp = 91*0.001;
Kidhai = 1 + (Kidhaadp*(1 + ((Cmatp*Pfatp/Patp)/Kidhiatp))/(Cmadp*Pfadp/Padp));

% Alpha-ketogluterate dehydrogenase
Kakgdeq0 = 6.93*10^-3;
Kakgdeq = Kakgdeq0*Pco2*Pscoa*Pnadh/(CmH*Pakg*Pcoash*Pnad);
kakgd = 5*7.70*10^-6*3.6*10^6;
Kakgdakg = 120*0.001;
Kakgdcoa = 55*0.001;
Kakgdnad= 21*0.001;
Kakgdiq = 6.9*0.001;
Kakgdir = 0.60*0.001;
Kakgdiatp = 50*0.001;
Kakgdaadp = 100*0.001;
Kakgdai = 1 + (Kakgdaadp/(Cmadp*Pfadp/Padp))*(1 + (Cmatp*Pfatp/Patp)/Kakgdiatp);

% Succinyl-CoA synthetase
kscoas = 0.2*0.582*10^(-3)*3.6*10^6;
Kscoasgdp = 16*0.001;
Kscoasscoa = 55*0.001;
Kscoaspi = 660.0*0.001;
Kscoascoa = 20*0.001;
Kscoassuc = 880*0.001;
Kscoasgtp = 11.1*0.001;
Kscoasia = 5.5*0.001;
Kscoasib = 100*0.001;
Kscoasic = 2000*0.001;
Kscoasip = 20*0.001;
Kscoasiq = 3000*0.001;
Kscoasir = 11.1*0.001;

% Succinate dehydrogenase
ksdh = 0.3*6.23*10^-3*3.6*10^6;
Ksdheq0 = 1.69;
Ksdheq = Ksdheq0*Pqh2*Pfum/(Psuc*Pcoq);
Ksdhiq = 1275*0.001;
Ksdhsuc = 467*0.001;
Ksdhia = 120*0.001;
Ksdhcoq = 480*0.001;
Ksdhqh2 = 2.45*0.001;
Ksdhfum = 1200*0.001;
Ksdhioaa = 1.5*0.001;
Ksdhasuc = 450*0.001;
Ksdhafum = 375*0.001;
Ksdhai = (1 + Cmoaa/Ksdhioaa + Cmsuc/Ksdhasuc + Cmfum/Ksdhafum)/(1 + Cmsuc/Ksdhasuc + Cmfum/Ksdhafum);

% Fumarase
kfum = 7.12*10^-3*3.6*10^6;
Kfumfum = 44.7*0.001;
Kfummal = 197.7*0.001;
Kfumicit = 3500*0.001;
Kfumiatp = 40*0.001;
Kfumiadp = 400*0.001;
Kfumigtp = 80*0.001;
Kfumigdp = 330*0.001;
Kfumeq0 = 4.04;
Kfumeq = Kfumeq0*Pmal/Pfum;
Kfumai = 1 + Cmcit/Kfumicit + (Cmatp*Pfatp/Patp)/Kfumiatp + (Cmadp*Pfadp/Padp)/Kfumiadp + (Cmgtp*Pfgtp/Pgtp)/Kfumigtp + (Cmgdp*Pfgdp/Pgdp)/Kfumigdp;

% Malate dehydrogenase 2
Kmdh2nad = 90.55*0.001;
Kmdh2mal = 250*0.001;
Kmdh2oaa = 6.128*0.001;
Kmdh2nadh = 2.580*0.001;
Kmdh2ia = 279*0.001;
Kmdh2ib = 360*0.001;
Kmdh2eq0 = 2*2.27*10^-12; %-12
Kmdh2eq = Kmdh2eq0*Poaa*Pnadh/(CmH*Pnad*Pmal);
Kmdh2ip = 5.5*0.001;
Kmdh2iq = 3.18*0.001;
Kmdh2iatp = 183.2*0.001;
Kmdh2iamp = 420.0*0.001;
Kmdh2iadp = 394.4*0.001;
kmdh2 = 6.94*10^-2*3.6*10^6;
Kmdh2ai= 1 + ((Cmadp*Pfadp/Padp)/Kmdh2iadp) + ((Cmatp*Pfatp/Patp)/Kmdh2iatp) + ((Cmamp*Pfamp/Pamp)/Kmdh2iamp);

% Glutamate-oxaloacetate transaminase 2
Kgot2asp = 0.89;
Kgot2akg = 3.22;
Kgot2oaa = 88*0.001;
Kgot2glu = 32.5;
Kgot2ia = 3900*0.001;
Kgot2iq = 10700*0.001;
Kgot2iakg = 26500*0.001;
Kgot2eq0 = 1.77;
Kgot2eq = Kgot2eq0*Poaa*Pglu/(Pasp*Pakg);
kgot2 = 7.96*3.6*10^6;
Kgot2ai = 1 + (Cmakg/Kgot2iakg);

% Mitochondrial Malic Enzyme
Keqmmalic= 34.4;
kcatmalic= 0.42;
Kmmalmalic=1.7; 0.34;
Kmnadmalic= 0.16;0.009;
Kiatpmalic= 0.5;

% Cytosolic Malic Enzyme
Kmnadcmalic= 0.0014;
Kmmalcmalic= 0.120;0.086;
KmCo2cmalic= 0.013;
Kmpyrcmalic= 0.0064;
Kmnadhcmalic= 0.0021;
Kinadcmalic= 0.00096;
Kimalcmalic= 0.22;
KiCo2cmalic= 0.0117;
Kipyrcmalic= 7.8; 
Kinadhcmalic= 0.002;
Keqcmalic= 34.4;
kcatcmalic= 0.336;

% Glutamate Dehydrogenase
Kmnadgdh=1.1;
Kmglugdh=3.5;
Kmnadhgdh=0.04;
Kmakggdh=1.1;
Kmnh3gdh=6;
Kiakggdh=0.25;
Kinh3gdh=6;
Kiglugdh= 3.5;
Kinadhgdh=0.004;
Kinadgdh=1;
Keqgdh= 0.003;
kcatgdh=53;

% ATP-Citrate Lyase
Kmcitcly= 0.0493;
Kicitcly= 0.0475;
Kmcoashcly= 0.0044;
Kicoashcly= 0.0061;
Kmoaacly= 0.177;
Kioaacly= 0.177;
Kmaccoacly= 0.0098;
Kiaccoacly= 0.0098;
Vcly= 49800;

% Pyruvate-hydrogen co-transporter flux
Kpyrh = 4.12*10^8*3.6;

% Glutamate-hydrogen co-transporter flux
Kgluh = 3.26*10^8*3.6;

% Citrate-malate exchange flux
Kcitmal = 73.1*3.6;

% a-ketoglutarate-malate exchange flux
Kakgmal = 0.346*3.6*10^6;
Kmmali = 0.4;
Kmmalm = 10;
Kmakgi =1.3; 0.3;
Kmakgm = 0.17;

% Malate-phosphate exchange flux
Kmalpi = 15.8*3.6;

% Aspartate-glutamate exchange flux
Keqaspglu= 0.6;
kaspglu = 7.48*10^-5*3.6*10^6;
Kiaspi = 0.028;
Kiaspm = 2.8;
Kiglui = 0.18;
Kiglum = 1.6;
Khaspglu = 10^-6.5;
m = 1.8;

% Glutaminase
Ccgln=60;
Kmglsgln= 12;
Kiglsglu=55;
Vfgls= 4.2;
Kglseq=1;

% Pyruvate carboxylase
Kmpyrpc=0.22;
Kmhco3pc=3.2;
Vfpc=0.000045*3800*60*70;

% Malate Dehydrogenase 1
Kmdh1nad = 0.114;
Kmdh1mal = 1.100;
Kmdh1oaa = 0.088;
Kmdh1nadh = 0.026;
Kmdh1ia = 0.0049;
Kmdh1ib = 0.063;
Kmdh1eq = 26644;
Kmdh1ip = 7.100;
Kmdh1iq = 0.9400;
kmdh1 = 1.6*10^4;

% Glutamate-oxaloacetate transaminase 1
Kgot1asp = 4400*0.001;
Kgot1akg = 380*0.001;
Kgot1oaa = 95*0.001;
Kgot1glu = 9600*0.001;
Kgot1ia = 3900*0.001;
Kgot1iq = 8400*0.001;
Kgot1iakg = 26500*0.001;
Kgot1eq0 = 1.77;
Kgot1eq = Kgot1eq0*Poaa*Pglu/(Pasp*Pakg);
kgot1 = 7.96*3.6*10^6;
Kgot1ai = 1 + (Ccakg/Kgot1iakg);

% Glutamate-pyruvate transaminase 1
kgpt1f=4*10^7;
kgpt1r=3.97*10^7;
Kgpt1ala=3;
Kgpt1akg=0.12;
Kgpt1pyr=0.23;
Kgpt1glu=8.1;
Kgpt1ipyr=0.23;
Kgpt1iglu=2.8;
Kgpt1IA=470;
Kgpt1IG=96;
Kgpt1RG=79.157;
Kgpt1eq=2.2;

% NAD Regeneration
fct4=0.8;
fct41=1.5;
fct6=1;

% Rate Equations
% Glycolysis

% Glucose transporters (GLUT1)
rperm = E0perm*0.5*100*2*Mperm*Kx11*(rmaxperm*Ceglc/Kglc - rmaxperm2*Ccglc/(Kglc/10))/(1+Ceglc/Kglc + Ccglc/(Kglc/10) );

% Hexokinase
rhk=E0hk*0.21*Kx12*((khkf*Ccglc*MgAtp/(Khkiglc*Khkmgatp))-(khkr*Ccg6p*MgAdp/(Khkig6p*Khkmgadp)))/(1+MgAtp/Khkimgatp +Ccglc/Khkiglc+Ccglc*MgAtp/(Khkiglc*Khkmgatp)+ MgAdp/Khkimgadp +Ccg6p/Khkig6p+Ccg6p*MgAdp/(Khkig6p*Khkmgadp)+...
    Ccglc*Ccg6p/(Khkiglc*Khkig6p1)+ Ccglc*Cc23p2g/(Khkiglc*Khki23p2g1)+ Ccglc*Ccgsh/(Khkiglc*KhkiGSH1)+ Ccglc*Ccg16p/(Khkiglc*Khkig16p1));

% Glucose Phosphate Isomerase
rpgi = E0pgi*0.005*0.8/10*(0.8*rfmgpi*(Ccg6p/Kfmgpi) - 1.1*rrmgpi*(Ccf6p/Krmgpi))/(1+(Ccg6p/Kfmgpi)+(Ccf6p/Krmgpi));


% 6-phosphofructo-2-kinase/fructose-2,6-biphosphatase 3
rpfk2a= KBP*AMPK*Vfpfk2*(Ccatp*Ccf6p-Ccf26p*Ccadp/Keqpfk2)/((Kipfk2atp*Kmpfk2f6p+Kmpfk2f6p*Ccatp+Kmpfk2atp*Ccf6p+Kmpfk2adp/Keqpfk2*Ccf26p+Kmpfk2f26p/Keqpfk2*Ccadp+...
       Ccatp*Ccf6p+Kmpfk2adp*Ccatp*Ccf26p/(Keqpfk2*Kipfk2atp)+Ccf26p*Ccadp/Keqpfk2+Kmpfk2atp*Ccf6p*Ccadp/Kipfk2adp+Ccatp*Ccf6p*Ccf26p/Kipfk2f26p+...
       Ccf6p*Ccf26p*Ccadp/(Kipfk2f6p*Keqpfk2))*(1+Ccpep/Kipfk2pep));

rf26bpase = Vff26bpase*Ccf26p/((1+Ccf6p/Kif26bpasef6p)*(Kmf26bpasef26p+Ccf26p));
rpfk2= E0pfkfb*(rpfk2a-0.85*Kx116*rf26bpase);

KBPcomb = KBP*AMPK*Vfpfk2/(0.85*Kx116*Vff26bpase);

% Phosphofructokinase
Npfk = 1 + Lpfk*(1+Ccatp/Kpfkatp)^4*(1+Mg/Kpfkmg)^4*(1+Cc23p2g/Kpfk23p2g)^4*(1+Cclac/Kpfkilac)^4/((1+Ccf6p/Kpfkf6p+Ccfbp/Kpfkfbp)^4*(1+Ccamp/Kpfkamp)^4*(1+Ccg16p/Kpfkg16p)^4*(1+pi/Kpfkpi)^4*(1+Ccf26p/Kpfkf26p)^4);
rpfk=E0pfk*2.38*1.1*Mpfk*Kx14*((kpfkf*Ccf6p*MgAtp/(Kpfkf6p*Kpfkmgatp))-(kpfkr*Ccfbp*MgAdp/(Kpfkfbp*Kpfkmgadp)))/(((1+Ccf6p/Kpfkf6p)*(1+MgAtp/Kpfkmgatp)+(1+Ccfbp/Kpfkfbp)*(1+MgAdp/Kpfkmgadp)-1)*Npfk);

% Aldolase
rald=E0ald*0.52*Kx15*(kaldf*Ccfbp/Kaldfbp-kaldr*Ccgap*Ccdhap/(Kaldgap*Kaldidhap))/(1+(Cc23p2g/Kaldi23p2g)+(Ccfbp/Kaldfbp)+(Kalddhap*Ccgap/(Kaldgap*Kaldidhap))*(1+(Cc23p2g/Kaldi23p2g))+(Ccdhap/Kaldidhap)+(Kalddhap*Ccfbp*Ccgap/(Kaldifbp*Kaldgap*Kaldidhap))+(Ccgap*Ccdhap/(Kaldgap*Kaldidhap)));

% Triose-phosphate Isomerase
rtpi = E0tpi*0.001*(rmftpi*Ccdhap/Kmftpi-50*Kx16*rmrtpi*Ccgap/Kmrtpi)/(1+Ccdhap/Kmftpi+Ccgap/Kmrtpi);

% Glyceraldehyde Phosphate Dehydrogenase
rgapd=E0gapd*1.4*Kx17*((kgapdf*pi*Ccgap*Cnad/(Kgapdnad*Kgapdipi*Kgapdigap))-(kgapdr*Cc13p2g*Cnadh*(CcH*1000)/(Kgapdi13p2g*Kgapdnadh)))/((1+Ccgap/Kgapdigap1)*(Ccgap/Kgapdigap+Cc13p2g/Kgapdi13p2g+Ccgap*pi/(Kgapdigap*Kgapdipi))+...
    Kgapd13p2g*Cnadh*(CcH*1000)/(Kgapdi13p2g*Kgapdnadh)+Kgapdgap*Cnad*pi/(Kgapdnad*Kgapdipi*Kgapdigap)+Ccgap*Cnad/(Kgapdigap*Kgapdinad)+Cc13p2g*Cnad/(Kgapdi13p2g*Kgapdinad)+...
    Kgapd13p2g*pi*Cnadh*(CcH*1000)/(Kgapdi13p2g*Kgapdnadh*Kgapdipi)+Ccgap*Cnadh*(CcH*1000)/(Kgapdigap*Kgapdinadh)+Cc13p2g*Cnadh*(CcH*1000)/(Kgapdi13p2g*Kgapdinadh)+Ccgap*Cnad*pi/(Kgapdnad*Kgapdipi*Kgapdigap)+...
    Kgapdgap*Cnad*pi*Cc13p2g/(Kgapdnad*Kgapdipi*Kgapdigap*Kgapdi13p2g1)+pi*Ccgap*Cnadh*(CcH*1000)/(Kgapdipi*Kgapdigap*Kgapdinadh)+Kgapd13p2g*pi*Cc13p2g*Cnadh*(CcH*1000)/(Kgapdipi*Kgapdigap*Kgapdnadh*Kgapdi13p2g1));

% Phosphoglycerate Kinase
rpgk=E0pgk*100*0.00256*Kx18*((kpgkf*Cc13p2g*MgAdp/(Kpgkimgadp*Kpgk13p2g))-((kpgkr*Cc3pg*MgAtp)/(Kpgkimgatp*Kpgk3pg)))/(1+(Cc13p2g/Kpgki13p2g)+(MgAdp/Kpgkimgadp)+(Cc13p2g*MgAdp/(Kpgkimgadp*Kpgk13p2g))+(Cc3pg/Kpgki3pg)+(MgAtp/Kpgkimgatp)+(Cc3pg*MgAdp/(Kpgkimgatp*Kpgk3pg)));

% Bisphosphoglycerate Mutase
rdpgm = k13dpg*Cc13p2g/(1 + Cc23p2g/k23dpg);

% Bisphosphoglycerate Phosphatase
rdpgase = rmdpgase*Cc23p2g/(Kmdpgase + Cc23p2g);

% Phosphoglycerate Mutase
rpgm=E0pgm*Mpgm*Kx115*((kpgmf*Cc3pg/Kpgm3pg)-(kpgmr*Cc2pg/Kpgm2pg))/(1+(Cc3pg/Kpgm3pg)+(Cc2pg/Kpgm2pg));

% Enolase
ren=E0en*14*Kx19*((kenf*Mg*Cc2pg/(Kenimg*Ken2pg))-(kenr*Mg*Ccpep/(Kenimg*Kenpep)))/(1+(Cc2pg/Keni2pg)+(Mg/Kenimg)+(Ccpep/Kenipep)+(Cc2pg*Mg/(Kenimg*Ken2pg))+(Ccpep*Mg/(Kenimg*Kenpep)));

% Pyruvate kinase
Npk = 1 + Lpk*((1+Ccatp/Kpkatp)^4*(1+Ccala/Kpkala)^4)/((1+Ccpep/Kpkpep+Ccpyr/Kpkpyr)^4*(1 +Ccfbp/Kpkfdp+Ccg16p/Kpkg16p)^4);
rpk = E0pk*0.15*Kx110*(rmpkf*(Ccpep/Kpkpep)*(MgAdp/Kpkmgadp)-rmpkr*(Ccpyr/Kpkpyr)*(MgAtp/Kpkmgatp))/(((1+Ccpep/Kpkpep)*(1+MgAdp/Kpkmgadp)+(1+Ccpyr/Kpkpyr)*(1+MgAtp/Kpkmgatp)-1)*Npk);

% Lactate dehydrogenase
rldh = E0ldh*2.0*0.75*Kx111*((kldhf*Cnadh*Ccpyr/(Kldhinadh*Kldhpyr)) - (kldhr*Cnad*Cclac/(Kldhinad*Kldhlac)))/((1+(Kldhnadh*Ccpyr/(Kldhinadh*Kldhpyr))+(Kldhnad*Cclac/(Kldhinad*Kldhlac)))*(1+Ccpyr/Kldhi2pyr)+...
    (Cnadh/Kldhinadh)+(Cnad/Kldhinad)+(Cnadh*Ccpyr/(Kldhinadh*Kldhpyr))+(Kldhnad*Cnadh*Cclac/(Kldhinadh*Kldhlac*Kldhinad))+(Kldhnadh*Cnad*Ccpyr/(Kldhinadh*Kldhpyr*Kldhinad))+...
    (Cnad*Cclac/(Kldhinad*Kldhlac))+(Cnadh*Ccpyr*Cclac/(Kldhinadh*Kldhpyr*Kldhilac))+(Cnad*Cclac*Ccpyr/(Kldhinad*Kldhlac*Kldhipyr)));

% Lactate Export
rmct= E0mct*2.73*1000*1000*(CcH*Cclac-CeHexpt*Celac)/(Kmclacmct*KicHmct+KmcHmct*Cclac+Kmclacmct*CcH*1000+KmcHmct*Celac+Kmclacmct*CeHexpt*1000+1000*CcH*Cclac+...
      1000*CeHexpt*Celac+(KmcHmct*Cclac*CeHexpt*1000/KicHmct)+(KmcHmct*CcH*1000*Celac/KicHmct)+(CcH*10^6*Cclac*CeHexpt/KieHmct)+(CcH*10^6*Celac*CeHexpt/(KicHmct)));

% Malate dehydrogenase 1
rmdh1 = 1*1*fct5*3.737*kmdh1*(Ccoaa*Cnadh - Cnad*Ccmal/Kmdh1eq)/(Kmdh1ia*Kmdh1oaa+ Kmdh1oaa*Cnadh + Kmdh1nadh*Ccoaa + Cnadh*Ccoaa +...
    (Kmdh1nadh*Ccoaa*Cnad/Kmdh1iq) + (Cnadh*Ccoaa*Ccmal/Kmdh1ip) + (Kmdh1ia*Kmdh1oaa/(Kmdh1iq*Kmdh1mal))*(Kmdh1nad*Ccmal + Kmdh1mal*Cnad +...
    (Kmdh1nad*Cnadh*Ccmal/Kmdh1ia) + Cnad*Ccmal + (Ccmal*Ccoaa*Cnad/Kmdh1ib)));

% glutamate-oxaloacetate transaminase 1
rgot1 =  1*fct5*1.37*10^(-7)*kgot1*(Ccasp*Ccakg - Ccoaa*Ccglu/Kgot1eq)/(Kgot1akg*Ccasp + Kgot1asp*Kgot1ai*Ccakg + Ccasp*Ccakg + (Kgot1asp*Ccakg*Ccglu/Kgot1iq) +...
    (Kgot1ia*Kgot1akg/(Kgot1iq*Kgot1oaa))*(Kgot1glu*Kgot1ai*Ccoaa + Kgot1oaa*Ccglu + (Kgot1glu*Ccasp*Ccoaa/Kgot1ia) + Ccoaa*Ccglu));

% glutamate-pyruvate transaminase 2
rgpt1= 1*1.66*10^(-6)*kgpt1f*((Ccala*0.01)*Cmakg-(Cmpyr*Cmglu/Kgpt1eq))/(Kgpt1ala*Cmakg+Kgpt1akg*(Ccala*0.01)+Cmakg*(Ccala*0.01)+(Kgpt1ala*Cmakg*Cmglu/Kgpt1iglu)+(Kgpt1akg*(Ccala*0.01)*(Ccala*0.01)/Kgpt1IA)+...
    (Kgpt1akg*(Ccala*0.01)*Cmglu/Kgpt1RG)+(kgpt1f/(kgpt1r*Kgpt1eq))*(Kgpt1pyr*Cmglu+Kgpt1glu*Cmpyr+Cmpyr*Cmglu+(Kgpt1akg*(Ccala*0.01)*Cmpyr/Kgpt1ipyr)+(Kgpt1pyr*Cmglu*Cmglu/Kgpt1IG)));


% Pentose Phosphate Pathway

% Glutathione Oxidation
rgshox =60000*kgshox*Ccgsh;

% Glutathione Reductase
rgssgr = 42*Mgssgr*Egssgr*(N1gssgr*Cnadph*Gssg - N2gssgr*Ccgsh*Ccgsh*Cnadp)/Detgssgr;

% Glucose 6-Phosphate Dehydrogenase
rg6pd = E0g6pd*45*Mg6pd*Eg6pd*(N1g6pd*Cnadp*Ccg6p- N2g6pd*Cc6pg*Cnadph)/Detg6pd;

% 6-Phosphogluconate Dehydrogenase
r6pgd = 40*4.2*M6pgd*E6pgd*(N16pgd*Cnadp*Cc6pg - N26pgd*Co2*Ccru5p*Cnadph)/Det6pgd;

% Ribose Phosphate Epimerase
rep = 10*Mep*(rmfep*(Ccru5p/Kfmep) - rmrep*(Ccxyl5p/Krmep))/(1 + (Ccru5p/Kfmep) + (Ccxyl5p/Krmep));
          
% Ribose Phosphate Isomerase
rrpi = 15*Mrpi*(rmfki*(Ccru5p/Kfmki) - rmrki*(Ccr5p/Krmki))/(1 + (Ccru5p/Kfmki) + (Ccr5p/Krmki));

% Phosphorybosylpyrophosphate Synthase
rprpps = 23*rmprpps*MgAtp*Ccr5p/((Kprppsatp + MgAtp)*(Kprppsr5p + Ccr5p));

% Transaldolose
rta = 15*(N1ta*Ccsh7p*Ccgap - N2ta*Cce4p*Ccf6p)*Eta/(D1ta*Ccsh7p + D2ta*Ccgap + D3ta*Cce4p + D4ta*Ccf6p + D5ta*Ccsh7p*Ccgap + D6ta*Cce4p*Ccf6p + D7ta*Ccgap*Ccf6p + D8ta*Ccsh7p*Cce4p);

% Transketolase
factor1=(k2tk1+k3tk1);
factor2=(k6tk1+k7tk1);
factor3=(k6tk2+k7tk2);

D1= k1tk1*k3tk1*factor2*factor3;
D2= k5tk1*k7tk1*factor1*factor3;
D3= k2tk1*k4tk1*factor2*factor3;
D4= k6tk1*k8tk1*factor1*factor3;
D5= k5tk2*k7tk2*factor1*factor2;
D6= k6tk2*k8tk2*factor1*factor2;
D7= k1tk1*k5tk1*(k3tk1+k7tk1)*factor3;
D8= k1tk1*k4tk1*factor2*factor3;
D9= k5tk2*k8tk2*factor1*factor2;
D10= k4tk1*k8tk1*(k2tk1+k6tk1)*factor3;
D11= k4tk1*k8tk2*factor2*(k2tk1+k6tk2);
D12= k1tk1*k5tk2*factor2*(k3tk1+k7tk2);
D13= k5tk1*k8tk1*factor1*factor3;
D14= k5tk2*k8tk1*factor1*(k6tk1+k7tk2);
D15= k5tk1*k8tk2*factor1*(k7tk1+k6tk2);


Denominator= D1*Ccxyl5p+ D2*Ccr5p+ D3*Ccgap+ D4*Ccsh7p+ D5*Cce4p+ D6*Ccf6p+ D7*Ccxyl5p*Ccr5p+ D8*Ccxyl5p*Ccgap+...
             D9*Cce4p*Ccf6p+ D10*Ccgap*Ccsh7p+ D11*Ccgap*Ccf6p+ D12*Ccxyl5p*Cce4p+ D13*Ccr5p*Ccsh7p+ D14*Cce4p*Ccsh7p+...
             D15*Ccf6p*Ccr5p;

R1= k1tk1*k3tk1*k5tk1*k7tk1*(k6tk2+k7tk2)*Ccxyl5p*Ccr5p;
R2= k1tk1*k3tk1*k5tk2*k7tk2*(k6tk1+k7tk1)*Ccxyl5p*Cce4p;
R3= k2tk1*k4tk1*k6tk1*k8tk1*(k6tk2+k7tk2)*Ccgap*Ccsh7p;
R4= k2tk1*k4tk1*k6tk2*k8tk2*(k6tk1+k7tk1)*Ccgap*Ccf6p;
R5= k5tk1*k6tk2*k7tk1*k8tk2*(k2tk1+k3tk1)*Ccr5p*Ccf6p;
R6= k5tk2*k6tk1*k7tk2*k8tk1*(k2tk1+k3tk1)*Cce4p*Ccsh7p;

rtkxyl5p= (-R1-R2+R3+R4)*Etk1/Denominator;
rtkgap= (R1+R2-R3-R4)*Etk1/Denominator;
rtkr5p= (-R1-R5+R3+R6)*Etk1/Denominator;
rtksh7p= (R1+R5-R3-R6)*Etk1/Denominator;
rtke4p= (-R2-R6+R4+R5)*Etk1/Denominator;
rtkf6p= (R2+R6-R4-R5)*Etk1/Denominator;


% TCA Cycle

% Pyruvate Dehydrogense
rpdhc = E0pdhc*fct4*2.4*10^(-5)*fct2*kpdhc*Cmpyr*Cmcoash*Cmnad*(1 - (Cmco2*Cmaccoa*Cmnadh)/(Kpdhceq*Cmpyr*Cmcoash*Cmnad))/(Kpdhcnad*Kpdhcai2*Cmpyr*Cmcoash +...
    Kpdhccoa*Kpdhcai1*Cmpyr*Cmnad + Kpdhcpyr*Cmcoash*Cmnad + Cmpyr*Cmcoash*Cmnad);
 
% Citrate Synthase
rcs = fct4*0.00001*5*fct3*fct2*kcs*(Cmoaa*Cmaccoa - (Cmcoash*Cmcit/Kcseq))/(Kcsia*Kcsaccoa*Kcsai1 + Kcsoaa*Kcsai1*Cmaccoa + Kcsaccoa*Kcsai2*Cmoaa + Cmoaa*Cmaccoa);

% Aconitase
racon =  fct4*fct3*fct2*kacon*(Cmcit - Cmicit/Kaconeq)/(Kaconcit + Cmcit + Kaconcit*Cmicit/Kaconicit);

% Isocitrate Dehydrogenase
ridh = fct4*0.001*fct2*kidh*(1 - (Cmakg*Cmnadh*Cmco2/(Kidheq*Cmnad*Cmicit)))/(1 + ((Kidhicit/Cmicit)^nh)*Kidhai + (Kidhnad/Cmnad)*(1 +...
    ((Kidhib/Cmicit)^nh)*Kidhai + (Cmnadh/Kidhiq)*Kidhai));

% Alpha Ketoglutarate Dehydrogenase
rakgd =fct4*fct2*kakgd*(1 - (Cmco2*Cmscoa*Cmnadh/(Kakgdeq*Cmakg*Cmcoash*Cmnad)))/(1 + (Kakgdakg*Kakgdai/Cmakg) + (Kakgdcoa*(1 + Cmscoa/Kakgdiq)/Cmcoash) +...
    (Kakgdnad/Cmnad)*(1 + Cmnadh/Kakgdir));

% Succinate COA Synthase
rscoas = fct4*fct3*fct2*kscoas*(Cmgdp*Cmscoa*Cmpi - Cmcoash*Cmsuc*Cmgtp)/(Kscoasia*Kscoasib*Kscoaspi + Kscoasib*Kscoaspi*Cmgdp + Kscoasia*Kscoasscoa*Cmpi + Kscoaspi*Cmgdp*Cmscoa +...
    Kscoasscoa*Cmgdp*Cmpi + Kscoasgdp*Cmscoa*Cmpi + Cmgdp*Cmscoa*Cmpi + (Kscoasia*Kscoasib*Kscoaspi/(Kscoascoa*Kscoasiq*Kscoasir))*(Kscoasir*Kscoassuc*Cmcoash +...
    Kscoasiq*Kscoascoa*Cmgtp + Kscoasgtp*Cmcoash*Cmsuc + Kscoassuc*Cmcoash*Cmgtp + Kscoascoa*Cmgtp*Cmsuc + Cmcoash*Cmsuc*Cmgtp + (Kscoassuc*Kscoasir*Cmgdp*Cmcoash/Kscoasia) +...
    (Kscoassuc*Kscoasir*Cmgdp*Cmscoa*Cmsuc/(Kscoasia*Kscoasib)) + (Kscoasgtp*Cmgdp*Cmcoash*Cmsuc/Kscoasia) + (Kscoasir*Kscoassuc*Cmgdp*Cmscoa*Cmpi*Cmcoash/(Kscoasia*Kscoasib*Kscoasic)) +...
    (Kscoasip*Kscoasgtp*Cmgdp*Cmscoa*Cmpi*Cmsuc/(Kscoasia*Kscoasib*Kscoasic)) + (Kscoasgtp*Cmgdp*Cmscoa*Cmcoash*Cmsuc/(Kscoasia*Kscoasib)) + (Kscoasgtp*Cmgdp*Cmscoa*Cmpi*Cmcoash*Cmsuc/(Kscoasia*Kscoasib*Kscoasic))) +...
    (Kscoasgdp*Cmscoa*Cmpi*Cmgtp/Kscoasir) + (Kscoasia*Kscoasscoa*Cmpi*Cmgtp/Kscoasir) + (Kscoasia*Kscoasscoa*Cmpi*Cmsuc*Cmgtp/(Kscoasiq*Kscoasir)) + (Kscoasgdp*Cmscoa*Cmpi*Cmgtp*Cmsuc/(Kscoasiq*Kscoasir)) +...
    (Kscoasgdp*Kscoasic*Cmscoa*Cmcoash*Cmsuc*Cmgtp/(Kscoasip*Kscoasiq*Kscoasir)) + (Kscoasia*Kscoasscoa*Cmpi*Cmcoash*Cmsuc*Cmgtp/(Kscoasip*Kscoasiq*Kscoasir)) +...
    (Kscoasgtp*Cmscoa*Cmpi*Cmcoash*Cmsuc*Cmgtp/(Kscoasip*Kscoasiq*Kscoasir)));

% Succinate Dehydrogenase
rsdh = fct4*10*fct3*fct2*ksdh*(Cmsuc*Cmcoq - Cmqh2*Cmfum/Ksdheq)/(Ksdhia*Ksdhcoq*Ksdhai + Ksdhcoq*Cmsuc + Ksdhsuc*Ksdhai*Cmcoq + Cmsuc*Cmcoq + (Ksdhsuc*Cmcoq*Cmfum/Ksdhiq) +...
    (Ksdhia*Ksdhcoq/(Ksdhiq*Ksdhqh2))*(Ksdhfum*Cmqh2*Ksdhai + Ksdhqh2*Cmfum + (Ksdhfum*Cmsuc*Cmqh2/Ksdhia) + Cmqh2*Cmfum));

% Fumarase
rfum = fct4*1*fct3*fct2*kfum*(Cmfum - Cmmal/Kfumeq)/(Kfumfum*Kfumai + Cmfum + Cmmal*Kfumfum/Kfummal);

% Glutamate Dehydrogenase
rgdh= 1*104.67*kcatgdh*(Cmnad*Cmglu-Cmnh3*Cmakg*Cmnadh/Keqgdh)/(Kinadgdh*Kmglugdh+Kmglugdh*Cmnad+Kmnadgdh*Cmglu+...
           Cmglu*Cmnad+Kinadgdh*Kmglugdh*Kmakggdh*Cmnh3/(Kmnh3gdh*Kiakggdh)+Kinadgdh*Kmglugdh*Cmnadh/Kinadhgdh+Kmglugdh*Cmnad*Cmnh3/Kinh3gdh+...
           Kinadgdh*Kmglugdh*Kmnadhgdh*Cmnh3*Cmakg/(Kmnh3gdh*Kiakggdh*Kinadhgdh)+Kmnadgdh*Cmglu*Cmnadh/Kinadhgdh+Kinadgdh*Kmglugdh*Cmakg*Cmnadh/(Kiakggdh*Kinadhgdh)+...
           Kinadgdh*Kmglugdh*Kmakggdh*Cmnh3*Cmnadh/(Kmnh3gdh*Kiakggdh*Kinadhgdh)+Cmnad*Cmglu*Cmnh3/Kiakggdh+Kinadgdh*Kmglugdh*Kmakggdh*Cmnh3*Cmakg*Cmnadh/(Kmnh3gdh*Kiakggdh*Kinadhgdh)+...
           Kmnadhgdh*Kmglugdh*Cmnad*Cmnh3*Cmakg/(Kmnh3gdh*Kiakggdh*Kinadhgdh)+Cmnad*Cmglu*Cmakg/Kiakggdh+Kinadgdh*Kmglugdh*Cmglu*Cmakg*Cmnadh/(Kiglugdh*Kiakggdh*Kinadhgdh)+...
           Cmnad*Cmglu*Cmnh3*Cmakg/(Kinh3gdh*Kiakggdh)+Kinadgdh*Kmglugdh*Cmglu*Cmnh3*Cmakg*Cmnadh/(Kmnh3gdh*Kiglugdh*Kiakggdh*Kinadhgdh));
       

% Malate Dehydrogenase 2
rmdh2 =  1.0*0.1*fct5*1.0454*fct2*kmdh2*(Cmnad*Cmmal - Cmoaa*Cmnadh/Kmdh2eq)/(Kmdh2ia*Kmdh2mal*Kmdh2ai + Kmdh2mal*Cmnad + Kmdh2nad*Kmdh2ai*Cmmal + Cmnad*Cmmal +...
    (Kmdh2nad*Cmmal*Cmnadh/Kmdh2iq) + (Cmnad*Cmmal*Cmoaa/Kmdh2ip) + (Kmdh2ia*Kmdh2mal/(Kmdh2iq*Kmdh2oaa))*(Kmdh2nadh*Kmdh2ai*Cmoaa + Kmdh2oaa*Cmnadh +...
    (Kmdh2nadh*Cmnad*Cmoaa/Kmdh2ia) + Cmnadh*Cmoaa + (Cmmal*Cmoaa*Cmnadh/Kmdh2ib)));

% Glutamate-oxaloacetate transaminase 2
rgot2 = 0.01*fct5*10^(-3)*fct2*kgot2*(Cmasp*Cmakg - Cmoaa*Cmglu/Kgot2eq)/(Kgot2akg*Cmasp + Kgot2asp*Kgot2ai*Cmakg + Cmasp*Cmakg + (Kgot2asp*Cmakg*Cmglu/Kgot2iq) +...
    (Kgot2ia*Kgot2akg/(Kgot2iq*Kgot2oaa))*(Kgot2glu*Kgot2ai*Cmoaa + Kgot2oaa*Cmglu + (Kgot2glu*Cmasp*Cmoaa/Kgot2ia) + Cmoaa*Cmglu));

% Mitochondrial Malic enzyme
rmmalic=5*0.9*1.088*fct2*kcatmalic*(1-(Cmpyr*Co2*Cmnadh)/(Cmmal*Cmnad*Keqmmalic))/(Kmmalmalic/Cmmal*(1+Cmatp/Kiatpmalic)+Kmnadmalic/Cmnad+Kmmalmalic*Kmnadmalic/(Cmmal*Cmnad));

% Cytosolic Malic Enzyme
rcmalic= fct6*70*7.39455*kcatcmalic*(Cnadp*Ccmal-Co2*Ccpyr*Cnadph/Keqcmalic)/(Kinadcmalic*Kmmalcmalic+Kmmalcmalic*Cnadp+Kmnadcmalic*Ccmal+...
           Cnadp*Ccmal+Kinadcmalic*Kmmalcmalic*Kmpyrcmalic*Co2/(KmCo2cmalic*Kipyrcmalic)+Kinadcmalic*Kmmalcmalic*Cnadph/Kinadhcmalic+Kmmalcmalic*Kmpyrcmalic*Cnadp*Co2/(KmCo2cmalic*Kipyrcmalic)+...
           Kmnadcmalic*Ccmal*Cnadph/Kinadhcmalic+Kmmalcmalic*Kmnadhcmalic*Cnadp*Co2*Ccpyr/(KmCo2cmalic*Kipyrcmalic*Kinadhcmalic)+...
           Kmmalcmalic*Kmpyrcmalic*Cnadp*Ccmal*Co2/(Kimalcmalic*KmCo2cmalic*Kipyrcmalic)+Kmmalcmalic*KiCo2cmalic*Kmnadhcmalic*Cnadp*Ccmal*Ccpyr/(Kimalcmalic*KmCo2cmalic*Kipyrcmalic*Kinadhcmalic)+...
           Kmnadcmalic*Ccmal*Ccpyr*Cnadph/(Kipyrcmalic*Kinadhcmalic)+Kmnadcmalic*Ccmal*Co2*Ccpyr*Cnadph/(KiCo2cmalic*Kipyrcmalic*Kinadhcmalic)+Kmmalcmalic*Kinadcmalic*Ccpyr*Cnadph/(Kipyrcmalic*Kinadhcmalic)+...
           Kmmalcmalic*Kmnadhcmalic*Cnadp*Ccmal*Co2*Ccpyr/(Kimalcmalic*KmCo2cmalic*Kipyrcmalic*Kinadhcmalic)+Kinadcmalic*Kmmalcmalic*Kmnadhcmalic*Co2*Ccpyr/(KmCo2cmalic*Kipyrcmalic*Kinadhcmalic)+...
           Kinadcmalic*Kmmalcmalic*Co2*Ccpyr*Cnadph/(KmCo2cmalic*Kipyrcmalic*Kinadhcmalic)+Kinadcmalic*Kmmalcmalic*Kmpyrcmalic*Co2*Cnadph/(KmCo2cmalic*Kipyrcmalic*Kinadhcmalic));
  
% ATP-Citrate Lyase           
rcly=0.7*0.5025*0.001*Vcly*(Cccit*Cccoash) /(Kicitcly*Kmcoashcly + Kmcitcly*Cccoash + Kmaccoacly*Cccit + Cccit*Cccoash +...
    (Kmcitcly*Cccoash*Ccaccoa/Kiaccoacly) + Cccit*Cccoash*Ccoaa/Kioaacly + (Kicitcly*Kmcoashcly/(Kmoaacly*Kiaccoacly))*(Kmaccoacly*Ccoaa + Kmoaacly*Ccaccoa +...
    (Kmaccoacly*Cccit*Ccoaa/Kicitcly) + Ccoaa*Ccaccoa + (Cccoash*Ccoaa*Ccaccoa/Kicoashcly)));        

% Pyruvate-hydrogen co-transporter flux
rpyrh =300*fct2*Kpyrh*(Ccpyr*CcH - Cmpyr*CmH*10);

% Glutamate-hydrogen co-transporter flux
rgluh = 1*0.01465*fct2*Kgluh*(Ccglu*CcH - Cmglu*CmH);

% Citrate-malate exchange flux 
rcitmal =6*fct2*0.008350*Kcitmal*(Cccit*Cmmal - Cmcit*Ccmal); %0.0000000015* 0.0000080*

% Alpha-ketoglutarate-malate exchange flux
rakgmal = 1*0.01*1*fct5*0.01897*fct2*Kakgmal*(Ccakg*Cmmal - Cmakg*Ccmal)/(Kmakgi*Kmmalm*(2 + Ccmal/Kmmali + Cmmal/Kmmalm + Ccakg/Kmakgi + Cmakg/Kmakgm +...
    Ccmal*Cmakg/(Kmmali*Kmakgm) + Cmmal*Ccakg/(Kmmalm*Kmakgi)));

% Malate-phosphate exchange flux
rmalpi = 1*0.0135*fct2*Kmalpi*(Ccmal*Cmpi - Cmmal*Ccpi);

% Aspartate-glutamate exchange flux
raspglu =  fct5*6.85*10^(-3)*fct2*kaspglu*(Keqaspglu*Ccasp*Cmglu*(CmH*1000)*10^(0) - Cmasp*Ccglu*(CcH*1000))/(Keqaspglu*Kiaspi*Kiglum*Khaspglu*(2*m + m*Ccasp/Kiaspi +...
    Ccasp*Cmglu*(CmH*1000)/(Kiaspi*Kiglum*Khaspglu) + m*Cmasp*(CcH*1000)/(Kiaspm*Khaspglu) + Cmasp*Ccglu*(CcH*1000)/(Kiaspm*Kiglui*Khaspglu) +...
    m*Cmasp/Kiaspm + m*Ccasp*(CmH*1000)/(Kiaspi*Khaspglu) + m*(CmH*1000)/Khaspglu + m*Ccglu*(CcH*1000)/(Kiglui*Khaspglu) + m*(CcH*1000)/Khaspglu +...
    m*Cmglu*(CmH*1000)/(Kisglum*Khaspglu)));

% Glutaminase
rgls= 1*fct41*fct4*0.8*0.9*75/7*(Vm/Vc)*Vfgls*(Ccgln-Cmglu/Kglseq)/(Kmglsgln*(1+Cmglu/Kiglsglu)+Ccgln);

% Pyruvate carboxylase
rpc= Vfpc*(Cmpyr*Co2-Cmoaa)/(Kmpyrpc*Kmhco3pc+Kmpyrpc*Co2+Kmhco3pc*Cmpyr+Cmpyr*Co2);


% Rate of ATP generation
ratpgen = 2.5*(rpdhc + ridh + rgdh + rakgd + rmdh2 + rmmalic)*(Vm/Vc)  ...
            + 1.5*(rsdh)*(Vm/Vc) ...
            + (rpgk - rcly) ...
            - (rpc + raspglu)*(Vm/Vc);

% ------------------------------------------------------------------------
%             Extracellular metabolites
% ------------------------------------------------------------------------

relac = qcell * rmct * Ccell;
reglc = qcell * rperm * Ccell;
rI = qI * Ccell;
rco2cell = qcell * Ccell * (-rpc + rmmalic*(Vm/Vc) + rcmalic + r6pgd + rpdhc*(Vm/Vc) + ridh*(Vm/Vc) + rakgd*(Vm/Vc)); % aqueous CO2 accumulation from metabolism

% ------------------------------------------------------------------------
%             pH controller
% ------------------------------------------------------------------------

Flow = 2*120;
pHdead = 0;
pHset = 6.95;
controller = 1;
LGfrac = 20;

Kp = 0.1;
KI = 0.1;

pHdelta = pHset + log10(CeH);

dpHI = KI*pHdelta; % For integral term
pHC = Kp*pHdelta + pHI;
if controller == 1
    % Control action: CO2 for acid, Na2CO3 for base, with deadband
    if pHC < 0 && abs(pHdelta) > pHdead
        Cgco2in = -1*pHC;
        basecon = 0;
    elseif pHC > 0 && abs(pHdelta) > pHdead
        Cgco2in = 0;
        basecon = 1*pHC;
    else
        Cgco2in = 0;
        basecon = 0;
        dpHI = 0; % Reset controller so integral action does not accumulate when in deadband
    end
else
    % Constant CO2 sparge rate
    Cgco2in = Cin;
    basecon = 0;
    dpHI = 0;
end

% ------------------------------------------------------------------------
%             CO2 balance kinetic equations
% ------------------------------------------------------------------------

% Kinetic/equilibrium form

Ceh2o = 55.5;

% CO2 G/L exchange
Hco2 = 2.09e3;
Ceqbco2 = Cgco2/Hco2*Ceh2o;
rco2ex = kla*(Ceqbco2 - Ceco2);

kco2eqb = 1e8;

% CO2 hydration, equilibrium
Kh = 1.5e-3;
kco2hyd = kco2eqb;
rco2hyd = kco2hyd*(Kh*Ceco2 - Ceh2co3);

% Dissolution of bicarbonate, equilibrium
Kah2dis = 2.5e-4;
kh2dis = kco2eqb;
rh2co3dis = kh2dis*(Kah2dis*Ceh2co3 - CeH*Cehco3);

% Carbonate reaction, equilibrium
Kaco3 = 4.69e-11;
kco3 = kco2eqb;
rco3 = kco3*(Kaco3*Cehco3 - CeH*Ceco3);

% ------------------------------------------------------------------------
%             Return rates and mass balance equations
% ------------------------------------------------------------------------

if rateFlag==1
    rate(1,1)=rperm;
    rate(1,2)=rhk;
    rate(1,3)=rpgi;
    rate(1,4)=rpfk;
    rate(1,5)=rald;
    rate(1,6)=rtpi;
    rate(1,7)=rgapd;
    rate(1,8)=rpgk;
    rate(1,9)=rdpgm; 
    rate(1,10)=rdpgase;
    rate(1,11)=rpgm;
    rate(1,12)=ren;
    rate(1,13)=rpk; 
    rate(1,14)=rldh;
    rate(1,15)=rmct;
    rate(1,16)=rgshox; 
    rate(1,17)=rgssgr; 
    rate(1,18)=rg6pd;
    rate(1,19)=r6pgd;
    rate(1,20)=rep;
    rate(1,21)=rrpi;
    rate(1,22)=rprpps;
    rate(1,23)=rtkxyl5p;
    rate(1,24)=rta;
    rate(1,25)=rtke4p;
    rate(1,26)=rpyrh; 
    rate(1,27)=rpdhc*(Vm/Vc);
    rate(1,28)=rcs*(Vm/Vc);
    rate(1,29)=racon*(Vm/Vc);
    rate(1,30)=ridh*(Vm/Vc);
    rate(1,31)=rakgd*(Vm/Vc);
    rate(1,32)=rscoas*(Vm/Vc);
    rate(1,33)=rsdh*(Vm/Vc);
    rate(1,34)=rfum*(Vm/Vc);
    rate(1,35)=rmdh2*(Vm/Vc); 
    rate(1,36)=rgot2*(Vm/Vc);
    rate(1,37)=rgluh*(Vm/Vc);
    rate(1,38)=rcitmal*(Vm/Vc);
    rate(1,39)=rakgmal*(Vm/Vc);
    rate(1,40)=rmalpi*(Vm/Vc);
    rate(1,41)=raspglu*(Vm/Vc);
    rate(1,42)=rmdh1;
    rate(1,43)=rgot1;
    rate(1,44)=rcmalic;
    rate(1,45)=rmmalic*(Vm/Vc);
    rate(1,46)=rgls;
    rate(1,47)=rgdh*(Vm/Vc);
    rate(1,48)=rcly;
    rate(1,49)=rpc*(Vm/Vc);
    rate(1,50)=rgpt1*(Vm/Vc);
    rate(1,51)=rcell;
    rate(1,52)=ratpgen;
    rate(1,53)=Akt;
    rate(1,54)=AMPK;
    rate(1,55)=KBPcomb;
    rate(1,56)=basecon;
    rate(1,57)=Cgco2in;
    dc = rate;
else
    dc = [-reglc % ex Gluc
        rperm-rhk % In Gluc
        rhk-rpgi-rg6pd % In G6P
        rpgi-rpfk-rpfk2+rta+rtkf6p %In F6P
        rpfk-rald % In F126
        rald-rtpi % In Dihydroxyacetone phosphate
        rald+rtpi-rgapd+rtkgap-rta % In Glyceraldehyde 3-Phosphate
        rgapd-rpgk % In 1,3-Bisphosphate glycerate
        rpgk-rpgm % In 3-Phosphoglycerate
        rpgm-ren % In 2-Phosphoglycerate
        ren-rpk % In Phosphoenol Pyruavte
        rpk-rldh-rpyrh+rcmalic % In Pyruvate
        -rpdhc+rmmalic+rpyrh*(Vc/Vm)-rpc+rgpt1 %  m Pyruvate 
        rldh-rmct % In Lactate
        relac % Ex Lactate 
        rg6pd-r6pgd % In 6 Phospoglycerate
        rgssgr-rg6pd-r6pgd-rcmalic % In NADP
        rgssgr-rgshox % In Glutathione
        r6pgd-rep-rrpi % In ru5p
        rep+rtkxyl5p % In xylp
        rrpi-rprpps+rtkr5p % In r5p
        rtksh7p-rta % In sh7p
        rta+rtke4p % In e4p
        rpdhc-rcs % m accoa
        racon-ridh % m icit
        ridh-rakgd+rakgmal-rgot2+rgdh-rgpt1 %m akg
        rldh-rgapd+rmdh1 % In NAD 
        rakgd-rscoas % m scoa
        rscoas-rsdh % m suc
        rsdh-rfum % m fum
        rfum-rmmalic-rmdh2-rakgmal-rcitmal+rmalpi % m mal
        -rcs+rmdh2+rgot2+rpc % m oaa 
        raspglu-rgot2 % m asp
        rgot2-raspglu-rgdh+rgluh+rgpt1 % m glu
        rgot1-rmdh1+rcly % In oaa 
        rmdh1+rakgmal*(Vm/Vc)+rcitmal*(Vm/Vc)-rmalpi*(Vm/Vc)-rcmalic % In mal
        -rgot1-raspglu*(Vm/Vc) % In asp
        -rakgmal*(Vm/Vc)-rgot1 % In akg
        rgot1+raspglu*(Vm/Vc)-rgluh*(Vm/Vc)+rgls % In glu
        rcs-racon+rcitmal % m cit 
        -rcitmal*(Vm/Vc)-rcly % In cit
        rpfk2
        rcell
        rI % inhibitor
        - reglc + 2*relac + 1000*(rco2cell + rco2ex + rh2co3dis + rco3 + 3*basecon) - Ccell/cellosmo*(mu/mumax*0.8+0.2) % osmo balance [metabolism + feed]
        Flow*Cgco2in - Flow*Cgco2 - LGfrac*rco2ex % co2 gas phase
        rco2ex - rco2hyd + rco2cell % co2 extracellular
        rco2hyd - rh2co3dis % h2co3 extracellular
        rh2co3dis + rco3 + relac/1000 % H+ extracellular
        rh2co3dis - rco3 % hco3 extracellular
        rco3 + basecon % co3 extracellular
        dpHI % pH controller integral action
        ]; 
end
end