search

sys-bio / roadrunner

libRoadRunner: A high-performance SBML simulator
http://libroadrunner.org/
Other
39 stars 25 forks source link

Difference in solution between copasi and libRoadRunner #420

Closed hsauro closed 6 years ago

hsauro commented 7 years ago

Run the SBML file from

http://seek.virtual-liver.de/models/5 (also enclosed) Koenig2014_Hepatic_Glucose_Model_annotated.zip

print r.getFloatingSpeciesIds()
m = r.simulate (0, 2000000, 100)
print "After simulating:"
print r.dv()
print r.getReducedJacobian()

Compare the reduced jacobian produced by COPASI and libRoadRunner, they are different in 13 places.

Screen-shotrs of both enclosed. Posisible problem with compartments?

t1 t2

0u812 commented 7 years ago

I think this is a case of COPASI using an amount-based Jacobian by default, whereas we apparently use a concentration-based Jacobian. Consider the following simple model:

model *mymodel()
  // Compartments and Species:
  compartment c0, c1;
  species A in c0, B in c1;
  // Reactions:
  J0: A -> ; A-5;
  J1: B -> ; B-5;
  // Species initializations:
  A = 10;
  B = 10;
  // Compartment initializations:
  c0 = 1;
  c1 = 10;
end

If we compute the steady state and get the Jacobian, it is:

A B
A -1 0
B 0 -1

However, if we load this model into COPASI, compute the steady state, and get the Jacobian, we get:

A B
A -1 0
B 0 -0.1

This is because the concentration of B is 5, but its amount is ten times that, hence the factor of 1/10. We can make roadrunner compute an amount-based Jacobian:

import roadrunner
roadrunner.Config.setValue(roadrunner.Config.ROADRUNNER_JACOBIAN_MODE, roadrunner.Config.ROADRUNNER_JACOBIAN_MODE_AMOUNTS)

Now the Jacobian matches COPASI's:

A B
A -1 0
B 0 -0.1

I haven't checked extensively to see if that's what's going on with the liver model, but I'm 99% sure that's what it is.

hsauro commented 7 years ago

I tried your fix and it doesn't make any different to Matthias's model, Still can't find stady state

hsauro commented 7 years ago

I can confirm that the time course simulations match Copasi. The problem is in computing the steady state using steadyState().

PS It appears to make no difference to the time course whether I compute the Jacobian in amounts of concentrations.

Hereis the model that breaks the steady state solver (Running 2.0.3)

import tellurium as te
import roadrunner

r = te.loada('''
// Created by libAntimony v2.9.4
function maximum(x, y)
  piecewise(x, x > y, y);
end

function minimum(x, y)
  piecewise(x, x < y, y);
end

model *Koenig2014_Hepatic_Glucose_Model_annotated()

  // Compartments and Species:
  compartment extern, cyto, mito;
  species $atp in cyto, $adp in cyto, $amp in cyto, utp in cyto, udp in cyto;
  species gtp in cyto, gdp in cyto, $nad in cyto, $nadh in cyto, $phos in cyto;
  species pp in cyto, $co2 in cyto, $h2o in cyto, $h in cyto, glc1p in cyto;
  species udpglc in cyto, glyglc in cyto, glc in cyto, glc6p in cyto, fru6p in cyto;
  species fru16bp in cyto, fru26bp in cyto, grap in cyto, dhap in cyto, bpg13 in cyto;
  species pg3 in cyto, pg2 in cyto, pep in cyto, pyr in cyto, oaa in cyto;
  species lac in cyto, $glc_ext in extern, $lac_ext in extern, $co2_mito in mito;
  species $phos_mito in mito, oaa_mito in mito, pep_mito in mito, $acoa_mito in mito;
  species pyr_mito in mito, $cit_mito in mito, $atp_mito in mito, $adp_mito in mito;
  species gtp_mito in mito, gdp_mito in mito, $coa_mito in mito, $nadh_mito in mito;
  species $nad_mito in mito, $h2o_mito in mito, $h_mito in mito;

  // Assignment Rules:
  nadh_tot := nadh + nad;
  atp_tot := atp + adp + amp;
  utp_tot := utp + udp + udpglc;
  gtp_tot := gtp + gdp;
  nadh_mito_tot := nadh_mito + nad_mito;
  atp_mito_tot := atp_mito + adp_mito;
  gtp_mito_tot := gtp_mito + gdp_mito;
  ins := x_ins2 + (x_ins1 - x_ins2)*glc_ext^x_ins4/(glc_ext^x_ins4 + x_ins3^x_ins4);
  ins_norm := maximum(0 pmol_per_l, ins - x_ins2);
  glu := x_glu2 + (x_glu1 - x_glu2)*(1 dimensionless - glc_ext^x_glu4/(glc_ext^x_glu4 + x_glu3^x_glu4));
  glu_norm := maximum(0 pmol_per_l, glu - x_glu2);
  epi := x_epi2 + (x_epi1 - x_epi2)*(1 dimensionless - glc_ext^x_epi4/(glc_ext^x_epi4 + x_epi3^x_epi4));
  epi_norm := maximum(0 pmol_per_l, epi - x_epi2);
  gamma := 0.5 dimensionless*((1 dimensionless - ins_norm/(ins_norm + K_ins)) + maximum(glu_norm/(glu_norm + K_glu), epi_f*epi_norm/(epi_norm + K_epi)));
  GK_gc_free := (glc^GK_n_gkrp/(glc^GK_n_gkrp + GK_km_glc1^GK_n_gkrp))*(1 dimensionless - GK_b*fru6p/(fru6p + GK_km_fru6p));
  GS_storage_factor := (1 dimensionless + GS_k1_max)*(GS_C - glyglc)/((GS_C - glyglc) + GS_k1_max*GS_C);
  GS_k_udpglc_native := GS_k1_nat/(glc6p + GS_k2_nat);
  GS_k_udpglc_phospho := GS_k1_phospho/(glc6p + GS_k2_phospho);
  GS_native := scale_glyglc*GS_Vmax*GS_storage_factor*udpglc/(GS_k_udpglc_native + udpglc);
  GS_phospho := scale_glyglc*GS_Vmax*GS_storage_factor*udpglc/(GS_k_udpglc_phospho + udpglc);
  GP_fmax := (1 dimensionless + GP_k1_max)*glyglc/(glyglc + GP_k1_max*GP_C);
  GP_vmax_native := scale_glyglc*GP_Vmax*GP_fmax*(GP_base_amp_native + (GP_max_amp_native - GP_base_amp_native)*amp/(amp + GP_ka_amp_native));
  GP_native := (GP_vmax_native/(GP_k_glyc_native*GP_k_p_native))*(glyglc*phos - glc1p/GP_keq)/((1 dimensionless + glyglc/GP_k_glyc_native)*(1 dimensionless + phos/GP_k_p_native) + 1 dimensionless + glc1p/GP_k_glc1p_native - 1 dimensionless);
  GP_vmax_phospho := scale_glyglc*GP_Vmax*GP_fmax*exp((-ln(2 dimensionless)/GP_ki_glc_phospho)*glc);
  GP_phospho := (GP_vmax_phospho/(GP_k_glyc_phospho*GP_k_p_phospho))*(glyglc*phos - glc1p/GP_keq)/((1 dimensionless + glyglc/GP_k_glyc_phospho)*(1 dimensionless + phos/GP_k_p_phospho) + 1 dimensionless + glc1p/GP_k_glc1p_phospho - 1 dimensionless);
  PFK2_native := (scale_gly*PFK2_Vmax*fru6p^PFK2_n_native/(fru6p^PFK2_n_native + PFK2_k_fru6p_native^PFK2_n_native))*atp/(atp + PFK2_k_atp_native);
  PFK2_phospho := (scale_gly*PFK2_Vmax*fru6p^PFK2_n_phospho/(fru6p^PFK2_n_phospho + PFK2_k_fru6p_phospho^PFK2_n_phospho))*atp/(atp + PFK2_k_atp_phospho);
  FBP2_native := (scale_gly*FBP2_Vmax/(1 dimensionless + fru6p/FBP2_ki_fru6p_native))*fru26bp/(FBP2_km_fru26bp_native + fru26bp);
  FBP2_phospho := (scale_gly*FBP2_Vmax/(1 dimensionless + fru6p/FBP2_ki_fru6p_phospho))*fru26bp/(FBP2_km_fru26bp_phospho + fru26bp);
  PK_f := fru16bp^PK_n_fbp/(PK_k_fbp^PK_n_fbp + fru16bp^PK_n_fbp);
  PK_f_p := fru16bp^PK_n_fbp_p/(PK_k_fbp_p^PK_n_fbp_p + fru16bp^PK_n_fbp_p);
  PK_alpha_inp := (1 dimensionless - PK_f)*(PK_alpha - PK_alpha_end) + PK_alpha_end;
  PK_alpha_p_inp := (1 dimensionless - PK_f_p)*(PK_alpha_p - PK_alpha_end) + PK_alpha_end;
  PK_pep_inp := (1 dimensionless - PK_f)*(PK_k_pep - PK_k_pep_end) + PK_k_pep_end;
  PK_pep_p_inp := (1 dimensionless - PK_f_p)*(PK_k_pep_p - PK_k_pep_end) + PK_k_pep_end;
  PK_native := ((scale_gly*PK_Vmax*PK_alpha_inp*pep^PK_n/(PK_pep_inp^PK_n + pep^PK_n))*adp/(adp + PK_k_adp))*(PK_base_act + (1 - PK_base_act)*PK_f);
  PK_phospho := ((scale_gly*PK_Vmax*PK_alpha_p_inp*pep^PK_n_p/(PK_pep_p_inp^PK_n_p + pep^PK_n_p))*adp/(adp + PK_k_adp))*(PK_base_act_p + (1 - PK_base_act_p)*PK_f_p);
  PDH_base := ((scale_gly*PDH_Vmax*pyr_mito/(pyr_mito + PDH_k_pyr))*nad_mito/(nad_mito + PDH_k_nad*(1 dimensionless + nadh_mito/PDH_ki_nadh)))*coa_mito/(coa_mito + PDH_k_coa*(1 dimensionless + acoa_mito/PDH_ki_acoa));
  PDH_nat := PDH_base*PDH_alpha_nat;
  PDH_p := PDH_base*PDH_alpha_p;
  HGP := GLUT2*conversion_factor;
  GNG := GPI*conversion_factor;
  GLY := -G16PI*conversion_factor;

  // Reactions:
  GLUT2: $glc_ext -> glc; (scale_gly*GLUT2_Vmax/GLUT2_km)*(glc_ext - glc/GLUT2_keq)/(1 dimensionless + glc_ext/GLUT2_km + glc/GLUT2_km);
  GPI: glc6p -> fru6p; (scale_gly*GPI_Vmax/GPI_km_glc6p)*(glc6p - fru6p/GPI_keq)/(1 dimensionless + glc6p/GPI_km_glc6p + fru6p/GPI_km_fru6p);
  G16PI: glc1p -> glc6p; (scale_glyglc*G16PI_Vmax/G16PI_km_glc1p)*(glc1p - glc6p/G16PI_keq)/(1 dimensionless + glc1p/G16PI_km_glc1p + glc6p/G16PI_km_glc6p);
  GK: glc + $atp => glc6p + $adp; (scale_gly*GK_Vmax*GK_gc_free*atp/(GK_km_atp + atp))*glc^GK_n/(glc^GK_n + GK_km_glc^GK_n);
  G6PASE: glc6p + $h2o => glc + $phos; scale_gly*G6PASE_Vmax*glc6p/(G6PASE_km_glc6p + glc6p);
  UPGASE: utp + glc1p -> udpglc + pp; (scale_glyglc*UPGASE_Vmax/(UPGASE_km_utp*UPGASE_km_glc1p))*(utp*glc1p - udpglc*pp/UPGASE_keq)/((1 dimensionless + utp/UPGASE_km_utp)*(1 dimensionless + glc1p/UPGASE_km_glc1p) + (1 dimensionless + udpglc/UPGASE_km_udpglc)*(1 dimensionless + pp/UPGASE_km_pp) - 1 dimensionless);
  PPASE: pp + $h2o => 2 $phos; scale_glyglc*PPASE_Vmax*pp/(pp + PPASE_km_pp);
  GS: udpglc => udp + glyglc; (1 dimensionless - gamma)*GS_native + gamma*GS_phospho;
  GP: glyglc + $phos -> glc1p; (1 dimensionless - gamma)*GP_native + gamma*GP_phospho;
  NDKGTP: $atp + gdp -> $adp + gtp; (scale_gly*NDKGTP_Vmax/(NDKGTP_km_atp*NDKGTP_km_gdp))*(atp*gdp - adp*gtp/NDKGTP_keq)/((1 dimensionless + atp/NDKGTP_km_atp)*(1 dimensionless + gdp/NDKGTP_km_gdp) + (1 dimensionless + adp/NDKGTP_km_adp)*(1 dimensionless + gtp/NDKGTP_km_gtp) - 1 dimensionless);
  NDKUTP: $atp + udp -> $adp + utp; (scale_glyglc*NDKUTP_Vmax/(NDKUTP_km_atp*NDKUTP_km_udp))*(atp*udp - adp*utp/NDKUTP_keq)/((1 dimensionless + atp/NDKUTP_km_atp)*(1 dimensionless + udp/NDKUTP_km_udp) + (1 dimensionless + adp/NDKUTP_km_adp)*(1 dimensionless + utp/NDKUTP_km_utp) - 1 dimensionless);
  AK: $atp + $amp -> 2 $adp; (scale_gly*AK_Vmax/(AK_km_atp*AK_km_amp))*(atp*amp - adp*adp/AK_keq)/((1 dimensionless + atp/AK_km_atp)*(1 dimensionless + amp/AK_km_amp) + (1 dimensionless + adp/AK_km_adp)*(1 dimensionless + adp/AK_km_adp) - 1 dimensionless);
  PFK2: fru6p + $atp => fru26bp + $adp; (1 dimensionless - gamma)*PFK2_native + gamma*PFK2_phospho;
  FBP2: fru26bp => fru6p + $phos; (1 dimensionless - gamma)*FBP2_native + gamma*FBP2_phospho;
  PFK1: fru6p + $atp => fru16bp + $adp; scale_gly*PFK1_Vmax*(1 dimensionless - 1 dimensionless/(1 dimensionless + fru26bp/PFK1_ka_fru26bp))*fru6p*atp/(PFK1_ki_fru6p*PFK1_km_atp + PFK1_km_fru6p*atp + PFK1_km_atp*fru6p + atp*fru6p);
  FBP1: fru16bp + $h2o => fru6p + $phos; (scale_gly*FBP1_Vmax/(1 dimensionless + fru26bp/FBP1_ki_fru26bp))*fru16bp/(fru16bp + FBP1_km_fru16bp);
  ALD: fru16bp -> grap + dhap; (scale_gly*ALD_Vmax/ALD_km_fru16bp)*(fru16bp - grap*dhap/ALD_keq)/(1 dimensionless + fru16bp/ALD_km_fru16bp + grap/ALD_ki1_grap + dhap*(grap + ALD_km_grap)/(ALD_km_dhap*ALD_ki1_grap) + fru16bp*grap/(ALD_km_fru16bp*ALD_ki2_grap));
  TPI: dhap -> grap; (scale_gly*TPI_Vmax/TPI_km_dhap)*(dhap - grap/TPI_keq)/(1 dimensionless + dhap/TPI_km_dhap + grap/TPI_km_grap);
  GAPDH: grap + $phos + $nad -> bpg13 + $nadh + $h; (scale_gly*GAPDH_Vmax/(GAPDH_k_nad*GAPDH_k_grap*GAPDH_k_p))*(nad*grap*phos - bpg13*nadh/GAPDH_keq)/((1 dimensionless + nad/GAPDH_k_nad)*(1 dimensionless + grap/GAPDH_k_grap)*(1 dimensionless + phos/GAPDH_k_p) + (1 dimensionless + nadh/GAPDH_k_nadh)*(1 dimensionless + bpg13/GAPDH_k_bpg13) - 1 dimensionless);
  PGK: $adp + bpg13 -> $atp + pg3; (scale_gly*PGK_Vmax/(PGK_k_adp*PGK_k_bpg13))*(adp*bpg13 - atp*pg3/PGK_keq)/((1 dimensionless + adp/PGK_k_adp)*(1 dimensionless + bpg13/PGK_k_bpg13) + (1 dimensionless + atp/PGK_k_atp)*(1 dimensionless + pg3/PGK_k_pg3) - 1 dimensionless);
  PGM: pg3 -> pg2; scale_gly*PGM_Vmax*(pg3 - pg2/PGM_keq)/(pg3 + PGM_k_pg3*(1 dimensionless + pg2/PGM_k_pg2));
  EN: pg2 -> pep; scale_gly*EN_Vmax*(pg2 - pep/EN_keq)/(pg2 + EN_k_pg2*(1 dimensionless + pep/EN_k_pep));
  PK: pep + $adp => pyr + $atp; (1 dimensionless - gamma)*PK_native + gamma*PK_phospho;
  PEPCK: oaa + gtp -> pep + gdp + $co2; (scale_gly*PEPCK_Vmax/(PEPCK_k_oaa*PEPCK_k_gtp))*(oaa*gtp - pep*gdp*co2/PEPCK_keq)/((1 dimensionless + oaa/PEPCK_k_oaa)*(1 dimensionless + gtp/PEPCK_k_gtp) + (1 dimensionless + pep/PEPCK_k_pep)*(1 dimensionless + gdp/PEPCK_k_gdp)*(1 dimensionless + co2/PEPCK_k_co2) - 1 dimensionless);
  PEPCKM: oaa_mito + gtp_mito -> pep_mito + gdp_mito + $co2_mito; (scale_gly*PEPCKM_Vmax/(PEPCK_k_oaa*PEPCK_k_gtp))*(oaa_mito*gtp_mito - pep_mito*gdp_mito*co2_mito/PEPCK_keq)/((1 dimensionless + oaa_mito/PEPCK_k_oaa)*(1 dimensionless + gtp_mito/PEPCK_k_gtp) + (1 dimensionless + pep_mito/PEPCK_k_pep)*(1 dimensionless + gdp_mito/PEPCK_k_gdp)*(1 dimensionless + co2_mito/PEPCK_k_co2) - 1 dimensionless);
  PC: $atp_mito + pyr_mito + $co2_mito => oaa_mito + $adp_mito + $phos_mito; (((scale_gly*PC_Vmax*atp_mito/(PC_k_atp + atp_mito))*pyr_mito/(PC_k_pyr + pyr_mito))*co2_mito/(PC_k_co2 + co2_mito))*acoa_mito^PC_n/(acoa_mito^PC_n + PC_k_acoa^PC_n);
  LDH: pyr + $nadh -> lac + $nad; (scale_gly*LDH_Vmax/(LDH_k_pyr*LDH_k_nadh))*(pyr*nadh - lac*nad/LDH_keq)/((1 dimensionless + nadh/LDH_k_nadh)*(1 dimensionless + pyr/LDH_k_pyr) + (1 dimensionless + lac/LDH_k_lac)*(1 dimensionless + nad/LDH_k_nad) - 1 dimensionless);
  LACT: $lac_ext -> lac; (scale_gly*LACT_Vmax/LACT_k_lac)*(lac_ext - lac/LACT_keq)/(1 dimensionless + lac_ext/LACT_k_lac + lac/LACT_k_lac);
  PYRTM: pyr -> pyr_mito; (scale_gly*PYRTM_Vmax/PYRTM_k_pyr)*(pyr - pyr_mito/PYRTM_keq)/(1 dimensionless + pyr/PYRTM_k_pyr + pyr_mito/PYRTM_k_pyr);
  PEPTM: pep_mito -> pep; (scale_gly*PEPTM_Vmax/PEPTM_k_pep)*(pep_mito - pep/PEPTM_keq)/(1 dimensionless + pep/PEPTM_k_pep + pep_mito/PEPTM_k_pep);
  PDH: pyr_mito + $coa_mito + $nad_mito => $acoa_mito + $co2_mito + $nadh_mito + $h_mito; (1 dimensionless - gamma)*PDH_nat + gamma*PDH_p;
  CS: $acoa_mito + oaa_mito + $h2o_mito -> $cit_mito + $coa_mito; (scale_gly*CS_Vmax/(CS_k_oaa*CS_k_acoa))*(acoa_mito*oaa_mito - cit_mito*coa_mito/CS_keq)/((1 dimensionless + acoa_mito/CS_k_acoa)*(1 dimensionless + oaa_mito/CS_k_oaa) + (1 dimensionless + cit_mito/CS_k_cit)*(1 dimensionless + coa_mito/CS_k_coa) - 1 dimensionless);
  NDKGTPM: $atp_mito + gdp_mito -> $adp_mito + gtp_mito; (scale_gly*NDKGTPM_Vmax/(NDKGTPM_k_atp*NDKGTPM_k_gdp))*(atp_mito*gdp_mito - adp_mito*gtp_mito/NDKGTPM_keq)/((1 dimensionless + atp_mito/NDKGTPM_k_atp)*(1 dimensionless + gdp_mito/NDKGTPM_k_gdp) + (1 dimensionless + adp_mito/NDKGTPM_k_adp)*(1 dimensionless + gtp_mito/NDKGTPM_k_gtp) - 1 dimensionless);
  OAAFLX:  => oaa_mito; scale_gly*OAAFLX_Vmax;
  ACOAFLX: $acoa_mito => ; scale_gly*ACOAFLX_Vmax;
  CITFLX: $cit_mito => ; scale_gly*CITFLX_Vmax;

  // Species initializations:
  atp = 2.8;
  atp has mM;
  adp = 0.8;
  adp has mM;
  amp = 0.16;
  amp has mM;
  utp = 0.27;
  utp has mM;
  udp = 0.09;
  udp has mM;
  gtp = 0.29;
  gtp has mM;
  gdp = 0.1;
  gdp has mM;
  nad = 1.22;
  nad has mM;
  nadh = 0.00056;
  nadh has mM;
  phos = 5;
  phos has mM;
  pp = 0.008;
  pp has mM;
  co2 = 5;
  co2 has mM;
  h2o = 0;
  h2o has mM;
  h = 0;
  h has mM;
  glc1p = 0.012;
  glc1p has mM;
  udpglc = 0.38;
  udpglc has mM;
  glyglc = 250;
  glyglc has mM;
  glc = 5;
  glc has mM;
  glc6p = 0.12;
  glc6p has mM;
  fru6p = 0.05;
  fru6p has mM;
  fru16bp = 0.02;
  fru16bp has mM;
  fru26bp = 0.004;
  fru26bp has mM;
  grap = 0.1;
  grap has mM;
  dhap = 0.03;
  dhap has mM;
  bpg13 = 0.3;
  bpg13 has mM;
  pg3 = 0.27;
  pg3 has mM;
  pg2 = 0.03;
  pg2 has mM;
  pep = 0.15;
  pep has mM;
  pyr = 0.1;
  pyr has mM;
  oaa = 0.01;
  oaa has mM;
  lac = 0.5;
  lac has mM;
  glc_ext = 3;
  glc_ext has mM;
  lac_ext = 1.2;
  lac_ext has mM;
  co2_mito = 5;
  co2_mito has mM;
  phos_mito = 5;
  phos_mito has mM;
  oaa_mito = 0.01;
  oaa_mito has mM;
  pep_mito = 0.15;
  pep_mito has mM;
  acoa_mito = 0.04;
  acoa_mito has mM;
  pyr_mito = 0.1;
  pyr_mito has mM;
  cit_mito = 0.32;
  cit_mito has mM;
  atp_mito = 2.8;
  atp_mito has mM;
  adp_mito = 0.8;
  adp_mito has mM;
  gtp_mito = 0.29;
  gtp_mito has mM;
  gdp_mito = 0.1;
  gdp_mito has mM;
  coa_mito = 0.055;
  coa_mito has mM;
  nadh_mito = 0.24;
  nadh_mito has mM;
  nad_mito = 0.98;
  nad_mito has mM;
  h2o_mito = 0;
  h2o_mito has mM;
  h_mito = 0;
  h_mito has mM;

  // Compartment initializations:
  extern = 10 dimensionless*Vcyto;
  extern has litre;
  cyto = Vcyto;
  cyto has litre;
  mito = 0.2 dimensionless*Vcyto;
  mito has litre;

  // Variable initializations:
  Vcyto = 1;
  Vcyto has litre;
  Vliver = 1.5;
  Vliver has litre;
  fliver = 0.583333333333334;
  fliver has dimensionless;
  bodyweight = 70;
  bodyweight has kilogram;
  sec_per_min = 60;
  sec_per_min has s_per_min;
  conversion_factor = (fliver*Vliver/Vcyto)*sec_per_min*1000 dimensionless/bodyweight;
  conversion_factor has s_per_min_kg;
  nadh_tot has mM;
  atp_tot has mM;
  utp_tot has mM;
  gtp_tot has mM;
  nadh_mito_tot has mM;
  atp_mito_tot has mM;
  gtp_mito_tot has mM;
  x_ins1 = 818.9;
  x_ins1 has pmol_per_l;
  x_ins2 = 0;
  x_ins2 has pmol_per_l;
  x_ins3 = 8.6;
  x_ins3 has mM;
  x_ins4 = 4.2;
  x_ins4 has dimensionless;
  ins has pmol_per_l;
  ins_norm has pmol_per_l;
  x_glu1 = 190;
  x_glu1 has pmol_per_l;
  x_glu2 = 37.9;
  x_glu2 has pmol_per_l;
  x_glu3 = 3.01;
  x_glu3 has mM;
  x_glu4 = 6.4;
  x_glu4 has dimensionless;
  glu has pmol_per_l;
  glu_norm has pmol_per_l;
  x_epi1 = 6090;
  x_epi1 has pmol_per_l;
  x_epi2 = 100;
  x_epi2 has pmol_per_l;
  x_epi3 = 3.1;
  x_epi3 has mM;
  x_epi4 = 8.4;
  x_epi4 has dimensionless;
  epi has pmol_per_l;
  epi_norm has pmol_per_l;
  K_val = 0.1;
  K_val has dimensionless;
  epi_f = 0.8;
  epi_f has dimensionless;
  K_ins = (x_ins1 - x_ins2)*K_val;
  K_ins has pmol_per_l;
  K_glu = (x_glu1 - x_glu2)*K_val;
  K_glu has pmol_per_l;
  K_epi = (x_epi1 - x_epi2)*K_val;
  K_epi has pmol_per_l;
  gamma has dimensionless;
  scale = 1 dimensionless/60 dimensionless;
  scale has dimensionless;
  scale_gly = scale;
  scale_gly has dimensionless;
  scale_glyglc = scale;
  scale_glyglc has dimensionless;
  GLUT2_keq = 1;
  GLUT2_keq has dimensionless;
  GLUT2_km = 42;
  GLUT2_km has mM;
  GLUT2_Vmax = 420;
  GLUT2_Vmax has mmol_per_s;
  GK_n_gkrp = 2;
  GK_n_gkrp has dimensionless;
  GK_km_glc1 = 15;
  GK_km_glc1 has mM;
  GK_km_fru6p = 0.01;
  GK_km_fru6p has mM;
  GK_b = 0.7;
  GK_b has dimensionless;
  GK_n = 1.6;
  GK_n has dimensionless;
  GK_km_glc = 7.5;
  GK_km_glc has mM;
  GK_km_atp = 0.26;
  GK_km_atp has mM;
  GK_Vmax = 25.2;
  GK_Vmax has mmol_per_s;
  GK_gc_free has dimensionless;
  G6PASE_km_glc6p = 2;
  G6PASE_km_glc6p has mM;
  G6PASE_Vmax = 18.9;
  G6PASE_Vmax has mmol_per_s;
  GPI_keq = 0.517060817492925;
  GPI_keq has dimensionless;
  GPI_km_glc6p = 0.182;
  GPI_km_glc6p has mM;
  GPI_km_fru6p = 0.071;
  GPI_km_fru6p has mM;
  GPI_Vmax = 420;
  GPI_Vmax has mmol_per_s;
  G16PI_keq = 15.7175540821514;
  G16PI_keq has dimensionless;
  G16PI_km_glc6p = 0.67;
  G16PI_km_glc6p has mM;
  G16PI_km_glc1p = 0.045;
  G16PI_km_glc1p has mM;
  G16PI_Vmax = 100;
  G16PI_Vmax has mmol_per_s;
  UPGASE_keq = 0.312237619153088;
  UPGASE_keq has dimensionless;
  UPGASE_km_utp = 0.563;
  UPGASE_km_utp has mM;
  UPGASE_km_glc1p = 0.172;
  UPGASE_km_glc1p has mM;
  UPGASE_km_udpglc = 0.049;
  UPGASE_km_udpglc has mM;
  UPGASE_km_pp = 0.166;
  UPGASE_km_pp has mM;
  UPGASE_Vmax = 80;
  UPGASE_Vmax has mmol_per_s;
  PPASE_km_pp = 0.005;
  PPASE_km_pp has mM;
  PPASE_Vmax = 2.4;
  PPASE_Vmax has mmol_per_s;
  GS_C = 500;
  GS_C has mM;
  GS_k1_max = 0.2;
  GS_k1_max has dimensionless;
  GS_k1_nat = 0.224;
  GS_k1_nat has mMmM;
  GS_k2_nat = 0.1504;
  GS_k2_nat has mM;
  GS_k1_phospho = 3.003;
  GS_k1_phospho has mMmM;
  GS_k2_phospho = 0.09029;
  GS_k2_phospho has mM;
  GS_Vmax = 13.2;
  GS_Vmax has mmol_per_s;
  GS_storage_factor has dimensionless;
  GS_k_udpglc_native has mM;
  GS_k_udpglc_phospho has mM;
  GS_native has mmol_per_s;
  GS_phospho has mmol_per_s;
  GP_keq = 0.211826505793075;
  GP_keq has per_mM;
  GP_k_glyc_native = 4.8;
  GP_k_glyc_native has mM;
  GP_k_glyc_phospho = 2.7;
  GP_k_glyc_phospho has mM;
  GP_k_glc1p_native = 120;
  GP_k_glc1p_native has mM;
  GP_k_glc1p_phospho = 2;
  GP_k_glc1p_phospho has mM;
  GP_k_p_native = 300;
  GP_k_p_native has mM;
  GP_k_p_phospho = 5;
  GP_k_p_phospho has mM;
  GP_ki_glc_phospho = 5;
  GP_ki_glc_phospho has mM;
  GP_ka_amp_native = 1;
  GP_ka_amp_native has mM;
  GP_base_amp_native = 0.03;
  GP_base_amp_native has dimensionless;
  GP_max_amp_native = 0.3;
  GP_max_amp_native has dimensionless;
  GP_Vmax = 6.8;
  GP_Vmax has mmol_per_s;
  GP_C = GS_C;
  GP_C has mM;
  GP_k1_max = GS_k1_max;
  GP_k1_max has dimensionless;
  GP_fmax has dimensionless;
  GP_vmax_native has mmol_per_s;
  GP_native has mmol_per_s;
  GP_vmax_phospho has mmol_per_s;
  GP_phospho has mmol_per_s;
  NDKGTP_keq = 1;
  NDKGTP_keq has dimensionless;
  NDKGTP_km_atp = 1.33;
  NDKGTP_km_atp has mM;
  NDKGTP_km_adp = 0.042;
  NDKGTP_km_adp has mM;
  NDKGTP_km_gtp = 0.15;
  NDKGTP_km_gtp has mM;
  NDKGTP_km_gdp = 0.031;
  NDKGTP_km_gdp has mM;
  NDKGTP_Vmax = 0;
  NDKGTP_Vmax has mmol_per_s;
  NDKUTP_keq = 1;
  NDKUTP_keq has dimensionless;
  NDKUTP_km_atp = 1.33;
  NDKUTP_km_atp has mM;
  NDKUTP_km_adp = 0.042;
  NDKUTP_km_adp has mM;
  NDKUTP_km_utp = 16;
  NDKUTP_km_utp has mM;
  NDKUTP_km_udp = 0.19;
  NDKUTP_km_udp has mM;
  NDKUTP_Vmax = 2940;
  NDKUTP_Vmax has mmol_per_s;
  AK_keq = 0.247390074904985;
  AK_keq has dimensionless;
  AK_km_atp = 0.09;
  AK_km_atp has mM;
  AK_km_amp = 0.08;
  AK_km_amp has mM;
  AK_km_adp = 0.11;
  AK_km_adp has mM;
  AK_Vmax = 0;
  AK_Vmax has mmol_per_s;
  PFK2_n_native = 1.3;
  PFK2_n_native has dimensionless;
  PFK2_n_phospho = 2.1;
  PFK2_n_phospho has dimensionless;
  PFK2_k_fru6p_native = 0.016;
  PFK2_k_fru6p_native has mM;
  PFK2_k_fru6p_phospho = 0.05;
  PFK2_k_fru6p_phospho has mM;
  PFK2_k_atp_native = 0.28;
  PFK2_k_atp_native has mM;
  PFK2_k_atp_phospho = 0.65;
  PFK2_k_atp_phospho has mM;
  PFK2_Vmax = 0.0042;
  PFK2_Vmax has mmol_per_s;
  PFK2_native has mmol_per_s;
  PFK2_phospho has mmol_per_s;
  FBP2_km_fru26bp_native = 0.01;
  FBP2_km_fru26bp_native has mM;
  FBP2_ki_fru6p_native = 0.0035;
  FBP2_ki_fru6p_native has mM;
  FBP2_km_fru26bp_phospho = 0.0005;
  FBP2_km_fru26bp_phospho has mM;
  FBP2_ki_fru6p_phospho = 0.01;
  FBP2_ki_fru6p_phospho has mM;
  FBP2_Vmax = 0.126;
  FBP2_Vmax has mmol_per_s;
  FBP2_native has mmol_per_s;
  FBP2_phospho has mmol_per_s;
  PFK1_km_atp = 0.111;
  PFK1_km_atp has mM;
  PFK1_km_fru6p = 0.077;
  PFK1_km_fru6p has mM;
  PFK1_ki_fru6p = 0.012;
  PFK1_ki_fru6p has mM;
  PFK1_ka_fru26bp = 0.001;
  PFK1_ka_fru26bp has mM;
  PFK1_Vmax = 7.182;
  PFK1_Vmax has mmol_per_s;
  FBP1_ki_fru26bp = 0.001;
  FBP1_ki_fru26bp has mM;
  FBP1_km_fru16bp = 0.0013;
  FBP1_km_fru16bp has mM;
  FBP1_Vmax = 4.326;
  FBP1_Vmax has mmol_per_s;
  ALD_keq = 9.76298897362969e-05;
  ALD_keq has mM;
  ALD_km_fru16bp = 0.0071;
  ALD_km_fru16bp has mM;
  ALD_km_dhap = 0.0364;
  ALD_km_dhap has mM;
  ALD_km_grap = 0.0071;
  ALD_km_grap has mM;
  ALD_ki1_grap = 0.0572;
  ALD_ki1_grap has mM;
  ALD_ki2_grap = 0.176;
  ALD_ki2_grap has mM;
  ALD_Vmax = 420;
  ALD_Vmax has mmol_per_s;
  TPI_keq = 0.054476985386756;
  TPI_keq has dimensionless;
  TPI_km_dhap = 0.59;
  TPI_km_dhap has mM;
  TPI_km_grap = 0.42;
  TPI_km_grap has mM;
  TPI_Vmax = 420;
  TPI_Vmax has mmol_per_s;
  GAPDH_keq = 0.086779866194594;
  GAPDH_keq has per_mM;
  GAPDH_k_nad = 0.05;
  GAPDH_k_nad has mM;
  GAPDH_k_grap = 0.005;
  GAPDH_k_grap has mM;
  GAPDH_k_p = 3.9;
  GAPDH_k_p has mM;
  GAPDH_k_nadh = 0.0083;
  GAPDH_k_nadh has mM;
  GAPDH_k_bpg13 = 0.0035;
  GAPDH_k_bpg13 has mM;
  GAPDH_Vmax = 420;
  GAPDH_Vmax has mmol_per_s;
  PGK_keq = 6.95864405248854;
  PGK_keq has dimensionless;
  PGK_k_adp = 0.35;
  PGK_k_adp has mM;
  PGK_k_atp = 0.48;
  PGK_k_atp has mM;
  PGK_k_bpg13 = 0.002;
  PGK_k_bpg13 has mM;
  PGK_k_pg3 = 1.2;
  PGK_k_pg3 has mM;
  PGK_Vmax = 420;
  PGK_Vmax has mmol_per_s;
  PGM_keq = 0.181375378837397;
  PGM_keq has dimensionless;
  PGM_k_pg3 = 5;
  PGM_k_pg3 has mM;
  PGM_k_pg2 = 1;
  PGM_k_pg2 has mM;
  PGM_Vmax = 420;
  PGM_Vmax has mmol_per_s;
  EN_keq = 0.054476985386756;
  EN_keq has dimensionless;
  EN_k_pep = 1;
  EN_k_pep has mM;
  EN_k_pg2 = 1;
  EN_k_pg2 has mM;
  EN_Vmax = 35.994;
  EN_Vmax has mmol_per_s;
  PK_n = 3.5;
  PK_n has dimensionless;
  PK_n_p = 3.5;
  PK_n_p has dimensionless;
  PK_n_fbp = 1.8;
  PK_n_fbp has dimensionless;
  PK_n_fbp_p = 1.8;
  PK_n_fbp_p has dimensionless;
  PK_alpha = 1;
  PK_alpha has dimensionless;
  PK_alpha_p = 1.1;
  PK_alpha_p has dimensionless;
  PK_alpha_end = 1;
  PK_alpha_end has dimensionless;
  PK_k_fbp = 0.00016;
  PK_k_fbp has mM;
  PK_k_fbp_p = 0.00035;
  PK_k_fbp_p has mM;
  PK_k_pep = 0.58;
  PK_k_pep has mM;
  PK_k_pep_p = 1.1;
  PK_k_pep_p has mM;
  PK_k_pep_end = 0.08;
  PK_k_pep_end has mM;
  PK_k_adp = 2.3;
  PK_k_adp has mM;
  PK_base_act = 0.08;
  PK_base_act has mM;
  PK_base_act_p = 0.04;
  PK_base_act_p has mM;
  PK_Vmax = 46.2;
  PK_Vmax has mmol_per_s;
  PK_f has dimensionless;
  PK_f_p has dimensionless;
  PK_alpha_inp has dimensionless;
  PK_alpha_p_inp has dimensionless;
  PK_pep_inp has mM;
  PK_pep_p_inp has mM;
  PK_native has mmol_per_s;
  PK_phospho has mmol_per_s;
  PEPCK_keq = 336.956521586429;
  PEPCK_keq has mM;
  PEPCK_k_pep = 0.237;
  PEPCK_k_pep has mM;
  PEPCK_k_gdp = 0.0921;
  PEPCK_k_gdp has mM;
  PEPCK_k_co2 = 25.5;
  PEPCK_k_co2 has mM;
  PEPCK_k_oaa = 0.0055;
  PEPCK_k_oaa has mM;
  PEPCK_k_gtp = 0.0222;
  PEPCK_k_gtp has mM;
  PEPCK_Vmax = 0;
  PEPCK_Vmax has mmol_per_s;
  PEPCKM_Vmax = 546;
  PEPCKM_Vmax has mmol_per_s;
  PC_k_atp = 0.22;
  PC_k_atp has mM;
  PC_k_pyr = 0.22;
  PC_k_pyr has mM;
  PC_k_co2 = 3.2;
  PC_k_co2 has mM;
  PC_k_acoa = 0.015;
  PC_k_acoa has mM;
  PC_n = 2.5;
  PC_n has dimensionless;
  PC_Vmax = 168;
  PC_Vmax has mmol_per_s;
  LDH_keq = 0.000278321076004752;
  LDH_keq has dimensionless;
  LDH_k_pyr = 0.495;
  LDH_k_pyr has mM;
  LDH_k_lac = 31.98;
  LDH_k_lac has mM;
  LDH_k_nad = 0.984;
  LDH_k_nad has mM;
  LDH_k_nadh = 0.027;
  LDH_k_nadh has mM;
  LDH_Vmax = 12.6;
  LDH_Vmax has mmol_per_s;
  LACT_keq = 1;
  LACT_keq has dimensionless;
  LACT_k_lac = 0.8;
  LACT_k_lac has mM;
  LACT_Vmax = 5.418;
  LACT_Vmax has mmol_per_s;
  PYRTM_keq = 1;
  PYRTM_keq has dimensionless;
  PYRTM_k_pyr = 0.1;
  PYRTM_k_pyr has mM;
  PYRTM_Vmax = 42;
  PYRTM_Vmax has mmol_per_s;
  PEPTM_keq = 1;
  PEPTM_keq has dimensionless;
  PEPTM_k_pep = 0.1;
  PEPTM_k_pep has mM;
  PEPTM_Vmax = 33.6;
  PEPTM_Vmax has mmol_per_s;
  PDH_k_pyr = 0.025;
  PDH_k_pyr has mM;
  PDH_k_coa = 0.013;
  PDH_k_coa has mM;
  PDH_k_nad = 0.05;
  PDH_k_nad has mM;
  PDH_ki_acoa = 0.035;
  PDH_ki_acoa has mM;
  PDH_ki_nadh = 0.036;
  PDH_ki_nadh has mM;
  PDH_alpha_nat = 5;
  PDH_alpha_nat has dimensionless;
  PDH_alpha_p = 1;
  PDH_alpha_p has dimensionless;
  PDH_Vmax = 13.44;
  PDH_Vmax has mmol_per_s;
  PDH_base has mmol_per_s;
  PDH_nat has mmol_per_s;
  PDH_p has mmol_per_s;
  CS_keq = 266599.030842759;
  CS_keq has dimensionless;
  CS_k_oaa = 0.002;
  CS_k_oaa has mM;
  CS_k_acoa = 0.016;
  CS_k_acoa has mM;
  CS_k_cit = 0.42;
  CS_k_cit has mM;
  CS_k_coa = 0.07;
  CS_k_coa has mM;
  CS_Vmax = 4.2;
  CS_Vmax has mmol_per_s;
  NDKGTPM_keq = 1;
  NDKGTPM_keq has dimensionless;
  NDKGTPM_k_atp = 1.33;
  NDKGTPM_k_atp has mM;
  NDKGTPM_k_adp = 0.042;
  NDKGTPM_k_adp has mM;
  NDKGTPM_k_gtp = 0.15;
  NDKGTPM_k_gtp has mM;
  NDKGTPM_k_gdp = 0.031;
  NDKGTPM_k_gdp has mM;
  NDKGTPM_Vmax = 420;
  NDKGTPM_Vmax has mmol_per_s;
  OAAFLX_Vmax = 0;
  OAAFLX_Vmax has mmol_per_s;
  ACOAFLX_Vmax = 0;
  ACOAFLX_Vmax has mmol_per_s;
  CITFLX_Vmax = 0;
  CITFLX_Vmax has mmol_per_s;
  HGP has mumol_per_min_kg;
  GNG has mumol_per_min_kg;
  GLY has mumol_per_min_kg;

  // Other declarations:
  var nadh_tot, atp_tot, utp_tot, gtp_tot, nadh_mito_tot, atp_mito_tot, gtp_mito_tot;
  var ins, ins_norm, glu, glu_norm, epi, epi_norm, gamma, GK_gc_free, GS_storage_factor;
  var GS_k_udpglc_native, GS_k_udpglc_phospho, GS_native, GS_phospho, GP_fmax;
  var GP_vmax_native, GP_native, GP_vmax_phospho, GP_phospho, PFK2_native;
  var PFK2_phospho, FBP2_native, FBP2_phospho, PK_f, PK_f_p, PK_alpha_inp;
  var PK_alpha_p_inp, PK_pep_inp, PK_pep_p_inp, PK_native, PK_phospho, PDH_base;
  var PDH_nat, PDH_p, HGP, GNG, GLY;
  const extern, Vcyto, cyto, mito, Vliver, fliver, bodyweight, sec_per_min;
  const conversion_factor, x_ins1, x_ins2, x_ins3, x_ins4, x_glu1, x_glu2;
  const x_glu3, x_glu4, x_epi1, x_epi2, x_epi3, x_epi4, K_val, epi_f, K_ins;
  const K_glu, K_epi, scale, scale_gly, scale_glyglc, GLUT2_keq, GLUT2_km;
  const GLUT2_Vmax, GK_n_gkrp, GK_km_glc1, GK_km_fru6p, GK_b, GK_n, GK_km_glc;
  const GK_km_atp, GK_Vmax, G6PASE_km_glc6p, G6PASE_Vmax, GPI_keq, GPI_km_glc6p;
  const GPI_km_fru6p, GPI_Vmax, G16PI_keq, G16PI_km_glc6p, G16PI_km_glc1p;
  const G16PI_Vmax, UPGASE_keq, UPGASE_km_utp, UPGASE_km_glc1p, UPGASE_km_udpglc;
  const UPGASE_km_pp, UPGASE_Vmax, PPASE_km_pp, PPASE_Vmax, GS_C, GS_k1_max;
  const GS_k1_nat, GS_k2_nat, GS_k1_phospho, GS_k2_phospho, GS_Vmax, GP_keq;
  const GP_k_glyc_native, GP_k_glyc_phospho, GP_k_glc1p_native, GP_k_glc1p_phospho;
  const GP_k_p_native, GP_k_p_phospho, GP_ki_glc_phospho, GP_ka_amp_native;
  const GP_base_amp_native, GP_max_amp_native, GP_Vmax, GP_C, GP_k1_max, NDKGTP_keq;
  const NDKGTP_km_atp, NDKGTP_km_adp, NDKGTP_km_gtp, NDKGTP_km_gdp, NDKGTP_Vmax;
  const NDKUTP_keq, NDKUTP_km_atp, NDKUTP_km_adp, NDKUTP_km_utp, NDKUTP_km_udp;
  const NDKUTP_Vmax, AK_keq, AK_km_atp, AK_km_amp, AK_km_adp, AK_Vmax, PFK2_n_native;
  const PFK2_n_phospho, PFK2_k_fru6p_native, PFK2_k_fru6p_phospho, PFK2_k_atp_native;
  const PFK2_k_atp_phospho, PFK2_Vmax, FBP2_km_fru26bp_native, FBP2_ki_fru6p_native;
  const FBP2_km_fru26bp_phospho, FBP2_ki_fru6p_phospho, FBP2_Vmax, PFK1_km_atp;
  const PFK1_km_fru6p, PFK1_ki_fru6p, PFK1_ka_fru26bp, PFK1_Vmax, FBP1_ki_fru26bp;
  const FBP1_km_fru16bp, FBP1_Vmax, ALD_keq, ALD_km_fru16bp, ALD_km_dhap;
  const ALD_km_grap, ALD_ki1_grap, ALD_ki2_grap, ALD_Vmax, TPI_keq, TPI_km_dhap;
  const TPI_km_grap, TPI_Vmax, GAPDH_keq, GAPDH_k_nad, GAPDH_k_grap, GAPDH_k_p;
  const GAPDH_k_nadh, GAPDH_k_bpg13, GAPDH_Vmax, PGK_keq, PGK_k_adp, PGK_k_atp;
  const PGK_k_bpg13, PGK_k_pg3, PGK_Vmax, PGM_keq, PGM_k_pg3, PGM_k_pg2, PGM_Vmax;
  const EN_keq, EN_k_pep, EN_k_pg2, EN_Vmax, PK_n, PK_n_p, PK_n_fbp, PK_n_fbp_p;
  const PK_alpha, PK_alpha_p, PK_alpha_end, PK_k_fbp, PK_k_fbp_p, PK_k_pep;
  const PK_k_pep_p, PK_k_pep_end, PK_k_adp, PK_base_act, PK_base_act_p, PK_Vmax;
  const PEPCK_keq, PEPCK_k_pep, PEPCK_k_gdp, PEPCK_k_co2, PEPCK_k_oaa, PEPCK_k_gtp;
  const PEPCK_Vmax, PEPCKM_Vmax, PC_k_atp, PC_k_pyr, PC_k_co2, PC_k_acoa;
  const PC_n, PC_Vmax, LDH_keq, LDH_k_pyr, LDH_k_lac, LDH_k_nad, LDH_k_nadh;
  const LDH_Vmax, LACT_keq, LACT_k_lac, LACT_Vmax, PYRTM_keq, PYRTM_k_pyr;
  const PYRTM_Vmax, PEPTM_keq, PEPTM_k_pep, PEPTM_Vmax, PDH_k_pyr, PDH_k_coa;
  const PDH_k_nad, PDH_ki_acoa, PDH_ki_nadh, PDH_alpha_nat, PDH_alpha_p, PDH_Vmax;
  const CS_keq, CS_k_oaa, CS_k_acoa, CS_k_cit, CS_k_coa, CS_Vmax, NDKGTPM_keq;
  const NDKGTPM_k_atp, NDKGTPM_k_adp, NDKGTPM_k_gtp, NDKGTPM_k_gdp, NDKGTPM_Vmax;
  const OAAFLX_Vmax, ACOAFLX_Vmax, CITFLX_Vmax;
end
''')

