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

<p><code>pH=-a*Acid</code> instead of <code>pH=pH-a*Acid</code> or <code>pH=7.4-a*Acid</code>. </p>

<p>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.</p>

<pre><code>library("sysBio")

glu &lt;- 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 &lt;- simulateModel(glu, seq(0,3,0.01))
plotResults(simResults, title="glu")
</code></pre>

<p>The head of the resulting data.frame is:</p>

<pre><code>  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
</code></pre>

<p>and resulting graph looks like:</p>

<p><a href="https://i.sstatic.net/kg6bG.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/kg6bG.png" alt="correct plot"></a></p>

CC BY-SA 4.0
2019-02-08T22:26:03.303Z
2019-02-08T22:26:03.303Z
barerd
<p>I would like to model a biological compartment with <a href="https://github.com/Vessy/sysBio" rel="nofollow noreferrer">sysBio R package</a>. It is fine but I'm having a difficult time with modeling the change in pH due to acidosis. I tried several things like:</p>

<pre><code>addMAreaction(glu, "pH = null", "acidosis", "base excess")
addMAreactRate(glu, "acidosis", "assigned", "ka*pH")
addMAreactRate(glu, "base excess", "assigned", "ka*pH*-1")
</code></pre>

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

<pre><code>library("sysBio")

glu &lt;- 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 &lt;- simulateModel(glu, times = seq(0,10,0.1))
plotResults(simResults, title="glu")
</code></pre>

<p>And the resulting plot is this where pH is increasing forever despite <code>Acid</code> is at a steady state concentration. </p>

<p><a href="https://i.sstatic.net/qUych.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/qUych.jpg" alt="Plot of the glu.df"></a></p>

<p><strong>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.</strong></p>

<p>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?</p>

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

CC BY-SA 4.0
2019-02-02T09:53:52.697Z
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.se].
WYSIWYG
<r><modelling>
How to model the effect of acid/base on pH
# How to model the effect of acid/base on pH

(Asked by: username_0 on 2019-02-01)

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. 

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.

**Tags:** r, modelling
> **Comments:**
> **username_1**: 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.se].  

---

## 1 answers

### Answer by username_0 on 2019-02-08. Score: 2 | Accepted answer

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:


---

