Team:Toulouse/Modeling
Modeling
Content
Metabolic networks
The aim of our project is to create a biological system able to produce two molecules: butyric acid and formic acid.
To achieve this, we need to modify the existing balance between genes and environmental condition in E.Coli, which is used because it's one of more well known strains by scientifics.
Indeed, we want to optimize butyrate and formate productions in our bacterium
by adjusting environmental conditions in order to obtain the desired concentrations of the associated acids.
Thus, before starting modeling, it can be helpful to get a panoramic view of both all metabolic ways and specific ways in E.coli. The following metabolic network represents all of the known metabolites and metabolic pathways for Escherichia coli K12 MG1655 (best known model) as of today.
It was obtained from the KEGG database.
Our first step was to identify the pathways in which our molecules take part, in order to have a clear understanding of their role and effect.
Figure 1: Kegg Metabolic pathways - Escherichia coli K-12 MG1655
Formate network
Formate is naturally produced by E. coli but at a level that is quite low. Our project requires that Apicoli produces out highter yield.
Hence we had to optimize its biosynthesis by studying the genes coding for the enzymes involved in the pathway. We decided to focus our efforts on the Pyruvate Formate Lyase (PFL), the enzyme that causes degradation of pyruvate, thus yielding formate.
The enzyme Pyruvate Formate Lyase catalyzes the transfer of coenzyme-A on pyruvate, leading to the conversion into formate and acetyl-CoA.
Figure 2: Reaction catalyzed by the enzyme Pyruvate Formate Lyase (PFL), EC 2.3.1.54
The subnetwork presented below was obtained from MetExplore platform
and presents all reactions from the KEGG and ByoCyc databases involved in the production or consumption of formate. This map will help us predict the likely consequences of a PFL-induced formate overproduction in Apicoli.
In fact, formate is harmful to our bacterium and is normally metabolized to other products to prevent toxicity effects and redox imbalance. We thus have to find the balance between producing enough formate to kill the varroa without killing Apicoli.
Figure 3: Metabolic network of all reactions involving formate happening in E. coli
Butyrate network
Contrary to formate, butyrate is not naturally produced by E. Coli. E. coli possesses an enzyme called (Butyryl-coA transferase (or Acetate-coA transferase) that yields butyrate, but this reaction cannot happen spontaneously in the organism due to the lack of Butanoyl-coA, its substrate (Fig. 5). Indeed, study of the biosynthesis pathway shows that the enzymes responsible for Butanoyl-coA production (EC.2.1.3.19 phosphate butyryltransferase , EC.1.3.1.44 trans-2-enoyl-CoA reductase (NAD+) , etc.) cannot be found in our strain. Hence, in order to obtain butyrate, we chose to introduce a complete production pathway relying on genes coming from different organisms in Apicoli (see Attract).
Figure 4: Metabolic network of all reactions involving butyrate happening in E. coli
Figure 5: Kegg Metabolic pathways - Butanoate metabolism - Escherichia coli K-12 MG1655
If green background, enzyme described by its EC number is naturally present in E.coli. In return if blank background.
Flux Balance Analysis (FBA)
Presentation
To go further in the development of our project, we decided to use a method called Flux Balance Analysis (FBA) and Flux Variability Analysis (FVA).
It is based on the model EC_iJO1366 [1]. This is the most recent model concerning E. coli K12 MG1655.
It is a stoichiometric model defining all metabolic ways known to this day in this particular strain. It has been modificated by introducing the following pathways. Here you will find the associated XML file.
The modifications we made are such that the described strain is now capable of producing butyrate.
Our modeling aims at determining the maximum butyrate and formate quantities our strain would be able to produce, depending on two initial conditions: oxygen and glucose flux. We set them as follows:
- Maximal oxygen entry flux: 5 mmol.gDW-1.h-1
- Maximal glucose entry flux: 0.3998 mmol.gDW-1.h-1
- pH: 7
It is interesting to note that FBA provides the produced and consumed metabolites in a quantitative way.
The setpoint for oxygen flux was choosen to simulate microaerobic conditions. The glucose flux setpoint was chosen according to the results of our tests with the Biosilta kit (see “Preliminary part”).
This kit was used to ensure a stable availability of glucose over time for our bacteria. Indeed the medium contains enzymes capable of catabolizing starch, thus gradually releasing glucose in the culture.
For known initial quantities of starch and enzymes, we are able to deduce the glucose release flux. Thus, we chose the appropriate enzyme concentration and polymer quantity in order to have a glucose rate of 0.3998 mmol.gDW-1.h-1.
The model provides results for the metabolites flux in mmol.gDW-1.h-1 but it is difficult to get an idea of the actual quantity this represents so we will convert it to mmol/L. To do this we chose a length of time of 13 hours for the day and 7 hours for the night, since our solution will primarily be deployed during summer. This means that butyric acid production time is estimated to 13 hours and formic acid production time to 7 hours.
Concerning biomass it is more complicated since the bacterium grows all the time. Thus in order to ensure a minimum production, biomass concentration (X) at the beginning is defined to X = 0.56 gDW.L-1.
$$ \textrm{Real unit} (mmol.L^{-1})= \textrm{Model unit} (mmol.gDW^{-1}\cdot h^{-1})\times X \times time $$
- X: 0,56 gDW.L-1
- time (formic acid): 13h
- time butyric acid: 7h
Acid/Base Balance
Another parameter has to be taken into account, and this is the acid/base balance. Indeed, our bacteria will produce the base but we are actually interested in the acid concentration. The formula below is used:
$$pH = pKa + log\frac{C_b}{C_a}$$
- pH: the medium used is buffered so for low acid concentrations pH = 7 is considered
- pKa: 3,7 for formic acid and 4,81 for butyric acid
- Cb: base concentration
- Ca: acid concentration
Finally, it should not be forgotten that FBA methods rely on a stoichiometric model. This implies that some biological realities might be overlooked.
For example, it has been demonstrated that PFL (Pyruvate Formate Lyase), one of the enzymes involved in the production of formate, is inhibited in aerobic conditions [3], a fact that is not taken into account by the Flux Balance Analysis. This can result in the model predicting a higher formate production than what will actually be observed.
BUTYRATE
Thus under the described conditions, our first object was to try to optimize butyrate production flux. As was said before, the support file of the stoichiometric model was modified in such a way as to add the lacking butyrate biosynthesis enzymes. This was done to have a model as close to our in vivo system as possible.
As expected, when we optimize the objective function
(butyrate production), it shows an optimum when all
the carbon available is used for butyrate production,
and none goes into biomass growth. So when there is no
growth, butyrate production is equal to ~0.016 mmol.L-1, see Fig.1 (~ Flux value of 0.352 mmol.gDW-1.h-1).
Figure 1: Butyric acid produced (mmol/L) depending on growth rate
Thus under the described conditions, our first object was to try to optimize butyrate production flux. As was said before, the support file of the stoichiometric model was modified in such a way as to add the lacking butyrate biosynthesis enzymes. This was done to have a model as close to our in vivo system as possible.
If we set a minimal value for growth rate, butyrate production drops (Fig. 2), as it can be seen on the graph below obtained for different FBA simulations. We have tried different minimal growth rates up to 0.025 which is the maximal value obtained when the objective function under FBA simulation is biomass production.
Figure 2: Butyrate flux depending on minimal growth rate
Finally, to understand the effect of the initial glucose
concentration on the maximum butyrate quantity we can expect,
we tested different glucose concentrations between 0.4 and 15
mmol.gDW-1.h-1 using FVA method (Flux Variability Analysis). Indeed, 0.4 mmol.gDW-1.h-1 corresponds to the glucose flux obtained with the Biosilta Kit for a culture time of 14 days, while 15 mmol.gDW-1.h-1 is the flux needed to reach maximum growth rate (same rates are applicable to formate, Data not shown. FBA, objective function = biomass production, Result = 0.0268).
As expected, the higher the initial glucose concentration is, the higher the level of produced metabolites will be.
Figure 3: Butyrate flux depending on growth rate
FORMATE
The maximal production implies no growth and is equal to ~ 0.006 mmol.L-1, Fig. 4 (~ flux value of 2,917 mmol.gDW-1.h-1).
Figure 4: Formic acid produced for different growth rate
As with butyrate, formate biosynthesis
is unfavorably altered when a high level
of constraint is applied to the minimal
growth rate (Fig. 5). And formate production
will be more important when glucose availability is higher (Fig. 6).
Figure 5: Formate flux (mmol.gDW-1.h-1) depending on growth level
Figure 6: Formate flux (mmol.gDW-1.h-1) depending on growth rate
DIFFUSION - By KU LEUVEN Team
We designed our device basing our reflexion on practical considerations: how would it cling to the hive? How to prevent trapped varroas from getting out?...
But an important question is that of the diffusion of our molecules in this trap. Indeed it is particularly critical when referring to butyric acid since we need to have a defined concentration of it at the top of the trap to attract the varroas in.
To model this, we've collaborated with KU Leuven Team. Below are presented the results of their work, for which we are very thankful.
To know the effect our molecules, butyric and formic acid, will have on the bees, we need to calculate their concentration in the environment.
The molecules will be transported in the air by diffusion and convection.
Diffusion is always present, whereas convection depends on the presence of an external source (for example wind). To make a realistic model, we need the following parameters: diffusion coefficients and wind velocity (when we want to take it into account). The production rates of the molecules were calculated earlier (you can find them here).
The diffusion and convection can be obtained by solving the convection-diffusion equation, where D is the diffusion coefficient and u is the velocity of the solvent, in this case air, and ∇ is a gradient and ∇2 is divergence:
To solve the equations, we used COMSOL Multiphysics. It is a finite element analysis, solver and Simulation software/FEA Software package for various physics and engineering applications, especially coupled phenomena, or multiphysics. COMSOL is a very user-friendly software, which guides you step by step through the program.
The following steps need to be taken in COMSOL to calculate the transport of the molecules.
First we chose the physics. To calculate diffusion, we used Chemical Species Transport, more precisely Transport of Diluted Species in 3D. We did a time dependent study to see how the molecules spread over time. Because of the constant flux, it is also possible to use a stationary solver.
Next, we implemented the geometry of the trap. (See Fig.1) To know the effect on the environment, we can also add a block of air next to or around the trap.
Figure 1: Geometry of the trap
Diffusion coefficients can be found in literature. We used an estimation based on a calculator [4]. This resulted in a diffusion coefficient of 0.148 cm²/s for formic acid and 0.0912 cm²/s for butyric acid.
We can also take convection into account by accounting for the wind, therefore as an estimation, we checked the average velocity and direction of the wind in the proper region.
To solve the differential equation, we also need to specify boundary conditions in the calculation domain.
If we take the wind into account, the faces perpendicular to the wind are considered as in- and outflow.
The face representing the ground in the air block does not have any flux through it, because the diffusion coefficient of our molecules is much larger in soil than in air.
The three remaining faces of the air block have open boundary conditions, the flux is not known there.
All the faces of the trap have a zero flux boundary condition, except for two which will be explained below.
The face of the trap that is connected to the outside has an open boundary condition.
The lower face of the trap, where the molecules are released, has a specified constant flux. This flux can be calculated through the given production rates [4]. For formic acid, at the end of 7h we expect to have 50μmol/L, so an average rate of 7,15μmol/L/h. For butyric acid, at the end of 13h we expect to have 150μmol/L. So this corresponds to 11,5μmol/L/h.
We assume that the bacteria are in a bag with a volume of 15ml and an area of 60cm². The flux goes through a face with an area of 0.005 m². So this gives a flux of 0.00497 µmol/m²/s for formic acid and 0.00799 µmol/m²/s for butyric acid.
The initial concentration are 0 mol/m³ everywhere except for the face with the flux. Note that we put an initial concentration of 0.001 mol/m³ because we are working in a logarithmic scale.
After meshing the geometry, we can compute the results (see below) and interpret them.
Diffusion of formic acid in our trap over time.
Diffusion of butyric acid in our trap over time.
We can see on these two movies that diffusion is effective throughout our trap, for both acids.
The weak link of our trap is the funnel. It is designed to keep varroas that fell in the trap from climbing out (which would definitely be counterproductive) but we also observe that it restricts diffusion.
This is not really an issue concerning formic acid since the varroas to kill will be at the bottom of the trap. It actually would be very interesting to be able to keep formic acid in the trap to make sure that the bees will absolutely not be affected by it.
Regarding butyric acid, a diffusion that is too low could prove to be a problem. Here it seems to be acceptable still, and butyric acid is a molecule that is detectable at very low concentrations.
If we were to slightly adjust our trap's design, we would need to take this duality into account: trying to facilitate butyric acid diffusion while still making sure that formic acid does not go out in the environment.
Discussion & Optimisation
As a reminder, required acid concentrations for our system to be effective are:
- butyric acid: 0.436 mol/L (436 mmol/L) (see results in vivo)
- formic acid: 0.010 mol/L (10 mmol/L) [2],[3]
The results presented above show that we won’t be able to obtain these concentrations without optimizing our system.
In particular, we can see that to optimize either butyrate or formate production, growth level should be as low as possible. Our first solution would thus be to exert a control on this growth rate by limiting the availability of the nitrogen source in our bacteria’s environment. To test the relevance of this method we used FVA (Fig 7.).
Figure 7: NH4 flux (mmol.gDW-1.h-1) depending on growth rates
By combining this graph with the precedent results, we are able to define the nitrogen concentration required to get the wanted growth rate, which in its turn depends on the associated butyrate and formate quantities.
Another solution to optimize the metabolites
production would be to change the initial carbon source.
That’s why we tested different possible carbon sources and
optimized either formate or butyrate production as the objective
function, keeping initial glucose
flux at 0.3998 mmol.gDW-1.h-1.
Results (Fig. 8) thus show the maximal concentrations
of either butyric or formic acid we can expect
(let’s note that the method used suggests that
there should be no cellular growth).
Figure 8: Butyric and formic acid productions depending on the carbon source
According to this modeling, the use of lactose, maltose or melto-triose/-pentose/-hexaose could enable us to obtain the required formic acid concentration. More generally and with no suprise, production of formic and butyric acids will be as higher as sugar is a high carbon source.
As for butyric acid, even though the aforementioned sugars yield a higher final concentration, it is still not enough. To meet our expectations, it would be interesting to work on the parameters (size, volume…) of our device and/or our trap.
Finally, we also can adjust the pH value to increase final acids concentrations. Indeed, so far, pH has been set at 7. However, for example, if we set the pH value at 6, we can increase formic and butyric acids concentrations by a factor of 10 (Fig. 9).
Figure 9: Formic (A) and Butyric (B) acids [mmol/L] depending on pH for different growth rates.
Concentrations will be higher with low growth rates and pH decreasing.
Annex: Produced and consumed metabolites
Objective function: butyrate production
| Consommation | Production | ||||
|---|---|---|---|---|---|
| M_cl_e | Chloride | 2.476E-5 | M_5mtr_e | 5-Methylthio-D-ribose | 3.372E-5 |
| M_glc_DASH_D_e | D-Glucose | 0.3998 | M_but_e | Butyrate (n-C4:0) | 0.28610836 |
| M_mg2_e | Cu2+ | 3.37E-6 | M_co2_e | CO2 | 1.0498941 |
| M_cu2_e | Cu2+ | 3.37E-6 | M_h2o_e | H2O | 1.4649394 |
| M_cbl1_e | Cob(I)alamin | 1.115E-6 | M_amob_c | S-Adenosyl-4-methylthio-2-oxobutanoate | 1.0E-8 |
| M_pi_e | Phosphate | 0.00463573 | M_h_e | H+ | 0.0443632 |
| M_mobd_e | Molybdate | 6.85E-7 | |||
| M_zn2_e | Zinc | 1.62E-6 | |||
| M_ni2_e | nickel | 1.535E-6 | |||
| M_mn2_e | Mn2+ | 3.29E-6 | |||
| M_k_e | potassium | 9.2845E-4 | |||
| M_cobalt2_e | Co2+ | 1.2E-7 | |||
| M_so4_e | Sulfate | 0.001259915 | |||
| M_fe2_e | Fe2+ | 7.7635E-5 | |||
| M_nh4_e | Ammonium | 0.052531395 | |||
| M_o2_e | O2 | 0.75361237 | |||
| M_ca2_e | Calcium | 2.476E-5 | |||
Objective function: formate production
| Consommation | Production | ||||
|---|---|---|---|---|---|
| M_cl_e | Chloride | 2.476E-5 | M_5mtr_e | 5-Methylthio-D-ribose | 3.372E-5 |
| M_h2o_e | H2O | 0.063322984 | M_amob_c | S-Adenosyl-4-methylthio-2-oxobutanoate | 1.0E-8 |
| M_co2_e | CO2 | 0.19225993 | M_h_e | H+ | 2.4309506 |
| M_glc_DASH_D_e | D-Glucose | 0.3998 | M_5drib_c | 5'-deoxyribose | 1.155E-6 |
| M_mg2_e | magnesium | 4.1265E-5 | M_4crsol_c | p-Cresol | 1.115E-6 |
| M_cu2_e | Cu2+ | 3.37E-6 | M_mththf_c | (2R,4S)-2-methyl-2,3,3,4-tetrahydroxytetrahydrofuran | 6.7E-6 |
| M_cbl1_e | Cob(I)alamin | 1.115E-6 | M_for_e | Formate | 2.3865874 |
| M_pi_e | Phosphate | 0.00463573 | |||
| M_mobd_e | Molybdate | 6.85E-7 | |||
| M_ni2_e | nickel | 1.535E-6 | |||
| M_zn2_e | Zinc | 1.62E-6 | |||
| M_mn2_e | Mn2+ | 3.29E-6 | |||
| M_k_e | potassium | 9.2845E-4 | |||
| M_cobalt2_e | Co2+ | 1.2E-7 | |||
| M_so4_e | Sulfate | 0.001259915 | |||
| M_fe2_e | Fe2+ | 7.7635E-5 | |||
| M_nh4_e | Ammonium | 0.052531395 | |||
| M_o2_e | O2 | 0.99086045 | |||
| M_ca2_e | Calcium | 2.476E-5 | |||
References
- [1]Orth JD, Conrad TM, Na J, Lerman JA, Nam H, Feist AM & Palsson BØ (2011) A comprehensive genome-scale reconstruction of Escherichia coli metabolism--2011. Mol. Syst. Biol. 7: 535
- [2] Methods for attracting honey bee parasitic mites US 8647615 B1 See more
- [3] Imdorf, A; Charriere, J; Rosenkranz, P (1999). Varroa control with formic acid. Coordination in Europe of research on integrated control of varroa mite in honey bee colonies, Agriculture Research Centre, Merelbeke, Belgium, Commission of the European Communites. pp. 18-26. See more
- [4] Lyman, W. J., Reehl, W. F., & Rosenblatt, D. H. (1982). Handbook of chemical property estimation methods: environmental behavior of organic compounds. McGraw-Hill; http://www.epa.gov/athens/learn2model/part-two/onsite/estdiffusion-ext.html.