roadrunner.Config.setValue(roadrunner.Config.ROADRUNNER_JACOBIAN_MODE, roadrunner.Config.ROADRUNNER_JACOBIAN_MODE_AMOUNTS)

print r.getFloatingSpeciesIds()
m = r.simulate (0, 2000, 100)

ids = r.getFloatingSpeciesIds()
concs = r.getFloatingSpeciesConcentrations()
print len(ids)

count = 0
for i in ids:
    print i, ": ", concs[count]
    count = count + 1

print r.steadyState()
0u812 commented 6 years ago

udp is one of the problematic variables, and is used in GS and NDKUTP reactions.

kirichoi commented 6 years ago

Update to this issue: I have tried to debug using C API but the values are correct here. And C API supposedly uses C++ functions. Any thoughts?

luciansmith commented 6 years ago

If the C API can calculate the steady state directly, that probably means that the python API is setting some other parameter value that interferes. Can you try re-running the C API with the exact same structure as Herbert's program? (Alternately, you might try the python API directly to make sure the problem still exists--it might have been fixed by some other change.)

hsauro commented 6 years ago

Kiri, have you tried calling the steady state function from C?

Herbert

On Tue, Nov 14, 2017 at 6:00 PM Lucian Smith notifications@github.com wrote:

If the C API can calculate the steady state directly, that probably means that the python API is setting some other parameter value that interferes. Can you try re-running the C API with the exact same structure as Herbert's program? (Alternately, you might try the python API directly to make sure the problem still exists--it might have been fixed by some other change.)

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/sys-bio/roadrunner/issues/420#issuecomment-344461762, or mute the thread https://github.com/notifications/unsubscribe-auth/ABAZDoG3gYtEX6VAX21wdZBKYemad6puks5s2kWNgaJpZM4Pq4_p .

kirichoi commented 6 years ago

Hmm, I still cannot find steady state solution in C.

hsauro commented 6 years ago

That's not the problem then. You'll have to create a stable version that you can debug through Python. I still don't quite understand why the Jacobian is ordered differently though when you called it from C.

Herbert

On Tue, Nov 14, 2017 at 7:41 PM Kiri Choi notifications@github.com wrote:

Hmm, I still cannot find steady state solution in C.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/sys-bio/roadrunner/issues/420#issuecomment-344476773, or mute the thread https://github.com/notifications/unsubscribe-auth/ABAZDk9CEDN0oxfIEzgPET878yKD1QwNks5s2l1KgaJpZM4Pq4_p .

kirichoi commented 6 years ago

Well, this is interesting. The roadrunner I compiled to test against Python debug library seems to be fine, at least in terms of reduced jacobian output and unscaled elasticity. The ordering is identical to C API. Any thoughts?

It could be something to do with dependencies or the way roadrunner was built.

luciansmith commented 6 years ago

I would recompile in non-debug mode just to make sure that doesn't make a difference. If that also works, you can either just close this issue and say 'it works now', or you can try to track it down. If you want to track it down:

hsauro commented 6 years ago

Are you saying the ordering from Python is now different from copasi?

Herbert

On Wed, Nov 15, 2017 at 2:19 PM Kiri Choi notifications@github.com wrote:

Well, this is interesting. The roadrunner I compiled to test against Python debug library seems to be fine, at least in terms of reduced jacobian output and unscaled elasticity. The ordering is identical to C API. Any thoughts?

It could be something to do with dependencies or the way roadrunner was built.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/sys-bio/roadrunner/issues/420#issuecomment-344747567, or mute the thread https://github.com/notifications/unsubscribe-auth/ABAZDsjXK-QcfE-JT8atii91EkCGbcZtks5s22M0gaJpZM4Pq4_p .

luciansmith commented 6 years ago

My reading was that it is now the same as copasi; i.e. he is getting different (correct) results in his setup than you did when you initially reported the bug.

kirichoi commented 6 years ago

Sorry for confusion. The ordering IS different from copasi.

hsauro commented 6 years ago

That concerns me. It should be the same, although it's possible for one or two row/column differences.

H

On Wed, Nov 15, 2017 at 3:32 PM Kiri Choi notifications@github.com wrote:

Sorry for confusion. The ordering IS different from copasi.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/sys-bio/roadrunner/issues/420#issuecomment-344764412, or mute the thread https://github.com/notifications/unsubscribe-auth/ABAZDlYjaT4g_sgneq0HbMW4Y_VM8BYpks5s23RBgaJpZM4Pq4_p .

