How to model the effect of acid/base on pH
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 &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")


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: