Using PEST for Bacteria Growth Parameters

[mek001]

Dear All,
I am trying to model denitrifying conditions. I have NO3 and NO2 concentrations over time. I found the source code, it looks like it works but it does not fit exactly. I searched on the internet, and up to this point I only saw PEST usage on log_k values in equilibrium. What I want is to PEST to fit values in the RATES section. Can I use PEST for such a fitting? If I can, how should I write my observation data in a .pst file? I am planning to optimize  Kn, Ka, and kv1 values.
Regards

Code:

Code: [Select]

SOLUTION 1 Eau DNS2
    temp      25
    pH        5.55
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    Acetate   0.008475
    C         0.0077 as HCO3
    Ca        0.0003818
    Cl        0.0007786
    Fe        0.0001255
    K         3.325e-05
    Mg        0.0003661
    N_v_      0.039 mmol/kgw   as N_v_O3-
    Na        0.00228
    S         0.0003363 as SO4
    Si        0.0005427
    -water    2.085 # kg
end

RATES
Denit
-start
1  REM Vitesse de dénitrification NO3/NO2 sur Acetate
2  SR1 = SR("E_N(3)/N(5)")
       3  SR2 = SR("E_HAcetate/HCO3-")
6  dg_min = -50 
7  ksi = 2       # turn-over of limited step / mol Nitrates (from Jin, Personal Communication)
8  n = 2/3 # mol ATP produced / mol Nitrates (from Jin, Personal Communication)
      10  Kn = 0.66 / 62000 #related to no3
      20  Ka = 1.20 / 59000 #related to acetate
30  kv1 = 1.60e-5/86400
50  KMono1 =(Mol("N_v_O3-")/(Mol("N_v_O3-")+ Kn))*(Mol("Acetate-")/(Mol("Acetate-")+ Ka))
55  rate1 = kv1 * KIN("biomass")*KMono1*(1-(SR1/SR2)^(1/ksi)*exp(-n*dg_min*1e3/(ksi*8.314*(t+273.15))))
# 56  IF (SR1/SR2)*exp(-n*dg_min*1e3/(8.314*(t+273.15)))>10^(-0.01) THEN rate1 = 0     # pe calculée à 0.01 près
56  if (Mol("N_v_O3-")< 1e-16) then rate1 = 0
      60  mole1 = rate1 * TIME
80  PUT(rate1, 1)
90  PUT(mole1, 2)
     100  SAVE mole1
-end

Denit2
-start
1  REM Vitesse de dénitrification NO2/N2O sur Acetate
2  SR3 = SR("E_N(1)/N(3)")
       3  SR4 = SR("E_HAcetate/HCO3-")
6  dg_min = -50 
7  ksi = 2       # turn-over of limited step / mol Nitrates (from Jin, Personal Communication)
8  n = 2/3 # mol ATP produced / mol Nitrates (from Jin, Personal Communication)
      10  Km = 0.66 / 46000 #related to no2
      20  Kb = 1.20 / 59000 #related to acetate
30  kv2 = 9.000e-6/86400
50  mole1 = GET(2)
      60  MO1 = MOL("Acetate-") - 0.25*mole1
70  KMono2 = (Mol("N_iii_O2-")/(Mol("N_iii_O2-")+ Km))*(MO1/(MO1 + Kb))
80  rate2 = kv2*KIN("biomass")*KMono2*(1-(SR3/SR4)^(1/ksi)*exp(-n*dg_min*1e3/(ksi*8.314*(t+273.15))))
85  if (Mol("N_iii_O2-") < 1e-16) then rate2 = 0
      90  mole2 = rate2 * TIME
     100  PUT(rate2, 3)   
     110  PUT(mole2, 4) 
     120  SAVE mole2
-end

Denit3
-start
1  REM Vitesse de dénitrification N2O/N2 sur Acetate
2  SR5 = SR("E_N(0)/N(1)")
       3  SR6 = SR("E_HAcetate/HCO3-")