luciansmith commented 6 years ago

Yeah, at this point you've basically confirmed the bug that Herbert saw: the Jacobian is different, which (presumably) makes the steady state impossible to find. (There is definitely a way in the roadrunner C API to call the 'steadystate' function, but I don't remember what it might be--it still might be nice to confirm this bit.) At this point, it sounds like you can just use the C API to debug the problem, figuring out why the Jacobian is wrong.

hsauro commented 6 years ago

If you have a stable debug version using the Python interface we can try to debug it tomorrow afternoon.

Or if you can get the C API to also show the bug then we can try that tomorrow afternoon.

Herbert

On Wed, Nov 15, 2017 at 3:46 PM Lucian Smith notifications@github.com wrote:

Yeah, at this point you've basically confirmed the bug that Herbert saw: the Jacobian is different, which (presumably) makes the steady state impossible to find. (There is definitely a way in the roadrunner C API to call the 'steadystate' function, but I don't remember what it might be--it still might be nice to confirm this bit.) At this point, it sounds like you can just use the C API to debug the problem, figuring out why the Jacobian is wrong.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/sys-bio/roadrunner/issues/420#issuecomment-344767153, or mute the thread https://github.com/notifications/unsubscribe-auth/ABAZDnWrVem1sglQcMn7Zdm04ObZ8v3Nks5s23eEgaJpZM4Pq4_p .

kirichoi commented 6 years ago

