Negative Values During the Simulation

[gmgmgm]

Hello everyone, I am trying to model the growth of sulfate-reducing bacteria in an environment with lactate and sulfate as substrates. During this process, lactate sometimes appears as a negative value, which leads to the error "Bad RK steps > 500 in cell". Below is my model. I appreciate your time and help.

Code: [Select]

SOLUTION_MASTER_SPECIES
    Lactate       Lactate          0     Lactate         90
    Acetate       Acetate          0     Acetate         60
SOLUTION_SPECIES
Acetate = Acetate
    log_k     0
Lactate = Lactate
    log_k     0
RATES
    Lactate_degradation
-start
 10 L = TOT("Lactate")
 20 muL = 0.0000417
 40 KL = 0.0045
 60 YL = 5.65
 80 if ( L<=0) then goto 180
150 rate = muL * (KIN("Biomass") + 0.002) / YL * [L /(KL + L)]
160 moles = rate * TIME
170 PUT(moles, 1)
180 SAVE moles
-end
    Sulfate_reduction
-start
 10 S = TOT("S(6)")
 30 muS = 0.000022   
 50 KS = 0.00005
 70 YS = 4.45       
 80 if ( S<=0 ) then goto 180 
150 rate = muS * (KIN("Biomass") + 0.002)/ YS * [S /(KS + S)]
160 moles = rate * TIME
170 PUT(moles, 2)
180 SAVE moles
-end
    Biomass
-start
 60 YL = 5.65   
 70 YS = 4.45       
200 rate = -YL * get(1) -YS * get(2) + 2.315e-7 * (KIN("Biomass") + 0.002)
210 dB = rate * time
220 save dB
-end

END

SOLUTION 1
    temp      25
    pH        3.6
    pe        4
    redox     pe
    units     mmol/l
    density   1
    Acetate   0.5
    Lactate   11.2
    S(6)      10.7
    -water    1 # kg
END
USE solution 1
KINETICS 1
Lactate_degradation
    -formula  Acetate  1 H2O  -1 Lactate  -1
    -m        1
    -m0       1
    -tol      1e-08
Sulfate_reduction
    -formula  SO4-2  -1 H+  -2 H2S  1 H2O  4
    -m        1
    -m0       1
    -tol      1e-08
Biomass
    -formula  H2O  1
    -m        1e-06
    -m0       1e-06
    -tol      1e-08
-steps       21600 in 60 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500

USER_GRAPH 1
    -headings               Time
    -axis_titles            "Time(h)" "Growth (mg/L)" ""
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 graph_x TOTAL_TIME/3600
20 graph_y kin("Biomass") * 1000
  -end
    -active                 true


USER_GRAPH 3
    -headings               Time
    -axis_titles            "Time(h)" "Sulphate (mM)" ""
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 graph_x TOTAL_TIME/3600
20 graph_y 1000 * tot("S(6)")
  -end
    -active                 true

USER_GRAPH 2
    -headings               Time Lactate Acetate
    -axis_titles            "Time(h)" "Concentration, mM" "Concentration, mM"
    -initial_solutions      true
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 graph_x TOTAL_TIME/3600
20 graph_y 1000 * tot("Lactate")
30 graph_sy 1000 * tot("Acetate")
  -end
    -active                 true

END

[dlparkhurst]

Your kinetic reaction formulas need to be different. -formula in KINETICS adds the net change in elements defined by the reaction. Charge is not considered, so your sulfate reduction reaction actually adds 8 H to the solution. You should consider the C,nH, and O change; when lactate is converted to acetate, CH2O is released. Below is a data block that assumes sulfate reduction occurs at the same rate as lactate conversion. You will have to do some more work if you want lactate to react at a different rate than sulfate. Possibly putting the carbon from lactate in a pool that determines the sulfate reduction reaction. Ultimately, all the carbon from lactate should be consumed.

Code: [Select]

KINETICS 1
Lactate_degradation
# C3H6O3  Lactic acid
# C2H4O2  Acetic acid
    -formula  Acetate  1 CH2O 1  Lactate  -1
    -m        1
    -m0       1
    -tol      1e-08
#Sulfate_reduction
#    -formula  SO4-2  -1 H+  -2 H2S  1 H2O  4
#    -m        1
#    -m0       1
#    -tol      1e-08
Biomass
    -formula  H2O  1
    -m        1e-06
    -m0       1e-06
    -tol      1e-08
-steps       21600 in 60 steps # seconds

[dlparkhurst]

I think there is also a logic error in your rates definitions. When you "SAVE moles" it is the product of rate times TIME. You multiply by TIME again in the Biomass definition. You should base each rate on the solution at the current time (as incremented automatically by KINETICS). You should not predict the concentration at the end of the time step. The kinetics algorithms are not straightforward extrapolations.

I found this problem by running the -cvode option. It leads to a different answer than the RK. I think -cvode is closer to the correct answer, but you need to fix the rate definitions.

Code: [Select]

RATES
    Lactate_degradation
-start
 10 L = TOT("Lactate")
 20 muL = 0.0000417
 40 KL = 0.0045
 60 YL = 5.65
 80 if ( L<=0) then goto 180
150 rate = muL * (KIN("Biomass") + 0.002) / YL * [L /(KL + L)]
160 moles = rate * TIME
170 PUT(moles, 1)
180 SAVE moles
-end
    Biomass
-start
 60 YL = 5.65   
 70 YS = 4.45       
200 rate = -YL * get(1) -YS * get(2) + 2.315e-7 * (KIN("Biomass") + 0.002)
210 dB = rate * time
220 save dB
-end