6  dg_min = -50 
7  ksi = 2       # turn-over of limited step / mol Nitrates (from Jin, Personal Communication)
8  n = 2/3 # mol ATP produced / mol Nitrates (from Jin, Personal Communication)
      10  Kp = 0.66 / 44000 #related to N2o?
      20  Kc = 1.20 / 59000 #related to acetate
50  kv3 = 0.00
25  mole1 = GET(2)
30  mole2 = GET(4)
      40  MO2 = MOL("Acetate-") - 0.25*mole1 - 0.25*mole2
50  KMono3 = (Mol("N_i_2O")/(Mol("N_i_2O")+ Kp)) * (MO2/(MO2 + Kc))
60  if (MOL("N_v_O3-") < 1e-03) then kv3 = 8.500e-6/86400
#70  if (MOL("N_v_O3-") < 1e-08 and MOL("N_iii_O2-") < 1e-07) then kv3 = 0.150e-6/86400
80  rate3 = kv3 * KIN("biomass")*KMono3 *(1-(SR5/SR6)^(1/ksi)*exp(-n*dg_min*1e3/(ksi*8.314*(t+273.15))))
85  if (Mol("N_i_2O")<1e-16) then rate3 = 0
      90  mole3 = rate3 * TIME
     100  PUT(rate, 5)   
     110  SAVE mole3
-end

Biomass
      -start
10 Y1 = 0.10*62000 # mg biomasse/moles Nitra
15 Y2 = 0.40*46000 # mg biomasse/moles Nitrites
16 Y3 = 0.10*44000 # mg biomasse/moles Protoxyde d'azote
20 b = 0.060/86400 # Facteur de décroissance (s-1)
30 rate1 = GET(1)
35 rate2 = GET(3)
36 rate3 = GET(5)
40 rate = Y1*rate1 + Y2*rate2 + Y3*rate3 - b * M
50 moles = - rate * TIME
60 SAVE moles
        -end
use solution 1

KINETICS 1
      Denit
-formula N_v_O3- -1.0 Acetate- -0.25 H+ -0.25 N_iii_O2- 1.0 CO2 0.5 H2O 0.5     
Denit2
-formula N_iii_O2- -1.0 Acetate- -0.25 H+ -1.25 N_i_2O 0.5 CO2 0.5 H2O 1.0
      Denit3
-formula N_i_2O -1.0 Acetate- -0.25 H+ -0.25 N_2 1.0 CO2 0.5  H2O 0.5
Biomass
            -formula H 0.0
            -m 1.913 # Concentration initiale en biomasse
-steps 561600 in 65 steps

INCREMENTAL_REACTIONS True


[dlparkhurst]

Pest does not care where an optimization parameter occurs in the input file. You just have to replace the parameter with a string that is defined to PEST. You can replace 0.66, 1.2, and 1.60e-5 with appropriate strings for the optimization.

Note that the -formula of KINETICS does not recognize the charge on the species. What you have will work because the charge is ignored and your equations are balanced. However, you are not really changing species concentrations, only element totals in solution (Acetate, N_v_, N_iii_, N_i_ and N_ are considered elements because they are defined in SOLUTION_MASTER_SPECIES). So, only the net change in elements is important. The following is the way PHREEQC actually works:

Code: [Select]

      Denit
-formula N_v_ -1.0  Acetate -0.25 N_iii_ +1.0 C2H3OO +1

or to my aesthetic for charge balance

      Denit
-formula N_v_ -1.0  HAcetate -0.25 N_iii_ +1.0 C2H3OOH +1

The other formulas are similar.

[mek001]

Dear dlparkhurst,
Thank you for your kind answer. I made the changes according to your advice. Now, PEST gives me an error. I created a template file. I have the initial NO3 value, and measurements on day 2 (172800 seconds) and day 9 (777600 seconds) and I wrote a simple code, if time is equal to these two, I am getting mmol/l values of NO3 and I will use these values for parameter estimation. I thought it needs to do something in the "ELSE" case, so I added Q=1 to just fill.

Code: [Select]

ptf@
SELECTED_OUTPUT 1
    -file                 Results.txt
    -high_precision       true
    -reset                false
USER_PUNCH 1
    -headings
    -start
10 if (TOTAL_TIME=172800) THEN GOTO 20 ELSE GOTO 30
20 PUNCH MOL("N_v_O3-")*1000
30 Q=1
40 if (TOTAL_TIME=777600) THEN GOTO 50 ELSE GOTO 60
50 PUNCH MOL("N_v_O3-")*1000
60 Q=1
    -end
SOLUTION 1 Eau DNS2
    temp      25
    pH        5.55
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    Acetate   0.008475
    C         0.0077 as HCO3
    Ca        0.0003818
    Cl        0.0007786
    Fe        0.0001255
    K         3.325e-05
    Mg        0.0003661
    N_v_      0.039 mmol/kgw   as N_v_O3-
    Na        0.00228
    S         0.0003363 as SO4
    Si        0.0005427
    -water    2.085 # kg

RATES
Denit
-start
1  REM Vitesse de dénitrification NO3/NO2 sur Acetate
2  SR1 = SR("E_N(3)/N(5)")
       3  SR2 = SR("E_HAcetate/HCO3-")
6  dg_min = -50 
7  ksi = 2       # turn-over of limited step / mol Nitrates (from Jin, Personal Communication)
8  n = 2/3 # mol ATP produced / mol Nitrates (from Jin, Personal Communication)
      10  Kn = 0.66 / 62000 #related to no3
      20  Ka = 1.20 / 59000 #related to acetate
30  kv1 = 1.60e-5/86400
50  KMono1 =(Mol("N_v_O3-")/(Mol("N_v_O3-")+ Kn))*(Mol("Acetate-")/(Mol("Acetate-")+ Ka))
55  rate1 = kv1 * KIN("biomass")*KMono1*(1-(SR1/SR2)^(1/ksi)*exp(-n*dg_min*1e3/(ksi*8.314*(t+273.15))))
# 56  IF (SR1/SR2)*exp(-n*dg_min*1e3/(8.314*(t+273.15)))>10^(-0.01) THEN rate1 = 0     # pe calculée à 0.01 près
56  if (Mol("N_v_O3-")< 1e-16) then rate1 = 0
      60  mole1 = rate1 * TIME
80  PUT(rate1, 1)
90  PUT(mole1, 2)
     100  SAVE mole1
-end

Denit2
-start
1  REM Vitesse de dénitrification NO2/N2O sur Acetate
2  SR3 = SR("E_N(1)/N(3)")
       3  SR4 = SR("E_HAcetate/HCO3-")
6  dg_min = -50 
7  ksi = 2       # turn-over of limited step / mol Nitrates (from Jin, Personal Communication)
8  n = 2/3 # mol ATP produced / mol Nitrates (from Jin, Personal Communication)
      10  Km = 0.66 / 46000 #related to no2
      20  Kb = 1.20 / 59000 #related to acetate
30  kv2 = 9.000e-6/86400
50  mole1 = GET(2)
      60  MO1 = MOL("Acetate-") - 0.25*mole1
70  KMono2 = (Mol("N_iii_O2-")/(Mol("N_iii_O2-")+ Km))*(MO1/(MO1 + Kb))
80  rate2 = kv2*KIN("biomass")*KMono2*(1-(SR3/SR4)^(1/ksi)*exp(-n*dg_min*1e3/(ksi*8.314*(t+273.15))))
85  if (Mol("N_iii_O2-") < 1e-16) then rate2 = 0
      90  mole2 = rate2 * TIME
     100  PUT(rate2, 3)   
     110  PUT(mole2, 4) 
     120  SAVE mole2
-end

Denit3
-start
1  REM Vitesse de dénitrification N2O/N2 sur Acetate
2  SR5 = SR("E_N(0)/N(1)")
       3  SR6 = SR("E_HAcetate/HCO3-")
6  dg_min = -50 
7  ksi = 2       # turn-over of limited step / mol Nitrates (from Jin, Personal Communication)
8  n = 2/3 # mol ATP produced / mol Nitrates (from Jin, Personal Communication)
      10  Kp = 0.66 / 44000 #related to N2o?
      20  Kc = 1.20 / 59000 #related to acetate
50  kv3 = 0.00
25  mole1 = GET(2)
30  mole2 = GET(4)
      40  MO2 = MOL("Acetate-") - 0.25*mole1 - 0.25*mole2
50  KMono3 = (Mol("N_i_2O")/(Mol("N_i_2O")+ Kp)) * (MO2/(MO2 + Kc))
60  if (MOL("N_v_O3-") < 1e-03) then kv3 = 8.500e-6/86400
#70  if (MOL("N_v_O3-") < 1e-08 and MOL("N_iii_O2-") < 1e-07) then kv3 = 0.150e-6/86400
80  rate3 = kv3 * KIN("biomass")*KMono3 *(1-(SR5/SR6)^(1/ksi)*exp(-n*dg_min*1e3/(ksi*8.314*(t+273.15))))
85  if (Mol("N_i_2O")<1e-16) then rate3 = 0
      90  mole3 = rate3 * TIME
     100  PUT(rate, 5)   
     110  SAVE mole3
-end

Biomass
      -start
10 Y1 = 0.10*62000 # mg biomasse/moles Nitra
15 Y2 = 0.40*46000 # mg biomasse/moles Nitrites
16 Y3 = 0.10*44000 # mg biomasse/moles Protoxyde d'azote
20 b = 0.060/86400 # Facteur de décroissance (s-1)
30 rate1 = GET(1)
35 rate2 = GET(3)
36 rate3 = GET(5)
40 rate = Y1*rate1 + Y2*rate2 + Y3*rate3 - b * M
50 moles = - rate * TIME
60 SAVE moles
        -end
use solution 1

KINETICS 1
Denit
    -formula  Acetate-  -0.25 C2H3OO  1 N_iii_  1 N_v_O3-  -1
    -m        1
    -m0       1
    -tol      1e-08
Denit2
    -formula  Acetate-  -0.25 CO2  0.5 H+  -1.25 H2O  1 N_i_2O  0.5 N_iii_O2-  -1
    -m        1
    -m0       1
    -tol      1e-08
Denit3
    -formula  Acetate-  -0.25 CO2  0.5 H+  -0.25 H2O  0.5 N_2  1 N_i_2O  -1
    -m        1
    -m0       1
    -tol      1e-08
Biomass
    -formula  H  0
    -m        1.913
    -m0       1.913
    -tol      1e-08
-steps       864000 in 240 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500

INCREMENTAL_REACTIONS True


control file: I have initial values for k1,k2, and k3 and I wrote them. The batch file tells me that "Insufficient entries on line 20 of file case.pst."

Code: [Select]

pcf
* control data
restart estimation
8 33 1 0 3
# no of parameters,observations, group parameters, prior infos , groups
1 1 single point 1 0 0
# no of template files, no of instruction files, ….
10.0 2.0 0.3 0.03 10
3.0 3.0 0.001
0.1
30 0.001 3 3 0.01 3
1 1 1
* parameter groups
k1 relative 1e-2 0.000001 switch 2 parabolic
* parameter data
k1 none relative 1.06e-5 1.0e-7 1.0e-3 k1 1.0 0.0 1
#param_name,type, initial, min, max, ..….
k2 none relative 2.03e-5 1.0e-7 1.0e-3 k1 1.0 0.0 1
k3 none relative 1.85e-10 1.0e-12 1.0e-5 k1 1.0 0.0 1
* observation groups
group_1
group_2
group_3
* observation data
c1 0.0295 1 group_1
# obs_name, value, weight, group no
c2 0.00418 1 group_1
* model command line
phreeqc.bat case.phrq case.out c:/phreeqc/database/new.dat scr.out
* model input/output
case.tpl case.phrq
case.ins case.prn
* prior information
Instruction file, I have only 2 measurements