Okay, so few things I found so far: 1) Elasticities are practically identical, except few columns are swapped. 2) Still cannot get steady-state. 2) Full Jacobian is identical except the labeling is wrong. I wonder if this is a simple issue with the way we add labels. When conservedMoietyAnalysis = True, I believe the ordering of the floating species ids changes. I wonder if this has something to do with that because the ordering of unchanged floating species ids matches perfectly. 3) Reduced Jacobian has issues with setting conserved moiety analysis as True. I think this might be something to be with the way we handle dependent species and how the algorithm keeps track of it? Keep in mind we use two different algorithms when calculating reduced and full jacobians.

This is 100% reproducible in C API. So we can debug on C API only.

I will attach few csv files for unscaled elasticity (ue.csv), reduced jacobian (rj.csv), and full jacobian (fj.csv) I got from recompiled roadrunner.

rrcsv.zip

hsauro commented 6 years ago

Let's look at this tomorrow afternoon.

Herbert

On Wed, Nov 15, 2017 at 4:37 PM Kiri Choi notifications@github.com wrote:

Okay, so few things I found so far:

  1. Elasticities are practically identical, except few columns are swapped.
  2. Still cannot get steady-state.
  3. Full Jacobian is identical except the labeling is wrong. I wonder if this is a simple issue with the way we add labels. When conservedMoietyAnalysis = True, I believe the ordering of the floating species ids changes. I wonder if this has something to do with that because the ordering of unchanged floating species ids matches perfectly.
  4. Reduced Jacobian has issues with setting conserved moiety analysis as True. I think this might be something to be with the way we handle dependent species and how the algorithm keeps track of it? Keep in mind we use two different algorithms when calculating reduced and full jacobians.

This is 100% reproducible in C API. So we can debug on C API only.

I will attach few csv files for unscaled elasticity (ue.csv), reduced jacobian (rj.csv), and full jacobian (fj.csv) I got from recompiled roadrunner.

rrcsv.zip https://github.com/sys-bio/roadrunner/files/1476853/rrcsv.zip

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/sys-bio/roadrunner/issues/420#issuecomment-344776552, or mute the thread https://github.com/notifications/unsubscribe-auth/ABAZDuECq1ln9foAn3s3Wq-8Kf802bN6ks5s24OJgaJpZM4Pq4_p .

kirichoi commented 6 years ago

Okay, I have looked into this further and found something quite disturbing. So one thing Herbert suggested was that maybe the reduced jacobian calculation (which uses rates of change) was having issue with conserved totals (unintentionally changing other values). So I tried to revert the calculation back to use unscaled elasticity (Nr x uE x Link I believe?) and while doing so, found out that our list of dependent species seems to differ from COPASI's (and thus our Nr as well).

RoadRunner: ['gtp', 'utp', 'gdp_mito'] COPASI: ['gtp', 'udpglc', 'gdp_mito']

This was confirmed in C API. This is an issue with libStruct. And I think this is causing all the other issues with Jacobian and etc.

kirichoi commented 6 years ago

On libStructural.cpp, the code simply puts the list of dependent species ids based on the difference in index (floating species[number of independent floating species + i]) I need to double check whether the ordering of floating species are based on inpdependent and dependent floating species.

hsauro commented 6 years ago

Let's look at this more closely tomorrow. One thing I can say, I believe libstruct is also used by the c# roadrunner and I believe I confirmed that it can simulate matthias's model. I'll check again tonight.

If we don't fix it tomorrow's can we can create a debug setup for me so that I try figuring it out at the weekend. I know the background theory better.

Herbert

On Thu, Nov 16, 2017 at 6:15 PM Kiri Choi notifications@github.com wrote:

On libStructural.cpp, the code simply puts the list of dependent species ids based on the difference in index (floating species[number of independent floating species + i]) I need to double check whether the ordering of floating species are based on inpdependent and dependent floating species.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/sys-bio/roadrunner/issues/420#issuecomment-345124941, or mute the thread https://github.com/notifications/unsubscribe-auth/ABAZDiczeG90RK77wZPmJk0OsVaO8zKTks5s3OwXgaJpZM4Pq4_p .

kirichoi commented 6 years ago

So apparently the ordering of species is very important for many different calculations. We use QR decomposition for conservation analysis and here's where the code swaps around the species ordering:

// reorder species //Todo: This causes c++ generated code be different from c#??
 for (unsigned int i = 0; i < P->numRows(); i++)
        {
            for (unsigned int j = 0; j < P->numCols(); j++)
            {
                double aNUm = (*P)(i,j);
                if ((*P)(i,j) == 1)
                {
                    spVec[j]=i;        //here it is!!!
                    break;
                }
            }
        }

P is permutation matrix. That todo comment is very interesting btw... I would really appreciate any comments.

hsauro commented 6 years ago

Actually I was wrong, it seems the C# roadrunenr can't solve the steady state either.

Herbert

On Thu, Nov 16, 2017 at 7:59 PM, Herbert Sauro uw.hsauro@gmail.com wrote:

Let's look at this more closely tomorrow. One thing I can say, I believe libstruct is also used by the c# roadrunner and I believe I confirmed that it can simulate matthias's model. I'll check again tonight.

If we don't fix it tomorrow's can we can create a debug setup for me so that I try figuring it out at the weekend. I know the background theory better.

Herbert

On Thu, Nov 16, 2017 at 6:15 PM Kiri Choi notifications@github.com wrote:

On libStructural.cpp, the code simply puts the list of dependent species ids based on the difference in index (floating species[number of independent floating species + i]) I need to double check whether the ordering of floating species are based on inpdependent and dependent floating species.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/sys-bio/roadrunner/issues/420#issuecomment-345124941, or mute the thread https://github.com/notifications/unsubscribe-auth/ABAZDiczeG90RK77wZPmJk0OsVaO8zKTks5s3OwXgaJpZM4Pq4_p .

hsauro commented 6 years ago

That doesn't look good, looks like its a left over from something. Yosef's tests on libStruct however pass. I don't think he has tried this model though.

Herbert

On Thu, Nov 16, 2017 at 8:02 PM, Kiri Choi notifications@github.com wrote:

So apparently the ordering of species is very important for many different calculations. We use QR decomposition for conservation analysis and here's where the code swaps around the species ordering:

// reorder species //Todo: This causes c++ generated code be different from c#?? for (unsigned int i = 0; i < P->numRows(); i++) { for (unsigned int j = 0; j < P->numCols(); j++) { double aNUm = (P)(i,j); if ((P)(i,j) == 1) { spVec[j]=i; //here it is!!! break; } } }

P is permutation matrix. That todo comment is very interesting btw... I would really appreciate any comments.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/sys-bio/roadrunner/issues/420#issuecomment-345139604, or mute the thread https://github.com/notifications/unsubscribe-auth/ABAZDipEwaY207xgjCn-xQfc4Yl_H-aFks5s3QVpgaJpZM4Pq4_p .

hsauro commented 6 years ago

I just ran the model on libstruct that Yosef has and it passed the 6 internal tests.

Herbert

On Thu, Nov 16, 2017 at 8:11 PM, Herbert Sauro uw.hsauro@gmail.com wrote:

That doesn't look good, looks like its a left over from something. Yosef's tests on libStruct however pass. I don't think he has tried this model though.

Herbert

On Thu, Nov 16, 2017 at 8:02 PM, Kiri Choi notifications@github.com wrote:

So apparently the ordering of species is very important for many different calculations. We use QR decomposition for conservation analysis and here's where the code swaps around the species ordering:

// reorder species //Todo: This causes c++ generated code be different from c#?? for (unsigned int i = 0; i < P->numRows(); i++) { for (unsigned int j = 0; j < P->numCols(); j++) { double aNUm = (P)(i,j); if ((P)(i,j) == 1) { spVec[j]=i; //here it is!!! break; } } }

P is permutation matrix. That todo comment is very interesting btw... I would really appreciate any comments.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/sys-bio/roadrunner/issues/420#issuecomment-345139604, or mute the thread https://github.com/notifications/unsubscribe-auth/ABAZDipEwaY207xgjCn-xQfc4Yl_H-aFks5s3QVpgaJpZM4Pq4_p .

kirichoi commented 6 years ago

Unless a test is run against the original model, i wouldn't trust the results. Could anyone confirm that utp is indeed not a dependent species in this case?

hsauro commented 6 years ago

In libStruct, utp is a DEPENDENT species

ls.getDependentSpecies() Out[8]: ('gtp', 'utp', 'gdp_mito')

ls.getIndependentSpecies() Out[9]: ('fru6p', 'pep', 'oaa_mito', 'glc6p', 'grap', 'pyr', 'glc1p', 'pyr_mito', 'dhap', 'pg3', 'udp', 'glc', 'gdp', 'gtp_mito', 'pp', 'lac', 'fru16bp', 'pg2', 'glyglc', 'pep_mito', 'bpg13', 'fru26bp', 'oaa', 'udpglc')

On Thu, Nov 16, 2017 at 9:08 PM, Kiri Choi notifications@github.com wrote:

Unless a test is run against the original model, i wouldn't trust the results. Could anyone confirm that utp is indeed not a dependent species in this case?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/sys-bio/roadrunner/issues/420#issuecomment-345146837, or mute the thread https://github.com/notifications/unsubscribe-auth/ABAZDsAYbM7y8bWbmE0xBwy7W5EYA65Pks5s3RTBgaJpZM4Pq4_p .

kirichoi commented 6 years ago

Well, according to COPASI, it is not suppose to.

kirichoi commented 6 years ago

