Jump to content
Geochemist's Workbench Support Forum
TBush

Including new chemicals with custom rate laws

Recommended Posts

Hi everyone,

I have been using GWB to model populations of microbes cycling sulfur in sediments.

I am interested in including a population of fermentative microbes that degrade cellulose and provide the acetate that is used as an electron donor in reduction reactions (like sulfate reduction).

I understand that GWB has some functionality for specifying custom rate laws (such as for microbial respiration) but I am unsure as to how to include a new chemical species, because the GWB database does not currently contain cellulose.

Also, will it be a problem that the production of acetate will be governed by a custom rate law and the consumption of acetate will be governed by the normal GWB functions for microbial respiration for the populations I have included.

Thanks very much!

Regards,

Tim Bush

Share this post


Link to post
Share on other sites

Hi Tim,

 

The GWB's thermodynamic databases are fully editable, so you can add new species as you see fit. You'll need some properties of the species, such as molecular weight and charge. You'll also need to write a balanced reaction to form the new species from existing basis or redox species in the thermo dataset. Finally, you'll need a log K for that reaction at one or more temperatures. If you don't know the log K, a value of 500 tells the programs that no data is available.

 

You can edit the datasets in a text editor like Notepad, or if you have GWB10 you can use TEdit, the graphical thermo data editor. You should check out the Thermo Datasets Appendix in the GWB Reference Manual and the Using TEdit section in the GWB Essentials Guide for more info.

 

Custom rate laws are a different matter. If you're interested in the kinetics of redox transformations (microbially catalyzed or not), for example, you could use the GWB's built in rate-law to specify the reaction and describe how quickly it proceeds. With the built-in rate law, the reaction proceeds at a rate proportional to a rate constant times the concentration of one or more species. The custom rate law doesn't change the reaction, it simply allows you to describe the reaction rate in a different form, if necessary. Perhaps you include an if statement to use a different rate constant below a certain pH, or the rate constant is multiplied by the sum of the concentrations of species.

 

It's not a problem to have one reaction described by the built-in rate law and another described with a custom rate law.

 

Hope this helps,

 

Brian Farrell

Aqueous Solutions LLC

Share this post


Link to post
Share on other sites

Hi Brian,

That helps alot, thanks.

So with regard to cellulose (which is insoluble) what category would you recommend I add it under in thermo.tdat? Also, what does the program actually do if you put 500 to indicate that the log K is not known, does this mean that it just doesn't make it dissolve in the fluid at all?

Regards,

Tim

Share this post


Link to post
Share on other sites

Hi Brian,

I tried adding cellulose as a "mineral" and it got added to the list ok, but it is not showing up within react so I can't include it in the model.

Cellulose is C6H12O6 so I included its formation as C6H12O6 -> 3*CH3COO- + 3*H+ and set the log Ks to all be 500.

Do you have any idea why this is not working?

Regards,

Tim

Share this post


Link to post
Share on other sites

Hi Brian,

So having played around a bit, the problem seems to be when all of the log Ks at set to 500. If I set only one of the log Ks at some arbitrary temperature to some arbitrary number then it works, and cellulose can be included. However, if they are all set to 500 then it no longer recognises cellulose as being in thermo.tdat.

Does this make any sense?

Regards,

Tim

Share this post


Link to post
Share on other sites

Hi Tim,

 

Yes, this makes sense. If you set the log K for all 8 principal temperatures to 500, then there is no data for the program to use.

 

If the temperature of your run is set to one of the principal temperatures in the thermo dataset, the GWB applications take the log K corresponding to that temperature from the data table. Where temperature differs from a principal temperature, the application fits each table to a polynomial with respect to temperature. Keep in mind, though, that a log K of 500 indicates to the program that there is no data. The program will not load species, minerals, or gases whose log K values do not span the temperature range of the calculation. So if all of your log Ks are 500, the program will not consider your mineral.

 

I think adding cellulose as a mineral makes sense. But since you're interested in redox transformations involving cellulose and acetate, you should probably also make a fictive redox species, like cellulose(aq), so that you can decouple cellulose from acetate. For example, I might make a redox couple between cellulose(aq) and acetate such that cellulose is much more stable. The reaction would look like what you have now for the mineral cellulose. Then, I would change the reaction for the mineral cellulose so that cellulose = cellulose(aq), with a log K such that cellulose is much more stable (less soluble). Then, you could decouple the cellulose(aq)/ acetate pair and define a kinetic rate law for the redox transformation of cellulose to acetate. As long as it's a kinetically dominated system, the arbitrary log Ks you specify should be okay.

 

For more on redox kinetics and fictive redox species, see Chapter 17 in the Geochemical and Biogeochemical Reaction Modeling text.

 

Regards,

Brian

Share this post


Link to post
Share on other sites

Hi Brian,

Thanks for that, that was really helpful.

I have tried to do as you suggested, but I think maybe I have done it a bit wrong.

Cellulose -> Cellulose (aq) has a log K of -100 for low temperatures

Cellulose(aq) -> 3*Acetate + 3*H+ has a log K of -100 for low temperatures

Does this mean that Cellulose -> Cellulose (aq) will tend to keep Cellulose as Cellulose and not Cellulose(aq), and Cellulose(aq) -> 3*Acetate + 3*H+ will tend

to stop Cellulose(aq) from dissolving into acetate? If Cellulose is more stable that Cellulose(aq) then surely this means that there will never be any available Cellulose(aq) to be transformed into acetate?

Also, how do i specify the rate law, is it as follows? Or should the rate law be in terms of Cellulose(aq)?

kinetic Cellulose rate_law = Cellulose.bas rate constant = 10^-6

Finally, in order to decouple Cellulose(aq) from acetate, should I use "decouple acetate"?

Regards,

Tim


Regards,

Tim

Share this post


Link to post
Share on other sites

Hi Brian,

I just had a quick additional question about the motivation behind this fictive redox species Cellulose(aq).

Is the idea that this allows me to vary the Cellulose without varying the acetate (because without the fictive intermediate redox species some of the cellulose might dissolve into acetate) and then all of the acetate can be generated by the kinetics of microbial respiration (which is what I want)? If so, then does this mean that the concentration of Cellulose(aq) in the system should basically always be zero?

Regards,
Tim

Share this post


Link to post
Share on other sites

Hi Brian,

Another question. Is it possible to specify the reaction for a custom rate law? What I want is for microbes to be doing the reaction Cellulose -> 3*Acetate + 3*H+ but with only one Monod term. Do I even need a custom rate law for that, and if so, is it possible to specify that this is the reaction?

Thanks for all of your help.

Regards,

Tim

Share this post


Link to post
Share on other sites

Hi Brian,

I just realised as well, would it not make more sense to set the log Ks as follows:

Cellulose(aq) -> 3*Acetate + 3*H+ highly predominates towards acetate

Cellulose -> Cellulose (aq) highly predominates towards Cellulose

 

then decouple acetate from Cellulose(aq).

 

This way, Cellulose will be decoupled from acetate by the fictive species cellulose(aq). The way I have the log Ks at the moment (as described above, very negative for both) means that all acetate produced

ends up precipitating out as Cellulose.

Maybe this is what you meant in the first place and I just got confused?


Sorry for the multiple posts all at once.

Thanks.

Regards,

Tim

Share this post


Link to post
Share on other sites

Hi Tim,

 

You're correct about the log K for the reaction between the cellulose and acetate. Because you're observing the cellulose reacting to form acetate, acetate needs to be the more stable form. Sorry about the mixup.

 

The command "decouple acetate" would decouple acetate from HCO3-. You want to decouple cellulose(aq) from acetate, so you'd need the "decouple cellulose(aq)" command. Using the GUI (Config - Redox couples), you'll be able to see the redox pair.

 

Whether you're using the Built-in rate law or a custom rate law, you specify the transformation reaction for two types of kinetic reactions: redox reactions and microbially catalyzed redox reactions. I'm not sure you need a custom rate law, however. Try looking at section 4.7.2 in the GWB Reaction Modeling Guide for a comparison of various rate laws used in microbial ecology and how to set them up. If you do end up using a custom rate law, a simple equation is all that's needed for a single-Modod type rate law, not a basic script.

 

A few more things to think about. Are you sure that the stoichiometry of the degradation reaction is correct? It doesn't look like a redox reaction. And just remember, because the log Ks are made up, this setup is basically assuming that thermodynamic limitations are not important and the cellulose degradation reaction will go to completion. You don't want to go using cellulose in any equilibrium calculations or activity diagrams.

 

Regards,

Brian

Share this post


Link to post
Share on other sites

Hi Brian,

Another question.

At the moment the model is working fine, but the size of the population of cellulose degrading microbes doesn't seem to change. I tried to

get it to update the biomass concentration using:

setgwbvar("biomass", newbiom)

but it didn't work.

Is it because that command only works if you actually write a .bas script?

If so is there some way to get it to update the biomass without writing a separate script?


Thanks for all of your help.

Regards,

Tim

Share this post


Link to post
Share on other sites
data = thermo.tdat

decouple Cellulose(aq)


# Properties of the sediment column


time end = 1 yr


pH = 8.1

Na+ = 3000 umolal

Cl- = 3000 umolal

SO4-- = 20.0 mmolal

CH3COO- = 0.1 umolal

O2(aq) = 100 umolal

CO2(aq) = 100 umolal

HCO3- = .8 umolal

HS- = .1 umolal

CH4(aq) = .1 umolal

Fe++ = 400.0 umolal

Cellulose = 0.00001 umolal

Porosity = 0.9

balance on Na+


react 10.0 g of Cellulose


# (1) Sulfate reducers

kinetic microbe-1 \

rxn = "CH3COO- + SO4-- -> 2*HCO3- + HS-", \

biomass = 0.001, rate_con = 1e-9, KA = 6.8e-5, KD = 5e-6, \

mpower(CH3COO-) = 1, mpowerD(CH3COO-) = 1, \

mpower(SO4--) = 1, mpowerA(SO4--) = 1, \

ATP_energy = -46.75, ATP_number = 1, order1 = 1/5, \

growth_yield = 4300, decay_con = 1e-8


# (2) Aerobic acetate oxidizers

kinetic microbe-2 \

rxn = "CH3COO- + 2*O2(aq) -> 2*HCO3- + H+", \

biomass = 0.001, rate_con = 1e-9, KA = 1e-7, KD = 5e-5, \

mpower(CH3COO-) = 1, mpowerD(CH3COO-) = 1, \

mpower(O2(aq)) = 1, mpowerA(O2(aq)) = 1, \

ATP_energy = -45.00, ATP_number = 8, order1 = 1/8, \

growth_yield = 20000, decay_con = 1e-8


# (3) Methanogens

kinetic microbe-3 \

rxn = "CH3COO- + H2O -> CH4(aq) + HCO3-", \

biomass = 0.001, rate_con = 1e-9, KA = 0, KD = 1e-3, \

mpower(CH3COO-) = 1, mpowerD(CH3COO-) = 1, \

ATP_energy = -45.00, ATP_number = 0.5, order1 = 1/2, \

growth_yield = 2000, decay_con = 1e-8



# (4) Cellulose degraders

kinetic microbe-4 rate_law = 'rate_con*molality("Cellulose")/(molality("Cellulose")+KD)', \

rate_con = 1e-9, KD = 1e-6, \

rxn = "Cellulose -> 3*CH3COO- + 3*H+", \

biomass = 1000.001, \


# (5) Green Sulfur Bacteria

kinetic microbe-5 rate_law = 'rate_con*molality("H2S(aq)")/(molality("H2S(aq)")+KD)', \

rate_con = 1e-9, KD = 1e-6, \

rxn = "H2S(aq) -> S-- + 2*H+", \

biomass = 0.001, \

Share this post


Link to post
Share on other sites

Hi Brian,

Another quick question.

The geochemists workbench documentation says that the default units of the rate constant (rate_con) are mol/cm2 sec, but when i input a rate_con React seems to default to the units mol/mg sec. What does this mean and why is it different? Does it mean moles reacted per mg of sediment per second?

Regards,

Tim

Share this post


Link to post
Share on other sites

Hi Tim,

 

