Flux control coefficients that are zero calculated as near-zero

[luciansmith]

(From Herbert):

"This model is giving the wrong flux control coefficient for reaction vAK,. I think its because the flux through the reaction at steady state is virtually zero (1E-15) and so what we're seeing is numerical noise. I tried the model on copasi and jwsonline and they get the correct answer which is zero for the flux control coefficients. I suspect they are zeroing out the values if the flux is below some threshold."

import tellurium as te
import roadrunner
import teUtils as tu
r = te.loada("""
// Created by libAntimony v2.13.2
function L_PFK(L, CiATP, KiATP, CAMP, KAMP, CF26BP, KF26BP, CF16BP, KF16BP, AT, AM, F16, F26)
  L*((1 + CiATP*(AT/KiATP))/(1 + AT/KiATP))^2*((1 + CAMP*(AM/KAMP))/(1 + AM/KAMP))^2*((1 + CF26BP*F26/KF26BP + CF16BP*F16/KF16BP)/(1 + F26/KF26BP + F16/KF16BP))^2;
end
function R_PFK(KmF6P, KmATP, g, AT, F6)
  1 + F6/KmF6P + AT/KmATP + g*(F6/KmF6P)*(AT/KmATP);
end
function T_PFK(CATP, KmATP, AT)
  1 + CATP*(AT/KmATP);
end
  // Compartments and Species:
  compartment extracellular, cytosol;
  species GLCi in cytosol, G6P in cytosol, F6P in cytosol, F16P in cytosol;
  species TRIO in cytosol, BPG in cytosol, P3G in cytosol, P2G in cytosol;
  species PEP in cytosol, PYR in cytosol, ACE in cytosol, NAD in cytosol;
  species NADH in cytosol, $CO2 in cytosol; # $Glyc in cytosol, $Trh in cytosol,
  species $GLCo in extracellular, $ETOH in cytosol, $GLY in cytosol; # $SUCC in cytosol,
  species ATP in cytosol, ADP in cytosol, AMP in cytosol;
  //species $F26BP in cytosol;
  // Assignment Rules:
  //ATP := (P - ADP)/2;
  //ADP := (SUM_P - (P^2*(1 - 4*KeqAK) + 2*SUM_P*P*(4*KeqAK - 1) + SUM_P^2)^0.5)/(1 - 4*KeqAK);
  //AMP := SUM_P - ATP - ADP;
  // Reactions:
  vGLK: GLCi + ATP -> G6P + ADP; e_vGLK * (cytosol*vGLK_VmGLK/(vGLK_KmGLKGLCi*vGLK_KmGLKATP))*(GLCi*ATP - G6P*ADP/vGLK_KeqGLK)/((1 + GLCi/vGLK_KmGLKGLCi + G6P/vGLK_KmGLKG6P)*(1 + ATP/vGLK_KmGLKATP + ADP/vGLK_KmGLKADP));
  vPGI: G6P -> F6P; e_vPGI * (cytosol*vPGI_VmPGI_2/vPGI_KmPGIG6P_2)*(G6P - F6P/vPGI_KeqPGI_2)/(1 + G6P/vPGI_KmPGIG6P_2 + F6P/vPGI_KmPGIF6P_2);
  vGLYCO: G6P + ATP => ADP; e_vGLYCO * cytosol*vGLYCO_KGLYCOGEN_3;
  vTreha: 2 G6P + ATP => ADP; e_vTreha * cytosol*vTreha_KTREHALOSE;
  vPFK: F6P + ATP => F16P + ADP; e_vPFK * cytosol*vPFK_VmPFK*gR*(F6P/KmPFKF6P)*(ATP/KmPFKATP)*R_PFK(KmPFKF6P, KmPFKATP, gR, ATP, F6P)/(R_PFK(KmPFKF6P, KmPFKATP, gR, ATP, F6P)^2 + L_PFK(Lzero, CiPFKATP, KiPFKATP, CPFKAMP, KPFKAMP, CPFKF26BP, KPFKF26BP, CPFKF16BP, KPFKF16BP, ATP, AMP, F16P, F26BP)*T_PFK(CPFKATP, KmPFKATP, ATP)^2);
  vALD: F16P -> 2 TRIO; e_vALD * (cytosol*vALD_VmALD/vALD_KmALDF16P)*(F16P - (KeqTPI/(1 + KeqTPI))*TRIO*(1/(1 + KeqTPI))*TRIO/vALD_KeqALD)/(1 + F16P/vALD_KmALDF16P + (KeqTPI/(1 + KeqTPI))*TRIO/vALD_KmALDGAP + (1/(1 + KeqTPI))*TRIO/vALD_KmALDDHAP + (KeqTPI/(1 + KeqTPI))*TRIO*(1/(1 + KeqTPI))*TRIO/(vALD_KmALDGAP*vALD_KmALDDHAP) + F16P*(KeqTPI/(1 + KeqTPI))*TRIO/(vALD_KmALDGAPi*vALD_KmALDF16P));
  vGAPDH: TRIO + NAD -> BPG + NADH; e_vGAPDH * cytosol*(vGAPDH_VmGAPDHf*(KeqTPI/(1 + KeqTPI))*TRIO*NAD/(vGAPDH_KmGAPDHGAP*vGAPDH_KmGAPDHNAD) - vGAPDH_VmGAPDHr*BPG*NADH/(vGAPDH_KmGAPDHBPG*vGAPDH_KmGAPDHNADH))/((1 + (KeqTPI/(1 + KeqTPI))*TRIO/vGAPDH_KmGAPDHGAP + BPG/vGAPDH_KmGAPDHBPG)*(1 + NAD/vGAPDH_KmGAPDHNAD + NADH/vGAPDH_KmGAPDHNADH));
  vPGK: BPG + ADP -> P3G + ATP; e_vPGK * (cytosol*vPGK_VmPGK/(vPGK_KmPGKP3G*vPGK_KmPGKATP))*(vPGK_KeqPGK*BPG*ADP - P3G*ATP)/((1 + BPG/vPGK_KmPGKBPG + P3G/vPGK_KmPGKP3G)*(1 + ATP/vPGK_KmPGKATP + ADP/vPGK_KmPGKADP));
  vPGM: P3G -> P2G; e_vPGM * (cytosol*vPGM_VmPGM/vPGM_KmPGMP3G)*(P3G - P2G/vPGM_KeqPGM)/(1 + P3G/vPGM_KmPGMP3G + P2G/vPGM_KmPGMP2G);
  vENO: P2G -> PEP; e_vENO * (cytosol*vENO_VmENO/vENO_KmENOP2G)*(P2G - PEP/vENO_KeqENO)/(1 + P2G/vENO_KmENOP2G + PEP/vENO_KmENOPEP);
  vPYK: PEP + ADP -> PYR + ATP; e_vPYK *(cytosol*vPYK_VmPYK/(vPYK_KmPYKPEP*vPYK_KmPYKADP))*(PEP*ADP - PYR*ATP/vPYK_KeqPYK)/((1 + PEP/vPYK_KmPYKPEP + PYR/vPYK_KmPYKPYR)*(1 + ATP/vPYK_KmPYKATP + ADP/vPYK_KmPYKADP));
  vPDC: PYR => ACE + $CO2; e_vPDC * cytosol*vPDC_VmPDC*(PYR^vPDC_nPDC/vPDC_KmPDCPYR^vPDC_nPDC)/(1 + PYR^vPDC_nPDC/vPDC_KmPDCPYR^vPDC_nPDC);
  vSUC: 2 ACE + 3 NAD + 4 ATP => 3 NADH + 4 ADP; e_vSUC * cytosol*vSUC_KSUCC*ACE;
  vGLT: $GLCo -> GLCi; e_vGLT * (vGLT_VmGLT/vGLT_KmGLTGLCo)*(GLCo - GLCi/vGLT_KeqGLT)/(1 + GLCo/vGLT_KmGLTGLCo + GLCi/vGLT_KmGLTGLCi + 0.91*GLCo*GLCi/(vGLT_KmGLTGLCo*vGLT_KmGLTGLCi));
  vADH: ACE + NADH -> NAD + $ETOH; e_vADH * -cytosol*((vADH_VmADH/(vADH_KiADHNAD*vADH_KmADHETOH))*(NAD*ETOH - NADH*ACE/vADH_KeqADH)/(1 + NAD/vADH_KiADHNAD + vADH_KmADHNAD*ETOH/(vADH_KiADHNAD*vADH_KmADHETOH) + vADH_KmADHNADH*ACE/(vADH_KiADHNADH*vADH_KmADHACE) + NADH/vADH_KiADHNADH + NAD*ETOH/(vADH_KiADHNAD*vADH_KmADHETOH) + vADH_KmADHNADH*NAD*ACE/(vADH_KiADHNAD*vADH_KiADHNADH*vADH_KmADHACE) + vADH_KmADHNAD*ETOH*NADH/(vADH_KiADHNAD*vADH_KmADHETOH*vADH_KiADHNADH) + NADH*ACE/(vADH_KiADHNADH*vADH_KmADHACE) + NAD*ETOH*ACE/(vADH_KiADHNAD*vADH_KmADHETOH*vADH_KiADHACE) + ETOH*NADH*ACE/(vADH_KiADHETOH*vADH_KiADHNADH*vADH_KmADHACE)));
  vG3PDH: TRIO + NADH => NAD + $GLY; e_vG3PDH * (cytosol*vG3PDH_VmG3PDH/(vG3PDH_KmG3PDHDHAP*vG3PDH_KmG3PDHNADH))*((1/(1 + KeqTPI))*TRIO*NADH - GLY*NAD/vG3PDH_KeqG3PDH)/((1 + (1/(1 + KeqTPI))*TRIO/vG3PDH_KmG3PDHDHAP + GLY/vG3PDH_KmG3PDHGLY)*(1 + NADH/vG3PDH_KmG3PDHNADH + NAD/vG3PDH_KmG3PDHNAD));
  vATP: ATP -> ADP; e_vATP * cytosol*vATP_KATPASE*ATP;
  // New reaction to model adelylate kinase
  vAK: 2 ADP -> AMP + ATP; EAK*(k1*ADP^2 - k2*AMP*ATP)
  k1 = 0.45; k2 = 1; EAK = 1;
  // Species initializations:
  GLCi = 0.087;
  G6P = 2.45;
  F6P = 0.62;
  F16P = 5.51;
  TRIO = 0.96;
  BPG = 0;
  P3G = 0.9;
  P2G = 0.12;
  PEP = 0.07;
  PYR = 1.85;
  ACE = 0.17;
  ATP = 4.1//6.31;
  ADP = 0
  AMP = 0
  NAD = 1.2;
  NADH = 0.39;
  # Glyc = 0;
  # Trh = 0;
  CO2 = 1;
  # SUCC = 0;
  GLCo = 50;
  ETOH = 50;
  GLY = 0.15;
  SUM_P = 4.1;
  F26BP = 0.02;
  e_vGLK=1;
  e_vPGI=1;
  e_vGLYCO=1;
  e_vTreha=1;
  e_vPFK=1;
  e_vALD=1;
  e_vGAPDH=1;
  e_vPGK=1;
  e_vPGM=1;
  e_vENO=1;
  e_vPYK=1;
  e_vPDC=1;
  e_vSUC=1;
  e_vGLT=1;
  e_vADH=1;
  e_vG3PDH=1;
  e_vATP=1;
  // Compartment initializations:
  extracellular = 1;
  cytosol = 1;
  // Variable initializations:
  KeqAK = 0.45;
  KeqAK has dimensionless;
  gR = 5.12;
  gR has dimensionless;
  KmPFKF6P = 0.1;
  KmPFKF6P has mM;
  KmPFKATP = 0.71;
  KmPFKATP has mM;
  Lzero = 0.66;
  Lzero has dimensionless;
  CiPFKATP = 100;
  CiPFKATP has dimensionless;
  KiPFKATP = 0.65;
  KiPFKATP has mM;
  CPFKAMP = 0.0845;
  CPFKAMP has dimensionless;
  KPFKAMP = 0.0995;
  KPFKAMP has mM;
  CPFKF26BP = 0.0174;
  CPFKF26BP has dimensionless;
  KPFKF26BP = 0.000682;
  KPFKF26BP has mM;
  CPFKF16BP = 0.397;
  CPFKF16BP has dimensionless;
  KPFKF16BP = 0.111;
  KPFKF16BP has mM;
  CPFKATP = 3;
  CPFKATP has dimensionless;
  KeqTPI = 0.045;
  KeqTPI has dimensionless;
  vGLK_VmGLK = 226.452;
  vGLK_VmGLK has mMpermin;
  vGLK_KmGLKGLCi = 0.08;
  vGLK_KmGLKGLCi has mM;
  vGLK_KmGLKATP = 0.15;
  vGLK_KmGLKATP has mM;
  vGLK_KeqGLK = 3800;
  vGLK_KeqGLK has dimensionless;
  vGLK_KmGLKG6P = 30;
  vGLK_KmGLKG6P has mM;
  vGLK_KmGLKADP = 0.23;
  vGLK_KmGLKADP has mM;
  vPGI_VmPGI_2 = 339.677;
  vPGI_VmPGI_2 has mMpermin;
  vPGI_KmPGIG6P_2 = 1.4;
  vPGI_KmPGIG6P_2 has mM;
  vPGI_KeqPGI_2 = 0.314;
  vPGI_KeqPGI_2 has dimensionless;
  vPGI_KmPGIF6P_2 = 0.3;
  vPGI_KmPGIF6P_2 has mM;
  vGLYCO_KGLYCOGEN_3 = 6;
  vGLYCO_KGLYCOGEN_3 has mMpermin;
  vTreha_KTREHALOSE = 2.4;
  vTreha_KTREHALOSE has mMpermin;
  vPFK_VmPFK = 182.903;
  vPFK_VmPFK has mMpermin;
  vALD_VmALD = 322.258;
  vALD_VmALD has mMpermin;
  vALD_KmALDF16P = 0.3;
  vALD_KmALDF16P has mM;
  vALD_KeqALD = 0.069;
  vALD_KeqALD has dimensionless;
  vALD_KmALDGAP = 2;
  vALD_KmALDGAP has mM;
  vALD_KmALDDHAP = 2.4;
  vALD_KmALDDHAP has mM;
  vALD_KmALDGAPi = 10;
  vALD_KmALDGAPi has mM;
  vGAPDH_VmGAPDHf = 1184.52;
  vGAPDH_VmGAPDHf has mMpermin;
  vGAPDH_KmGAPDHGAP = 0.21;
  vGAPDH_KmGAPDHGAP has mM;
  vGAPDH_KmGAPDHNAD = 0.09;
  vGAPDH_KmGAPDHNAD has mM;
  vGAPDH_VmGAPDHr = 6549.8;
  vGAPDH_VmGAPDHr has mMpermin;
  vGAPDH_KmGAPDHBPG = 0.0098;
  vGAPDH_KmGAPDHBPG has mM;
  vGAPDH_KmGAPDHNADH = 0.06;
  vGAPDH_KmGAPDHNADH has mM;
  vPGK_VmPGK = 1306.45;
  vPGK_VmPGK has mMpermin;
  vPGK_KmPGKP3G = 0.53;
  vPGK_KmPGKP3G has mM;
  vPGK_KmPGKATP = 0.3;
  vPGK_KmPGKATP has mM;
  vPGK_KeqPGK = 3200;
  vPGK_KeqPGK has dimensionless;
  vPGK_KmPGKBPG = 0.003;
  vPGK_KmPGKBPG has mM;
  vPGK_KmPGKADP = 0.2;
  vPGK_KmPGKADP has mM;
  vPGM_VmPGM = 2525.81;
  vPGM_VmPGM has mMpermin;
  vPGM_KmPGMP3G = 1.2;
  vPGM_KmPGMP3G has mM;
  vPGM_KeqPGM = 0.19;
  vPGM_KeqPGM has dimensionless;
  vPGM_KmPGMP2G = 0.08;
  vPGM_KmPGMP2G has mM;
  vENO_VmENO = 365.806;
  vENO_VmENO has mMpermin;
  vENO_KmENOP2G = 0.04;
  vENO_KmENOP2G has mM;
  vENO_KeqENO = 6.7;
  vENO_KeqENO has dimensionless;
  vENO_KmENOPEP = 0.5;
  vENO_KmENOPEP has mM;
  vPYK_VmPYK = 1088.71;
  vPYK_VmPYK has mMpermin;
  vPYK_KmPYKPEP = 0.14;
  vPYK_KmPYKPEP has mM;
  vPYK_KmPYKADP = 0.53;
  vPYK_KmPYKADP has mM;
  vPYK_KeqPYK = 6500;
  vPYK_KeqPYK has dimensionless;
  vPYK_KmPYKPYR = 21;
  vPYK_KmPYKPYR has mM;
  vPYK_KmPYKATP = 1.5;
  vPYK_KmPYKATP has mM;
  vPDC_VmPDC = 174.194;
  vPDC_VmPDC has mMpermin;
  vPDC_nPDC = 1.9;
  vPDC_nPDC has dimensionless;
  vPDC_KmPDCPYR = 4.33;
  vPDC_KmPDCPYR has mM;
  vSUC_KSUCC = 21.4;
  vGLT_VmGLT = 97.264;
  vGLT_VmGLT has mmolepermin;
  vGLT_KmGLTGLCo = 1.1918;
  vGLT_KmGLTGLCo has mM;
  vGLT_KeqGLT = 1;
  vGLT_KeqGLT has mM;
  vGLT_KmGLTGLCi = 1.1918;
  vGLT_KmGLTGLCi has mM;
  vADH_VmADH = 810;
  vADH_VmADH has mMpermin;
  vADH_KiADHNAD = 0.92;
  vADH_KiADHNAD has mM;
  vADH_KmADHETOH = 17;
  vADH_KmADHETOH has mM;
  vADH_KeqADH = 6.9e-05;
  vADH_KeqADH has dimensionless;
  vADH_KmADHNAD = 0.17;
  vADH_KmADHNAD has mM;
  vADH_KmADHNADH = 0.11;
  vADH_KmADHNADH has mM;
  vADH_KiADHNADH = 0.031;
  vADH_KiADHNADH has mM;
  vADH_KmADHACE = 1.11;
  vADH_KiADHACE = 1.1;
  vADH_KiADHETOH = 90;
  vG3PDH_VmG3PDH = 70.15;
  vG3PDH_KmG3PDHDHAP = 0.4;
  vG3PDH_KmG3PDHDHAP has mM;
  vG3PDH_KmG3PDHNADH = 0.023;
  vG3PDH_KmG3PDHNADH has mM;
  vG3PDH_KeqG3PDH = 4300;
  vG3PDH_KeqG3PDH has dimensionless;
  vG3PDH_KmG3PDHGLY = 1;
  vG3PDH_KmG3PDHGLY has mM;
  vG3PDH_KmG3PDHNAD = 0.93;
  vG3PDH_KmG3PDHNAD has mM;
  vATP_KATPASE = 33.7;
  vATP_KATPASE has permin;
  // Other declarations:
  const extracellular, cytosol, KeqAK, gR, KmPFKF6P, KmPFKATP, Lzero, CiPFKATP;
  const KiPFKATP, CPFKAMP, KPFKAMP, CPFKF26BP, KPFKF26BP, CPFKF16BP, KPFKF16BP;
  const CPFKATP, KeqTPI;
""")
m = r.simulate (0, 10, 100)
#r.plot()
r.setDiffStepSize(0.05)
r.steadyState()
tu.plotting.plotFloatingSpecies(r)
tu.plotting.plotReactionRates(r)
print (r.ATP + r.ADP + r.AMP)
print ('2ATP + ADP = ', 2*r.ATP + r.ADP)
print ('Pyr = ', r.PYR)
print ('NAD = ', r.NAD)
print ('TRIO = ', r.TRIO)
print ('F16P = ', r.F16P)
print ('P3G = ', r.P3G)
print ('PEP = ', r.PEP)
print ('BPG = ', r.BPG)
print ('F6P = ', r.F6P)
print ('P2G = ', r.P2G)
print ('ACE = ', r.ACE)
print ('G6P = ', r.G6P)
print ('vPFK = ', r.vPFK)
# refAK = r.vAK
# r.EAK = r.vAK*1.02
# r.steadyState()
# cc = ((r.vAK - refAK)/refAK)/0.02
# print (cc)

[hsauro]

Starting at line 3636 in getCC in rrRoadRunner.cpp, we see

    double variableValue = getVariableValue(variableType, variableIndex);
    double parameterValue = impl->getParameterValue(parameterType, parameterIndex);
    return getuCC(variableNameMod, parameterName) * parameterValue / variableValue;

if the variableType == vtFlux and if abs (variableValue) < fluxThreshold then return zero.

The fluxThreshold should be settable just like we can set other constants like r.setDiffStepSize(0.05)

The default fluxThreshold can be set to 1E-9 which I think COPASI might be using. One thing I noticed is that copasi are not doing any zeroing when they compute the unscaled coefficients. It seems to happen when the scaling occurs. Their code isn't that transparent.

Use something like the following code to replace the above:

if (variableType == vtFlux) {
        double variableValue = getVariableValue(variableType, variableIndex);
         if (abs (variableValue) > fluxThreshold) {
           double parameterValue = impl->getParameterValue(parameterType, parameterIndex);
           return getuCC(variableNameMod, parameterName) * parameterValue / variableValue;
         } else
           return 0.0;
}

[luciansmith]

Fixed!