Code: [Select]

pif @ # read line 1 or 2 (l1 or l2); data to be read between position 1 and 25 i
l2 [c1]1:25
l1 [c2]1:25

My batch file:

Code: [Select]

@echo on
pest case
pause
I will be very happy if you can help me.
Regards

[dlparkhurst]

It's been a while since I used Pest, but you need to use unique strings in the file for the parameters you want to fit. Something like

Code: [Select]

      10  Kn = xxxK1xxx/ 62000 #related to no3
      20  Ka = xxxK2xxx/ 59000 #related to acetate
      30  kv1 = xxxK3xxx/86400

If you only have data at two points, then run the calculations to those two points. BTW, you will not be able to fit three parameters with only two points. With INCREMENTAL_REACTIONS true, you would use the following times:

Code: [Select]

-steps   172800 604800

This is not complete, but your input file would look more like this:

Code: [Select]

ptf@
SELECTED_OUTPUT 1
    -file                 Results.txt
    -high_precision       true
    -reset                false
USER_PUNCH 1
    -headings
    -start
#10 if (TOTAL_TIME=172800) THEN GOTO 20 ELSE GOTO 30
20 PUNCH MOL("N_v_O3-")*1000
#30 Q=1
#40 if (TOTAL_TIME=777600) THEN GOTO 50 ELSE GOTO 60
#50 PUNCH MOL("N_v_O3-")*1000
#60 Q=1
    -end
SOLUTION 1 Eau DNS2
    temp      25
    pH        5.55
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    Acetate   0.008475
    C         0.0077 as HCO3
    Ca        0.0003818
    Cl        0.0007786
    Fe        0.0001255
    K         3.325e-05
    Mg        0.0003661
    N_v_      0.039 mmol/kgw   as N_v_O3-
    Na        0.00228
    S         0.0003363 as SO4
    Si        0.0005427
    -water    2.085 # kg
END
RATES
Denit
-start
1  REM Vitesse de dénitrification NO3/NO2 sur Acetate
2  SR1 = SR("E_N(3)/N(5)")
       3  SR2 = SR("E_HAcetate/HCO3-")
6  dg_min = -50
7  ksi = 2       # turn-over of limited step / mol Nitrates (from Jin, Personal Communication)
8  n = 2/3 # mol ATP produced / mol Nitrates (from Jin, Personal Communication)
      10  Kn = 0.66 / 62000 #related to no3
      20  Ka = 1.20 / 59000 #related to acetate
30  kv1 = 1.60e-5/86400
      10  Kn = xxxK1xxx/ 62000 #related to no3
      20  Ka = xxxK2xxx/ 59000 #related to acetate
30  kv1 = xxxK3xxx/86400
50  KMono1 =(Mol("N_v_O3-")/(Mol("N_v_O3-")+ Kn))*(Mol("Acetate-")/(Mol("Acetate-")+ Ka))
55  rate1 = kv1 * KIN("biomass")*KMono1*(1-(SR1/SR2)^(1/ksi)*exp(-n*dg_min*1e3/(ksi*8.314*(t+273.15))))
# 56  IF (SR1/SR2)*exp(-n*dg_min*1e3/(8.314*(t+273.15)))>10^(-0.01) THEN rate1 = 0     # pe calculée à 0.01 près
56  if (Mol("N_v_O3-")< 1e-16) then rate1 = 0
      60  mole1 = rate1 * TIME
80  PUT(rate1, 1)
90  PUT(mole1, 2)
     100  SAVE mole1
-end

Denit2
-start
1  REM Vitesse de dénitrification NO2/N2O sur Acetate
2  SR3 = SR("E_N(1)/N(3)")
       3  SR4 = SR("E_HAcetate/HCO3-")