Could this be something to do with the model itself? It might be the case where udp, utp, and udpglc forms a conserved cycle.

hsauro commented 6 years ago

That depends on the algorithm that is used to reduce the stoichiometry matrix.

Ask Ysosef to run the model thorugh PySCeS to see what it gets. I can try running it through Jarnac to see what it gets.

Herbert

On Thu, Nov 16, 2017 at 9:11 PM, Kiri Choi notifications@github.com wrote:

Well, according to COPASI, it is not suppose to.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/sys-bio/roadrunner/issues/420#issuecomment-345147136, or mute the thread https://github.com/notifications/unsubscribe-auth/ABAZDuIU9aK-mcBN6iJCuMM-YWQjNXBEks5s3RVjgaJpZM4Pq4_p .

hsauro commented 6 years ago

There is a conserved cycle involving udp and utp, give the three reactions:

  UPGASE: utp + glc1p -> udpglc + pp;  v;
  GS: udpglc -> udp + glyglc;  v;
  NDKUTP: atp + udp -> adp + utp;   v;

there is a cyle from utp to udpglc to udp and back to utp.

There are three conserved cycles in all.

Herbert

On Thu, Nov 16, 2017 at 9:16 PM, Kiri Choi notifications@github.com wrote:

Could this be something to do with the model itself? It might be the case where udp, utp, and udpglc forms a conserved cycle.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/sys-bio/roadrunner/issues/420#issuecomment-345147744, or mute the thread https://github.com/notifications/unsubscribe-auth/ABAZDh8iszE745xWYTpguMwbAnT6XEl3ks5s3RavgaJpZM4Pq4_p .

hsauro commented 6 years ago

Results from Jarnac, for the conserved cycle:

udp + utp + glyudp

Jarnac sets udp and utp as independent variables and glyudp as the dependent variable, that because in the stoichiometry matrix udp and adp are in the first rows. If I were to change the order of the declared species, eg to

udpgly, then utp then udp, now Jarnac says:

udpglc and utp are the independent variables and idp the dependent variable. Also Jarnac uses a simpler algorithm and we may get different reductions anyway.

Pysces is the more similar one to compare to.

The order of the species depends on 1) The particular algorithm used and 2) the order of species in the SBML file.

In the SBML file I used utp was declared first followed by udp, that explains the Jarnac result.

If I load this sbml file into copasi it gives a completely different order from Jarnac but udp and utp are still both independent variables. However the order from copasi is more similar to what libstruct gives except for the last entry in the reduced stoichiometry matrix which has utp in copasi and udpglc in libstruct.

I think both copasi and libstruct and doing the correct calculation.

Herbert

On Thu, Nov 16, 2017 at 9:31 PM, Herbert Sauro uw.hsauro@gmail.com wrote:

There is a conserved cycle involving udp and utp, give the three reactions:

  UPGASE: utp + glc1p -> udpglc + pp;  v;
  GS: udpglc -> udp + glyglc;  v;
  NDKUTP: atp + udp -> adp + utp;   v;

there is a cyle from utp to udpglc to udp and back to utp.

There are three conserved cycles in all.

Herbert

On Thu, Nov 16, 2017 at 9:16 PM, Kiri Choi notifications@github.com wrote:

Could this be something to do with the model itself? It might be the case where udp, utp, and udpglc forms a conserved cycle.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/sys-bio/roadrunner/issues/420#issuecomment-345147744, or mute the thread https://github.com/notifications/unsubscribe-auth/ABAZDh8iszE745xWYTpguMwbAnT6XEl3ks5s3RavgaJpZM4Pq4_p .

kirichoi commented 6 years ago

So the reason why I am focused on ordering is because it determines dependent and independent species. (e.g. the last 3 species of reordered list will be treated as dependent species) The problem with this approach is that it affects everything from reduced stoichiometry matrix to L/L0 matrix. (Full stoichiometry matrix is identical with the output of COPASI but with swapped rows. Reduced stoichiometry however, is not the same with COPASI because the code takes what it assumes to be independent species.) Everything is based on for loop that goes through either list of independent species or dependent species. Sometimes this affects the labels only, but most of the time it affects the contents.

hsauro commented 6 years ago

It shouldn't matter that the labels in copasi and libroadrunner are different, what's important is that any matrix arithmetic is done using compatible row and column orderings.

Herbert

On Thu, Nov 16, 2017 at 11:00 PM Kiri Choi notifications@github.com wrote:

So the reason why I am focused on ordering is because it determines dependent and independent species. (e.g. the last 3 species of reordered list will be treated as dependent species) The problem with this approach is that it affects everything from reduced stoichiometry matrix to L/L0 matrix. (Full stoichiometry matrix is identical with the output of COPASI but with swapped rows. Reduced stoichiometry however, is not the same with COPASI because the code takes what it assumes to be independent species.) Everything is based on for loop that goes through either list of independent species or dependent species. Sometimes this affects the labels only, but most of the time it affects the contents.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/sys-bio/roadrunner/issues/420#issuecomment-345161912, or mute the thread https://github.com/notifications/unsubscribe-auth/ABAZDoGxl6_hdiRwtPKo-OeBVOdIiafDks5s3S7tgaJpZM4Pq4_p .

kirichoi commented 6 years ago

That's what I am saying. I think the ordering affects more than labeling right now. When things like reduced stoichiometry matrix is affect by this, everything else that depends on it will be incorrect. Right now, I don't think there are any issues with the jacobian calculation. It's doing what it's supposed to do, but with incorrect matrices to began with. And that is caused by ordering of species and how it determines independent and dependent species.

matthiaskoenig commented 6 years ago

@all I can't help with the libroadrunner debugging, but if you need any information about the model/constructs I used please ask M

On Fri, Nov 17, 2017 at 8:10 AM, Kiri Choi notifications@github.com wrote:

That's what I am saying. I think the ordering affects more than labeling right now. When things like reduced stoichiometry matrix is affect by this, everything else that depends on it will be incorrect. Right now, I don't think there are any issues with the jacobian calculation. It's doing what it's supposed to do, but with incorrect matrices to began with. And that is caused by ordering of species and how it determines independent and dependent species.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/sys-bio/roadrunner/issues/420#issuecomment-345163730, or mute the thread https://github.com/notifications/unsubscribe-auth/AA29ukTN1AtIiR5Zmxsdt8vOPtZOFtPcks5s3TGDgaJpZM4Pq4_p .

-- Dr. Matthias König Junior Group Leader LiSyM - Systems Medicine of the Liver Humboldt-University Berlin, Institute of Biology, Institute for Theoretical Biology https://www.livermetabolism.com konigmatt@googlemail.com Tel: +49 30 20938450 Tel: +49 176 81168480

hsauro commented 6 years ago

This morning I'm going to try reducing the model size, maybe it will be easier to figure out what is going on.

Herbert

On Thu, Nov 16, 2017 at 11:11 PM Kiri Choi notifications@github.com wrote:

That's what I am saying. I think the ordering affects more than labeling right now. When things like reduced stoichiometry matrix is affect by this, everything else that depends on it will be incorrect. Right now, I don't think there are any issues with the jacobian calculation. It's doing what it's supposed to do, but with incorrect matrices to began with. And that is caused by ordering of species and how it determines independent and dependent species.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/sys-bio/roadrunner/issues/420#issuecomment-345163730, or mute the thread https://github.com/notifications/unsubscribe-auth/ABAZDrjvtWn0KJu792I6uQkns4N8uYpsks5s3TGDgaJpZM4Pq4_p .

kirichoi commented 6 years ago

Update on this issue:

As suggested, this has everything to do with the ordering and dependent species. The problematic part is the link matrix. Reduced stoichiometry matrix is simply generated by putting 1 in the diagonal for independent species. This is problematic because in Nr and uEE, the values are swapped when the ordering changes. In Link matrix, that is not the case; ordering is swapped but the values are not. Thus we see that one incorrect value for udp. Also, the last row/column differ because values for udpglc, instead of utp, are used for calculation.

hsauro commented 6 years ago

ok, gettting closer.

Herbert

On Wed, Nov 22, 2017 at 9:57 PM, Kiri Choi notifications@github.com wrote:

Update on this issue:

As suggested, this has everything to do with the ordering and dependent species. The problematic part is the link matrix. Reduced stoichiometry matrix is simply generated by putting 1 in diagonal for independent species and others for independent species. Now, that means that even though the ordering differs in link matrix, the values are the same as COPASI. This is problematic because in Nr and uEE, the values are swapped along with the ordering. In Link matrix, that is not the case; ordering is swapped but the values are not.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/sys-bio/roadrunner/issues/420#issuecomment-346536846, or mute the thread https://github.com/notifications/unsubscribe-auth/ABAZDsGrZIjbVPpdl2R-addsdLw9q1loks5s5QlRgaJpZM4Pq4_p .

kirichoi commented 6 years ago

We reached on a conclusion that our reduced Jacobian might be correct after all. Our link matrix is technically correct, I think, but simply picking different species as dependent species in udp, utp, and udpglc cycle.

As for the reason why we can't calculate steady-state, the reduced Jacobian is not invertible to begin with. We will try to implement nleq2 as suggested in the issue https://github.com/sys-bio/roadrunner/issues/426, while attempting to expose more options for steady-state solver. We also intend to implement algorithm similar to COPASI, where is looks for steady-state solution multiple time until it finds one.