2
$\begingroup$

I would like to model a biological compartment with sysBio R package. It is fine but I'm having a difficult time with modeling the change in pH due to acidosis. I tried several things like:

addMAreaction(glu, "pH = null", "acidosis", "base excess")
addMAreactRate(glu, "acidosis", "assigned", "ka*pH")
addMAreactRate(glu, "base excess", "assigned", "ka*pH*-1")

which caused continuously decreasing pH levels. I'll leave out the intermediate steps but the simplified version of the last code looks like this:

library("sysBio")

glu <- newModel("glu")
addSpecies(glu, "Lactate", 0)
addSpecies(glu, "Glucose", 1000)
addSpecies(glu, "pH", 7.4)
addSpecies(glu, "Acid", 0)

addMAreaction(glu, "Glucose = Lactate + 0.01*Acid", "forwardL", 
              "backwardL", name = "glycolysis, sort of")

addMAreactRate(glu, "forwardL", "assigned", "kgl*Glucose")  
addMAreactRate(glu, "backwardL", "assigned", "ka*pH+klg*Lactate")

addParameters(glu, "kgl", 0.01)
addParameters(glu, "klg", 0.01)
addParameters(glu, "ka", 0.05)
addParameters(glu, "a", 0.001)

addRule(glu, "pH rule", "ODEs", "pH=7.4-a*Acid")

makeModel(glu)
simResults <- simulateModel(glu, times = seq(0,10,0.1))
plotResults(simResults, title="glu")

And the resulting plot is this where pH is increasing forever despite Acid is at a steady state concentration.

Plot of the glu.df

Please note that the abovementioned code is just a simplification. It doesn't even look like aerobic glycolysis where 1 glucose produces lactate and 2 hydrogen atoms etc.

I intend to include several acid- and base-creating reactions and evaluate the resulting pH. So could someone point me in the right direction about how to model the relationship between pH and acids or bases?

The answer does not need to be explicit code. I'll appreciate any solid comment on my design errors.

Improve this question
$\endgroup$
3
  • $\begingroup$ I'm unfamiliar with sysBio, but it looks like you may just be plotting a relationship, and need to know the relationship. Is the Hendersen Hasselbalch equation running behind the scenes, or is that what you're looking for? I'd note, if you're trying to model acid base in human physiology, you're making a number of assumptions by leaving out the major systems for acid base regulation. $\endgroup$
    – De Novo
    Feb 1, 2019 at 22:14
  • $\begingroup$ You are right. But please keep in mind that this is just a minimal working code sufficient to reproduce the error. You may think of it like a reaction where A reversibly converts to 2*B + 2+H, and H reversibly causes P to sink to 0. At least this was how I modeled it. For example I tried to add rules to limit pH within 0 and 14 as well but these are not included in the example. Apparently some familiarity with either SBML, ordinary differential equations, or maybe sysBio is needed. I couldn't reach the original author of the package so I'm trying to see if I'm missing sth fundamental. $\endgroup$
    – barerd
    Feb 2, 2019 at 4:23
  • $\begingroup$ Since your question pertains to the use of a package and not really about the biology behind the model, your question would not be on-topic here. It should, however, be on-topic on Bioinformatics. $\endgroup$
    – WYSIWYG
    Feb 2, 2019 at 9:53

1 Answer 1

Reset to default
2
$\begingroup$

I figured that one has to include only the changes in the rules formula. Keeping the initial value of the Species or Parameter when defining circular or self-referencing formulas caused the error. So the correct rule becomes:

pH=-a*Acid instead of pH=pH-a*Acid or pH=7.4-a*Acid.

For convenience, I include the correct code crowded with several reactions producing both acids and bases, where pH gradually decreases as expected. The names are fictional, and just are better placeholders instead of x1 and x2, etc.

library("sysBio")

glu <- newModel("glu")
addSpecies(glu, "Glucose", 1000)
addSpecies(glu, "Lactate", 0)
addSpecies(glu, "Pyruvate", 0)
addSpecies(glu, "Oxalate", 0)
addSpecies(glu, "Acid", 0)
addSpecies(glu, "Base", 0)
addSpecies(glu, "pH", 7.4)

addMAreaction(glu, "Glucose = Lactate + Acid", "forwardL", 
              "backwardL", name = "glycolysis")
addMAreaction(glu, "Lactate = Pyruvate + Acid", "forwardP", 
              "backwardP", name = "lactolysis")
addMAreaction(glu, "Pyruvate = Oxalate + Base", "forwardO", 
              "backwardO", name = "oxalate synthesis")

addMAreactRate(glu, "forwardL", "assigned", 
               "k_glucose_to_lactate*Glucose")
addMAreactRate(glu, "backwardL", "assigned", 
               "k_acid*pH+k_lactate_to_glucose*Lactate")
addMAreactRate(glu, "forwardP", "assigned", 
               "k_lactate_to_pyruvate*Lactate")
addMAreactRate(glu, "backwardP", "assigned", 
               "k_acid*pH+k_pyruvate_to_lactate*Pyruvate")
addMAreactRate(glu, "forwardO", "assigned", 
               "k_pyruvate_to_oxalate*Pyruvate")
addMAreactRate(glu, "backwardO", "assigned", 
               "k_base*pH+k_oxalate_to_pyruvate*Oxalate")

addParameters(glu, "k_glucose_to_lactate", 0.02)
addParameters(glu, "k_lactate_to_glucose", 0.001)
addParameters(glu, "k_lactate_to_pyruvate", 0.01)
addParameters(glu, "k_pyruvate_to_lactate", 0.0001)
addParameters(glu, "k_pyruvate_to_oxalate", 0.01)
addParameters(glu, "k_oxalate_to_pyruvate", 0.001)
addParameters(glu, "k_acid", 0.0001)
addParameters(glu, "k_base", 0.0001)

addRule(glu, "pH rule", "ODEs", "pH=-k_acid*Acid+k_base*Base")

makeModel(glu)
simResults <- simulateModel(glu, seq(0,3,0.01))
plotResults(simResults, title="glu")

The head of the resulting data.frame is:

  time   Glucose  Lactate   Pyruvate      Oxalate     Acid         Base       pH
1 0.00 1000.0000   0.0000  0.0000000 0.000000e+00   0.0000 0.000000e+00 7.400000
2 0.01  844.0001 155.0615  0.9384259 1.402282e-05 156.9384 1.402282e-05 7.399914
3 0.02  785.0869 210.3646  4.5476263 7.892761e-04 219.4615 7.892761e-04 7.399719
4 0.03  772.1557 218.6462  9.1924839 5.655403e-03 237.0425 5.655403e-03 7.399489
5 0.04  768.8212 217.2720 13.8876293 1.917807e-02 245.0856 1.917807e-02 7.399248
6 0.05  767.0383 214.4883 18.4278053 4.550630e-02 251.4350 4.550630e-02 7.399000

and resulting graph looks like:

correct plot

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.