6  dg_min = -50
7  ksi = 2       # turn-over of limited step / mol Nitrates (from Jin, Personal Communication)
8  n = 2/3 # mol ATP produced / mol Nitrates (from Jin, Personal Communication)
      10  Km = 0.66 / 46000 #related to no2
      20  Kb = 1.20 / 59000 #related to acetate
30  kv2 = 9.000e-6/86400
50  mole1 = GET(2)
      60  MO1 = MOL("Acetate-") - 0.25*mole1
70  KMono2 = (Mol("N_iii_O2-")/(Mol("N_iii_O2-")+ Km))*(MO1/(MO1 + Kb))
80  rate2 = kv2*KIN("biomass")*KMono2*(1-(SR3/SR4)^(1/ksi)*exp(-n*dg_min*1e3/(ksi*8.314*(t+273.15))))
85  if (Mol("N_iii_O2-") < 1e-16) then rate2 = 0
      90  mole2 = rate2 * TIME
     100  PUT(rate2, 3)   
     110  PUT(mole2, 4)
     120  SAVE mole2
-end

Denit3
-start
1  REM Vitesse de dénitrification N2O/N2 sur Acetate
2  SR5 = SR("E_N(0)/N(1)")
       3  SR6 = SR("E_HAcetate/HCO3-")
6  dg_min = -50
7  ksi = 2       # turn-over of limited step / mol Nitrates (from Jin, Personal Communication)
8  n = 2/3 # mol ATP produced / mol Nitrates (from Jin, Personal Communication)
      10  Kp = 0.66 / 44000 #related to N2o?
      20  Kc = 1.20 / 59000 #related to acetate
50  kv3 = 0.00
25  mole1 = GET(2)
30  mole2 = GET(4)
      40  MO2 = MOL("Acetate-") - 0.25*mole1 - 0.25*mole2
50  KMono3 = (Mol("N_i_2O")/(Mol("N_i_2O")+ Kp)) * (MO2/(MO2 + Kc))
60  if (MOL("N_v_O3-") < 1e-03) then kv3 = 8.500e-6/86400
#70  if (MOL("N_v_O3-") < 1e-08 and MOL("N_iii_O2-") < 1e-07) then kv3 = 0.150e-6/86400
80  rate3 = kv3 * KIN("biomass")*KMono3 *(1-(SR5/SR6)^(1/ksi)*exp(-n*dg_min*1e3/(ksi*8.314*(t+273.15))))
85  if (Mol("N_i_2O")<1e-16) then rate3 = 0
      90  mole3 = rate3 * TIME
     100  PUT(rate, 5)   
     110  SAVE mole3
-end

Biomass
      -start
10 Y1 = 0.10*62000 # mg biomasse/moles Nitra
15 Y2 = 0.40*46000 # mg biomasse/moles Nitrites
16 Y3 = 0.10*44000 # mg biomasse/moles Protoxyde d'azote
20 b = 0.060/86400 # Facteur de décroissance (s-1)
30 rate1 = GET(1)
35 rate2 = GET(3)
36 rate3 = GET(5)
40 rate = Y1*rate1 + Y2*rate2 + Y3*rate3 - b * M
50 moles = - rate * TIME
60 SAVE moles
        -end
END
use solution 1

KINETICS 1
Denit
    -formula  Acetate-  -0.25 C2H3OO  1 N_iii_  1 N_v_O3-  -1
    -m        1
    -m0       1
    -tol      1e-08
Denit2
    -formula  Acetate-  -0.25 CO2  0.5 H+  -1.25 H2O  1 N_i_2O  0.5 N_iii_O2-  -1
    -m        1
    -m0       1
    -tol      1e-08
Denit3
    -formula  Acetate-  -0.25 CO2  0.5 H+  -0.25 H2O  0.5 N_2  1 N_i_2O  -1
    -m        1
    -m0       1
    -tol      1e-08
Biomass
    -formula  H  0
    -m        1.913
    -m0       1.913
    -tol      1e-08
-steps   172800 604800
-step_divide 1
-runge_kutta 3
-bad_step_max 500

INCREMENTAL_REACTIONS True