If you go to "more reply options" you have the option to add attachments. In the future, you can attach your .rea files and custom thermo datasets instead of pasting in the commands. By the way, whenever you modify a thermo dataset it's a good idea to give it a unique name instead of saving over the standard datasets that we distribute.

 

It takes more than just the one line setgwbvar("biomass", newbiom) to have the biomass updated in your script. The example in section 5.2 of the Reaction Modeling Guide includes the additional code you need to include in your custom rate law (not a separate script). Additionally, I noticed that you don't even carry the biomass concentration in the rate law for the cellulose degraders, nor is there a growth yield.

 

I would use the built-in rate law for now, then consider a custom rate law when it's necessary. To get a simple single-Monod type rate law, you should set the "capital omega" (order2) to 0 to set the entire thermodynamic potential factor to 1. Or, you can use the default settings to prevent the reaction from going past the point of thermodynamic equilibrium (which in this case won't make much of a difference). Either way, it should look something like this:

# (4) Cellulose degraders
kinetic microbe-4 \
rxn = "Cellulose -> 3*CH3COO- + 3*H+", \
biomass = 1000.001, rate_con = 1e-9, KD = 1e-6, \
mpower(Celulose) = 1, mpowerD(Cellulose) = 1, \
ATP_energy = 0, ATP_number = 0, order1 = 1.0, order2 = 0, \
growth_yield = some number, decay_con = some number

Hope this helps,
Brian

Share this post


Link to post
Share on other sites

Regarding the default units of the rate constant (in the absence of promoting or inhibiting species), they depend on the type of reaction. The rate constant for a mineral dissolution reaction will have units of mol/cm^2/sec (per square cm of the mineral), but for a microbially catalyzed reaction it's mol/mg/sec (per mg of catalyzing biomass). Table 4.1 in the Reaction Modeling Guide should help explain the difference.

 

Brian

Share this post


Link to post
Share on other sites

Hi Brian,

I tried what you suggested, but it produced the following error:

"Unknown command: mpower(Cellulose)"

and then doesn't read in the rest of the information like the growth yield

Regards,

Tim

Share this post


Link to post
Share on other sites

Hi Tim,

 

It's possible there's a typo in there. First make sure you have the correct thermo dataset loaded, then try using the GUI to select your species and enter your values.

 

Regards,

Brian

Share this post


Link to post
Share on other sites

Hi Brian,

So I have a new problem that is similar to the problem that we previously discussed with the cellulose degraders.

We wish to implement a population of bacteria carrying out the reaction "H2S(aq) -> S-- + 2*H+", and have done so as follows:



# (4)
kinetic microbe-4 \
rxn = "H2S(aq) -> S-- + 2*H+", \
biomass = 0.001, rate_con = 1e-9, KD = 1e-6, \
mpower(H2S(aq)) = 1, mpowerD(H2S(aq)) = 1, \
ATP_energy = 0, ATP_number = 0, order1 = 1.0, order2 = 0, \
growth_yield = some number, decay_con = some number

However, this results in the population becoming unrealistically big because there is no limitation. This is actually also a population of anoxygenic photosynthesizers so it doesn't really make sense to use the normal thermodynamic limitation (the TPF thing) that GWB normally uses, because of the fact that a large proportion of the energy to carry out the reaction comes from light.

We would like to limit the population by some logistic term of the form (1 - biomass/ MAX_biomass) instead. I tried specifying this in a .bas script but I couldn't get it to work because I couldn't get the biomass concentration to update properly. Otherwise, we want the growth to proceed as normal, single Monod kinetics, growth yield, decay, etc.

Could you please help me out with how to implement such a limitation term?

Thanks for all of your help.

Regards,

Tim

Share this post


Link to post
Share on other sites

Hi Tim,

 

I'm out of the office this week so I probably won't be able to respond very quickly. I think you'll want to set up a custom rate law that looks something like the last example in section 5.2 in the Reaction Modeling Guide. It includes some code for updating the biomass as the catabolic reaction is carried out. You'll want to replace or modify that part with some code of your own.

 

Regards,

Brian

Share this post


Link to post
Share on other sites

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.

Guest
Reply to this topic...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.

Loading...

×
×
  • Create New...